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
|
/* isnan()
* signbit()
* isfinite()
*
* Floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* double ceil(), floor(), frexp(), ldexp(); --- gone
* int signbit(), isnan(), isfinite();
* double x, y;
* int expnt, n;
*
* y = floor(x); -gone
* y = ceil(x); -gone
* y = frexp( x, &expnt ); -gone
* y = ldexp( x, n ); -gone
* n = signbit(x);
* n = isnan(x);
* n = isfinite(x);
*
*
*
* DESCRIPTION:
*
* All four routines return a double precision floating point
* result.
*
* floor() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* ceil() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* frexp() extracts the exponent from x. It returns an integer
* power of two to expnt and the significand between 0.5 and 1
* to y. Thus x = y * 2**expn.
*
* ldexp() multiplies x by 2**n.
*
* signbit(x) returns 1 if the sign bit of x is 1, else 0.
*
* These functions are part of the standard C run time library
* for many but not all C compilers. The ones supplied are
* written in C for either DEC or IEEE arithmetic. They should
* be used only if your compiler library does not already have
* them.
*
* The IEEE versions assume that denormal numbers are implemented
* in the arithmetic. Some modifications will be required if
* the arithmetic has abrupt rather than gradual underflow.
*/
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include <Python.h>
#include <numpy/ndarrayobject.h>
#include "mconf.h"
/* XXX: horrible hacks, but those cephes macros are buggy and just plain ugly anywa.
* We should use npy_* macros instead once npy_math can be used reliably by
* packages outside numpy
*/
#undef isnan
#undef signbit
#undef isfinite
#define isnan(x) ((x) != (x))
int cephes_isnan(double x)
{
return isnan(x);
}
int isfinite(double x)
{
return !isnan((x) + (-x));
}
static int isbigendian(void)
{
const union {
npy_uint32 i;
char c[4];
} bint = {0x01020304};
if (bint.c[0] == 1) {
return 1;
}
return 0;
}
int signbit(double x)
{
union
{
double d;
short s[4];
int i[2];
} u;
u.d = x;
/*
* Tuis is stupid, we test for endianness every time, but that the easiest
* way I can see without using platform checks - for scipy 0.8.0, we should
* use npy_math
*/
#if SIZEOF_INT == 4
if (isbigendian()) {
return u.i[1] < 0;
} else {
return u.i[0] < 0;
}
#else /* SIZEOF_INT != 4 */
if (isbigendian()) {
return u.s[3] < 0;
} else {
return u.s[0] < 0;
}
#endif /* SIZEOF_INT */
}
|
