← Back to branch summary

~jtaylor/ubuntu/precise/python-numpy/multiarch-fix-818867

« back to all changes in this revision

Viewing changes to numpy/core/src/multiarray/calculation.c

  • Committer: Bazaar Package Importer
  • Author(s): Sandro Tosi
  • Date: 2010-10-07 10:19:13 UTC
  • mfrom: (7.1.5 sid)
  • Revision ID: james.westby@ubuntu.com-20101007101913-8b1kmt8ho4upcl9s
Tags: 1:1.4.1-5
http://bugs.debian.org/581058
* debian/patches/10_use_local_python.org_object.inv_sphinx.diff
  - fixed small typo in description
* debian/patches/changeset_r8364.diff
  - fix memory corruption (double free); thanks to Joseph Barillari for the
    report and to Michael Gilbert for pushing resolution; Closes: #581058

Show diffs side-by-side

added added

removed removed

Lines of Context:
numpy/core/src/multiarray/calculation.c
 
1
#define PY_SSIZE_T_CLEAN
 
2
#include <Python.h>
 
3
#include "structmember.h"
 
4
 
 
5
#define _MULTIARRAYMODULE
 
6
#define NPY_NO_PREFIX
 
7
#include "numpy/arrayobject.h"
 
8
 
 
9
#include "npy_config.h"
 
10
 
 
11
#include "common.h"
 
12
#include "number.h"
 
13
 
 
14
#include "calculation.h"
 
15
 
 
16
/* FIXME: just remove _check_axis ? */
 
17
#define _check_axis PyArray_CheckAxis
 
18
#define PyAO PyArrayObject
 
19
 
 
20
static double
 
21
power_of_ten(int n)
 
22
{
 
23
    static const double p10[] = {1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8};
 
24
    double ret;
 
25
    if (n < 9) {
 
26
        ret = p10[n];
 
27
    }
 
28
    else {
 
29
        ret = 1e9;
 
30
        while (n-- > 9) {
 
31
            ret *= 10.;
 
32
        }
 
33
    }
 
34
    return ret;
 
35
}
 
36
 
 
37
/*NUMPY_API
 
38
 * ArgMax
 
39
 */
 
40
NPY_NO_EXPORT PyObject *
 
41
PyArray_ArgMax(PyArrayObject *op, int axis, PyArrayObject *out)
 
42
{
 
43
    PyArrayObject *ap = NULL, *rp = NULL;
 
44
    PyArray_ArgFunc* arg_func;
 
45
    char *ip;
 
46
    intp *rptr;
 
47
    intp i, n, m;
 
48
    int elsize;
 
49
    int copyret = 0;
 
50
    NPY_BEGIN_THREADS_DEF;
 
51
 
 
52
    if ((ap=(PyAO *)_check_axis(op, &axis, 0)) == NULL) {
 
53
        return NULL;
 
54
    }
 
55
    /*
 
56
     * We need to permute the array so that axis is placed at the end.
 
57
     * And all other dimensions are shifted left.
 
58
     */
 
59
    if (axis != ap->nd-1) {
 
60
        PyArray_Dims newaxes;
 
61
        intp dims[MAX_DIMS];
 
62
        int i;
 
63
 
 
64
        newaxes.ptr = dims;
 
65
        newaxes.len = ap->nd;
 
66
        for (i = 0; i < axis; i++) dims[i] = i;
 
67
        for (i = axis; i < ap->nd - 1; i++) dims[i] = i + 1;
 
68
        dims[ap->nd - 1] = axis;
 
69
        op = (PyAO *)PyArray_Transpose(ap, &newaxes);
 
70
        Py_DECREF(ap);
 
71
        if (op == NULL) {
 
72
            return NULL;
 
73
        }
 
74
    }
 
75
    else {
 
76
        op = ap;
 
77
    }
 
78
 
 
79
    /* Will get native-byte order contiguous copy. */
 
80
    ap = (PyArrayObject *)
 
81
        PyArray_ContiguousFromAny((PyObject *)op,
 
82
                                  op->descr->type_num, 1, 0);
 
83
    Py_DECREF(op);
 
84
    if (ap == NULL) {
 
85
        return NULL;
 
86
    }
 
87
    arg_func = ap->descr->f->argmax;
 
88
    if (arg_func == NULL) {
 
89
        PyErr_SetString(PyExc_TypeError, "data type not ordered");
 
90
        goto fail;
 
91
    }
 
92
    elsize = ap->descr->elsize;
 
93
    m = ap->dimensions[ap->nd-1];
 
94
    if (m == 0) {
 
95
        PyErr_SetString(PyExc_ValueError,
 
96
                        "attempt to get argmax/argmin "\
 
97
                        "of an empty sequence");
 
98
        goto fail;
 
99
    }
 
100
 
 
101
    if (!out) {
 
102
        rp = (PyArrayObject *)PyArray_New(ap->ob_type, ap->nd-1,
 
103
                                          ap->dimensions, PyArray_INTP,
 
104
                                          NULL, NULL, 0, 0,
 
105
                                          (PyObject *)ap);
 
106
        if (rp == NULL) {
 
107
            goto fail;
 
108
        }
 
109
    }
 
110
    else {
 
111
        if (PyArray_SIZE(out) !=
 
112
                PyArray_MultiplyList(ap->dimensions, ap->nd - 1)) {
 
113
            PyErr_SetString(PyExc_TypeError,
 
114
                            "invalid shape for output array.");
 
115
        }
 
116
        rp = (PyArrayObject *)\
 
117
            PyArray_FromArray(out,
 
118
                              PyArray_DescrFromType(PyArray_INTP),
 
119
                              NPY_CARRAY | NPY_UPDATEIFCOPY);
 
120
        if (rp == NULL) {
 
121
            goto fail;
 
122
        }
 
123
        if (rp != out) {
 
124
            copyret = 1;
 
125
        }
 
126
    }
 
127
 
 
128
    NPY_BEGIN_THREADS_DESCR(ap->descr);
 
129
    n = PyArray_SIZE(ap)/m;
 
130
    rptr = (intp *)rp->data;
 
131
    for (ip = ap->data, i = 0; i < n; i++, ip += elsize*m) {
 
132
        arg_func(ip, m, rptr, ap);
 
133
        rptr += 1;
 
134
    }
 
135
    NPY_END_THREADS_DESCR(ap->descr);
 
136
 
 
137
    Py_DECREF(ap);
 
138
    if (copyret) {
 
139
        PyArrayObject *obj;
 
140
        obj = (PyArrayObject *)rp->base;
 
141
        Py_INCREF(obj);
 
142
        Py_DECREF(rp);
 
143
        rp = obj;
 
144
    }
 
145
    return (PyObject *)rp;
 
146
 
 
147
 fail:
 
148
    Py_DECREF(ap);
 
149
    Py_XDECREF(rp);
 
150
    return NULL;
 
151
}
 
