← 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/npymath/npy_math_complex.c.src

  • 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/npymath/npy_math_complex.c.src
 
1
/*
 
2
 * Implement some C99-compatible complex math functions
 
3
 *
 
4
 * Most of the code is taken from the msun library in FreeBSD (HEAD @ 30th June
 
5
 * 2009), under the following license:
 
6
 *
 
7
 * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG>
 
8
 * All rights reserved.
 
9
 *
 
10
 * Redistribution and use in source and binary forms, with or without
 
11
 * modification, are permitted provided that the following conditions
 
12
 * are met:
 
13
 * 1. Redistributions of source code must retain the above copyright
 
14
 *    notice, this list of conditions and the following disclaimer.
 
15
 * 2. Redistributions in binary form must reproduce the above copyright
 
16
 *    notice, this list of conditions and the following disclaimer in the
 
17
 *    documentation and/or other materials provided with the distribution.
 
18
 *
 
19
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 
20
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
21
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
22
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 
23
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 
24
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 
25
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 
26
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 
27
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 
28
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 
29
 * SUCH DAMAGE.
 
30
 */
 
31
#include "npy_math_common.h"
 
32
#include "npy_math_private.h"
 
33
 
 
34
/*==========================================================
 
35
 * Custom implementation of missing complex C99 functions
 
36
 *=========================================================*/
 
37
 
 
38
/**begin repeat
 
39
 * #type = float,double,npy_longdouble#
 
40
 * #ctype = npy_cfloat,npy_cdouble,npy_clongdouble#
 
41
 * #c = f, , l#
 
42
 * #C = F, , L#
 
43
 * #TMAX = FLT_MAX, DBL_MAX, LDBL_MAX#
 
44
 */
 
45
#ifndef HAVE_CABS@C@
 
46
@type@ npy_cabs@c@(@ctype@ z)
 
47
{
 
48
    return npy_hypot@c@(npy_creal@c@(z), npy_cimag@c@(z));
 
49
}
 
50
#endif
 
51
 
 
52
#ifndef HAVE_CARG@C@
 
53
@type@ npy_carg@c@(@ctype@ z)
 
54
{
 
55
    return npy_atan2@c@(npy_cimag@c@(z), npy_creal@c@(z));
 
56
}
 
57
#endif
 
58
 
 
59
#ifndef HAVE_CEXP@C@
 
60
@ctype@ npy_cexp@c@(@ctype@ z)
 
61
{
 
62
    @type@ x, c, s;
 
63
    @type@ r, i;
 
64
    @ctype@ ret;
 
65
 
 
66
    r = npy_creal@c@(z);
 
67
    i = npy_cimag@c@(z);
 
68
 
 
69
    if (npy_isfinite(r)) {
 
70
        x = npy_exp@c@(r);
 
71
 
 
72
        c = npy_cos@c@(i);
 
73
        s = npy_sin@c@(i);
 
74
 
 
75
        if (npy_isfinite(i)) {
 
76
            ret = npy_cpack@c@(x * c, x * s);
 
77
        } else {
 
78
            ret = npy_cpack@c@(NPY_NAN, npy_copysign@c@(NPY_NAN, i));
 
79
        }
 
80
 
 
81
    } else  if (npy_isnan(r)) {
 
82
        /* r is nan */
 
83
        if (i == 0) {
 
84
            ret = npy_cpack@c@(r, 0);
 
85
        } else {
 
86
            ret = npy_cpack@c@(r, npy_copysign@c@(NPY_NAN, i));
 
87
        }
 
88
    } else {
 
89
        /* r is +- inf */
 
90
        if (r > 0) {
 
91
            if (i == 0) {
 
92
                ret = npy_cpack@c@(r, i);
 
93
            } else if (npy_isfinite(i)) {
 
94
                c = npy_cos@c@(i);
 
95
                s = npy_sin@c@(i);
 
96
 
 
97
                ret = npy_cpack@c@(r * c, r * s);
 
98
            } else {
 
99
                /* x = +inf, y = +-inf | nan */
 
100
                ret = npy_cpack@c@(r, NPY_NAN);
 
101
            }
 
102
        } else {
 
103
            if (npy_isfinite(i)) {
 
104
                x = npy_exp@c@(r);
 
105
                c = npy_cos@c@(i);
 
106
                s = npy_sin@c@(i);
 
107
 
 
108
                ret = npy_cpack@c@(x * c, x * s);
 
109
            } else {
 
110
                /* x = -inf, y = nan | +i inf */
 
111
                ret = npy_cpack@c@(0, 0);
 
112
            }
 
113
        }
 
114
    }
 
115
 
 
116
    return ret;
 
117
}
 
