1
#define NPY_NO_DEPRECATED_API NPY_API_VERSION
2
#include <Python.h>
3
#include <structmember.h>
4
#include <string.h>
5

6
#define _MULTIARRAYMODULE
7
#include "numpy/arrayobject.h"
8
#include "numpy/npy_3kcompat.h"
9
#include "numpy/npy_math.h"
10
#include "npy_config.h"
11
#include "templ_common.h" /* for npy_mul_with_overflow_intp */
12
#include "lowlevel_strided_loops.h" /* for npy_bswap8 */
13
#include "alloc.h"
14
#include "ctors.h"
15
#include "common.h"
16

17

18
/*
19
 * Returns -1 if the array is monotonic decreasing,
20
 * +1 if the array is monotonic increasing,
21
 * and 0 if the array is not monotonic.
22
 */
23
static int
24 1
check_array_monotonic(const double *a, npy_intp lena)
25
{
26
    npy_intp i;
27
    double next;
28
    double last;
29

30 1
    if (lena == 0) {
31
        /* all bin edges hold the same value */
32
        return 1;
33
    }
34 1
    last = a[0];
35

36
    /* Skip repeated values at the beginning of the array */
37 1
    for (i = 1; (i < lena) && (a[i] == last); i++);
38

39 1
    if (i == lena) {
40
        /* all bin edges hold the same value */
41
        return 1;
42
    }
43

44 1
    next = a[i];
45 1
    if (last < next) {
46
        /* Possibly monotonic increasing */
47 1
        for (i += 1; i < lena; i++) {
48 1
            last = next;
49 1
            next = a[i];
50 1
            if (last > next) {
51
                return 0;
52
            }
53
        }
54
        return 1;
55
    }
56
    else {
57
        /* last > next, possibly monotonic decreasing */
58 1
        for (i += 1; i < lena; i++) {
59 1
            last = next;
60 1
            next = a[i];
61 1
            if (last < next) {
62
                return 0;
63
            }
64
        }
65
        return -1;
66
    }
67
}
68

69
/* Find the minimum and maximum of an integer array */
70
static void
71
minmax(const npy_intp *data, npy_intp data_len, npy_intp *mn, npy_intp *mx)
72
{
73 1
    npy_intp min = *data;
74 1
    npy_intp max = *data;
75

76 1
    while (--data_len) {
77 1
        const npy_intp val = *(++data);
78 1
        if (val < min) {
79
            min = val;
80
        }
81 1
        else if (val > max) {
82 1
            max = val;
83
        }
84
    }
85

86 1
    *mn = min;
87 1
    *mx = max;
88
}
89

90
/*
91
 * arr_bincount is registered as bincount.
92
 *
93
 * bincount accepts one, two or three arguments. The first is an array of
94
 * non-negative integers The second, if present, is an array of weights,
95
 * which must be promotable to double. Call these arguments list and
96
 * weight. Both must be one-dimensional with len(weight) == len(list). If
97
 * weight is not present then bincount(list)[i] is the number of occurrences
98
 * of i in list.  If weight is present then bincount(self,list, weight)[i]
99
 * is the sum of all weight[j] where list [j] == i.  Self is not used.
100
 * The third argument, if present, is a minimum length desired for the
101
 * output array.
102
 */
103
NPY_NO_EXPORT PyObject *
104 1
arr_bincount(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
105
{
106 1
    PyObject *list = NULL, *weight = Py_None, *mlength = NULL;
107 1
    PyArrayObject *lst = NULL, *ans = NULL, *wts = NULL;
108
    npy_intp *numbers, *ians, len, mx, mn, ans_size;
109 1
    npy_intp minlength = 0;
110
    npy_intp i;
111
    double *weights , *dans;
112
    static char *kwlist[] = {"list", "weights", "minlength", NULL};
113

114 1
    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|OO:bincount",
115
                kwlist, &list, &weight, &mlength)) {
116
            goto fail;
117
    }
118

119 1
    lst = (PyArrayObject *)PyArray_ContiguousFromAny(list, NPY_INTP, 1, 1);
120 1
    if (lst == NULL) {
121
        goto fail;
122
    }
123 1
    len = PyArray_SIZE(lst);
124

125
    /*
126
     * This if/else if can be removed by changing the argspec to O|On above,
127
     * once we retire the deprecation
128
     */
129 1
    if (mlength == Py_None) {
130
        /* NumPy 1.14, 2017-06-01 */
131 1
        if (DEPRECATE("0 should be passed as minlength instead of None; "
132
                      "this will error in future.") < 0) {
133
            goto fail;
134
        }
135
    }
136 1
    else if (mlength != NULL) {
137 1
        minlength = PyArray_PyIntAsIntp(mlength);
138 1
        if (error_converting(minlength)) {
139
            goto fail;
140
        }
141
    }
142

143 1
    if (minlength < 0) {
144 1
        PyErr_SetString(PyExc_ValueError,
145
                        "'minlength' must not be negative");
146 1
        goto fail;
147
    }
148

149
    /* handle empty list */
150 1
    if (len == 0) {
151 1
        ans = (PyArrayObject *)PyArray_ZEROS(1, &minlength, NPY_INTP, 0);
152 1
        if (ans == NULL){
153
            goto fail;
154
        }
155 1
        Py_DECREF(lst);
156
        return (PyObject *)ans;
157
    }
158

159 1
    numbers = (npy_intp *)PyArray_DATA(lst);
160 1
    minmax(numbers, len, &mn, &mx);
161 1
    if (mn < 0) {
162 0
        PyErr_SetString(PyExc_ValueError,
163
                "'list' argument must have no negative elements");
164 0
        goto fail;
165
    }
166 1
    ans_size = mx + 1;
167 1
    if (mlength != Py_None) {
168 1
        if (ans_size < minlength) {
169 1
            ans_size = minlength;
170
        }
171
    }
172 1
    if (weight == Py_None) {
173 1
        ans = (PyArrayObject *)PyArray_ZEROS(1, &ans_size, NPY_INTP, 0);
174 1
        if (ans == NULL) {
175
            goto fail;
176
        }
177 1
        ians = (npy_intp *)PyArray_DATA(ans);
178 1
        NPY_BEGIN_ALLOW_THREADS;
179 1
        for (i = 0; i < len; i++)
180 1
            ians[numbers[i]] += 1;
181 1
        NPY_END_ALLOW_THREADS;
182 1
        Py_DECREF(lst);
183
    }
184
    else {
185 1
        wts = (PyArrayObject *)PyArray_ContiguousFromAny(
186
                                                weight, NPY_DOUBLE, 1, 1);
187 1
        if (wts == NULL) {
188
            goto fail;
189
        }
190 1
        weights = (double *)PyArray_DATA(wts);
191 1
        if (PyArray_SIZE(wts) != len) {
192 0
            PyErr_SetString(PyExc_ValueError,
193
                    "The weights and list don't have the same length.");
194 0
            goto fail;
195
        }
196 1
        ans = (PyArrayObject *)PyArray_ZEROS(1, &ans_size, NPY_DOUBLE, 0);
197 1
        if (ans == NULL) {
198
            goto fail;
199
        }
200 1
        dans = (double *)PyArray_DATA(ans);
201 1
        NPY_BEGIN_ALLOW_THREADS;
202 1
        for (i = 0; i < len; i++) {
203 1
            dans[numbers[i]] += weights[i];
204
        }
205 1
        NPY_END_ALLOW_THREADS;
206 1
        Py_DECREF(lst);
207 1
        Py_DECREF(wts);
208
    }
209
    return (PyObject *)ans;
210

211 1
fail:
212 1
    Py_XDECREF(lst);
213 1
    Py_XDECREF(wts);
214 1
    Py_XDECREF(ans);
215
    return NULL;
216
}
217

218
/* Internal function to expose check_array_monotonic to python */
219
NPY_NO_EXPORT PyObject *
220 1
arr__monotonicity(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
221
{
222
    static char *kwlist[] = {"x", NULL};
223 1
    PyObject *obj_x = NULL;
224 1
    PyArrayObject *arr_x = NULL;
225
    long monotonic;
226
    npy_intp len_x;
227 1
    NPY_BEGIN_THREADS_DEF;
228

229 1
    if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|_monotonicity", kwlist,
230
                                     &obj_x)) {
231
        return NULL;
232
    }
233

234
    /*
235
     * TODO:
236
     *  `x` could be strided, needs change to check_array_monotonic
237
     *  `x` is forced to double for this check
238
     */
239 1
    arr_x = (PyArrayObject *)PyArray_FROMANY(
240
        obj_x, NPY_DOUBLE, 1, 1, NPY_ARRAY_CARRAY_RO);
241 1
    if (arr_x == NULL) {
242
        return NULL;
243
    }
244

245 1
    len_x = PyArray_SIZE(arr_x);
246 1
    NPY_BEGIN_THREADS_THRESHOLDED(len_x)
247 1
    monotonic = check_array_monotonic(
248 1
        (const double *)PyArray_DATA(arr_x), len_x);
249 1
    NPY_END_THREADS
250 1
    Py_DECREF(arr_x);
251

252 1
    return PyLong_FromLong(monotonic);
253
}
254

255
/*
256
 * Returns input array with values inserted sequentially into places
257
 * indicated by the mask
258
 */