152
 
 
153
/*NUMPY_API
 
154
 * ArgMin
 
155
 */
 
156
NPY_NO_EXPORT PyObject *
 
157
PyArray_ArgMin(PyArrayObject *ap, int axis, PyArrayObject *out)
 
158
{
 
159
    PyObject *obj, *new, *ret;
 
160
 
 
161
    if (PyArray_ISFLEXIBLE(ap)) {
 
162
        PyErr_SetString(PyExc_TypeError,
 
163
                        "argmax is unsupported for this type");
 
164
        return NULL;
 
165
    }
 
166
    else if (PyArray_ISUNSIGNED(ap)) {
 
167
        obj = PyInt_FromLong((long) -1);
 
168
    }
 
169
    else if (PyArray_TYPE(ap) == PyArray_BOOL) {
 
170
        obj = PyInt_FromLong((long) 1);
 
171
    }
 
172
    else {
 
173
        obj = PyInt_FromLong((long) 0);
 
174
    }
 
175
    new = PyArray_EnsureAnyArray(PyNumber_Subtract(obj, (PyObject *)ap));
 
176
    Py_DECREF(obj);
 
177
    if (new == NULL) {
 
178
        return NULL;
 
179
    }
 
180
    ret = PyArray_ArgMax((PyArrayObject *)new, axis, out);
 
181
    Py_DECREF(new);
 
182
    return ret;
 
183
}
 
184
 
 
185
/*NUMPY_API
 
186
 * Max
 
187
 */
 
188
NPY_NO_EXPORT PyObject *
 
189
PyArray_Max(PyArrayObject *ap, int axis, PyArrayObject *out)
 
190
{
 
191
    PyArrayObject *arr;
 
192
    PyObject *ret;
 
193
 
 
194
    if ((arr=(PyArrayObject *)_check_axis(ap, &axis, 0)) == NULL) {
 
195
        return NULL;
 
196
    }
 
197
    ret = PyArray_GenericReduceFunction(arr, n_ops.maximum, axis,
 
198
                                        arr->descr->type_num, out);
 
199
    Py_DECREF(arr);
 
200
    return ret;
 
201
}
 
202
 
 
203
/*NUMPY_API
 
204
 * Min
 
205
 */
 
206
NPY_NO_EXPORT PyObject *
 
207
PyArray_Min(PyArrayObject *ap, int axis, PyArrayObject *out)
 
208
{
 
209
    PyArrayObject *arr;
 
210
    PyObject *ret;
 
211
 
 
212
    if ((arr=(PyArrayObject *)_check_axis(ap, &axis, 0)) == NULL) {
 
213
        return NULL;
 
214
    }
 
215
    ret = PyArray_GenericReduceFunction(arr, n_ops.minimum, axis,
 
216
                                        arr->descr->type_num, out);
 
217
    Py_DECREF(arr);
 
218
    return ret;
 
219
}
 
220
 
 
221
/*NUMPY_API
 
222
 * Ptp
 
223
 */
 
224
NPY_NO_EXPORT PyObject *
 
225
PyArray_Ptp(PyArrayObject *ap, int axis, PyArrayObject *out)
 
226
{
 
227
    PyArrayObject *arr;
 
228
    PyObject *ret;
 
229
    PyObject *obj1 = NULL, *obj2 = NULL;
 
230
 
 
231
    if ((arr=(PyArrayObject *)_check_axis(ap, &axis, 0)) == NULL) {
 
232
        return NULL;
 
233
    }
 
234
    obj1 = PyArray_Max(arr, axis, out);
 
235
    if (obj1 == NULL) {
 
236
        goto fail;
 
237
    }
 
238
    obj2 = PyArray_Min(arr, axis, NULL);
 
239
    if (obj2 == NULL) {
 
240
        goto fail;
 
241
    }
 
242
    Py_DECREF(arr);
 
243
    if (out) {
 
244
        ret = PyObject_CallFunction(n_ops.subtract, "OOO", out, obj2, out);
 
245
    }
 
246
    else {
 
247
        ret = PyNumber_Subtract(obj1, obj2);
 
248
    }
 
249
    Py_DECREF(obj1);
 
250
    Py_DECREF(obj2);
 
251
    return ret;
 
252
 
 
253
 fail:
 
254
    Py_XDECREF(arr);
 
255
    Py_XDECREF(obj1);
 
256
    Py_XDECREF(obj2);
 
257
    return NULL;
 
258
}
 
