|
1
|
|
/*
|
|
2
|
|
* vim:syntax=c
|
|
3
|
|
* A small module to implement missing C99 math capabilities required by numpy
|
|
4
|
|
*
|
|
5
|
|
* Please keep this independent of python ! Only basic types (npy_longdouble)
|
|
6
|
|
* can be used, otherwise, pure C, without any use of Python facilities
|
|
7
|
|
*
|
|
8
|
|
* How to add a function to this section
|
|
9
|
|
* -------------------------------------
|
|
10
|
|
*
|
|
11
|
|
* Say you want to add `foo`, these are the steps and the reasons for them.
|
|
12
|
|
*
|
|
13
|
|
* 1) Add foo to the appropriate list in the configuration system. The
|
|
14
|
|
* lists can be found in numpy/core/setup.py lines 63-105. Read the
|
|
15
|
|
* comments that come with them, they are very helpful.
|
|
16
|
|
*
|
|
17
|
|
* 2) The configuration system will define a macro HAVE_FOO if your function
|
|
18
|
|
* can be linked from the math library. The result can depend on the
|
|
19
|
|
* optimization flags as well as the compiler, so can't be known ahead of
|
|
20
|
|
* time. If the function can't be linked, then either it is absent, defined
|
|
21
|
|
* as a macro, or is an intrinsic (hardware) function.
|
|
22
|
|
*
|
|
23
|
|
* i) Undefine any possible macros:
|
|
24
|
|
*
|
|
25
|
|
* #ifdef foo
|
|
26
|
|
* #undef foo
|
|
27
|
|
* #endif
|
|
28
|
|
*
|
|
29
|
|
* ii) Avoid as much as possible to declare any function here. Declaring
|
|
30
|
|
* functions is not portable: some platforms define some function inline
|
|
31
|
|
* with a non standard identifier, for example, or may put another
|
|
32
|
|
* identifier which changes the calling convention of the function. If you
|
|
33
|
|
* really have to, ALWAYS declare it for the one platform you are dealing
|
|
34
|
|
* with:
|
|
35
|
|
*
|
|
36
|
|
* Not ok:
|
|
37
|
|
* double exp(double a);
|
|
38
|
|
*
|
|
39
|
|
* Ok:
|
|
40
|
|
* #ifdef SYMBOL_DEFINED_WEIRD_PLATFORM
|
|
41
|
|
* double exp(double);
|
|
42
|
|
* #endif
|
|
43
|
|
*
|
|
44
|
|
* Some of the code is taken from msun library in FreeBSD, with the following
|
|
45
|
|
* notice:
|
|
46
|
|
*
|
|
47
|
|
* ====================================================
|
|
48
|
|
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
49
|
|
*
|
|
50
|
|
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
51
|
|
* Permission to use, copy, modify, and distribute this
|
|
52
|
|
* software is freely granted, provided that this notice
|
|
53
|
|
* is preserved.
|
|
54
|
|
* ====================================================
|
|
55
|
|
*/
|
|
56
|
|
#include "npy_math_private.h"
|
|
57
|
|
|
|
58
|
|
/*
|
|
59
|
|
*****************************************************************************
|
|
60
|
|
** BASIC MATH FUNCTIONS **
|
|
61
|
|
*****************************************************************************
|
|
62
|
|
*/
|
|
63
|
|
|
|
64
|
|
/* Original code by Konrad Hinsen. */
|
|
65
|
|
#ifndef HAVE_EXPM1
|
|
66
|
|
NPY_INPLACE double npy_expm1(double x)
|
|
67
|
|
{
|
|
68
|
|
if (npy_isinf(x) && x > 0) {
|
|
69
|
|
return x;
|
|
70
|
|
}
|
|
71
|
|
else {
|
|
72
|
|
const double u = npy_exp(x);
|
|
73
|
|
|
|
74
|
|
if (u == 1.0) {
|
|
75
|
|
return x;
|
|
76
|
|
} else if (u - 1.0 == -1.0) {
|
|
77
|
|
return -1;
|
|
78
|
|
} else {
|
|
79
|
|
return (u - 1.0) * x/npy_log(u);
|
|
80
|
|
}
|
|
81
|
|
}
|
|
82
|
|
}
|
|
83
|
|
#endif
|
|
84
|
|
|
|
85
|
|
#ifndef HAVE_LOG1P
|
|
86
|
|
NPY_INPLACE double npy_log1p(double x)
|
|
87
|
|
{
|
|
88
|
|
if (npy_isinf(x) && x > 0) {
|
|
89
|
|
return x;
|
|
90
|
|
}
|
|
91
|
|
else {
|
|
92
|
|
const double u = 1. + x;
|
|
93
|
|
const double d = u - 1.;
|
|
94
|
|
|
|
95
|
|
if (d == 0) {
|
|
96
|
|
return x;
|
|
97
|
|
} else {
|
|
98
|
|
return npy_log(u) * x / d;
|
|
99
|
|
}
|
|
100
|
|
}
|
|
101
|
|
}
|
|
102
|
|
#endif
|
|
103
|
|
|
|
104
|
|
/* Taken from FreeBSD mlib, adapted for numpy
|
|
105
|
|
*
|
|
106
|
|
* XXX: we could be a bit faster by reusing high/low words for inf/nan
|
|
107
|
|
* classification instead of calling npy_isinf/npy_isnan: we should have some
|
|
108
|
|
* macros for this, though, instead of doing it manually
|
|
109
|
|
*/
|
|
110
|
|
#ifndef HAVE_ATAN2
|
|
111
|
|
/* XXX: we should have this in npy_math.h */
|
|
112
|
|
#define NPY_DBL_EPSILON 1.2246467991473531772E-16
|
|
113
|
|
NPY_INPLACE double npy_atan2(double y, double x)
|
|
114
|
|
{
|
|
115
|
|
npy_int32 k, m, iy, ix, hx, hy;
|
|
116
|
|
npy_uint32 lx,ly;
|
|
117
|
|
double z;
|
|
118
|
|
|
|
119
|
|
EXTRACT_WORDS(hx, lx, x);
|
|
120
|
|
ix = hx & 0x7fffffff;
|
|
121
|
|
EXTRACT_WORDS(hy, ly, y);
|
|
122
|
|
iy = hy & 0x7fffffff;
|
|
123
|
|
|
|
124
|
|
/* if x or y is nan, return nan */
|
|
125
|
|
if (npy_isnan(x * y)) {
|
|
126
|
|
return x + y;
|
|
127
|
|
}
|
|
128
|
|
|
|
129
|
|
if (x == 1.0) {
|
|
130
|
|
return npy_atan(y);
|
|
131
|
|
}
|
|
132
|
|
|
|
133
|
|
m = 2 * (npy_signbit((x)) != 0) + (npy_signbit((y)) != 0);
|
|
134
|
|
if (y == 0.0) {
|
|
135
|
|
switch(m) {
|
|
136
|
|
case 0:
|
|
137
|
|
case 1: return y; /* atan(+-0,+anything)=+-0 */
|
|
138
|
|
case 2: return NPY_PI;/* atan(+0,-anything) = pi */
|
|
139
|
|
case 3: return -NPY_PI;/* atan(-0,-anything) =-pi */
|
|
140
|
|
}
|
|
141
|
|
}
|
|
142
|
|
|
|
143
|
|
if (x == 0.0) {
|
|
144
|
|
return y > 0 ? NPY_PI_2 : -NPY_PI_2;
|
|
145
|
|
}
|
|
146
|
|
|
|
147
|
|
if (npy_isinf(x)) {
|
|
148
|
|
if (npy_isinf(y)) {
|
|
149
|
|
switch(m) {
|
|
150
|
|
case 0: return NPY_PI_4;/* atan(+INF,+INF) */
|
|
151
|
|
case 1: return -NPY_PI_4;/* atan(-INF,+INF) */
|
|
152
|
|
case 2: return 3.0*NPY_PI_4;/*atan(+INF,-INF)*/
|
|
153
|
|
case 3: return -3.0*NPY_PI_4;/*atan(-INF,-INF)*/
|
|
154
|
|
}
|
|
155
|
|
} else {
|
|
156
|
|
switch(m) {
|
|
157
|
|
case 0: return NPY_PZERO; /* atan(+...,+INF) */
|
|
158
|
|
case 1: return NPY_NZERO; /* atan(-...,+INF) */
|
|
159
|
|
case 2: return NPY_PI; /* atan(+...,-INF) */
|
|
160
|
|
case 3: return -NPY_PI; /* atan(-...,-INF) */
|
|
161
|
|
}
|
|
162
|
|
}
|
|
163
|
|
}
|
|
164
|
|
|
|
165
|
|
if (npy_isinf(y)) {
|
|
166
|
|
return y > 0 ? NPY_PI_2 : -NPY_PI_2;
|
|
167
|
|
}
|
|
168
|
|
|
|
169
|
|
/* compute y/x */
|
|
170
|
|
k = (iy - ix) >> 20;
|
|
171
|
|
if (k > 60) { /* |y/x| > 2**60 */
|
|
172
|
|
z = NPY_PI_2 + 0.5 * NPY_DBL_EPSILON;
|
|
173
|
|
m &= 1;
|
|
174
|
|
} else if (hx < 0 && k < -60) {
|
|
175
|
|
z = 0.0; /* 0 > |y|/x > -2**-60 */
|
|
176
|
|
} else {
|
|
177
|
|
z = npy_atan(npy_fabs(y/x)); /* safe to do y/x */
|
|
178
|
|
}
|
|
179
|
|
|
|
180
|
|
switch (m) {
|
|
181
|
|
case 0: return z ; /* atan(+,+) */
|
|
182
|
|
case 1: return -z ; /* atan(-,+) */
|
|
183
|
|
case 2: return NPY_PI - (z - NPY_DBL_EPSILON);/* atan(+,-) */
|
|
184
|
|
default: /* case 3 */
|
|
185
|
|
return (z - NPY_DBL_EPSILON) - NPY_PI;/* atan(-,-) */
|
|
186
|
|
}
|
|
187
|
|
}
|
|
188
|
|
|
|
189
|
|
#endif
|
|
190
|
|
|
|
191
|
|
#ifndef HAVE_HYPOT
|
|
192
|
|
NPY_INPLACE double npy_hypot(double x, double y)
|
|
193
|
|
{
|
|
194
|
|
double yx;
|
|
195
|
|
|
|
196
|
|
if (npy_isinf(x) || npy_isinf(y)) {
|
|
197
|
|
return NPY_INFINITY;
|
|
198
|
|
}
|
|
199
|
|
|
|
200
|
|
if (npy_isnan(x) || npy_isnan(y)) {
|
|
201
|
|
return NPY_NAN;
|
|
202
|
|
}
|
|
203
|
|
|
|
204
|
|
x = npy_fabs(x);
|
|
205
|
|
y = npy_fabs(y);
|
|
206
|
|
if (x < y) {
|
|
207
|
|
double temp = x;
|
|
208
|
|
x = y;
|
|
209
|
|
y = temp;
|
|
210
|
|
}
|
|
211
|
|
if (x == 0.) {
|
|
212
|
|
return 0.;
|
|
213
|
|
}
|
|
214
|
|
else {
|
|
215
|
|
yx = y/x;
|
|
216
|
|
return x*npy_sqrt(1.+yx*yx);
|
|
217
|
|
}
|
|
218
|
|
}
|
|
219
|
|
#endif
|
|
220
|
|
|
|
221
|
|
#ifndef HAVE_ACOSH
|
|
222
|
|
NPY_INPLACE double npy_acosh(double x)
|
|
223
|
|
{
|
|
224
|
|
if (x < 1.0) {
|
|
225
|
|
return NPY_NAN;
|
|
226
|
|
}
|
|
227
|
|
|
|
228
|
|
if (npy_isfinite(x)) {
|
|
229
|
|
if (x > 1e8) {
|
|
230
|
|
return npy_log(x) + NPY_LOGE2;
|
|
231
|
|
}
|
|
232
|
|
else {
|
|
233
|
|
double u = x - 1.0;
|
|
234
|
|
return npy_log1p(u + npy_sqrt(2*u + u*u));
|
|
235
|
|
}
|
|
236
|
|
}
|
|
237
|
|
return x;
|
|
238
|
|
}
|
|
239
|
|
#endif
|
|
240
|
|
|
|
241
|
|
#ifndef HAVE_ASINH
|
|
242
|
|
NPY_INPLACE double npy_asinh(double xx)
|
|
243
|
|
{
|
|
244
|
|
double x, d;
|
|
245
|
|
int sign;
|
|
246
|
|
if (xx < 0.0) {
|
|
247
|
|
sign = -1;
|
|
248
|
|
x = -xx;
|
|
249
|
|
}
|
|
250
|
|
else {
|
|
251
|
|
sign = 1;
|
|
252
|
|
x = xx;
|
|
253
|
|
}
|
|
254
|
|
if (x > 1e8) {
|
|
255
|
|
d = x;
|
|
256
|
|
} else {
|
|
257
|
|
d = npy_sqrt(x*x + 1);
|
|
258
|
|
}
|
|
259
|
|
return sign*npy_log1p(x*(1.0 + x/(d+1)));
|
|
260
|
|
}
|
|
261
|
|
#endif
|
|
262
|
|
|
|
263
|
|
#ifndef HAVE_ATANH
|
|
264
|
|
NPY_INPLACE double npy_atanh(double x)
|
|
265
|
|
{
|
|
266
|
|
if (x > 0) {
|
|
267
|
|
return -0.5*npy_log1p(-2.0*x/(1.0 + x));
|
|
268
|
|
}
|
|
269
|
|
else {
|
|
270
|
|
return 0.5*npy_log1p(2.0*x/(1.0 - x));
|
|
271
|
|
}
|
|
272
|
|
}
|
|
273
|
|
#endif
|
|
274
|
|
|
|
275
|
|
#ifndef HAVE_RINT
|
|
276
|
|
#if defined(_MSC_VER) && (_MSC_VER == 1500) && !defined(_WIN64)
|
|
277
|
|
#pragma optimize("", off)
|
|
278
|
|
#endif
|
|
279
|
|
NPY_INPLACE double npy_rint(double x)
|
|
280
|
|
{
|
|
281
|
|
double y, r;
|
|
282
|
|
|
|
283
|
|
y = npy_floor(x);
|
|
284
|
|
r = x - y;
|
|
285
|
|
|
|
286
|
|
if (r > 0.5) {
|
|
287
|
|
y += 1.0;
|
|
288
|
|
}
|
|
289
|
|
|
|
290
|
|
/* Round to nearest even */
|
|
291
|
|
if (r == 0.5) {
|
|
292
|
|
r = y - 2.0*npy_floor(0.5*y);
|
|
293
|
|
if (r == 1.0) {
|
|
294
|
|
y += 1.0;
|
|
295
|
|
}
|
|
296
|
|
}
|
|
297
|
|
return y;
|
|
298
|
|
}
|
|
299
|
|
#if defined(_MSC_VER) && (_MSC_VER == 1500) && !defined(_WIN64)
|
|
300
|
|
#pragma optimize("", on)
|
|
301
|
|
#endif
|
|
302
|
|
#endif
|
|
303
|
|
|
|
304
|
|
#ifndef HAVE_TRUNC
|
|
305
|
|
NPY_INPLACE double npy_trunc(double x)
|
|
306
|
|
{
|
|
307
|
|
return x < 0 ? npy_ceil(x) : npy_floor(x);
|
|
308
|
|
}
|
|
309
|
|
#endif
|
|
310
|
|
|
|
311
|
|
#ifndef HAVE_EXP2
|
|
312
|
|
NPY_INPLACE double npy_exp2(double x)
|
|
313
|
|
{
|
|
314
|
|
return npy_exp(NPY_LOGE2*x);
|
|
315
|
|
}
|
|
316
|
|
#endif
|
|
317
|
|
|
|
318
|
|
#ifndef HAVE_LOG2
|
|
319
|
|
NPY_INPLACE double npy_log2(double x)
|
|
320
|
|
{
|
|
321
|
|
#ifdef HAVE_FREXP
|
|
322
|
|
if (!npy_isfinite(x) || x <= 0.) {
|
|
323
|
|
/* special value result */
|
|
324
|
|
return npy_log(x);
|
|
325
|
|
}
|
|
326
|
|
else {
|
|
327
|
|
/*
|
|
328
|
|
* fallback implementation copied from python3.4 math.log2
|
|
329
|
|
* provides int(log(2**i)) == i for i 1-64 in default rounding mode.
|
|
330
|
|
*
|
|
331
|
|
* We want log2(m * 2**e) == log(m) / log(2) + e. Care is needed when
|
|
332
|
|
* x is just greater than 1.0: in that case e is 1, log(m) is negative,
|
|
333
|
|
* and we get significant cancellation error from the addition of
|
|
334
|
|
* log(m) / log(2) to e. The slight rewrite of the expression below
|
|
335
|
|
* avoids this problem.
|
|
336
|
|
*/
|
|
337
|
|
int e;
|
|
338
|
|
double m = frexp(x, &e);
|
|
339
|
|
if (x >= 1.0) {
|
|
340
|
|
return log(2.0 * m) / log(2.0) + (e - 1);
|
|
341
|
|
}
|
|
342
|
|
else {
|
|
343
|
|
return log(m) / log(2.0) + e;
|
|
344
|
|
}
|
|
345
|
|
}
|
|
346
|
|
#else
|
|
347
|
|
/* does not provide int(log(2**i)) == i */
|
|
348
|
|
return NPY_LOG2E * npy_log(x);
|
|
349
|
|
#endif
|
|
350
|
|
}
|
|
351
|
|
#endif
|
|
352
|
|
|
|
353
|
|
/*
|
|
354
|
|
* if C99 extensions not available then define dummy functions that use the
|
|
355
|
|
* double versions for
|
|
356
|
|
*
|
|
357
|
|
* sin, cos, tan
|
|
358
|
|
* sinh, cosh, tanh,
|
|
359
|
|
* fabs, floor, ceil, rint, trunc
|
|
360
|
|
* sqrt, log10, log, exp, expm1
|
|
361
|
|
* asin, acos, atan,
|
|
362
|
|
* asinh, acosh, atanh
|
|
363
|
|
*
|
|
364
|
|
* hypot, atan2, pow, fmod, modf
|
|
365
|
|
* ldexp, frexp
|
|
366
|
|
*
|
|
367
|
|
* We assume the above are always available in their double versions.
|
|
368
|
|
*
|
|
369
|
|
* NOTE: some facilities may be available as macro only instead of functions.
|
|
370
|
|
* For simplicity, we define our own functions and undef the macros. We could
|
|
371
|
|
* instead test for the macro, but I am lazy to do that for now.
|
|
372
|
|
*/
|
|
373
|
|
|
|
374
|
|
/**begin repeat
|
|
375
|
|
* #type = npy_longdouble, npy_float#
|
|
376
|
|
* #TYPE = NPY_LONGDOUBLE, FLOAT#
|
|
377
|
|
* #c = l,f#
|
|
378
|
|
* #C = L,F#
|
|
379
|
|
*/
|
|
380
|
|
|
|
381
|
|
/**begin repeat1
|
|
382
|
|
* #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
|
|
383
|
|
* log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2#
|
|
384
|
|
* #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10,
|
|
385
|
|
* LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2#
|
|
386
|
|
*/
|
|
387
|
|
|
|
388
|
|
#ifdef @kind@@c@
|
|
389
|
|
#undef @kind@@c@
|
|
390
|
|
#endif
|
|
391
|
|
#ifndef HAVE_@KIND@@C@
|
|
392
|
|
NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
|
|
393
|
|
{
|
|
394
|
|
return (@type@) npy_@kind@((double)x);
|
|
395
|
|
}
|
|
396
|
|
#endif
|
|
397
|
|
|
|
398
|
|
/**end repeat1**/
|
|
399
|
|
|
|
400
|
|
/**begin repeat1
|
|
401
|
|
* #kind = atan2,hypot,pow,fmod,copysign#
|
|
402
|
|
* #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN#
|
|
403
|
|
*/
|
|
404
|
|
#ifdef @kind@@c@
|
|
405
|
|
#undef @kind@@c@
|
|
406
|
|
#endif
|
|
407
|
|
#ifndef HAVE_@KIND@@C@
|
|
408
|
|
NPY_INPLACE @type@ npy_@kind@@c@(@type@ x, @type@ y)
|
|
409
|
|
{
|
|
410
|
|
return (@type@) npy_@kind@((double)x, (double) y);
|
|
411
|
|
}
|
|
412
|
|
#endif
|
|
413
|
|
/**end repeat1**/
|
|
414
|
|
|
|
415
|
|
#ifdef modf@c@
|
|
416
|
|
#undef modf@c@
|
|
417
|
|
#endif
|
|
418
|
|
#ifndef HAVE_MODF@C@
|
|
419
|
|
NPY_INPLACE @type@ npy_modf@c@(@type@ x, @type@ *iptr)
|
|
420
|
|
{
|
|
421
|
|
double niptr;
|
|
422
|
|
double y = npy_modf((double)x, &niptr);
|
|
423
|
|
*iptr = (@type@) niptr;
|
|
424
|
|
return (@type@) y;
|
|
425
|
|
}
|
|
426
|
|
#endif
|
|
427
|
|
|
|
428
|
|
#ifdef ldexp@c@
|
|
429
|
|
#undef ldexp@c@
|
|
430
|
|
#endif
|
|
431
|
|
#ifndef HAVE_LDEXP@C@
|
|
432
|
|
NPY_INPLACE @type@ npy_ldexp@c@(@type@ x, int exp)
|
|
433
|
|
{
|
|
434
|
|
return (@type@) npy_ldexp((double)x, exp);
|
|
435
|
|
}
|
|
436
|
|
#endif
|
|
437
|
|
|
|
438
|
|
#ifdef frexp@c@
|
|
439
|
|
#undef frexp@c@
|
|
440
|
|
#endif
|
|
441
|
|
#ifndef HAVE_FREXP@C@
|
|
442
|
|
NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
|
|
443
|
|
{
|
|
444
|
|
return (@type@) npy_frexp(x, exp);
|
|
445
|
|
}
|
|
446
|
|
#endif
|
|
447
|
|
|
|
448
|
|
/**end repeat**/
|
|
449
|
|
|
|
450
|
|
|
|
451
|
|
/*
|
|
452
|
|
* Decorate all the math functions which are available on the current platform
|
|
453
|
|
*/
|
|
454
|
|
|
|
455
|
|
/**begin repeat
|
|
456
|
|
* #type = npy_longdouble, npy_double, npy_float#
|
|
457
|
|
* #c = l,,f#
|
|
458
|
|
* #C = L,,F#
|
|
459
|
|
*/
|
|
460
|
|
/**begin repeat1
|
|
461
|
|
* #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
|
|
462
|
|
* log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2#
|
|
463
|
|
* #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10,
|
|
464
|
|
* LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2#
|
|
465
|
|
*/
|
|
466
|
|
#ifdef HAVE_@KIND@@C@
|
|
467
|
1
|
NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
|
|
468
|
|
{
|
|
469
|
1
|
return @kind@@c@(x);
|
|
470
|
|
}
|
|
471
|
|
#endif
|
|
472
|
|
|
|
473
|
|
/**end repeat1**/
|
|
474
|
|
|
|
475
|
|
/**begin repeat1
|
|
476
|
|
* #kind = atan2,hypot,pow,fmod,copysign#
|
|
477
|
|
* #KIND = ATAN2,HYPOT,POW,FMOD,COPYSIGN#
|
|
478
|
|
*/
|
|
479
|
|
#ifdef HAVE_@KIND@@C@
|
|
480
|
1
|
NPY_INPLACE @type@ npy_@kind@@c@(@type@ x, @type@ y)
|
|
481
|
|
{
|
|
482
|
1
|
return @kind@@c@(x, y);
|
|
483
|
|
}
|
|
484
|
|
#endif
|
|
485
|
|
/**end repeat1**/
|
|
486
|
|
|
|
487
|
|
#ifdef HAVE_MODF@C@
|
|
488
|
0
|
NPY_INPLACE @type@ npy_modf@c@(@type@ x, @type@ *iptr)
|
|
489
|
|
{
|
|
490
|
1
|
return modf@c@(x, iptr);
|
|
491
|
|
}
|
|
492
|
|
#endif
|
|
493
|
|
|
|
494
|
|
#ifdef HAVE_LDEXP@C@
|
|
495
|
0
|
NPY_INPLACE @type@ npy_ldexp@c@(@type@ x, int exp)
|
|
496
|
|
{
|
|
497
|
1
|
return ldexp@c@(x, exp);
|
|
498
|
|
}
|
|
499
|
|
#endif
|
|
500
|
|
|
|
501
|
|
#ifdef HAVE_FREXP@C@
|
|
502
|
0
|
NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
|
|
503
|
|
{
|
|
504
|
1
|
return frexp@c@(x, exp);
|
|
505
|
|
}
|
|
506
|
|
#endif
|
|
507
|
|
|
|
508
|
|
/* C99 but not mandatory */
|
|
509
|
|
|
|
510
|
|
#ifndef HAVE_CBRT@C@
|
|
511
|
|
NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
|
|
512
|
|
{
|
|
513
|
|
/* don't set invalid flag */
|
|
514
|
|
if (npy_isnan(x)) {
|
|
515
|
|
return NPY_NAN;
|
|
516
|
|
}
|
|
517
|
|
else if (x < 0) {
|
|
518
|
|
return -npy_pow@c@(-x, 1. / 3.);
|
|
519
|
|
}
|
|
520
|
|
else {
|
|
521
|
|
return npy_pow@c@(x, 1. / 3.);
|
|
522
|
|
}
|
|
523
|
|
}
|
|
524
|
|
#else
|
|
525
|
1
|
NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
|
|
526
|
|
{
|
|
527
|
1
|
return cbrt@c@(x);
|
|
528
|
|
}
|
|
529
|
|
#endif
|
|
530
|
|
|
|
531
|
|
/**end repeat**/
|
|
532
|
|
|
|
533
|
|
|
|
534
|
|
/*
|
|
535
|
|
* Non standard functions
|
|
536
|
|
*/
|
|
537
|
|
|
|
538
|
|
/**begin repeat
|
|
539
|
|
* #type = npy_float, npy_double, npy_longdouble#
|
|
540
|
|
* #c = f, ,l#
|
|
541
|
|
* #C = F, ,L#
|
|
542
|
|
*/
|
|
543
|
|
|
|
544
|
1
|
@type@ npy_heaviside@c@(@type@ x, @type@ h0)
|
|
545
|
|
{
|
|
546
|
1
|
if (npy_isnan(x)) {
|
|
547
|
|
return (@type@) NPY_NAN;
|
|
548
|
|
}
|
|
549
|
1
|
else if (x == 0) {
|
|
550
|
|
return h0;
|
|
551
|
|
}
|
|
552
|
1
|
else if (x < 0) {
|
|
553
|
|
return (@type@) 0.0;
|
|
554
|
|
}
|
|
555
|
|
else {
|
|
556
|
1
|
return (@type@) 1.0;
|
|
557
|
|
}
|
|
558
|
|
}
|
|
559
|
|
|
|
560
|
|
#define LOGE2 NPY_LOGE2@c@
|
|
561
|
|
#define LOG2E NPY_LOG2E@c@
|
|
562
|
|
#define RAD2DEG (180.0@c@/NPY_PI@c@)
|
|
563
|
|
#define DEG2RAD (NPY_PI@c@/180.0@c@)
|
|
564
|
|
|
|
565
|
1
|
NPY_INPLACE @type@ npy_rad2deg@c@(@type@ x)
|
|
566
|
|
{
|
|
567
|
1
|
return x*RAD2DEG;
|
|
568
|
|
}
|
|
569
|
|
|
|
570
|
1
|
NPY_INPLACE @type@ npy_deg2rad@c@(@type@ x)
|
|
571
|
|
{
|
|
572
|
1
|
return x*DEG2RAD;
|
|
573
|
|
}
|
|
574
|
|
|
|
575
|
0
|
NPY_INPLACE @type@ npy_log2_1p@c@(@type@ x)
|
|
576
|
|
{
|
|
577
|
1
|
return LOG2E*npy_log1p@c@(x);
|
|
578
|
|
}
|
|
579
|
|
|
|
580
|
0
|
NPY_INPLACE @type@ npy_exp2_m1@c@(@type@ x)
|
|
581
|
|
{
|
|
582
|
0
|
return npy_expm1@c@(LOGE2*x);
|
|
583
|
|
}
|
|
584
|
|
|
|
585
|
1
|
NPY_INPLACE @type@ npy_logaddexp@c@(@type@ x, @type@ y)
|
|
586
|
|
{
|
|
587
|
1
|
if (x == y) {
|
|
588
|
|
/* Handles infinities of the same sign without warnings */
|
|
589
|
1
|
return x + LOGE2;
|
|
590
|
|
}
|
|
591
|
|
else {
|
|
592
|
1
|
const @type@ tmp = x - y;
|
|
593
|
1
|
if (tmp > 0) {
|
|
594
|
1
|
return x + npy_log1p@c@(npy_exp@c@(-tmp));
|
|
595
|
|
}
|
|
596
|
1
|
else if (tmp <= 0) {
|
|
597
|
1
|
return y + npy_log1p@c@(npy_exp@c@(tmp));
|
|
598
|
|
}
|
|
599
|
|
else {
|
|
600
|
|
/* NaNs */
|
|
601
|
|
return tmp;
|
|
602
|
|
}
|
|
603
|
|
}
|
|
604
|
|
}
|
|
605
|
|
|
|
606
|
1
|
NPY_INPLACE @type@ npy_logaddexp2@c@(@type@ x, @type@ y)
|
|
607
|
|
{
|
|
608
|
1
|
if (x == y) {
|
|
609
|
|
/* Handles infinities of the same sign without warnings */
|
|
610
|
1
|
return x + 1;
|
|
611
|
|
}
|
|
612
|
|
else {
|
|
613
|
1
|
const @type@ tmp = x - y;
|
|
614
|
1
|
if (tmp > 0) {
|
|
615
|
1
|
return x + npy_log2_1p@c@(npy_exp2@c@(-tmp));
|
|
616
|
|
}
|
|
617
|
1
|
else if (tmp <= 0) {
|
|
618
|
1
|
return y + npy_log2_1p@c@(npy_exp2@c@(tmp));
|
|
619
|
|
}
|
|
620
|
|
else {
|
|
621
|
|
/* NaNs */
|
|
622
|
|
return tmp;
|
|
623
|
|
}
|
|
624
|
|
}
|
|
625
|
|
}
|
|
626
|
|
|
|
627
|
|
/*
|
|
628
|
|
* Python version of divmod.
|
|
629
|
|
*
|
|
630
|
|
* The implementation is mostly copied from cpython 3.5.
|
|
631
|
|
*/
|
|
632
|
|
NPY_INPLACE @type@
|
|
633
|
1
|
npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus)
|
|
634
|
|
{
|
|
635
|
|
@type@ div, mod, floordiv;
|
|
636
|
|
|
|
637
|
1
|
mod = npy_fmod@c@(a, b);
|
|
638
|
|
|
|
639
|
1
|
if (!b) {
|
|
640
|
|
/* If b == 0, return result of fmod. For IEEE is nan */
|
|
641
|
1
|
*modulus = mod;
|
|
642
|
1
|
return mod;
|
|
643
|
|
}
|
|
644
|
|
|
|
645
|
|
/* a - mod should be very nearly an integer multiple of b */
|
|
646
|
1
|
div = (a - mod) / b;
|
|
647
|
|
|
|
648
|
|
/* adjust fmod result to conform to Python convention of remainder */
|
|
649
|
1
|
if (mod) {
|
|
650
|
1
|
if ((b < 0) != (mod < 0)) {
|
|
651
|
1
|
mod += b;
|
|
652
|
1
|
div -= 1.0@c@;
|
|
653
|
|
}
|
|
654
|
|
}
|
|
655
|
|
else {
|
|
656
|
|
/* if mod is zero ensure correct sign */
|
|
657
|
1
|
mod = npy_copysign@c@(0, b);
|
|
658
|
|
}
|
|
659
|
|
|
|
660
|
|
/* snap quotient to nearest integral value */
|
|
661
|
1
|
if (div) {
|
|
662
|
1
|
floordiv = npy_floor@c@(div);
|
|
663
|
1
|
if (div - floordiv > 0.5@c@)
|
|
664
|
0
|
floordiv += 1.0@c@;
|
|
665
|
|
}
|
|
666
|
|
else {
|
|
667
|
|
/* if div is zero ensure correct sign */
|
|
668
|
1
|
floordiv = npy_copysign@c@(0, a/b);
|
|
669
|
|
}
|
|
670
|
|
|
|
671
|
1
|
*modulus = mod;
|
|
672
|
1
|
return floordiv;
|
|
673
|
|
}
|
|
674
|
|
|
|
675
|
|
#undef LOGE2
|
|
676
|
|
#undef LOG2E
|
|
677
|
|
#undef RAD2DEG
|
|
678
|
|
#undef DEG2RAD
|
|
679
|
|
|
|
680
|
|
/**end repeat**/
|
|
681
|
|
|
|
682
|
|
/**begin repeat
|
|
683
|
|
*
|
|
684
|
|
* #type = npy_uint, npy_ulong, npy_ulonglong#
|
|
685
|
|
* #c = u,ul,ull#
|
|
686
|
|
*/
|
|
687
|
|
NPY_INPLACE @type@
|
|
688
|
0
|
npy_gcd@c@(@type@ a, @type@ b)
|
|
689
|
|
{
|
|
690
|
|
@type@ c;
|
|
691
|
1
|
while (a != 0) {
|
|
692
|
1
|
c = a;
|
|
693
|
1
|
a = b%a;
|
|
694
|
1
|
b = c;
|
|
695
|
|
}
|
|
696
|
0
|
return b;
|
|
697
|
|
}
|
|
698
|
|
|
|
699
|
|
NPY_INPLACE @type@
|
|
700
|
0
|
npy_lcm@c@(@type@ a, @type@ b)
|
|
701
|
|
{
|
|
702
|
1
|
@type@ gcd = npy_gcd@c@(a, b);
|
|
703
|
1
|
return gcd == 0 ? 0 : a / gcd * b;
|
|
704
|
|
}
|
|
705
|
|
/**end repeat**/
|
|
706
|
|
|
|
707
|
|
/**begin repeat
|
|
708
|
|
*
|
|
709
|
|
* #type = (npy_int, npy_long, npy_longlong)*2#
|
|
710
|
|
* #c = (,l,ll)*2#
|
|
711
|
|
* #func=gcd*3,lcm*3#
|
|
712
|
|
*/
|
|
713
|
|
NPY_INPLACE @type@
|
|
714
|
0
|
npy_@func@@c@(@type@ a, @type@ b)
|
|
715
|
|
{
|
|
716
|
1
|
return npy_@func@u@c@(a < 0 ? -a : a, b < 0 ? -b : b);
|
|
717
|
|
}
|
|
718
|
|
/**end repeat**/
|
|
719
|
|
|
|
720
|
|
/* Unlike LCM and GCD, we need byte and short variants for the shift operators,
|
|
721
|
|
* since the result is dependent on the width of the type
|
|
722
|
|
*/
|
|
723
|
|
/**begin repeat
|
|
724
|
|
*
|
|
725
|
|
* #type = byte, short, int, long, longlong#
|
|
726
|
|
* #c = hh,h,,l,ll#
|
|
727
|
|
*/
|
|
728
|
|
/**begin repeat1
|
|
729
|
|
*
|
|
730
|
|
* #u = u,#
|
|
731
|
|
* #is_signed = 0,1#
|
|
732
|
|
*/
|
|
733
|
|
NPY_INPLACE npy_@u@@type@
|
|
734
|
0
|
npy_lshift@u@@c@(npy_@u@@type@ a, npy_@u@@type@ b)
|
|
735
|
|
{
|
|
736
|
1
|
if (NPY_LIKELY((size_t)b < sizeof(a) * CHAR_BIT)) {
|
|
737
|
1
|
return a << b;
|
|
738
|
|
}
|
|
739
|
|
else {
|
|
740
|
|
return 0;
|
|
741
|
|
}
|
|
742
|
|
}
|
|
743
|
|
NPY_INPLACE npy_@u@@type@
|
|
744
|
0
|
npy_rshift@u@@c@(npy_@u@@type@ a, npy_@u@@type@ b)
|
|
745
|
|
{
|
|
746
|
1
|
if (NPY_LIKELY((size_t)b < sizeof(a) * CHAR_BIT)) {
|
|
747
|
1
|
return a >> b;
|
|
748
|
|
}
|
|
749
|
|
#if @is_signed@
|
|
750
|
1
|
else if (a < 0) {
|
|
751
|
|
return (npy_@u@@type@)-1; /* preserve the sign bit */
|
|
752
|
|
}
|
|
753
|
|
#endif
|
|
754
|
|
else {
|
|
755
|
0
|
return 0;
|
|
756
|
|
}
|
|
757
|
|
}
|
|
758
|
|
/**end repeat1**/
|
|
759
|
|
/**end repeat**/
|
BvB93