259
NPY_NO_EXPORT PyObject *
260 1
arr_insert(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
261
{
262
    char *src, *dest;
263
    npy_bool *mask_data;
264
    PyArray_Descr *dtype;
265
    PyArray_CopySwapFunc *copyswap;
266
    PyObject *array0, *mask0, *values0;
267
    PyArrayObject *array, *mask, *values;
268
    npy_intp i, j, chunk, nm, ni, nv;
269

270
    static char *kwlist[] = {"input", "mask", "vals", NULL};
271 1
    NPY_BEGIN_THREADS_DEF;
272 1
    values = mask = NULL;
273

274 1
    if (!PyArg_ParseTupleAndKeywords(args, kwdict, "O!OO:place", kwlist,
275
                &PyArray_Type, &array0, &mask0, &values0)) {
276
        return NULL;
277
    }
278

279 1
    array = (PyArrayObject *)PyArray_FromArray((PyArrayObject *)array0, NULL,
280
                                    NPY_ARRAY_CARRAY | NPY_ARRAY_WRITEBACKIFCOPY);
281 1
    if (array == NULL) {
282
        goto fail;
283
    }
284

285 1
    ni = PyArray_SIZE(array);
286 1
    dest = PyArray_DATA(array);
287 1
    chunk = PyArray_DESCR(array)->elsize;
288 1
    mask = (PyArrayObject *)PyArray_FROM_OTF(mask0, NPY_BOOL,
289
                                NPY_ARRAY_CARRAY | NPY_ARRAY_FORCECAST);
290 1
    if (mask == NULL) {
291
        goto fail;
292
    }
293

294 1
    nm = PyArray_SIZE(mask);
295 1
    if (nm != ni) {
296 0
        PyErr_SetString(PyExc_ValueError,
297
                        "place: mask and data must be "
298
                        "the same size");
299 0
        goto fail;
300
    }
301

302 1
    mask_data = PyArray_DATA(mask);
303 1
    dtype = PyArray_DESCR(array);
304 1
    Py_INCREF(dtype);
305

306 1
    values = (PyArrayObject *)PyArray_FromAny(values0, dtype,
307
                                    0, 0, NPY_ARRAY_CARRAY, NULL);
308 1
    if (values == NULL) {
309
        goto fail;
310
    }
311

312 1
    nv = PyArray_SIZE(values); /* zero if null array */
313 1
    if (nv <= 0) {
314
        npy_bool allFalse = 1;
315
        i = 0;
316

317 1
        while (allFalse && i < ni) {
318 1
            if (mask_data[i]) {
319
                allFalse = 0;
320
            } else {
321 1
                i++;
322
            }
323
        }
324 1
        if (!allFalse) {
325 1
            PyErr_SetString(PyExc_ValueError,
326
                            "Cannot insert from an empty array!");
327 1
            goto fail;
328
        } else {
329 1
            Py_XDECREF(values);
330 1
            Py_XDECREF(mask);
331 1
            PyArray_ResolveWritebackIfCopy(array);
332 1
            Py_XDECREF(array);
333 1
            Py_RETURN_NONE;
334
        }
335
    }
336

337 1
    src = PyArray_DATA(values);
338 1
    j = 0;
339

340 1
    copyswap = PyArray_DESCR(array)->f->copyswap;
341 1
    NPY_BEGIN_THREADS_DESCR(PyArray_DESCR(array));
342 1
    for (i = 0; i < ni; i++) {
343 1
        if (mask_data[i]) {
344 1
            if (j >= nv) {
345 1
                j = 0;
346
            }
347

348 1
            copyswap(dest + i*chunk, src + j*chunk, 0, array);
349 1
            j++;
350
        }
351
    }
352 1
    NPY_END_THREADS;
353

354 1
    Py_XDECREF(values);
355 1
    Py_XDECREF(mask);
356 1
    PyArray_ResolveWritebackIfCopy(array);
357 1
    Py_DECREF(array);
358 1
    Py_RETURN_NONE;
359

360 0
 fail:
361 1
    Py_XDECREF(mask);
362 1
    PyArray_ResolveWritebackIfCopy(array);
363 1
    Py_XDECREF(array);
364 1
    Py_XDECREF(values);
365
    return NULL;
366
}
367

368
#define LIKELY_IN_CACHE_SIZE 8
369

370
#ifdef __INTEL_COMPILER
371
#pragma intel optimization_level 0
372
#endif
373
static NPY_INLINE npy_intp
374
_linear_search(const npy_double key, const npy_double *arr, const npy_intp len, const npy_intp i0)
375
{
376
    npy_intp i;
377

378 1
    for (i = i0; i < len && key >= arr[i]; i++);
379 1
    return i - 1;
380
}
381

382
/** @brief find index of a sorted array such that arr[i] <= key < arr[i + 1].
383
 *
384
 * If an starting index guess is in-range, the array values around this
385
 * index are first checked.  This allows for repeated calls for well-ordered
386
 * keys (a very common case) to use the previous index as a very good guess.
387
 *
388
 * If the guess value is not useful, bisection of the array is used to
389
 * find the index.  If there is no such index, the return values are:
390
 *     key < arr[0] -- -1
391
 *     key == arr[len - 1] -- len - 1
392
 *     key > arr[len - 1] -- len
393
 * The array is assumed contiguous and sorted in ascending order.
394
 *
395
 * @param key key value.
396
 * @param arr contiguous sorted array to be searched.
397
 * @param len length of the array.
398
 * @param guess initial guess of index
399
 * @return index
400
 */
401
static npy_intp
402 1
binary_search_with_guess(const npy_double key, const npy_double *arr,
403
                         npy_intp len, npy_intp guess)
404
{
405 1
    npy_intp imin = 0;
406 1
    npy_intp imax = len;
407

408
    /* Handle keys outside of the arr range first */
409 1
    if (key > arr[len - 1]) {
410
        return len;
411
    }
412 1
    else if (key < arr[0]) {
413
        return -1;
414
    }
415

416
    /*
417
     * If len <= 4 use linear search.
418
     * From above we know key >= arr[0] when we start.
419
     */
420 1
    if (len <= 4) {
421 1
        return _linear_search(key, arr, len, 1);
422
    }
423

424 1
    if (guess > len - 3) {
425 1
        guess = len - 3;
426
    }
427 1
    if (guess < 1)  {
428 1
        guess = 1;
429
    }
430

431
    /* check most likely values: guess - 1, guess, guess + 1 */
432 1
    if (key < arr[guess]) {
433 1
        if (key < arr[guess - 1]) {
434 1
            imax = guess - 1;
435
            /* last attempt to restrict search to items in cache */
436 1
            if (guess > LIKELY_IN_CACHE_SIZE &&
437 0
                        key >= arr[guess - LIKELY_IN_CACHE_SIZE]) {
438 0
                imin = guess - LIKELY_IN_CACHE_SIZE;
439
            }
440
        }
441
        else {
442
            /* key >= arr[guess - 1] */
443 1
            return guess - 1;
444
        }
445
    }
446
    else {
447
        /* key >= arr[guess] */
448 1
        if (key < arr[guess + 1]) {
449
            return guess;
450
        }
451
        else {
452
            /* key >= arr[guess + 1] */
453 1
            if (key < arr[guess + 2]) {
454 1
                return guess + 1;
455
            }
456
            else {
457
                /* key >= arr[guess + 2] */
458 1
                imin = guess + 2;
459
                /* last attempt to restrict search to items in cache */
460 1
                if (guess < len - LIKELY_IN_CACHE_SIZE - 1 &&
461 1
                            key < arr[guess + LIKELY_IN_CACHE_SIZE]) {
462 0
                    imax = guess + LIKELY_IN_CACHE_SIZE;
463
                }
464
            }
465
        }
466
    }
467

468
    /* finally, find index by bisection */
469 1
    while (imin < imax) {
470 1
        const npy_intp imid = imin + ((imax - imin) >> 1);
471 1
        if (key >= arr[imid]) {
472 1
            imin = imid + 1;
473
        }
474
        else {
475
            imax = imid;
476
        }
477
    }
478 1
    return imin - 1;
479
}
480

481
#undef LIKELY_IN_CACHE_SIZE
482

483
NPY_NO_EXPORT PyObject *
484 1
arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
485
{
486

487
    PyObject *fp, *xp, *x;
488 1
    PyObject *left = NULL, *right = NULL;
489 1
    PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL;
490
    npy_intp i, lenx, lenxp;
491
    npy_double lval, rval;
492
    const npy_double *dy, *dx, *dz;
493 1
    npy_double *dres, *slopes = NULL;
494

495
    static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL};
496

497 1
    NPY_BEGIN_THREADS_DEF;
498

499 1
    if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO:interp", kwlist,
500
                                     &x, &xp, &fp, &left, &right)) {
501
        return NULL;
502
    }
503

504 1
    afp = (PyArrayObject *)PyArray_ContiguousFromAny(fp, NPY_DOUBLE, 1, 1);
505 1
    if (afp == NULL) {
506
        return NULL;
507
    }
508 1
    axp = (PyArrayObject *)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1);
509 1
    if (axp == NULL) {
510
        goto fail;
511
    }
512 1
    ax = (PyArrayObject *)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 0, 0);
513 1
    if (ax == NULL) {
514
        goto fail;
515
    }
516 1
    lenxp = PyArray_SIZE(axp);
517 1
    if (lenxp == 0) {
518 1
        PyErr_SetString(PyExc_ValueError,
519
                "array of sample points is empty");
520 1
        goto fail;
521
    }
522 1
    if (PyArray_SIZE(afp) != lenxp) {
523 1
        PyErr_SetString(PyExc_ValueError,
524
                "fp and xp are not of the same length.");
525 1
        goto fail;
526
    }
527

528 1
    af = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ax),
529
                                            PyArray_DIMS(ax), NPY_DOUBLE);
530 1
    if (af == NULL) {
531
        goto fail;
532
    }
533 1
    lenx = PyArray_SIZE(ax);
534

535 1
    dy = (const npy_double *)PyArray_DATA(afp);
536 1
    dx = (const npy_double *)PyArray_DATA(axp);
537 1
    dz = (const npy_double *)PyArray_DATA(ax);
538 1
    dres = (npy_double *)PyArray_DATA(af);
539
    /* Get left and right fill values. */
540 1
    if ((left == NULL) || (left == Py_None)) {
541 1
        lval = dy[0];
542
    }
543
    else {
544 1
        lval = PyFloat_AsDouble(left);
545 1
        if (error_converting(lval)) {
546
            goto fail;
547
        }
548
    }
549 1
    if ((right == NULL) || (right == Py_None)) {
550 1
        rval = dy[lenxp - 1];
551
    }
552
    else {
553 1
        rval = PyFloat_AsDouble(right);
554 1
        if (error_converting(rval)) {
555
            goto fail;
556
        }
557
    }
