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**/

Read our documentation on viewing source code .

Loading