259
 
 
260
 
 
261
 
 
262
/*NUMPY_API
 
263
 * Set variance to 1 to by-pass square-root calculation and return variance
 
264
 * Std
 
265
 */
 
266
NPY_NO_EXPORT PyObject *
 
267
PyArray_Std(PyArrayObject *self, int axis, int rtype, PyArrayObject *out,
 
268
            int variance)
 
269
{
 
270
    return __New_PyArray_Std(self, axis, rtype, out, variance, 0);
 
271
}
 
272
 
 
273
NPY_NO_EXPORT PyObject *
 
274
__New_PyArray_Std(PyArrayObject *self, int axis, int rtype, PyArrayObject *out,
 
275
                  int variance, int num)
 
276
{
 
277
    PyObject *obj1 = NULL, *obj2 = NULL, *obj3 = NULL, *new = NULL;
 
278
    PyObject *ret = NULL, *newshape = NULL;
 
279
    int i, n;
 
280
    intp val;
 
281
 
 
282
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
283
        return NULL;
 
284
    }
 
285
    /* Compute and reshape mean */
 
286
    obj1 = PyArray_EnsureAnyArray(PyArray_Mean((PyAO *)new, axis, rtype, NULL));
 
287
    if (obj1 == NULL) {
 
288
        Py_DECREF(new);
 
289
        return NULL;
 
290
    }
 
291
    n = PyArray_NDIM(new);
 
292
    newshape = PyTuple_New(n);
 
293
    if (newshape == NULL) {
 
294
        Py_DECREF(obj1);
 
295
        Py_DECREF(new);
 
296
        return NULL;
 
297
    }
 
298
    for (i = 0; i < n; i++) {
 
299
        if (i == axis) {
 
300
            val = 1;
 
301
        }
 
302
        else {
 
303
            val = PyArray_DIM(new,i);
 
304
        }
 
305
        PyTuple_SET_ITEM(newshape, i, PyInt_FromLong((long)val));
 
306
    }
 
307
    obj2 = PyArray_Reshape((PyAO *)obj1, newshape);
 
308
    Py_DECREF(obj1);
 
309
    Py_DECREF(newshape);
 
310
    if (obj2 == NULL) {
 
311
        Py_DECREF(new);
 
312
        return NULL;
 
313
    }
 
314
 
 
315
    /* Compute x = x - mx */
 
316
    obj1 = PyArray_EnsureAnyArray(PyNumber_Subtract((PyObject *)new, obj2));
 
317
    Py_DECREF(obj2);
 
318
    if (obj1 == NULL) {
 
319
        Py_DECREF(new);
 
320
        return NULL;
 
321
    }
 
322
    /* Compute x * x */
 
323
    if (PyArray_ISCOMPLEX(obj1)) {
 
324
        obj3 = PyArray_Conjugate((PyAO *)obj1, NULL);
 
325
    }
 
326
    else {
 
327
        obj3 = obj1;
 
328
        Py_INCREF(obj1);
 
329
    }
 
330
    if (obj3 == NULL) {
 
331
        Py_DECREF(new);
 
332
        return NULL;
 
333
    }
 
334
    obj2 = PyArray_EnsureAnyArray                                      \
 
335
        (PyArray_GenericBinaryFunction((PyAO *)obj1, obj3, n_ops.multiply));
 
336
    Py_DECREF(obj1);
 
337
    Py_DECREF(obj3);
 
338
    if (obj2 == NULL) {
 
339
        Py_DECREF(new);
 
340
        return NULL;
 
341
    }
 
342
    if (PyArray_ISCOMPLEX(obj2)) {
 
343
        obj3 = PyObject_GetAttrString(obj2, "real");
 
344
        switch(rtype) {
 
345
        case NPY_CDOUBLE:
 
346
            rtype = NPY_DOUBLE;
 
347
            break;
 
348
        case NPY_CFLOAT:
 
349
            rtype = NPY_FLOAT;
 
350
            break;
 
351
        case NPY_CLONGDOUBLE:
 
352
            rtype = NPY_LONGDOUBLE;
 
353
            break;
 
354
        }
 
355
    }
 
356
    else {
 
357
        obj3 = obj2;
 
358
        Py_INCREF(obj2);
 
359
    }
 
360
    if (obj3 == NULL) {
 
361
        Py_DECREF(new);
 
362
        return NULL;
 
363
    }
 
364
    /* Compute add.reduce(x*x,axis) */
 
365
    obj1 = PyArray_GenericReduceFunction((PyAO *)obj3, n_ops.add,
 
366
                                         axis, rtype, NULL);
 
367
    Py_DECREF(obj3);
 
368
    Py_DECREF(obj2);
 
369
    if (obj1 == NULL) {
 
370
        Py_DECREF(new);
 
371
        return NULL;
 
372
    }
 
373
    n = PyArray_DIM(new,axis);
 