558

559
    /* binary_search_with_guess needs at least a 3 item long array */
560 1
    if (lenxp == 1) {
561 1
        const npy_double xp_val = dx[0];
562 1
        const npy_double fp_val = dy[0];
563

564 1
        NPY_BEGIN_THREADS_THRESHOLDED(lenx);
565 1
        for (i = 0; i < lenx; ++i) {
566 1
            const npy_double x_val = dz[i];
567 1
            dres[i] = (x_val < xp_val) ? lval :
568 1
                                         ((x_val > xp_val) ? rval : fp_val);
569
        }
570 1
        NPY_END_THREADS;
571
    }
572
    else {
573 1
        npy_intp j = 0;
574

575
        /* only pre-calculate slopes if there are relatively few of them. */
576 1
        if (lenxp <= lenx) {
577 1
            slopes = PyArray_malloc((lenxp - 1) * sizeof(npy_double));
578 1
            if (slopes == NULL) {
579 0
                PyErr_NoMemory();
580 0
                goto fail;
581
            }
582
        }
583

584 1
        NPY_BEGIN_THREADS;
585

586 1
        if (slopes != NULL) {
587 1
            for (i = 0; i < lenxp - 1; ++i) {
588 1
                slopes[i] = (dy[i+1] - dy[i]) / (dx[i+1] - dx[i]);
589
            }
590
        }
591

592 1
        for (i = 0; i < lenx; ++i) {
593 1
            const npy_double x_val = dz[i];
594

595 1
            if (npy_isnan(x_val)) {
596 1
                dres[i] = x_val;
597 1
                continue;
598
            }
599

600 1
            j = binary_search_with_guess(x_val, dx, lenxp, j);
601 1
            if (j == -1) {
602 1
                dres[i] = lval;
603
            }
604 1
            else if (j == lenxp) {
605 1
                dres[i] = rval;
606
            }
607 1
            else if (j == lenxp - 1) {
608 1
                dres[i] = dy[j];
609
            }
610 1
            else if (dx[j] == x_val) {
611
                /* Avoid potential non-finite interpolation */
612 1
                dres[i] = dy[j];
613
            }
614
            else {
615 1
                const npy_double slope =
616 1
                        (slopes != NULL) ? slopes[j] :
617 1
                        (dy[j+1] - dy[j]) / (dx[j+1] - dx[j]);
618

619
                /* If we get nan in one direction, try the other */
620 1
                dres[i] = slope*(x_val - dx[j]) + dy[j];
621 1
                if (NPY_UNLIKELY(npy_isnan(dres[i]))) {
622 1
                    dres[i] = slope*(x_val - dx[j+1]) + dy[j+1];
623 1
                    if (NPY_UNLIKELY(npy_isnan(dres[i])) && dy[j] == dy[j+1]) {
624 1
                        dres[i] = dy[j];
625
                    }
626
                }
627
            }
628
        }
629

630 1
        NPY_END_THREADS;
631
    }
632

633 1
    PyArray_free(slopes);
634 1
    Py_DECREF(afp);
635 1
    Py_DECREF(axp);
636 1
    Py_DECREF(ax);
637 1
    return PyArray_Return(af);
638

639 0
fail:
640 1
    Py_XDECREF(afp);
641 1
    Py_XDECREF(axp);
642 1
    Py_XDECREF(ax);
643 1
    Py_XDECREF(af);
644
    return NULL;
645
}
646

647
/* As for arr_interp but for complex fp values */
648
NPY_NO_EXPORT PyObject *
649 1
arr_interp_complex(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
650
{
651

652
    PyObject *fp, *xp, *x;
653 1
    PyObject *left = NULL, *right = NULL;
654 1
    PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL;
655
    npy_intp i, lenx, lenxp;
656

657
    const npy_double *dx, *dz;
658
    const npy_cdouble *dy;
659
    npy_cdouble lval, rval;
660 1
    npy_cdouble *dres, *slopes = NULL;
661

662
    static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL};
663

664 1
    NPY_BEGIN_THREADS_DEF;
665

666 1
    if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO:interp_complex",
667
                                     kwlist, &x, &xp, &fp, &left, &right)) {
668
        return NULL;
669
    }
670

671 1
    afp = (PyArrayObject *)PyArray_ContiguousFromAny(fp, NPY_CDOUBLE, 1, 1);
672

673 1
    if (afp == NULL) {
674
        return NULL;
675
    }
676

677 1
    axp = (PyArrayObject *)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1);
678 1
    if (axp == NULL) {
679
        goto fail;
680
    }
681 1
    ax = (PyArrayObject *)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 0, 0);
682 1
    if (ax == NULL) {
683
        goto fail;
684
    }
685 1
    lenxp = PyArray_SIZE(axp);
686 1
    if (lenxp == 0) {
687 0
        PyErr_SetString(PyExc_ValueError,
688
                "array of sample points is empty");
689 0
        goto fail;
690
    }
691 1
    if (PyArray_SIZE(afp) != lenxp) {
692 0
        PyErr_SetString(PyExc_ValueError,
693
                "fp and xp are not of the same length.");
694 0
        goto fail;
695
    }
696

697 1
    lenx = PyArray_SIZE(ax);
698 1
    dx = (const npy_double *)PyArray_DATA(axp);
699 1
    dz = (const npy_double *)PyArray_DATA(ax);
700

701 1
    af = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ax),
702
                                            PyArray_DIMS(ax), NPY_CDOUBLE);
703 1
    if (af == NULL) {
704
        goto fail;
705
    }
706

707 1
    dy = (const npy_cdouble *)PyArray_DATA(afp);
708 1
    dres = (npy_cdouble *)PyArray_DATA(af);
709
    /* Get left and right fill values. */
710 1
    if ((left == NULL) || (left == Py_None)) {
711 1
        lval = dy[0];
712
    }
713
    else {
714 1
        lval.real = PyComplex_RealAsDouble(left);
715 1
        if (error_converting(lval.real)) {
716
            goto fail;
717
        }
718 1
        lval.imag = PyComplex_ImagAsDouble(left);
719 1
        if (error_converting(lval.imag)) {
720
            goto fail;
721
        }
722
    }
723

724 1
    if ((right == NULL) || (right == Py_None)) {
725 1
        rval = dy[lenxp - 1];
726
    }
727
    else {
728 1
        rval.real = PyComplex_RealAsDouble(right);
729 1
        if (error_converting(rval.real)) {
730
            goto fail;
731
        }
732 1
        rval.imag = PyComplex_ImagAsDouble(right);
733 1
        if (error_converting(rval.imag)) {
734
            goto fail;
735
        }
736
    }
737

738
    /* binary_search_with_guess needs at least a 3 item long array */
739 1
    if (lenxp == 1) {
740 0
        const npy_double xp_val = dx[0];
741 0
        const npy_cdouble fp_val = dy[0];
742

743 0
        NPY_BEGIN_THREADS_THRESHOLDED(lenx);
744 0
        for (i = 0; i < lenx; ++i) {
745 0
            const npy_double x_val = dz[i];
746 0
            dres[i] = (x_val < xp_val) ? lval :
747
              ((x_val > xp_val) ? rval : fp_val);
748
        }
749 0
        NPY_END_THREADS;
750
    }
751
    else {
752 1
        npy_intp j = 0;
753

754
        /* only pre-calculate slopes if there are relatively few of them. */
755 1
        if (lenxp <= lenx) {
756 1
            slopes = PyArray_malloc((lenxp - 1) * sizeof(npy_cdouble));
757 1
            if (slopes == NULL) {
758 0
                PyErr_NoMemory();
759 0
                goto fail;
760
            }
761
        }
762

763 1
        NPY_BEGIN_THREADS;
764

765 1
        if (slopes != NULL) {
766 1
            for (i = 0; i < lenxp - 1; ++i) {
767 1
                const double inv_dx = 1.0 / (dx[i+1] - dx[i]);
768 1
                slopes[i].real = (dy[i+1].real - dy[i].real) * inv_dx;
769 1
                slopes[i].imag = (dy[i+1].imag - dy[i].imag) * inv_dx;
770
            }
771
        }
772

773 1
        for (i = 0; i < lenx; ++i) {
774 1
            const npy_double x_val = dz[i];
775

776 1
            if (npy_isnan(x_val)) {
777 0
                dres[i].real = x_val;
778 0
                dres[i].imag = 0.0;
779 0
                continue;
780
            }
781

782 1
            j = binary_search_with_guess(x_val, dx, lenxp, j);
783 1
            if (j == -1) {
784 1
                dres[i] = lval;
785
            }
786 1
            else if (j == lenxp) {
787 1
                dres[i] = rval;
788
            }
789 1
            else if (j == lenxp - 1) {
790 1
                dres[i] = dy[j];
791
            }
792 1
            else if (dx[j] == x_val) {
793
                /* Avoid potential non-finite interpolation */
794 1
                dres[i] = dy[j];
795
            }
796
            else {
797
                npy_cdouble slope;
798 1
                if (slopes != NULL) {
799 1
                    slope = slopes[j];
800
                }
801
                else {
802 1
                    const npy_double inv_dx = 1.0 / (dx[j+1] - dx[j]);
803 1
                    slope.real = (dy[j+1].real - dy[j].real) * inv_dx;
804 1
                    slope.imag = (dy[j+1].imag - dy[j].imag) * inv_dx;
805
                }
806

807
                /* If we get nan in one direction, try the other */
808 1
                dres[i].real = slope.real*(x_val - dx[j]) + dy[j].real;
809 1
                if (NPY_UNLIKELY(npy_isnan(dres[i].real))) {
810 1
                    dres[i].real = slope.real*(x_val - dx[j+1]) + dy[j+1].real;
811 1
                    if (NPY_UNLIKELY(npy_isnan(dres[i].real)) &&
812 1
                            dy[j].real == dy[j+1].real) {
813 1
                        dres[i].real = dy[j].real;
814
                    }
815
                }
816 1
                dres[i].imag = slope.imag*(x_val - dx[j]) + dy[j].imag;
817 1
                if (NPY_UNLIKELY(npy_isnan(dres[i].imag))) {
818 1
                    dres[i].imag = slope.imag*(x_val - dx[j+1]) + dy[j+1].imag;
819 1
                    if (NPY_UNLIKELY(npy_isnan(dres[i].imag)) &&
820 1
                            dy[j].imag == dy[j+1].imag) {
821 1
                        dres[i].imag = dy[j].imag;
822
                    }
823
                }
824
            }
825
        }
826

827 1
        NPY_END_THREADS;
828
    }