118
#endif
 
119
 
 
120
#ifndef HAVE_CLOG@C@
 
121
@ctype@ npy_clog@c@(@ctype@ z)
 
122
{
 
123
    return npy_cpack@c@(npy_log@c@ (npy_cabs@c@ (z)), npy_carg@c@ (z));
 
124
}
 
125
#endif
 
126
 
 
127
#ifndef HAVE_CSQRT@C@
 
128
 
 
129
/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
 
130
#define THRESH  (@TMAX@ / (1 + NPY_SQRT2@c@))
 
131
 
 
132
@ctype@ npy_csqrt@c@(@ctype@ z)
 
133
{
 
134
    @ctype@ result;
 
135
    @type@ a, b;
 
136
    @type@ t;
 
137
    int scale;
 
138
 
 
139
    a = npy_creal@c@(z);
 
140
    b = npy_cimag@c@(z);
 
141
 
 
142
    /* Handle special cases. */
 
143
    if (a == 0 && b == 0)
 
144
        return (npy_cpack@c@(0, b));
 
145
    if (npy_isinf(b))
 
146
        return (npy_cpack@c@(NPY_INFINITY, b));
 
147
    if (npy_isnan(a)) {
 
148
        t = (b - b) / (b - b);  /* raise invalid if b is not a NaN */
 
149
        return (npy_cpack@c@(a, t));    /* return NaN + NaN i */
 
150
    }
 
151
    if (npy_isinf(a)) {
 
152
        /*
 
153
         * csqrt(inf + NaN i)  = inf +  NaN i
 
154
         * csqrt(inf + y i)    = inf +  0 i
 
155
         * csqrt(-inf + NaN i) = NaN +- inf i
 
156
         * csqrt(-inf + y i)   = 0   +  inf i
 
157
         */
 
158
        if (npy_signbit(a))
 
159
            return (npy_cpack@c@(npy_fabs@c@(b - b), npy_copysign@c@(a, b)));
 
160
        else
 
161
            return (npy_cpack@c@(a, npy_copysign@c@(b - b, b)));
 
162
    }
 
163
    /*
 
164
     * The remaining special case (b is NaN) is handled just fine by
 
165
     * the normal code path below.
 
166
     */
 
167
 
 
168
    /* Scale to avoid overflow. */
 
169
    if (npy_fabs@c@(a) >= THRESH || npy_fabs@c@(b) >= THRESH) {
 
170
        a *= 0.25;
 
171
        b *= 0.25;
 
172
        scale = 1;
 
173
    } else {
 
174
        scale = 0;
 
175
    }
 
176
 
 
177
    /* Algorithm 312, CACM vol 10, Oct 1967. */
 
178
    if (a >= 0) {
 
179
        t = npy_sqrt@c@((a + npy_hypot@c@(a, b)) * 0.5);
 
180
        result = npy_cpack@c@(t, b / (2 * t));
 
181
    } else {
 
182
        t = npy_sqrt@c@((-a + npy_hypot@c@(a, b)) * 0.5);
 
183
        result = npy_cpack@c@(npy_fabs@c@(b) / (2 * t), npy_copysign@c@(t, b));
 
184
    }
 
185
 
 
186
    /* Rescale. */
 
187
    if (scale)
 
188
        return (npy_cpack@c@(npy_creal@c@(result) * 2, npy_cimag@c@(result)));
 
189
    else
 
190
        return (result);
 
191
}
 
192
#undef THRESH
 
193
#endif
 
