1 |
aw0a |
1 |
/*********************************************************************\ |
2 |
|
|
d1mach: Ascend Replacement for d1mach.f |
3 |
|
|
by Ben Allan |
4 |
|
|
Created: September, 1994 |
5 |
|
|
Version: $Revision: 1.3 $ |
6 |
|
|
Date last modified: $Date: 1998/07/06 10:56:12 $ |
7 |
|
|
|
8 |
|
|
This file is part of the Ascend fortran subroutine collection. |
9 |
|
|
|
10 |
|
|
Copyright (C) Benjamin Andrew Allan |
11 |
|
|
|
12 |
|
|
The ascend fortran subroutine collection is free software; you can redistribute |
13 |
|
|
it and/or modify it under the terms of the GNU General Public License as |
14 |
|
|
published by the Free Software Foundation; either version 2 of the |
15 |
|
|
License, or (at your option) any later version. |
16 |
|
|
Most of the sources in the ascend fortran subroutine collection are public |
17 |
|
|
domain and available from NETLIB. See newsgroup sci.math.numerical-analysis. |
18 |
|
|
Sources from netlib are not restricted by the GNU license and are marked as |
19 |
|
|
such. |
20 |
|
|
|
21 |
|
|
The Ascend fortran subroutine collection is distributed in hope that it will be |
22 |
|
|
useful, but WITHOUT ANY WARRANTY; without even the implied warranty of |
23 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
24 |
|
|
General Public License for more details. |
25 |
|
|
|
26 |
|
|
You should have received a copy of the GNU General Public License along with |
27 |
|
|
the program; if not, write to the Free Software Foundation, Inc., 675 |
28 |
|
|
Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING. |
29 |
|
|
COPYING is found in ../compiler. |
30 |
|
|
\*********************************************************************/ |
31 |
|
|
|
32 |
|
|
/* d1mach.c. Ben Allan |
33 |
|
|
C replacement for d1mach.f in terms of ANSI constants. |
34 |
|
|
The LINPACK d1mach.f is not such that f77 compilers pick |
35 |
|
|
the right set of constants automatically. |
36 |
|
|
|
37 |
|
|
Observed equivalences: |
38 |
|
|
F D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. |
39 |
|
|
C DBL_MIN |
40 |
|
|
F D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. |
41 |
|
|
C DBL_MAX |
42 |
|
|
F D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. |
43 |
|
|
C DBL_EPSILON/2 |
44 |
|
|
F D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. |
45 |
|
|
C DBL_EPSILON |
46 |
|
|
F D1MACH( 5) = LOG10(B) |
47 |
|
|
C NONE |
48 |
|
|
B: FLT_RADIX |
49 |
|
|
EMIN: DBL_MIN_EXP |
50 |
|
|
EMAX: DBL_MAX_EXP |
51 |
|
|
T: DBL_MANT_DIG |
52 |
|
|
|
53 |
|
|
On alphas d1mach(3)=DBL_EPSILON/2 for some reason. Returning |
54 |
|
|
DBL_EPSILON may result in 1 bit of conservatism in some codes, but |
55 |
|
|
this is the price of portability. |
56 |
|
|
|
57 |
|
|
*/ |
58 |
|
|
#include <stdlib.h> |
59 |
|
|
#include <stdio.h> |
60 |
|
|
#include <float.h> |
61 |
|
|
#include <math.h> |
62 |
|
|
|
63 |
|
|
|
64 |
|
|
/* Commentary and Apology: |
65 |
|
|
* We used to have a bunch of #ifdef's here trying to figure out which |
66 |
|
|
* platform we were on and then blessing our d1mach function with the |
67 |
|
|
* proper number of underbars (d1mach vs d1mach_) so that the linker |
68 |
|
|
* would not whine about missing symbols due to the insanity of |
69 |
|
|
* whether or not the f77 compiler puts an underbar on the symbols it |
70 |
|
|
* generates. Of course, just to make life fun, it's not strictly |
71 |
|
|
* platform dependent, since some f77 compilers accept flags that turn |
72 |
|
|
* the underbars on or off. Given this lunacy and the wasted time of |
73 |
|
|
* trying to tack down this bug every time it occurs, we've decided to |
74 |
|
|
* just duplicate the function and be finished with the underbar |
75 |
|
|
* madness. We realize this sucks and we apologize for it, but at |
76 |
|
|
* this point,``Frankly, my dears, we don't give a damn.'' |
77 |
|
|
*/ |
78 |
|
|
|
79 |
|
|
double d1mach(int *i) { |
80 |
|
|
switch (*i) { |
81 |
|
|
case 1: |
82 |
|
|
return DBL_MIN; |
83 |
|
|
case 2: |
84 |
|
|
return DBL_MAX; |
85 |
|
|
case 3: |
86 |
|
|
return DBL_EPSILON; |
87 |
|
|
case 4: |
88 |
|
|
return DBL_EPSILON; |
89 |
|
|
case 5: |
90 |
|
|
return log10((double)FLT_RADIX); |
91 |
|
|
default: |
92 |
|
|
fprintf(stderr," D1MACH - I OUT OF BOUNDS %d",*i); |
93 |
|
|
abort(); |
94 |
|
|
} |
95 |
|
|
} |
96 |
|
|
|
97 |
|
|
double d1mach_(int *i) { |
98 |
|
|
return d1mach(i); |
99 |
|
|
} |
100 |
|
|
|
101 |
|
|
double D1MACH(int *i) { |
102 |
|
|
return d1mach(i); |
103 |
|
|
} |
104 |
|
|
|
105 |
|
|
|
106 |
|
|
|