829 1
    PyArray_free(slopes);
830

831 1
    Py_DECREF(afp);
832 1
    Py_DECREF(axp);
833 1
    Py_DECREF(ax);
834 1
    return PyArray_Return(af);
835

836 0
fail:
837 0
    Py_XDECREF(afp);
838 0
    Py_XDECREF(axp);
839 0
    Py_XDECREF(ax);
840 0
    Py_XDECREF(af);
841
    return NULL;
842
}
843

844
static const char *EMPTY_SEQUENCE_ERR_MSG = "indices must be integral: the provided " \
845
    "empty sequence was inferred as float. Wrap it with " \
846
    "'np.array(indices, dtype=np.intp)'";
847

848
static const char *NON_INTEGRAL_ERROR_MSG = "only int indices permitted";
849

850
/* Convert obj to an ndarray with integer dtype or fail */
851
static PyArrayObject *
852 1
astype_anyint(PyObject *obj) {
853
    PyArrayObject *ret;
854

855 1
    if (!PyArray_Check(obj)) {
856
        /* prefer int dtype */
857 1
        PyArray_Descr *dtype_guess = NULL;
858 1
        if (PyArray_DTypeFromObject(obj, NPY_MAXDIMS, &dtype_guess) < 0) {
859 1
            return NULL;
860
        }
861 1
        if (dtype_guess == NULL) {
862 1
            if (PySequence_Check(obj) && PySequence_Size(obj) == 0) {
863 1
                PyErr_SetString(PyExc_TypeError, EMPTY_SEQUENCE_ERR_MSG);
864
            }
865
            return NULL;
866
        }
867 1
        ret = (PyArrayObject*)PyArray_FromAny(obj, dtype_guess, 0, 0, 0, NULL);
868 1
        if (ret == NULL) {
869
            return NULL;
870
        }
871
    }
872
    else {
873 1
        ret = (PyArrayObject *)obj;
874 1
        Py_INCREF(ret);
875
    }
876

877 1
    if (!(PyArray_ISINTEGER(ret) || PyArray_ISBOOL(ret))) {
878
        /* ensure dtype is int-based */
879 1
        PyErr_SetString(PyExc_TypeError, NON_INTEGRAL_ERROR_MSG);
880 1
        Py_DECREF(ret);
881
        return NULL;
882
    }
883

884
    return ret;
885
}
886

887
/*
888
 * Converts a Python sequence into 'count' PyArrayObjects
889
 *
890
 * seq         - Input Python object, usually a tuple but any sequence works.
891
 *               Must have integral content.
892
 * paramname   - The name of the parameter that produced 'seq'.
893
 * count       - How many arrays there should be (errors if it doesn't match).
894
 * op          - Where the arrays are placed.
895
 */
896 1
static int int_sequence_to_arrays(PyObject *seq,
897
                              char *paramname,
898
                              int count,
899
                              PyArrayObject **op
900
                              )
901
{
902
    int i;
903

904 1
    if (!PySequence_Check(seq) || PySequence_Size(seq) != count) {
905 0
        PyErr_Format(PyExc_ValueError,
906
                "parameter %s must be a sequence of length %d",
907
                paramname, count);
908 0
        return -1;
909
    }
910

911 1
    for (i = 0; i < count; ++i) {
912 1
        PyObject *item = PySequence_GetItem(seq, i);
913 1
        if (item == NULL) {
914
            goto fail;
915
        }
916 1
        op[i] = astype_anyint(item);
917 1
        Py_DECREF(item);
918 1
        if (op[i] == NULL) {
919
            goto fail;
920
        }
921
    }
922

923
    return 0;
924

925 1
fail:
926 1
    while (--i >= 0) {
927 0
        Py_XDECREF(op[i]);
928 0
        op[i] = NULL;
929
    }
930
    return -1;
931
}
932

933
/* Inner loop for ravel_multi_index */
934
static int
935 1
ravel_multi_index_loop(int ravel_ndim, npy_intp *ravel_dims,
936
                        npy_intp *ravel_strides,
937
                        npy_intp count,
938
                        NPY_CLIPMODE *modes,
939
                        char **coords, npy_intp *coords_strides)
940
{
941
    int i;
942
    char invalid;
943
    npy_intp j, m;
944

945
    /*
946
     * Check for 0-dimensional axes unless there is nothing to do.
947
     * An empty array/shape cannot be indexed at all.
948
     */
949 1
    if (count != 0) {
950 1
        for (i = 0; i < ravel_ndim; ++i) {
951 1
            if (ravel_dims[i] == 0) {
952 1
                PyErr_SetString(PyExc_ValueError,
953
                        "cannot unravel if shape has zero entries (is empty).");
954 1
                return NPY_FAIL;
955
            }
956
        }
957
    }
958

959 1
    NPY_BEGIN_ALLOW_THREADS;
960 1
    invalid = 0;
961 1
    while (count--) {
962
        npy_intp raveled = 0;
963 1
        for (i = 0; i < ravel_ndim; ++i) {
964 1
            m = ravel_dims[i];
965 1
            j = *(npy_intp *)coords[i];
966 1
            switch (modes[i]) {
967 1
                case NPY_RAISE:
968 1
                    if (j < 0 || j >= m) {
969
                        invalid = 1;
970
                        goto end_while;
971
                    }
972
                    break;
973 1
                case NPY_WRAP:
974 1
                    if (j < 0) {
975 1
                        j += m;
976 1
                        if (j < 0) {
977 0
                            j = j % m;
978 0
                            if (j != 0) {
979 0
                                j += m;
980
                            }
981
                        }
982
                    }
983 1
                    else if (j >= m) {
984 1
                        j -= m;
985 1
                        if (j >= m) {
986 0
                            j = j % m;
987
                        }
988
                    }
989
                    break;
990 1
                case NPY_CLIP:
991 1
                    if (j < 0) {
992
                        j = 0;
993
                    }
994 1
                    else if (j >= m) {
995 1
                        j = m - 1;
996
                    }
997
                    break;
998

999
            }
1000 1
            raveled += j * ravel_strides[i];
1001

1002 1
            coords[i] += coords_strides[i];
1003
        }
1004 1
        *(npy_intp *)coords[ravel_ndim] = raveled;
1005 1
        coords[ravel_ndim] += coords_strides[ravel_ndim];
1006
    }
1007 1
end_while:
1008 1
    NPY_END_ALLOW_THREADS;
1009 1
    if (invalid) {
1010 1
        PyErr_SetString(PyExc_ValueError,
1011
              "invalid entry in coordinates array");
1012 1
        return NPY_FAIL;
1013
    }
1014
    return NPY_SUCCEED;
1015
}
1016