374
    Py_DECREF(new);
 
375
    n = (n-num);
 
376
    if (n == 0) {
 
377
        n = 1;
 
378
    }
 
379
    obj2 = PyFloat_FromDouble(1.0/((double )n));
 
380
    if (obj2 == NULL) {
 
381
        Py_DECREF(obj1);
 
382
        return NULL;
 
383
    }
 
384
    ret = PyNumber_Multiply(obj1, obj2);
 
385
    Py_DECREF(obj1);
 
386
    Py_DECREF(obj2);
 
387
 
 
388
    if (!variance) {
 
389
        obj1 = PyArray_EnsureAnyArray(ret);
 
390
        /* sqrt() */
 
391
        ret = PyArray_GenericUnaryFunction((PyAO *)obj1, n_ops.sqrt);
 
392
        Py_DECREF(obj1);
 
393
    }
 
394
    if (ret == NULL || PyArray_CheckExact(self)) {
 
395
        return ret;
 
396
    }
 
397
    if (PyArray_Check(self) && self->ob_type == ret->ob_type) {
 
398
        return ret;
 
399
    }
 
400
    obj1 = PyArray_EnsureArray(ret);
 
401
    if (obj1 == NULL) {
 
402
        return NULL;
 
403
    }
 
404
    ret = PyArray_View((PyAO *)obj1, NULL, self->ob_type);
 
405
    Py_DECREF(obj1);
 
406
    if (out) {
 
407
        if (PyArray_CopyAnyInto(out, (PyArrayObject *)ret) < 0) {
 
408
            Py_DECREF(ret);
 
409
            return NULL;
 
410
        }
 
411
        Py_DECREF(ret);
 
412
        Py_INCREF(out);
 
413
        return (PyObject *)out;
 
414
    }
 
415
    return ret;
 
416
}
 
417
 
 
418
 
 
419
/*NUMPY_API
 
420
 *Sum
 
421
 */
 
422
NPY_NO_EXPORT PyObject *
 
423
PyArray_Sum(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 
424
{
 
425
    PyObject *new, *ret;
 
426
 
 
427
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
428
        return NULL;
 
429
    }
 
430
    ret = PyArray_GenericReduceFunction((PyAO *)new, n_ops.add, axis,
 
431
                                        rtype, out);
 
432
    Py_DECREF(new);
 
433
    return ret;
 
434
}
 
435
 
 
436
/*NUMPY_API
 
437
 * Prod
 
438
 */
 
439
NPY_NO_EXPORT PyObject *
 
440
PyArray_Prod(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 
441
{
 
442
    PyObject *new, *ret;
 
443
 
 
444
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
445
        return NULL;
 
446
    }
 
447
    ret = PyArray_GenericReduceFunction((PyAO *)new, n_ops.multiply, axis,
 
448
                                        rtype, out);
 
449
    Py_DECREF(new);
 
450
    return ret;
 
451
}
 
452
 
 
453
/*NUMPY_API
 
454
 *CumSum
 
455
 */
 
456
NPY_NO_EXPORT PyObject *
 
457
PyArray_CumSum(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 
458
{
 
459
    PyObject *new, *ret;
 
460
 
 
461
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
462
        return NULL;
 
463
    }
 
464
    ret = PyArray_GenericAccumulateFunction((PyAO *)new, n_ops.add, axis,
 
465
                                            rtype, out);
 
466
    Py_DECREF(new);
 
467
    return ret;
 
468
}
 
469
 
 
470
/*NUMPY_API
 
471
 * CumProd
 
472
 */
 
473
NPY_NO_EXPORT PyObject *
 
474
PyArray_CumProd(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 
475
{
 
476
    PyObject *new, *ret;
 
477
 
 
478
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
479
        return NULL;
 
480
    }
 
481
 
 
482
    ret = PyArray_GenericAccumulateFunction((PyAO *)new,
 
483
                                            n_ops.multiply, axis,
 
484
                                            rtype, out);
 
485
    Py_DECREF(new);
 
486
    return ret;
 
487
}
 
488
 
 
489
/*NUMPY_API
 
490
 * Round
 
491
 */
 
492
NPY_NO_EXPORT PyObject *
 
493
PyArray_Round(PyArrayObject *a, int decimals, PyArrayObject *out)
 