194
 
 
195
#ifndef HAVE_CPOW@C@
 
196
@ctype@ npy_cpow@c@ (@ctype@ x, @ctype@ y)
 
197
{
 
198
    @ctype@ b;
 
199
    @type@ br, bi, yr, yi;
 
200
 
 
201
    yr = npy_creal@c@(y);
 
202
    yi = npy_cimag@c@(y);
 
203
    b = npy_clog@c@(x);
 
204
    br = npy_creal@c@(b);
 
205
    bi = npy_cimag@c@(b);
 
206
 
 
207
    return npy_cexp@c@(npy_cpack@c@(br * yr - bi * yi, br * yi + bi * yr));
 
208
}
 
209
#endif
 
210
 
 
211
#ifndef HAVE_CCOS@C@
 
212
@ctype@ npy_ccos@c@(@ctype@ z)
 
213
{
 
214
    @type@ x, y;
 
215
    x = npy_creal@c@(z);
 
216
    y = npy_cimag@c@(z);
 
217
    return npy_cpack@c@(npy_cos@c@(x) * npy_cosh@c@(y), -(npy_sin@c@(x) * npy_sinh@c@(y)));
 
218
}
 
219
#endif
 
220
 
 
221
#ifndef HAVE_CSIN@C@
 
222
@ctype@ npy_csin@c@(@ctype@ z)
 
223
{
 
224
    @type@ x, y;
 
225
    x = npy_creal@c@(z);
 
226
    y = npy_cimag@c@(z);
 
227
    return npy_cpack@c@(npy_sin@c@(x) * npy_cosh@c@(y), npy_cos@c@(x) * npy_sinh@c@(y));
 
228
}
 
229
#endif
 
230
/**end repeat**/
 
231
 
 
232
/*==========================================================
 
233
 * Decorate all the functions which are available natively
 
234
 *=========================================================*/
 
235
 
 
236
/**begin repeat
 
237
 * #type = float, double, npy_longdouble#
 
238
 * #ctype = npy_cfloat, npy_cdouble, npy_clongdouble#
 
239
 * #c = f, , l#
 
240
 * #C = F, , L#
 
241
 */
 
242
 
 
243
/**begin repeat1
 
244
 * #kind = cabs,carg#
 
245
 * #KIND = CABS,carg#
 
246
 */
 
247
#ifdef HAVE_@KIND@@C@
 
248
@type@ npy_@kind@@c@(@ctype@ z)
 
249
{
 
250
    __@ctype@_to_c99_cast z1 = {z};
 
251
    return @kind@@c@(z1.c99_z);
 
252
}
 
253
#endif
 
254
/**end repeat1**/
 
255
 
 
256
/**begin repeat1
 
257
 * #kind = cexp,clog,csqrt,ccos,csin#
 
258
 * #KIND = CEXP,CLOG,CSQRT,CCOS,CSIN#
 
259
 */
 
260
#ifdef HAVE_@KIND@@C@
 
261
@ctype@ npy_@kind@@c@(@ctype@ z)
 
262
{
 
263
    __@ctype@_to_c99_cast z1 = {z};
 
264
    __@ctype@_to_c99_cast ret;
 
265
    ret.c99_z = @kind@@c@(z1.c99_z);
 
266
    return ret.npy_z;
 
267
}
 
268
#endif
 
269
/**end repeat1**/
 
270
 
 
271
/**begin repeat1
 
272
 * #kind = cpow#
 
273
 * #KIND = CPOW#
 
274
 */
 
275
#ifdef HAVE_@KIND@@C@
 
276
@ctype@ npy_@kind@@c@(@ctype@ x, @ctype@ y)
 
277
{
 
278
    __@ctype@_to_c99_cast x1 = {x};
 
279
    __@ctype@_to_c99_cast y1 = {y};
 
280
    __@ctype@_to_c99_cast ret;
 
281
    ret.c99_z = @kind@@c@(x1.c99_z, y1.c99_z);
 
282
    return ret.npy_z;
 
283
}
 
284
#endif
 
285
/**end repeat1**/
 
286
 
 
287
/**end repeat**/