1017
/* ravel_multi_index implementation - see add_newdocs.py */
1018
NPY_NO_EXPORT PyObject *
1019 1
arr_ravel_multi_index(PyObject *self, PyObject *args, PyObject *kwds)
1020
{
1021
    int i;
1022 1
    PyObject *mode0=NULL, *coords0=NULL;
1023 1
    PyArrayObject *ret = NULL;
1024 1
    PyArray_Dims dimensions={0,0};
1025
    npy_intp s, ravel_strides[NPY_MAXDIMS];
1026 1
    NPY_ORDER order = NPY_CORDER;
1027
    NPY_CLIPMODE modes[NPY_MAXDIMS];
1028

1029
    PyArrayObject *op[NPY_MAXARGS];
1030
    PyArray_Descr *dtype[NPY_MAXARGS];
1031
    npy_uint32 op_flags[NPY_MAXARGS];
1032

1033 1
    NpyIter *iter = NULL;
1034

1035 1
    char *kwlist[] = {"multi_index", "dims", "mode", "order", NULL};
1036

1037 1
    memset(op, 0, sizeof(op));
1038 1
    dtype[0] = NULL;
1039

1040 1
    if (!PyArg_ParseTupleAndKeywords(args, kwds,
1041
                        "OO&|OO&:ravel_multi_index", kwlist,
1042
                     &coords0,
1043
                     PyArray_IntpConverter, &dimensions,
1044
                     &mode0,
1045
                     PyArray_OrderConverter, &order)) {
1046
        goto fail;
1047
    }
1048

1049 1
    if (dimensions.len+1 > NPY_MAXARGS) {
1050 0
        PyErr_SetString(PyExc_ValueError,
1051
                    "too many dimensions passed to ravel_multi_index");
1052 0
        goto fail;
1053
    }
1054

1055 1
    if (!PyArray_ConvertClipmodeSequence(mode0, modes, dimensions.len)) {
1056
       goto fail;
1057
    }
1058

1059 1
    switch (order) {
1060 1
        case NPY_CORDER:
1061 1
            s = 1;
1062 1
            for (i = dimensions.len-1; i >= 0; --i) {
1063 1
                ravel_strides[i] = s;
1064 1
                if (npy_mul_with_overflow_intp(&s, s, dimensions.ptr[i])) {
1065 1
                    PyErr_SetString(PyExc_ValueError,
1066
                        "invalid dims: array size defined by dims is larger "
1067
                        "than the maximum possible size.");
1068 1
                    goto fail;
1069
                }
1070
            }
1071
            break;
1072
        case NPY_FORTRANORDER:
1073
            s = 1;
1074 1
            for (i = 0; i < dimensions.len; ++i) {
1075 1
                ravel_strides[i] = s;
1076 1
                if (npy_mul_with_overflow_intp(&s, s, dimensions.ptr[i])) {
1077 1
                    PyErr_SetString(PyExc_ValueError,
1078
                        "invalid dims: array size defined by dims is larger "
1079
                        "than the maximum possible size.");
1080 1
                    goto fail;
1081
                }
1082
            }
1083
            break;
1084 0
        default:
1085 0
            PyErr_SetString(PyExc_ValueError,
1086
                            "only 'C' or 'F' order is permitted");
1087 0
            goto fail;
1088
    }
1089

1090
    /* Get the multi_index into op */
1091 1
    if (int_sequence_to_arrays(coords0, "multi_index", dimensions.len, op) < 0) {
1092
        goto fail;
1093
    }
1094

1095 1
    for (i = 0; i < dimensions.len; ++i) {
1096 1
        op_flags[i] = NPY_ITER_READONLY|
1097
                      NPY_ITER_ALIGNED;
1098
    }
1099 1
    op_flags[dimensions.len] = NPY_ITER_WRITEONLY|
1100
                               NPY_ITER_ALIGNED|
1101
                               NPY_ITER_ALLOCATE;
1102 1
    dtype[0] = PyArray_DescrFromType(NPY_INTP);
1103 1
    for (i = 1; i <= dimensions.len; ++i) {
1104 1
        dtype[i] = dtype[0];
1105
    }
1106

1107 1
    iter = NpyIter_MultiNew(dimensions.len+1, op, NPY_ITER_BUFFERED|
1108
                                                  NPY_ITER_EXTERNAL_LOOP|
1109
                                                  NPY_ITER_ZEROSIZE_OK,
1110
                                                  NPY_KEEPORDER,
1111
                                                  NPY_SAME_KIND_CASTING,
1112
                                                  op_flags, dtype);
1113 1
    if (iter == NULL) {
1114
        goto fail;
1115
    }
1116

1117 1
    if (NpyIter_GetIterSize(iter) != 0) {
1118
        NpyIter_IterNextFunc *iternext;
1119
        char **dataptr;
1120
        npy_intp *strides;
1121
        npy_intp *countptr;
1122

1123 1
        iternext = NpyIter_GetIterNext(iter, NULL);
1124 1
        if (iternext == NULL) {
1125
            goto fail;
1126
        }
1127 1
        dataptr = NpyIter_GetDataPtrArray(iter);
1128 1
        strides = NpyIter_GetInnerStrideArray(iter);
1129 1
        countptr = NpyIter_GetInnerLoopSizePtr(iter);
1130

1131
        do {
1132 1
            if (ravel_multi_index_loop(dimensions.len, dimensions.ptr,
1133
                        ravel_strides, *countptr, modes,
1134
                        dataptr, strides) != NPY_SUCCEED) {
1135
                goto fail;
1136
            }
1137 1
        } while(iternext(iter));
1138
    }
1139

1140 1
    ret = NpyIter_GetOperandArray(iter)[dimensions.len];
1141 1
    Py_INCREF(ret);
1142

1143 1
    Py_DECREF(dtype[0]);
1144 1
    for (i = 0; i < dimensions.len; ++i) {
1145 1
        Py_XDECREF(op[i]);
1146
    }
1147 1
    npy_free_cache_dim_obj(dimensions);
1148 1
    NpyIter_Deallocate(iter);
1149 1
    return PyArray_Return(ret);
1150

1151 1
fail:
1152 1
    Py_XDECREF(dtype[0]);
1153 1
    for (i = 0; i < dimensions.len; ++i) {
1154 1
        Py_XDECREF(op[i]);
1155
    }
1156 1
    npy_free_cache_dim_obj(dimensions);
1157 1
    NpyIter_Deallocate(iter);
1158 1
    return NULL;
1159
}
1160

1161

1162
/*
1163
 * Inner loop for unravel_index
1164
 * order must be NPY_CORDER or NPY_FORTRANORDER
1165
 */
1166
static int
1167 1
unravel_index_loop(int unravel_ndim, npy_intp const *unravel_dims,
1168
                   npy_intp unravel_size, npy_intp count,
1169
                   char *indices, npy_intp indices_stride,
1170
                   npy_intp *coords, NPY_ORDER order)
1171
{
1172
    int i, idx;
1173 1
    int idx_start = (order == NPY_CORDER) ? unravel_ndim - 1: 0;
1174 1
    int idx_step = (order == NPY_CORDER) ? -1 : 1;
1175 1
    char invalid = 0;
1176 1
    npy_intp val = 0;
1177

1178 1
    NPY_BEGIN_ALLOW_THREADS;
1179
    /* NPY_KEEPORDER or NPY_ANYORDER have no meaning in this setting */
1180
    assert(order == NPY_CORDER || order == NPY_FORTRANORDER);
1181 1
    while (count--) {
1182 1
        val = *(npy_intp *)indices;
1183 1
        if (val < 0 || val >= unravel_size) {
1184
            invalid = 1;
1185
            break;
1186
        }
1187
        idx = idx_start;
1188 1
        for (i = 0; i < unravel_ndim; ++i) {
1189
            /*
1190
             * Using a local seems to enable single-divide optimization
1191
             * but only if the / precedes the %
1192
             */
1193 1
            npy_intp tmp = val / unravel_dims[idx];
1194 1
            coords[idx] = val % unravel_dims[idx];
1195 1
            val = tmp;
1196 1
            idx += idx_step;
1197
        }
1198 1
        coords += unravel_ndim;
1199 1
        indices += indices_stride;
1200
    }
1201 1
    NPY_END_ALLOW_THREADS;
1202 1
    if (invalid) {
1203 1
        PyErr_Format(PyExc_ValueError,
1204
            "index %" NPY_INTP_FMT " is out of bounds for array with size "
1205
            "%" NPY_INTP_FMT,
1206
            val, unravel_size
1207
        );
1208 1
        return NPY_FAIL;
1209
    }
1210
    return NPY_SUCCEED;
1211
}
1212

1213
/* unravel_index implementation - see add_newdocs.py */
1214
NPY_NO_EXPORT PyObject *
1215 1
arr_unravel_index(PyObject *self, PyObject *args, PyObject *kwds)
1216
{
1217 1
    PyObject *indices0 = NULL;
1218 1
    PyObject *ret_tuple = NULL;
1219 1
    PyArrayObject *ret_arr = NULL;
1220 1
    PyArrayObject *indices = NULL;
1221 1
    PyArray_Descr *dtype = NULL;
1222 1
    PyArray_Dims dimensions = {0, 0};
1223 1
    NPY_ORDER order = NPY_CORDER;
1224
    npy_intp unravel_size;
1225

1226 1
    NpyIter *iter = NULL;
1227
    int i, ret_ndim;
1228
    npy_intp ret_dims[NPY_MAXDIMS], ret_strides[NPY_MAXDIMS];
1229

1230 1
    char *kwlist[] = {"indices", "shape", "order", NULL};
1231

1232
    /*
1233
     * TODO: remove this in favor of warning raised in the dispatcher when
1234
     * __array_function__ is enabled by default.
1235
     */
1236

1237
    /*
1238
     * Continue to support the older "dims" argument in place
1239
     * of the "shape" argument. Issue an appropriate warning
1240
     * if "dims" is detected in keywords, then replace it with
1241
     * the new "shape" argument and continue processing as usual.
1242
     */
1243 1
    if (kwds) {
1244
        PyObject *dims_item, *shape_item;
1245 1
        dims_item = _PyDict_GetItemStringWithError(kwds, "dims");
1246 1
        if (dims_item == NULL && PyErr_Occurred()){
1247
            return NULL;
1248
        }
1249 1
        shape_item = _PyDict_GetItemStringWithError(kwds, "shape");
1250 1
        if (shape_item == NULL && PyErr_Occurred()){
1251
            return NULL;
1252
        }
1253 1
        if (dims_item != NULL && shape_item == NULL) {
1254 1
            if (DEPRECATE("'shape' argument should be"
1255
                          " used instead of 'dims'") < 0) {
1256
                return NULL;
1257
            }
1258 1
            if (PyDict_SetItemString(kwds, "shape", dims_item) < 0) {
1259
                return NULL;
1260
            }
1261 1
            if (PyDict_DelItemString(kwds, "dims") < 0) {
1262
                return NULL;
1263
            }
1264
        }
1265
    }
1266

1267 1
    if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO&|O&:unravel_index",
1268
                    kwlist,
1269
                    &indices0,
1270
                    PyArray_IntpConverter, &dimensions,
1271
                    PyArray_OrderConverter, &order)) {
1272
        goto fail;
1273
    }
1274

1275 1
    unravel_size = PyArray_OverflowMultiplyList(dimensions.ptr, dimensions.len);
1276 1
    if (unravel_size == -1) {
1277 1
        PyErr_SetString(PyExc_ValueError,
1278
                        "dimensions are too large; arrays and shapes with "
1279
                        "a total size greater than 'intp' are not supported.");
1280 1
        goto fail;
1281
    }
1282

1283 1
    indices = astype_anyint(indices0);
1284 1
    if (indices == NULL) {
1285
        goto fail;
1286
    }
1287

1288 1
    dtype = PyArray_DescrFromType(NPY_INTP);
1289 1
    if (dtype == NULL) {
1290
        goto fail;
1291
    }
1292

1293 1
    iter = NpyIter_New(indices, NPY_ITER_READONLY|
1294
                                NPY_ITER_ALIGNED|
1295
                                NPY_ITER_BUFFERED|
1296
                                NPY_ITER_ZEROSIZE_OK|
1297
                                NPY_ITER_DONT_NEGATE_STRIDES|
1298
                                NPY_ITER_MULTI_INDEX,
1299
                                NPY_KEEPORDER, NPY_SAME_KIND_CASTING,
1300
                                dtype);
1301 1
    if (iter == NULL) {
1302
        goto fail;
1303
    }
1304

1305
    /*
1306
     * Create the return array with a layout compatible with the indices
1307
     * and with a dimension added to the end for the multi-index
1308
     */
1309 1
    ret_ndim = PyArray_NDIM(indices) + 1;
1310 1
    if (NpyIter_GetShape(iter, ret_dims) != NPY_SUCCEED) {
1311
        goto fail;
1312
    }
1313 1
    ret_dims[ret_ndim-1] = dimensions.len;