494
{
 
495
    PyObject *f, *ret = NULL, *tmp, *op1, *op2;
 
496
    int ret_int=0;
 
497
    PyArray_Descr *my_descr;
 
498
    if (out && (PyArray_SIZE(out) != PyArray_SIZE(a))) {
 
499
        PyErr_SetString(PyExc_ValueError,
 
500
                        "invalid output shape");
 
501
        return NULL;
 
502
    }
 
503
    if (PyArray_ISCOMPLEX(a)) {
 
504
        PyObject *part;
 
505
        PyObject *round_part;
 
506
        PyObject *new;
 
507
        int res;
 
508
 
 
509
        if (out) {
 
510
            new = (PyObject *)out;
 
511
            Py_INCREF(new);
 
512
        }
 
513
        else {
 
514
            new = PyArray_Copy(a);
 
515
            if (new == NULL) {
 
516
                return NULL;
 
517
            }
 
518
        }
 
519
 
 
520
        /* new.real = a.real.round(decimals) */
 
521
        part = PyObject_GetAttrString(new, "real");
 
522
        if (part == NULL) {
 
523
            Py_DECREF(new);
 
524
            return NULL;
 
525
        }
 
526
        part = PyArray_EnsureAnyArray(part);
 
527
        round_part = PyArray_Round((PyArrayObject *)part,
 
528
                                   decimals, NULL);
 
529
        Py_DECREF(part);
 
530
        if (round_part == NULL) {
 
531
            Py_DECREF(new);
 
532
            return NULL;
 
533
        }
 
534
        res = PyObject_SetAttrString(new, "real", round_part);
 
535
        Py_DECREF(round_part);
 
536
        if (res < 0) {
 
537
            Py_DECREF(new);
 
538
            return NULL;
 
539
        }
 
540
 
 
541
        /* new.imag = a.imag.round(decimals) */
 
542
        part = PyObject_GetAttrString(new, "imag");
 
543
        if (part == NULL) {
 
544
            Py_DECREF(new);
 
545
            return NULL;
 
546
        }
 
547
        part = PyArray_EnsureAnyArray(part);
 
548
        round_part = PyArray_Round((PyArrayObject *)part,
 
549
                                   decimals, NULL);
 
550
        Py_DECREF(part);
 
551
        if (round_part == NULL) {
 
552
            Py_DECREF(new);
 
553
            return NULL;
 
554
        }
 
555
        res = PyObject_SetAttrString(new, "imag", round_part);
 
556
        Py_DECREF(round_part);
 
557
        if (res < 0) {
 
558
            Py_DECREF(new);
 
559
            return NULL;
 
560
        }
 
561
        return new;
 
562
    }
 
563
    /* do the most common case first */
 
564
    if (decimals >= 0) {
 
565
        if (PyArray_ISINTEGER(a)) {
 
566
            if (out) {
 
567
                if (PyArray_CopyAnyInto(out, a) < 0) {
 
568
                    return NULL;
 
569
                }
 
570
                Py_INCREF(out);
 
571
                return (PyObject *)out;
 
572
            }
 
573
            else {
 
574
                Py_INCREF(a);
 
575
                return (PyObject *)a;
 
576
            }
 
577
        }
 
578
        if (decimals == 0) {
 
579
            if (out) {
 
580
                return PyObject_CallFunction(n_ops.rint, "OO", a, out);
 
581
            }
 
582
            return PyObject_CallFunction(n_ops.rint, "O", a);
 
583
        }
 
584
        op1 = n_ops.multiply;
 
585
        op2 = n_ops.true_divide;
 
586
    }
 
587
    else {
 
588
        op1 = n_ops.true_divide;
 
589
        op2 = n_ops.multiply;
 
590
        decimals = -decimals;
 
591
    }
 
592
    if (!out) {
 
593
        if (PyArray_ISINTEGER(a)) {
 
594
            ret_int = 1;
 
595
            my_descr = PyArray_DescrFromType(NPY_DOUBLE);
 
596
        }
 
597
        else {
 
598
            Py_INCREF(a->descr);
 
599
            my_descr = a->descr;
 
600
        }
 
601
        out = (PyArrayObject *)PyArray_Empty(a->nd, a->dimensions,
 
602
                                             my_descr,
 
603
                                             PyArray_ISFORTRAN(a));
 
604
        if (out == NULL) {
 
605
            return NULL;
 
606
        }
 
607
    }
 
608
    else {
 
609
        Py_INCREF(out);
 
610
    }
 
611
    f = PyFloat_FromDouble(power_of_ten(decimals));
 
612
    if (f == NULL) {
 
613
        return NULL;
 
614
    }
 
615
    ret = PyObject_CallFunction(op1, "OOO", a, f, out);
 
616
    if (ret == NULL) {
 
617
        goto finish;
 
618
    }
 
619
    tmp = PyObject_CallFunction(n_ops.rint, "OO", ret, ret);
 
620
    if (tmp == NULL) {
 
621
        Py_DECREF(ret);
 
622
        ret = NULL;
 
623
        goto finish;
 
624
    }
 
625
    Py_DECREF(tmp);
 
626
    tmp = PyObject_CallFunction(op2, "OOO", ret, f, ret);
 
627
    if (tmp == NULL) {
 
628
        Py_DECREF(ret);
 
629
        ret = NULL;
 
630
        goto finish;
 
631
    }
 
632
    Py_DECREF(tmp);
 
633
 
 
634
 finish:
 
635
    Py_DECREF(f);
 
636
    Py_DECREF(out);
 
637
    if (ret_int) {
 
638
        Py_INCREF(a->descr);
 
639
        tmp = PyArray_CastToType((PyArrayObject *)ret,
 
640
                                 a->descr, PyArray_ISFORTRAN(a));
 
641
        Py_DECREF(ret);
 
642
        return tmp;
 
643
    }
 
644
    return ret;
 
645
}
 
646
 
 
647
 
 
648
/*NUMPY_API
 
649
 * Mean
 
650
 */
 
651
NPY_NO_EXPORT PyObject *
 
652
PyArray_Mean(PyArrayObject *self, int axis, int rtype, PyArrayObject *out)
 
653
{
 
654
    PyObject *obj1 = NULL, *obj2 = NULL;
 
655
    PyObject *new, *ret;
 
656
 
 
657
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
658
        return NULL;
 
659
    }
 
660
    obj1 = PyArray_GenericReduceFunction((PyAO *)new, n_ops.add, axis,
 
661
                                         rtype, out);
 
