1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
#include "f2c.h"
#ifdef KR_headers
double log(), f__cabs(), atan2();
#define ANSI(x) ()
#else
#define ANSI(x) x
#undef abs
#include "math.h"
#ifdef __cplusplus
extern "C" {
#endif
extern double f__cabs(double, double);
#endif
#ifndef NO_DOUBLE_EXTENDED
#ifndef GCC_COMPARE_BUG_FIXED
#ifndef Pre20000310
#ifdef Comment
Some versions of gcc, such as 2.95.3 and 3.0.4, are buggy under -O2 or -O3:
on IA32 (Intel 80x87) systems, they may do comparisons on values computed
in extended-precision registers. This can lead to the test "s > s0" that
was used below being carried out incorrectly. The fix below cannot be
spoiled by overzealous optimization, since the compiler cannot know
whether gcc_bug_bypass_diff_F2C will be nonzero. (We expect it always
to be zero. The weird name is unlikely to collide with anything.)
An example (provided by Ulrich Jakobus) where the bug fix matters is
double complex a, b
a = (.1099557428756427618354862829619, .9857360542953131909982289471372)
b = log(a)
An alternative to the fix below would be to use 53-bit rounding precision,
but the means of specifying this 80x87 feature are highly unportable.
#endif /*Comment*/
#define BYPASS_GCC_COMPARE_BUG
double (*gcc_bug_bypass_diff_F2C) ANSI((double*,double*));
static double
#ifdef KR_headers
diff1(a,b) double *a, *b;
#else
diff1(double *a, double *b)
#endif
{ return *a - *b; }
#endif /*Pre20000310*/
#endif /*GCC_COMPARE_BUG_FIXED*/
#endif /*NO_DOUBLE_EXTENDED*/
#ifdef KR_headers
VOID z_log(r, z) doublecomplex *r, *z;
#else
void z_log(doublecomplex *r, doublecomplex *z)
#endif
{
double s, s0, t, t2, u, v;
double zi = z->i, zr = z->r;
#ifdef BYPASS_GCC_COMPARE_BUG
double (*diff) ANSI((double*,double*));
#endif
r->i = atan2(zi, zr);
#ifdef Pre20000310
r->r = log( f__cabs( zr, zi ) );
#else
if (zi < 0)
zi = -zi;
if (zr < 0)
zr = -zr;
if (zr < zi) {
t = zi;
zi = zr;
zr = t;
}
t = zi/zr;
s = zr * sqrt(1 + t*t);
/* now s = f__cabs(zi,zr), and zr = |zr| >= |zi| = zi */
if ((t = s - 1) < 0)
t = -t;
if (t > .01)
r->r = log(s);
else {
#ifdef Comment
log(1+x) = x - x^2/2 + x^3/3 - x^4/4 + - ...
= x(1 - x/2 + x^2/3 -+...)
[sqrt(y^2 + z^2) - 1] * [sqrt(y^2 + z^2) + 1] = y^2 + z^2 - 1, so
sqrt(y^2 + z^2) - 1 = (y^2 + z^2 - 1) / [sqrt(y^2 + z^2) + 1]
#endif /*Comment*/
#ifdef BYPASS_GCC_COMPARE_BUG
if (!(diff = gcc_bug_bypass_diff_F2C))
diff = diff1;
#endif
t = ((zr*zr - 1.) + zi*zi) / (s + 1);
t2 = t*t;
s = 1. - 0.5*t;
u = v = 1;
do {
s0 = s;
u *= t2;
v += 2;
s += u/v - t*u/(v+1);
}
#ifdef BYPASS_GCC_COMPARE_BUG
while(s - s0 > 1e-18 || (*diff)(&s,&s0) > 0.);
#else
while(s > s0);
#endif
r->r = s*t;
}
#endif
}
#ifdef __cplusplus
}
#endif
|