1314 1
    if (NpyIter_CreateCompatibleStrides(iter,
1315 1
                dimensions.len*sizeof(npy_intp), ret_strides) != NPY_SUCCEED) {
1316
        goto fail;
1317
    }
1318 1
    ret_strides[ret_ndim-1] = sizeof(npy_intp);
1319

1320
    /* Remove the multi-index and inner loop */
1321 1
    if (NpyIter_RemoveMultiIndex(iter) != NPY_SUCCEED) {
1322
        goto fail;
1323
    }
1324 1
    if (NpyIter_EnableExternalLoop(iter) != NPY_SUCCEED) {
1325
        goto fail;
1326
    }
1327

1328 1
    ret_arr = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, dtype,
1329
                            ret_ndim, ret_dims, ret_strides, NULL, 0, NULL);
1330 1
    dtype = NULL;
1331 1
    if (ret_arr == NULL) {
1332
        goto fail;
1333
    }
1334

1335 1
    if (order != NPY_CORDER && order != NPY_FORTRANORDER) {
1336 0
        PyErr_SetString(PyExc_ValueError,
1337
                        "only 'C' or 'F' order is permitted");
1338 0
        goto fail;
1339
    }
1340 1
    if (NpyIter_GetIterSize(iter) != 0) {
1341
        NpyIter_IterNextFunc *iternext;
1342
        char **dataptr;
1343
        npy_intp *strides;
1344
        npy_intp *countptr, count;
1345 1
        npy_intp *coordsptr = (npy_intp *)PyArray_DATA(ret_arr);
1346

1347 1
        iternext = NpyIter_GetIterNext(iter, NULL);
1348 1
        if (iternext == NULL) {
1349
            goto fail;
1350
        }
1351 1
        dataptr = NpyIter_GetDataPtrArray(iter);
1352 1
        strides = NpyIter_GetInnerStrideArray(iter);
1353 1
        countptr = NpyIter_GetInnerLoopSizePtr(iter);
1354

1355
        do {
1356 1
            count = *countptr;
1357 1
            if (unravel_index_loop(dimensions.len, dimensions.ptr,
1358
                                   unravel_size, count, *dataptr, *strides,
1359
                                   coordsptr, order) != NPY_SUCCEED) {
1360
                goto fail;
1361
            }
1362 1
            coordsptr += count * dimensions.len;
1363 1
        } while (iternext(iter));
1364
    }
1365

1366

1367 1
    if (dimensions.len == 0 && PyArray_NDIM(indices) != 0) {
1368
        /*
1369
         * There's no index meaning "take the only element 10 times"
1370
         * on a zero-d array, so we have no choice but to error. (See gh-580)
1371
         *
1372
         * Do this check after iterating, so we give a better error message
1373
         * for invalid indices.
1374
         */
1375 1
        PyErr_SetString(PyExc_ValueError,
1376
                "multiple indices are not supported for 0d arrays");
1377 1
        goto fail;
1378
    }
1379

1380
    /* Now make a tuple of views, one per index */
1381 1
    ret_tuple = PyTuple_New(dimensions.len);
1382 1
    if (ret_tuple == NULL) {
1383
        goto fail;
1384
    }
1385 1
    for (i = 0; i < dimensions.len; ++i) {
1386
        PyArrayObject *view;
1387

1388 1
        view = (PyArrayObject *)PyArray_NewFromDescrAndBase(
1389
                &PyArray_Type, PyArray_DescrFromType(NPY_INTP),
1390
                ret_ndim - 1, ret_dims, ret_strides,
1391 1
                PyArray_BYTES(ret_arr) + i*sizeof(npy_intp),
1392
                NPY_ARRAY_WRITEABLE, NULL, (PyObject *)ret_arr);
1393 1
        if (view == NULL) {
1394
            goto fail;
1395
        }
1396 1
        PyTuple_SET_ITEM(ret_tuple, i, PyArray_Return(view));
1397
    }
1398

1399 1
    Py_DECREF(ret_arr);
1400 1
    Py_XDECREF(indices);
1401 1
    npy_free_cache_dim_obj(dimensions);
1402 1
    NpyIter_Deallocate(iter);
1403

1404 1
    return ret_tuple;
1405

1406 1
fail:
1407 1
    Py_XDECREF(ret_tuple);
1408 1
    Py_XDECREF(ret_arr);
1409 1
    Py_XDECREF(dtype);
1410 1
    Py_XDECREF(indices);
1411 1
    npy_free_cache_dim_obj(dimensions);
1412 1
    NpyIter_Deallocate(iter);
1413 1
    return NULL;
1414
}
1415

1416

1417
/* Can only be called if doc is currently NULL */
1418
NPY_NO_EXPORT PyObject *
1419 1
arr_add_docstring(PyObject *NPY_UNUSED(dummy), PyObject *args)
1420
{
1421
    PyObject *obj;
1422
    PyObject *str;
1423
    #if (PY_VERSION_HEX >= 0x030700A2)
1424
    const char *docstr;
1425
    #else
1426
    char *docstr;
1427
    #endif
1428
    static char *msg = "already has a different docstring";
1429

1430
    /* Don't add docstrings */
1431 1
    if (Py_OptimizeFlag > 1) {
1432 0
        Py_RETURN_NONE;
1433
    }
1434

1435 1
    if (!PyArg_ParseTuple(args, "OO!:add_docstring", &obj, &PyUnicode_Type, &str)) {
1436
        return NULL;
1437
    }
1438

1439 1
    docstr = PyUnicode_AsUTF8(str);
1440 1
    if (docstr == NULL) {
1441
        return NULL;
1442
    }
1443

1444
#define _ADDDOC(doc, name)                                              \
1445
        if (!(doc)) {                                                   \
1446
            doc = docstr;                                               \
1447
            Py_INCREF(str);  /* hold on to string (leaks reference) */  \
1448
        }                                                               \
1449
        else if (strcmp(doc, docstr) != 0) {                            \
1450
            PyErr_Format(PyExc_RuntimeError, "%s method %s", name, msg); \
1451
            return NULL;                                                \
1452
        }
1453

1454 1
    if (Py_TYPE(obj) == &PyCFunction_Type) {
1455 1
        PyCFunctionObject *new = (PyCFunctionObject *)obj;
1456 1
        _ADDDOC(new->m_ml->ml_doc, new->m_ml->ml_name);
1457
    }
1458 1
    else if (Py_TYPE(obj) == &PyType_Type) {
1459 1
        PyTypeObject *new = (PyTypeObject *)obj;
1460 1
        _ADDDOC(new->tp_doc, new->tp_name);
1461
    }
1462 1
    else if (Py_TYPE(obj) == &PyMemberDescr_Type) {
1463 1
        PyMemberDescrObject *new = (PyMemberDescrObject *)obj;
1464 1
        _ADDDOC(new->d_member->doc, new->d_member->name);
1465
    }
1466 1
    else if (Py_TYPE(obj) == &PyGetSetDescr_Type) {
1467 1
        PyGetSetDescrObject *new = (PyGetSetDescrObject *)obj;
1468 1
        _ADDDOC(new->d_getset->doc, new->d_getset->name);
1469
    }
1470 1
    else if (Py_TYPE(obj) == &PyMethodDescr_Type) {
1471 1
        PyMethodDescrObject *new = (PyMethodDescrObject *)obj;
1472 1
        _ADDDOC(new->d_method->ml_doc, new->d_method->ml_name);
1473
    }
1474
    else {
1475
        PyObject *doc_attr;
1476

1477 1
        doc_attr = PyObject_GetAttrString(obj, "__doc__");
1478 1
        if (doc_attr != NULL && doc_attr != Py_None &&
1479 1
                (PyUnicode_Compare(doc_attr, str) != 0)) {
1480 1
            Py_DECREF(doc_attr);
1481 1
            if (PyErr_Occurred()) {
1482
                /* error during PyUnicode_Compare */
1483
                return NULL;
1484
            }
1485 1
            PyErr_Format(PyExc_RuntimeError, "object %s", msg);
1486 1
            return NULL;
1487
        }
1488 1
        Py_XDECREF(doc_attr);
1489

1490 1
        if (PyObject_SetAttrString(obj, "__doc__", str) < 0) {
1491 1
            PyErr_SetString(PyExc_TypeError,
1492
                            "Cannot set a docstring for that object");
1493 1
            return NULL;
1494
        }
1495 1
        Py_RETURN_NONE;
1496
    }
1497

1498
#undef _ADDDOC
1499

1500 1
    Py_RETURN_NONE;
1501
}
1502

1503
#if defined NPY_HAVE_SSE2_INTRINSICS
1504
#include <emmintrin.h>
1505
#endif
1506

1507
#ifdef NPY_HAVE_NEON
1508
    typedef npy_uint64 uint64_unaligned __attribute__((aligned(16)));
1509
    static NPY_INLINE int32_t
1510
    sign_mask(uint8x16_t input)
1511
    {
1512
        int8x8_t m0 = vcreate_s8(0x0706050403020100ULL);
1513
        uint8x16_t v0 = vshlq_u8(vshrq_n_u8(input, 7), vcombine_s8(m0, m0));
1514
        uint64x2_t v1 = vpaddlq_u32(vpaddlq_u16(vpaddlq_u8(v0)));
1515
        return (int)vgetq_lane_u64(v1, 0) + ((int)vgetq_lane_u64(v1, 1) << 8);
1516
    }
1517
#endif
1518
/*
1519
 * This function packs boolean values in the input array into the bits of a
1520
 * byte array. Truth values are determined as usual: 0 is false, everything
1521
 * else is true.
1522
 */
1523
static NPY_INLINE void
1524 1
pack_inner(const char *inptr,
1525
           npy_intp element_size,   /* in bytes */
1526
           npy_intp n_in,
1527
           npy_intp in_stride,
1528
           char *outptr,
1529
           npy_intp n_out,
1530
           npy_intp out_stride,
1531
           char order)
