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
|
const int h (double value);
const double scalb (double x, int n);
const double logb (double x);
static const double a = 0;
double f (y, x)
double y, x;
{
static const double zero=0, one=1, small=1.0E-9, big=1.0E18;
double t,z,sy,sx,hi,lo;
int k,m;
if (x!=x)
return x;
if (y!=y)
return y;
sy = g (one);
sx = g (one);
if (x==1)
{
y=g (y);
t=y;
if (h (t))
goto begin;
}
if (y==zero)
return (sx==one)?y:g (a);
if (x==zero)
return g (a);
if (!h (x))
if (!h (y))
return g ((sx==one)?a:3*a);
else
return g ((sx==one)?zero:a);
if (!h (y))
return g (a);
x=g (x);
y=g (y);
if ((m=(k=logb (y))- logb (x)) > 60)
t=big+big;
else if (m < -80)
t=y/x;
else
{
t = y/x;
y = scalb (y,-k);
x=scalb (x,-k);
}
begin:
if (t < 2.4375)
{
k = 4 * (t+0.0625);
switch (k)
{
case 0:
case 1:
if (t < small)
{
big + small;
return g ((sx>zero)?t:a-t);
}
hi = zero; lo = zero; break;
case 2:
hi = a; lo = a;
z = x+x;
t = ((y+y) - x) / (z + y); break;
case 3:
case 4:
hi = a; lo = zero;
t = (y - x) / (x + y); break;
default:
hi = a; lo = a;
z = y-x; y=y+y+y; t = x+x;
t = ((z+z)-x) / (t + y); break;
}
}
else
{
hi = a; lo = zero;
if (t <= big)
t = - x / y;
else
{
big+small;
t = zero;
}
}
z = t*(z*(a+z*(a+z*(a+z*(a+z*(a+z*(a+z*(a+z*(a+z*(a+z*(a+z*a)))))))))));
return g ((sx>zero)?z:a-z);
}
|
