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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
static double numarray_zero = 0.0;
static double raiseDivByZero(void)
{
return 1.0/numarray_zero;
}
static double raiseNegDivByZero(void)
{
return -1.0/numarray_zero;
}
static double num_log(double x)
{
if (x == 0.0)
return raiseNegDivByZero();
else
return log(x);
}
static double num_log10(double x)
{
if (x == 0.0)
return raiseNegDivByZero();
else
return log10(x);
}
static double num_pow(double x, double y)
{
int z = (int) y;
if ((x < 0.0) && (y != z))
return raiseDivByZero();
else
return pow(x, y);
}
/* Inverse hyperbolic trig functions from Numeric */
static double num_acosh(double x)
{
return log(x + sqrt((x-1.0)*(x+1.0)));
}
static double num_asinh(double xx)
{
double x;
int sign;
if (xx < 0.0) {
sign = -1;
x = -xx;
}
else {
sign = 1;
x = xx;
}
return sign*log(x + sqrt(x*x+1.0));
}
static double num_atanh(double x)
{
return 0.5*log((1.0+x)/(1.0-x));
}
/* NUM_CROUND (in numcomplex.h) also calls num_round */
static double num_round(double x)
{
return (x >= 0) ? floor(x+0.5) : ceil(x-0.5);
}
/* The following routine is used in the event of a detected integer *
** divide by zero so that a floating divide by zero is generated. *
** This is done since numarray uses the floating point exception *
** sticky bits to detect errors. The last bit is an attempt to *
** prevent optimization of the divide by zero away, the input value *
** should always be 0 *
*/
static int int_dividebyzero_error(long value, long unused) {
double dummy;
dummy = 1./numarray_zero;
if (dummy) /* to prevent optimizer from eliminating expression */
return 0;
else
return 1;
}
/* Likewise for Integer overflows */
static int int_overflow_error(Float64 value) {
double dummy;
dummy = pow(1.e10, fabs(value/2));
if (dummy) /* to prevent optimizer from eliminating expression */
return (int) value;
else
return 1;
}
static int umult64_overflow(UInt64 a, UInt64 b)
{
UInt64 ah, al, bh, bl, w, x, y, z;
ah = (a >> 32);
al = (a & 0xFFFFFFFFL);
bh = (b >> 32);
bl = (b & 0xFFFFFFFFL);
/* 128-bit product: z*2**64 + (x+y)*2**32 + w */
w = al*bl;
x = bh*al;
y = ah*bl;
z = ah*bh;
/* *c = ((x + y)<<32) + w; */
return z || (x>>32) || (y>>32) ||
(((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 32);
}
static int smult64_overflow(Int64 a0, Int64 b0)
{
UInt64 a, b;
UInt64 ah, al, bh, bl, w, x, y, z;
/* Convert to non-negative quantities */
if (a0 < 0) { a = -a0; } else { a = a0; }
if (b0 < 0) { b = -b0; } else { b = b0; }
ah = (a >> 32);
al = (a & 0xFFFFFFFFL);
bh = (b >> 32);
bl = (b & 0xFFFFFFFFL);
w = al*bl;
x = bh*al;
y = ah*bl;
z = ah*bh;
/*
UInt64 c = ((x + y)<<32) + w;
if ((a0 < 0) ^ (b0 < 0))
*c = -c;
else
*c = c
*/
return z || (x>>31) || (y>>31) ||
(((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 31);
}
|