662
    obj2 = PyFloat_FromDouble((double) PyArray_DIM(new,axis));
 
663
    Py_DECREF(new);
 
664
    if (obj1 == NULL || obj2 == NULL) {
 
665
        Py_XDECREF(obj1);
 
666
        Py_XDECREF(obj2);
 
667
        return NULL;
 
668
    }
 
669
    if (!out) {
 
670
        ret = PyNumber_Divide(obj1, obj2);
 
671
    }
 
672
    else {
 
673
        ret = PyObject_CallFunction(n_ops.divide, "OOO", out, obj2, out);
 
674
    }
 
675
    Py_DECREF(obj1);
 
676
    Py_DECREF(obj2);
 
677
    return ret;
 
678
}
 
679
 
 
680
/*NUMPY_API
 
681
 * Any
 
682
 */
 
683
NPY_NO_EXPORT PyObject *
 
684
PyArray_Any(PyArrayObject *self, int axis, PyArrayObject *out)
 
685
{
 
686
    PyObject *new, *ret;
 
687
 
 
688
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
689
        return NULL;
 
690
    }
 
691
    ret = PyArray_GenericReduceFunction((PyAO *)new,
 
692
                                        n_ops.logical_or, axis,
 
693
                                        PyArray_BOOL, out);
 
694
    Py_DECREF(new);
 
695
    return ret;
 
696
}
 
697
 
 
698
/*NUMPY_API
 
699
 * All
 
700
 */
 
701
NPY_NO_EXPORT PyObject *
 
702
PyArray_All(PyArrayObject *self, int axis, PyArrayObject *out)
 
703
{
 
704
    PyObject *new, *ret;
 
705
 
 
706
    if ((new = _check_axis(self, &axis, 0)) == NULL) {
 
707
        return NULL;
 
708
    }
 
709
    ret = PyArray_GenericReduceFunction((PyAO *)new,
 
710
                                        n_ops.logical_and, axis,
 
711
                                        PyArray_BOOL, out);
 
712
    Py_DECREF(new);
 
713
    return ret;
 
714
}
 
715
 
 
716
 
 
717
static PyObject *
 
718
_GenericBinaryOutFunction(PyArrayObject *m1, PyObject *m2, PyArrayObject *out,
 
719
                          PyObject *op)
 
720
{
 
721
    if (out == NULL) {
 
722
        return PyObject_CallFunction(op, "OO", m1, m2);
 
723
    }
 
724
    else {
 
725
        return PyObject_CallFunction(op, "OOO", m1, m2, out);
 
726
    }
 
727
}
 
728
 
 
729
static PyObject *
 
730
_slow_array_clip(PyArrayObject *self, PyObject *min, PyObject *max, PyArrayObject *out)
 
731
{
 
732
    PyObject *res1=NULL, *res2=NULL;
 
733
 
 
734
    if (max != NULL) {
 
735
        res1 = _GenericBinaryOutFunction(self, max, out, n_ops.minimum);
 
736
        if (res1 == NULL) {
 
737
            return NULL;
 
738
        }
 
739
    }
 
740
    else {
 
741
        res1 = (PyObject *)self;
 
742
        Py_INCREF(res1);
 
743
    }
 
744
 
 
745
    if (min != NULL) {
 
746
        res2 = _GenericBinaryOutFunction((PyArrayObject *)res1,
 
747
                                         min, out, n_ops.maximum);
 
748
        if (res2 == NULL) {
 
749
            Py_XDECREF(res1);
 
750
            return NULL;
 
751
        }
 
752
    }
 
753
    else {
 
754
        res2 = res1;
 
755
        Py_INCREF(res2);
 
756
    }
 
757
    Py_DECREF(res1);
 
758
    return res2;
 
759
}
 
760
 
 
761
/*NUMPY_API
 
762
 * Clip
 
763
 */
 
764
NPY_NO_EXPORT PyObject *
 
765
PyArray_Clip(PyArrayObject *self, PyObject *min, PyObject *max, PyArrayObject *out)
 