1532
{
1533
    /*
1534
     * Loop through the elements of inptr.
1535
     * Determine whether or not it is nonzero.
1536
     *  Yes: set corresponding bit (and adjust build value)
1537
     *  No:  move on
1538
     * Every 8th value, set the value of build and increment the outptr
1539
     */
1540 1
    npy_intp index = 0;
1541 1
    int remain = n_in % 8;              /* uneven bits */
1542

1543
#if defined NPY_HAVE_SSE2_INTRINSICS && defined HAVE__M_FROM_INT64
1544 1
    if (in_stride == 1 && element_size == 1 && n_out > 2) {
1545 1
        __m128i zero = _mm_setzero_si128();
1546
        /* don't handle non-full 8-byte remainder */
1547 1
        npy_intp vn_out = n_out - (remain ? 1 : 0);
1548 1
        vn_out -= (vn_out & 1);
1549 1
        for (index = 0; index < vn_out; index += 2) {
1550
            unsigned int r;
1551 1
            npy_uint64 a = *(npy_uint64*)inptr;
1552 1
            npy_uint64 b = *(npy_uint64*)(inptr + 8);
1553 1
            if (order == 'b') {
1554 1
                a = npy_bswap8(a);
1555 1
                b = npy_bswap8(b);
1556
            }
1557
            
1558
            /* note x86 can load unaligned */
1559 1
            __m128i v = _mm_set_epi64(_m_from_int64(b), _m_from_int64(a));
1560
            /* false -> 0x00 and true -> 0xFF (there is no cmpneq) */
1561 1
            v = _mm_cmpeq_epi8(v, zero);
1562 1
            v = _mm_cmpeq_epi8(v, zero);
1563
            /* extract msb of 16 bytes and pack it into 16 bit */
1564 1
            r = _mm_movemask_epi8(v);
1565
            /* store result */
1566 1
            memcpy(outptr, &r, 1);
1567 1
            outptr += out_stride;
1568 1
            memcpy(outptr, (char*)&r + 1, 1);
1569 1
            outptr += out_stride;
1570 1
            inptr += 16;
1571
        }
1572
    }
1573
#elif defined NPY_HAVE_NEON
1574
    if (in_stride == 1 && element_size == 1 && n_out > 2) {
1575
        /* don't handle non-full 8-byte remainder */
1576
        npy_intp vn_out = n_out - (remain ? 1 : 0);
1577
        vn_out -= (vn_out & 1);
1578
        for (index = 0; index < vn_out; index += 2) {
1579
            unsigned int r;
1580
            npy_uint64 a = *((uint64_unaligned*)inptr);
1581
            npy_uint64 b = *((uint64_unaligned*)(inptr + 8));
1582
            if (order == 'b') {
1583
                a = npy_bswap8(a);
1584
                b = npy_bswap8(b);
1585
            }
1586
            uint64x2_t v = vcombine_u64(vcreate_u64(a), vcreate_u64(b));
1587
            uint64x2_t zero = vdupq_n_u64(0);
1588
            /* false -> 0x00 and true -> 0xFF */
1589
            v = vreinterpretq_u64_u8(vmvnq_u8(vceqq_u8(vreinterpretq_u8_u64(v), vreinterpretq_u8_u64(zero))));
1590
            /* extract msb of 16 bytes and pack it into 16 bit */
1591
            uint8x16_t input = vreinterpretq_u8_u64(v);
1592
            r = sign_mask(input);
1593
            /* store result */
1594
            memcpy(outptr, &r, 1);
1595
            outptr += out_stride;
1596
            memcpy(outptr, (char*)&r + 1, 1);
1597
            outptr += out_stride;
1598
            inptr += 16;
1599
        }
1600
    }
1601
#endif
1602

1603 1
    if (remain == 0) {                  /* assumes n_in > 0 */
1604 1
        remain = 8;
1605
    }
1606
    /* Don't reset index. Just handle remainder of above block */
1607 1
    for (; index < n_out; index++) {
1608 1
        unsigned char build = 0;
1609
        int i, maxi;
1610
        npy_intp j;
1611

1612 1
        maxi = (index == n_out - 1) ? remain : 8;
1613 1
        if (order == 'b') {
1614 1
            for (i = 0; i < maxi; i++) {
1615 1
                build <<= 1;
1616 1
                for (j = 0; j < element_size; j++) {
1617 1
                    build |= (inptr[j] != 0);
1618
                }
1619 1
                inptr += in_stride;
1620
            }
1621 1
            if (index == n_out - 1) {
1622 1
                build <<= 8 - remain;
1623
            }
1624
        }
1625
        else
1626
        {
1627 1
            for (i = 0; i < maxi; i++) {
1628 1
                build >>= 1;
1629 1
                for (j = 0; j < element_size; j++) {
1630 1
                    build |= (inptr[j] != 0) ? 128 : 0;
1631
                }
1632 1
                inptr += in_stride;
1633
            }
1634 1
            if (index == n_out - 1) {
1635 1
                build >>= 8 - remain;
1636
            }
1637
        }
1638 1
        *outptr = (char)build;
1639 1
        outptr += out_stride;
1640
    }
1641
}
1642

1643
static PyObject *
1644 1
pack_bits(PyObject *input, int axis, char order)
1645
{
1646
    PyArrayObject *inp;
1647 1
    PyArrayObject *new = NULL;
1648 1
    PyArrayObject *out = NULL;
1649
    npy_intp outdims[NPY_MAXDIMS];
1650
    int i;
1651
    PyArrayIterObject *it, *ot;
1652 1
    NPY_BEGIN_THREADS_DEF;
1653

1654 1
    inp = (PyArrayObject *)PyArray_FROM_O(input);
1655

1656 1
    if (inp == NULL) {
1657
        return NULL;
1658
    }
1659 1
    if (!PyArray_ISBOOL(inp) && !PyArray_ISINTEGER(inp)) {
1660 1
        PyErr_SetString(PyExc_TypeError,
1661
                "Expected an input array of integer or boolean data type");
1662 1
        Py_DECREF(inp);
1663
        goto fail;
1664
    }
1665

1666 1
    new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0);
1667 1
    Py_DECREF(inp);
1668 1
    if (new == NULL) {
1669
        return NULL;
1670
    }
1671

1672 1
    if (PyArray_NDIM(new) == 0) {
1673
        char *optr, *iptr;
1674

1675 0
        out = (PyArrayObject *)PyArray_NewFromDescr(
1676 0
                Py_TYPE(new), PyArray_DescrFromType(NPY_UBYTE),
1677
                0, NULL, NULL, NULL,
1678
                0, NULL);
1679 0
        if (out == NULL) {
1680
            goto fail;
1681
        }
1682 0
        optr = PyArray_DATA(out);
1683 0
        iptr = PyArray_DATA(new);
1684 0
        *optr = 0;
1685 0
        for (i = 0; i < PyArray_ITEMSIZE(new); i++) {
1686 0
            if (*iptr != 0) {
1687 0
                *optr = 1;
1688 0
                break;
1689
            }
1690 0
            iptr++;
1691
        }
1692
        goto finish;
1693
    }
1694

1695

1696
    /* Setup output shape */
1697 1
    for (i = 0; i < PyArray_NDIM(new); i++) {
1698 1
        outdims[i] = PyArray_DIM(new, i);
1699
    }
1700

1701
    /*
1702
     * Divide axis dimension by 8
1703
     * 8 -> 1, 9 -> 2, 16 -> 2, 17 -> 3 etc..
1704
     */
1705 1
    outdims[axis] = ((outdims[axis] - 1) >> 3) + 1;
1706

1707
    /* Create output array */
1708 1
    out = (PyArrayObject *)PyArray_NewFromDescr(
1709 1
            Py_TYPE(new), PyArray_DescrFromType(NPY_UBYTE),
1710
            PyArray_NDIM(new), outdims, NULL, NULL,
1711 1
            PyArray_ISFORTRAN(new), NULL);
1712 1
    if (out == NULL) {
1713
        goto fail;
1714
    }
1715
    /* Setup iterators to iterate over all but given axis */
1716 1
    it = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)new, &axis);
1717 1
    ot = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)out, &axis);
1718 1
    if (it == NULL || ot == NULL) {
1719 0
        Py_XDECREF(it);
1720 0
        Py_XDECREF(ot);
1721
        goto fail;
1722
    }
1723

1724 1
    NPY_BEGIN_THREADS_THRESHOLDED(PyArray_DIM(out, axis));
1725 1
    while (PyArray_ITER_NOTDONE(it)) {
1726 1
        pack_inner(PyArray_ITER_DATA(it), PyArray_ITEMSIZE(new),
1727
                   PyArray_DIM(new, axis), PyArray_STRIDE(new, axis),
1728
                   PyArray_ITER_DATA(ot), PyArray_DIM(out, axis),
1729
                   PyArray_STRIDE(out, axis), order);
1730 1
        PyArray_ITER_NEXT(it);
1731 1
        PyArray_ITER_NEXT(ot);
1732
    }
1733 1
    NPY_END_THREADS;
1734

1735 1
    Py_DECREF(it);
1736 1
    Py_DECREF(ot);
1737

1738 1
finish:
1739 1
    Py_DECREF(new);
1740
    return (PyObject *)out;
1741

1742 0
fail:
1743 0
    Py_XDECREF(new);
1744 1
    Py_XDECREF(out);
1745
    return NULL;
1746
}
1747

1748
static PyObject *
1749 1
unpack_bits(PyObject *input, int axis, PyObject *count_obj, char order)
1750
{
1751
    static int unpack_init = 0;
1752
    /*
1753
     * lookuptable for bitorder big as it has been around longer
1754
     * bitorder little is handled via byteswapping in the loop
1755
     */
1756
    static union {
1757
        npy_uint8  bytes[8];
1758
        npy_uint64 uint64;
1759
    } unpack_lookup_big[256];
1760
    PyArrayObject *inp;
1761 1
    PyArrayObject *new = NULL;
1762 1
    PyArrayObject *out = NULL;
1763
    npy_intp outdims[NPY_MAXDIMS];
1764
    int i;
1765
    PyArrayIterObject *it, *ot;
1766
    npy_intp count, in_n, in_tail, out_pad, in_stride, out_stride;
1767 1
    NPY_BEGIN_THREADS_DEF;
1768

1769 1
    inp = (PyArrayObject *)PyArray_FROM_O(input);
1770

1771 1
    if (inp == NULL) {
1772
        return NULL;
1773
    }
1774 1
    if (PyArray_TYPE(inp) != NPY_UBYTE) {
1775 0
        PyErr_SetString(PyExc_TypeError,
1776
                "Expected an input array of unsigned byte data type");
1777 0
        Py_DECREF(inp);
1778
        goto fail;
1779
    }
1780

1781 1
    new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0);