766
{
 
767
    PyArray_FastClipFunc *func;
 
768
    int outgood = 0, ingood = 0;
 
769
    PyArrayObject *maxa = NULL;
 
770
    PyArrayObject *mina = NULL;
 
771
    PyArrayObject *newout = NULL, *newin = NULL;
 
772
    PyArray_Descr *indescr, *newdescr;
 
773
    char *max_data, *min_data;
 
774
    PyObject *zero;
 
775
 
 
776
    if ((max == NULL) && (min == NULL)) {
 
777
        PyErr_SetString(PyExc_ValueError, "array_clip: must set either max "\
 
778
                        "or min");
 
779
        return NULL;
 
780
    }
 
781
 
 
782
    func = self->descr->f->fastclip;
 
783
    if (func == NULL || (min != NULL && !PyArray_CheckAnyScalar(min)) ||
 
784
        (max != NULL && !PyArray_CheckAnyScalar(max))) {
 
785
        return _slow_array_clip(self, min, max, out);
 
786
    }
 
787
    /* Use the fast scalar clip function */
 
788
 
 
789
    /* First we need to figure out the correct type */
 
790
    indescr = NULL;
 
791
    if (min != NULL) {
 
792
        indescr = PyArray_DescrFromObject(min, NULL);
 
793
        if (indescr == NULL) {
 
794
            return NULL;
 
795
        }
 
796
    }
 
797
    if (max != NULL) {
 
798
        newdescr = PyArray_DescrFromObject(max, indescr);
 
799
        Py_XDECREF(indescr);
 
800
        if (newdescr == NULL) {
 
801
            return NULL;
 
802
        }
 
803
    }
 
804
    else {
 
805
        /* Steal the reference */
 
806
        newdescr = indescr;
 
807
    }
 
808
 
 
809
 
 
810
    /*
 
811
     * Use the scalar descriptor only if it is of a bigger
 
812
     * KIND than the input array (and then find the
 
813
     * type that matches both).
 
814
     */
 
815
    if (PyArray_ScalarKind(newdescr->type_num, NULL) >
 
816
        PyArray_ScalarKind(self->descr->type_num, NULL)) {
 
817
        indescr = _array_small_type(newdescr, self->descr);
 
818
        func = indescr->f->fastclip;
 
819
        if (func == NULL) {
 
820
            return _slow_array_clip(self, min, max, out);
 
821
        }
 
822
    }
 
823
    else {
 
824
        indescr = self->descr;
 
825
        Py_INCREF(indescr);
 
826
    }
 
827
    Py_DECREF(newdescr);
 
828
 
 
829
    if (!PyDataType_ISNOTSWAPPED(indescr)) {
 
830
        PyArray_Descr *descr2;
 
831
        descr2 = PyArray_DescrNewByteorder(indescr, '=');
 
832
        Py_DECREF(indescr);
 
833
        if (descr2 == NULL) {
 
834
            goto fail;
 
835
        }
 
836
        indescr = descr2;
 
837
    }
 
838
 
 
839
    /* Convert max to an array */
 
840
    if (max != NULL) {
 
841
        maxa = (NPY_AO *)PyArray_FromAny(max, indescr, 0, 0,
 
842
                                         NPY_DEFAULT, NULL);
 
843
        if (maxa == NULL) {
 
844
            return NULL;
 
845
        }
 
846
    }
 
847
    else {
 
848
        /* Side-effect of PyArray_FromAny */
 
849
        Py_DECREF(indescr);
 
850
    }
 
851
 
 
852
    /*
 
853
     * If we are unsigned, then make sure min is not < 0
 
854
     * This is to match the behavior of _slow_array_clip
 
855
     *
 
856
     * We allow min and max to go beyond the limits
 
857
     * for other data-types in which case they
 
858
     * are interpreted as their modular counterparts.
 
859
    */
 
860
    if (min != NULL) {
 
861
        if (PyArray_ISUNSIGNED(self)) {
 
862
            int cmp;
 
863
            zero = PyInt_FromLong(0);
 
864
            cmp = PyObject_RichCompareBool(min, zero, Py_LT);
 
865
            if (cmp == -1) {
 
866
                Py_DECREF(zero);
 
867
                goto fail;
 
868
            }
 
869
            if (cmp == 1) {
 
870
                min = zero;
 
871
            }
 
872
            else {
 
873
                Py_DECREF(zero);
 
874
                Py_INCREF(min);
 
875
            }
 
876
        }
 
877
        else {
 
878
            Py_INCREF(min);
 
879
        }
 
880
 
 
881
        /* Convert min to an array */
 
882
        Py_INCREF(indescr);
 
883
        mina = (NPY_AO *)PyArray_FromAny(min, indescr, 0, 0,
 
884
                                         NPY_DEFAULT, NULL);
 
885
        Py_DECREF(min);
 
886
        if (mina == NULL) {
 
887
            goto fail;
 
888
        }
 
889
    }
 
890
 
 
891
 
 
892
    /*
 
893
     * Check to see if input is single-segment, aligned,
 
894
     * and in native byteorder
 
895
     */
 
896
    if (PyArray_ISONESEGMENT(self) && PyArray_CHKFLAGS(self, ALIGNED) &&
 
897
        PyArray_ISNOTSWAPPED(self) && (self->descr == indescr)) {
 
898
        ingood = 1;
 
899
    }
 
900
    if (!ingood) {
 
901
        int flags;
 
902
 
 
903
        if (PyArray_ISFORTRAN(self)) {
 
904
            flags = NPY_FARRAY;
 
905
        }
 
906
        else {
 
907
            flags = NPY_CARRAY;
 
908
        }
 
909
        Py_INCREF(indescr);
 
910
        newin = (NPY_AO *)PyArray_FromArray(self, indescr, flags);
 
911
        if (newin == NULL) {
 
912
            goto fail;
 
913
        }
 
914
    }
 
915
    else {
 
916
        newin = self;
 
917
        Py_INCREF(newin);
 
918
    }
 
919
 
 
920
    /*
 
921
     * At this point, newin is a single-segment, aligned, and correct
 
922
     * byte-order array of the correct type
 
923
     *
 
924
     * if ingood == 0, then it is a copy, otherwise,
 
925
     * it is the original input.
 
926
     */
 
927
 
 
928
    /*
 
929
     * If we have already made a copy of the data, then use
 
930
     * that as the output array
 
931
     */
 
932
    if (out == NULL && !ingood) {
 
933
        out = newin;
 
934
    }
 
935
 
 
936
    /*
 
937
     * Now, we know newin is a usable array for fastclip,
 
938
     * we need to make sure the output array is available
 
939
     * and usable
 
940
     */
 
941
    if (out == NULL) {
 
942
        Py_INCREF(indescr);
 
943
        out = (NPY_AO*)PyArray_NewFromDescr(self->ob_type,
 
944
                                            indescr, self->nd,
 
945
                                            self->dimensions,
 
946
                                            NULL, NULL,
 
947
                                            PyArray_ISFORTRAN(self),
 
948
                                            (PyObject *)self);
 
949
        if (out == NULL) {
 
950
            goto fail;
 
951
        }
 
952
        outgood = 1;
 
953
    }
 
954
    else Py_INCREF(out);
 
955
    /* Input is good at this point */
 
956
    if (out == newin) {
 
957
        outgood = 1;
 
958
    }
 
959
    if (!outgood && PyArray_ISONESEGMENT(out) &&
 
960
        PyArray_CHKFLAGS(out, ALIGNED) && PyArray_ISNOTSWAPPED(out) &&
 
961
        PyArray_EquivTypes(out->descr, indescr)) {
 
962
        outgood = 1;
 
963
    }
 
964
 
 
965
    /*
 
966
     * Do we still not have a suitable output array?
 
967
     * Create one, now
 
968
     */
 
969
    if (!outgood) {
 
970
        int oflags;
 
971
        if (PyArray_ISFORTRAN(out))
 
972
            oflags = NPY_FARRAY;
 
973
        else
 
974
            oflags = NPY_CARRAY;
 
975
        oflags |= NPY_UPDATEIFCOPY | NPY_FORCECAST;
 
976
        Py_INCREF(indescr);
 
977
        newout = (NPY_AO*)PyArray_FromArray(out, indescr, oflags);
 
978
        if (newout == NULL) {
 
979
            goto fail;
 
980
        }
 
981
    }
 
982
    else {
 
983
        newout = out;
 
984
        Py_INCREF(newout);
 
985
    }
 
986
 
 
987
    /* make sure the shape of the output array is the same */
 
988
    if (!PyArray_SAMESHAPE(newin, newout)) {
 
989
        PyErr_SetString(PyExc_ValueError, "clip: Output array must have the"
 
990
                        "same shape as the input.");
 
991
        goto fail;
 
992
    }
 
993
    if (newout->data != newin->data) {
 
994
        memcpy(newout->data, newin->data, PyArray_NBYTES(newin));
 
995
    }
 
996
 
 
997
    /* Now we can call the fast-clip function */
 
998
    min_data = max_data = NULL;
 
999
    if (mina != NULL) {
 
1000
        min_data = mina->data;
 
1001
    }
 
1002
    if (maxa != NULL) {
 
1003
        max_data = maxa->data;
 
1004
    }
 
1005
    func(newin->data, PyArray_SIZE(newin), min_data, max_data, newout->data);
 
1006
 
 
1007
    /* Clean up temporary variables */
 
1008
    Py_XDECREF(mina);
 
1009
    Py_XDECREF(maxa);
 
1010
    Py_DECREF(newin);
 
1011
    /* Copy back into out if out was not already a nice array. */
 
1012
    Py_DECREF(newout);
 
1013
    return (PyObject *)out;
 
1014
 
 
1015
 fail:
 
1016
    Py_XDECREF(maxa);
 
1017
    Py_XDECREF(mina);
 
1018
    Py_XDECREF(newin);
 
1019
    PyArray_XDECREF_ERR(newout);
 
1020
    return NULL;
 
1021
}
 