1782 1
    Py_DECREF(inp);
1783 1
    if (new == NULL) {
1784
        return NULL;
1785
    }
1786

1787 1
    if (PyArray_NDIM(new) == 0) {
1788
        /* Handle 0-d array by converting it to a 1-d array */
1789
        PyArrayObject *temp;
1790 0
        PyArray_Dims newdim = {NULL, 1};
1791 0
        npy_intp shape = 1;
1792

1793 0
        newdim.ptr = &shape;
1794 0
        temp = (PyArrayObject *)PyArray_Newshape(new, &newdim, NPY_CORDER);
1795 0
        Py_DECREF(new);
1796 0
        if (temp == NULL) {
1797 0
            return NULL;
1798
        }
1799 0
        new = temp;
1800
    }
1801

1802
    /* Setup output shape */
1803 1
    for (i = 0; i < PyArray_NDIM(new); i++) {
1804 1
        outdims[i] = PyArray_DIM(new, i);
1805
    }
1806

1807
    /* Multiply axis dimension by 8 */
1808 1
    outdims[axis] *= 8;
1809 1
    if (count_obj != Py_None) {
1810 1
        count = PyArray_PyIntAsIntp(count_obj);
1811 1
        if (error_converting(count)) {
1812
            goto fail;
1813
        }
1814 1
        if (count < 0) {
1815 1
            outdims[axis] += count;
1816 1
            if (outdims[axis] < 0) {
1817 1
                PyErr_Format(PyExc_ValueError,
1818
                             "-count larger than number of elements");
1819 1
                goto fail;
1820
            }
1821
        }
1822
        else {
1823 1
            outdims[axis] = count;
1824
        }
1825
    }
1826

1827
    /* Create output array */
1828 1
    out = (PyArrayObject *)PyArray_NewFromDescr(
1829 1
            Py_TYPE(new), PyArray_DescrFromType(NPY_UBYTE),
1830
            PyArray_NDIM(new), outdims, NULL, NULL,
1831 1
            PyArray_ISFORTRAN(new), NULL);
1832 1
    if (out == NULL) {
1833
        goto fail;
1834
    }
1835

1836
    /* Setup iterators to iterate over all but given axis */
1837 1
    it = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)new, &axis);
1838 1
    ot = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)out, &axis);
1839 1
    if (it == NULL || ot == NULL) {
1840 0
        Py_XDECREF(it);
1841 0
        Py_XDECREF(ot);
1842
        goto fail;
1843
    }
1844

1845
    /*
1846
     * setup lookup table under GIL, 256 8 byte blocks representing 8 bits
1847
     * expanded to 1/0 bytes
1848
     */
1849 1
    if (unpack_init == 0) {
1850
        npy_intp j;
1851 1
        for (j=0; j < 256; j++) {
1852
            npy_intp k;
1853 1
            for (k=0; k < 8; k++) {
1854 1
                npy_uint8 v = (j & (1 << k)) == (1 << k);
1855 1
                unpack_lookup_big[j].bytes[7 - k] = v;
1856
            }
1857
        }
1858 1
        unpack_init = 1;
1859
    }
1860

1861 1
    count = PyArray_DIM(new, axis) * 8;
1862 1
    if (outdims[axis] > count) {
1863 1
        in_n = count / 8;
1864 1
        in_tail = 0;
1865 1
        out_pad = outdims[axis] - count;
1866
    }
1867
    else {
1868 1
        in_n = outdims[axis] / 8;
1869 1
        in_tail = outdims[axis] % 8;
1870 1
        out_pad = 0;
1871
    }
1872

1873 1
    in_stride = PyArray_STRIDE(new, axis);
1874 1
    out_stride = PyArray_STRIDE(out, axis);
1875

1876 1
    NPY_BEGIN_THREADS_THRESHOLDED(PyArray_Size((PyObject *)out) / 8);
1877

1878 1
    while (PyArray_ITER_NOTDONE(it)) {
1879
        npy_intp index;
1880 1
        unsigned const char *inptr = PyArray_ITER_DATA(it);
1881 1
        char *outptr = PyArray_ITER_DATA(ot);
1882

1883 1
        if (out_stride == 1) {
1884
            /* for unity stride we can just copy out of the lookup table */
1885 1
            if (order == 'b') {
1886 1
                for (index = 0; index < in_n; index++) {
1887 1
                    npy_uint64 v = unpack_lookup_big[*inptr].uint64;
1888 1
                    memcpy(outptr, &v, 8);
1889 1
                    outptr += 8;
1890 1
                    inptr += in_stride;
1891
                }
1892
            }
1893
            else {
1894 1
                for (index = 0; index < in_n; index++) {
1895 1
                    npy_uint64 v = unpack_lookup_big[*inptr].uint64;
1896
                    if (order != 'b') {
1897 1
                        v = npy_bswap8(v);
1898
                    }
1899 1
                    memcpy(outptr, &v, 8);
1900 1
                    outptr += 8;
1901 1
                    inptr += in_stride;
1902
                }
1903
            }
1904
            /* Clean up the tail portion */
1905 1
            if (in_tail) {
1906 1
                npy_uint64 v = unpack_lookup_big[*inptr].uint64;
1907 1
                if (order != 'b') {
1908 1
                    v = npy_bswap8(v);
1909
                }
1910 1
                memcpy(outptr, &v, in_tail);
1911
            }
1912
            /* Add padding */
1913 1
            else if (out_pad) {
1914 1
                memset(outptr, 0, out_pad);
1915
            }
1916
        }
1917
        else {
1918 1
            if (order == 'b') {
1919 1
                for (index = 0; index < in_n; index++) {
1920 1
                    for (i = 0; i < 8; i++) {
1921 1
                        *outptr = ((*inptr & (128 >> i)) != 0);
1922 1
                        outptr += out_stride;
1923
                    }
1924 1
                    inptr += in_stride;
1925
                }
1926
                /* Clean up the tail portion */
1927 1
                for (i = 0; i < in_tail; i++) {
1928 1
                    *outptr = ((*inptr & (128 >> i)) != 0);
1929 1
                    outptr += out_stride;
1930
                }
1931
            }
1932
            else {
1933 1
                for (index = 0; index < in_n; index++) {
1934 1
                    for (i = 0; i < 8; i++) {
1935 1
                        *outptr = ((*inptr & (1 << i)) != 0);
1936 1
                        outptr += out_stride;
1937
                    }
1938 1
                    inptr += in_stride;
1939
                }
1940
                /* Clean up the tail portion */
1941 1
                for (i = 0; i < in_tail; i++) {
1942 1
                    *outptr = ((*inptr & (1 << i)) != 0);
1943 1
                    outptr += out_stride;
1944
                }
1945
            }
1946
            /* Add padding */
1947 0
            for (index = 0; index < out_pad; index++) {
1948 0
                *outptr = 0;
1949 0
                outptr += out_stride;
1950
            }
1951
        }
1952

1953 1
        PyArray_ITER_NEXT(it);
1954 1
        PyArray_ITER_NEXT(ot);
1955
    }
1956 1
    NPY_END_THREADS;
1957

1958 1
    Py_DECREF(it);
1959 1
    Py_DECREF(ot);
1960

1961 1
    Py_DECREF(new);
1962
    return (PyObject *)out;
1963

1964 0
fail:
1965 1
    Py_XDECREF(new);
1966 1
    Py_XDECREF(out);
1967
    return NULL;
1968
}
1969

1970

1971
NPY_NO_EXPORT PyObject *
1972 1
io_pack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
1973
{
1974
    PyObject *obj;
1975 1
    int axis = NPY_MAXDIMS;
1976
    static char *kwlist[] = {"in", "axis", "bitorder", NULL};
1977 1
    char c = 'b';
1978 1
    const char * order_str = NULL;
1979

1980 1
    if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&s:pack" , kwlist,
1981
                &obj, PyArray_AxisConverter, &axis, &order_str)) {
1982
        return NULL;
1983
    }
1984 1
    if (order_str != NULL) {
1985 1
        if (strncmp(order_str, "little", 6) == 0)
1986
            c = 'l';
1987 1
        else if (strncmp(order_str, "big", 3) == 0)
1988
            c = 'b';
1989
        else {
1990 0
            PyErr_SetString(PyExc_ValueError,
1991
                    "'order' must be either 'little' or 'big'");
1992 0
            return NULL;
1993
        }
1994
    }
1995 1
    return pack_bits(obj, axis, c);
1996
}
1997

1998

1999
NPY_NO_EXPORT PyObject *
2000 1
io_unpack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
2001
{
2002
    PyObject *obj;
2003 1
    int axis = NPY_MAXDIMS;
2004 1
    PyObject *count = Py_None;
2005
    static char *kwlist[] = {"in", "axis", "count", "bitorder", NULL};
2006 1
    const char * c = NULL;
2007

2008 1
    if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&Os:unpack" , kwlist,
2009
                &obj, PyArray_AxisConverter, &axis, &count, &c)) {
2010
        return NULL;
2011
    }
2012 1
    if (c == NULL) {
2013 1
        c = "b";
2014
    }
2015 1
    if (c[0] != 'l' && c[0] != 'b') {
2016 1
        PyErr_SetString(PyExc_ValueError,
2017
                    "'order' must begin with 'l' or 'b'");
2018 1
        return NULL;
2019
    }
2020 1
    return unpack_bits(obj, axis, count, c[0]);
2021
}

Read our documentation on viewing source code .

Loading