1022
 
 
1023
 
 
1024
/*NUMPY_API
 
1025
 * Conjugate
 
1026
 */
 
1027
NPY_NO_EXPORT PyObject *
 
1028
PyArray_Conjugate(PyArrayObject *self, PyArrayObject *out)
 
1029
{
 
1030
    if (PyArray_ISCOMPLEX(self)) {
 
1031
        if (out == NULL) {
 
1032
            return PyArray_GenericUnaryFunction(self,
 
1033
                                                n_ops.conjugate);
 
1034
        }
 
1035
        else {
 
1036
            return PyArray_GenericBinaryFunction(self,
 
1037
                                                 (PyObject *)out,
 
1038
                                                 n_ops.conjugate);
 
1039
        }
 
1040
    }
 
1041
    else {
 
1042
        PyArrayObject *ret;
 
1043
        if (out) {
 
1044
            if (PyArray_CopyAnyInto(out, self) < 0) {
 
1045
                return NULL;
 
1046
            }
 
1047
            ret = out;
 
1048
        }
 
1049
        else {
 
1050
            ret = self;
 
1051
        }
 
1052
        Py_INCREF(ret);
 
1053
        return (PyObject *)ret;
 
1054
    }
 
1055
}
 
1056
 
 
1057
/*NUMPY_API
 
1058
 * Trace
 
1059
 */
 
1060
NPY_NO_EXPORT PyObject *
 
1061
PyArray_Trace(PyArrayObject *self, int offset, int axis1, int axis2,
 
1062
              int rtype, PyArrayObject *out)
 
1063
{
 
1064
    PyObject *diag = NULL, *ret = NULL;
 
1065
 
 
1066
    diag = PyArray_Diagonal(self, offset, axis1, axis2);
 
1067
    if (diag == NULL) {
 
1068
        return NULL;
 
1069
    }
 
1070
    ret = PyArray_GenericReduceFunction((PyAO *)diag, n_ops.add, -1, rtype, out);
 
1071
    Py_DECREF(diag);
 
1072
    return ret;
 
1073
}
 
1074