← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/sandbox/numexpr/complex_functions.inc

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/numexpr/complex_functions.inc
 
1
 
 
2
 
 
3
/* constants */
 
4
static cdouble nc_1 = {1., 0.};
 
5
static cdouble nc_half = {0.5, 0.};
 
6
static cdouble nc_i = {0., 1.};
 
7
static cdouble nc_i2 = {0., 0.5};
 
8
/*
 
9
static cdouble nc_mi = {0., -1.};
 
10
static cdouble nc_pi2 = {M_PI/2., 0.};
 
11
*/
 
12
 
 
13
static void
 
14
nc_sum(cdouble *a, cdouble *b, cdouble *r)
 
15
{
 
16
        r->real = a->real + b->real;
 
17
        r->imag = a->imag + b->imag;
 
18
        return;
 
19
}
 
20
 
 
21
static void
 
22
nc_diff(cdouble *a, cdouble *b, cdouble *r)
 
23
{
 
24
        r->real = a->real - b->real;
 
25
        r->imag = a->imag - b->imag;
 
26
        return;
 
27
}
 
28
 
 
29
static void
 
30
nc_neg(cdouble *a, cdouble *r)
 
31
{
 
32
        r->real = -a->real;
 
33
        r->imag = -a->imag;
 
34
        return;
 
35
}
 
36
 
 
37
static void
 
38
nc_prod(cdouble *a, cdouble *b, cdouble *r)
 
39
{
 
40
        double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
 
41
        r->real = ar*br - ai*bi;
 
42
        r->imag = ar*bi + ai*br;
 
43
        return;
 
44
}
 
45
 
 
46
static void
 
47
nc_quot(cdouble *a, cdouble *b, cdouble *r)
 
48
{
 
49
 
 
50
        double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
 
51
        double d = br*br + bi*bi;
 
52
        r->real = (ar*br + ai*bi)/d;
 
53
        r->imag = (ai*br - ar*bi)/d;
 
54
        return;
 
55
}
 
56
 
 
57
static void
 
58
nc_floor_quot(cdouble *a, cdouble *b, cdouble *r)
 
59
{
 
60
        double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
 
61
        double d = br*br + bi*bi;
 
62
        r->real = floor((ar*br + ai*bi)/d);
 
63
        r->imag = 0;
 
64
        return;
 
65
}
 
66
 
 
67
static void
 
68
nc_sqrt(cdouble *x, cdouble *r)
 
69
{
 
70
        double s,d;
 
71
        if (x->real == 0. && x->imag == 0.)
 
72
                *r = *x;
 
73
        else {
 
74
                s = sqrt(0.5*(fabs(x->real) + hypot(x->real,x->imag)));
 
75
                d = 0.5*x->imag/s;
 
76
                if (x->real > 0.) {
 
77
                        r->real = s;
 
78
                        r->imag = d;
 
79
                }
 
80
                else if (x->imag >= 0.) {
 
81
                        r->real = d;
 
82
                        r->imag = s;
 
83
                }
 
84
                else {
 
85
                        r->real = -d;
 
86
                        r->imag = -s;
 
87
                }
 
88
        }
 
89
        return;
 
90
}
 
91
 
 
92
static void
 
93
nc_log(cdouble *x, cdouble *r)
 
94
{
 
95
        double l = hypot(x->real,x->imag);
 
96
        r->imag = atan2(x->imag, x->real);
 
97
        r->real = log(l);
 
98
        return;
 
99
}
 
100
 
 
101
static void
 
102
nc_log1p(cdouble *x, cdouble *r)
 
103
{
 
104
        double l = hypot(x->real + 1.0,x->imag);
 
105
        r->imag = atan2(x->imag, x->real + 1.0);
 
106
        r->real = log(l);
 
107
        return;
 
108
}
 
109
 
 
110
static void
 
111
nc_exp(cdouble *x, cdouble *r)
 
112
{
 
113
        double a = exp(x->real);
 
114
        r->real = a*cos(x->imag);
 
115
        r->imag = a*sin(x->imag);
 
116
        return;
 
117
}
 
118
 
 
119
static void
 
120
nc_expm1(cdouble *x, cdouble *r)
 
121
{
 
122
        double a = exp(x->real);
 
123
        r->real = a*cos(x->imag) - 1.0;
 
124
        r->imag = a*sin(x->imag);
 
125
        return;
 
126
}
 
127
 
 
128
static void
 
129
nc_pow(cdouble *a, cdouble *b, cdouble *r)
 
130
{
 
131
        intp n;
 
132
        double ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
 
133
 
 
134
        if (br == 0. && bi == 0.) {
 
135
                r->real = 1.;
 
136
                r->imag = 0.;
 
137
                return;
 
138
        }
 
139
        if (ar == 0. && ai == 0.) {
 
140
                r->real = 0.;
 
141
                r->imag = 0.;
 
142
                return;
 
143
        }
 
144
        if (bi == 0 && (n=(intp)br) == br) {
 
145
            if (n > -100 && n < 100) {
 
146
                cdouble p, aa;
 
147
                intp mask = 1;
 
148
                if (n < 0) n = -n;
 
149
                aa = nc_1;
 
150
                p.real = ar; p.imag = ai;
 
151
                while (1) {
 
152
                        if (n & mask)
 
153
                            nc_prod(&aa,&p,&aa);
 
154
                        mask <<= 1;
 
155
                        if (n < mask || mask <= 0) break;
 
156
                        nc_prod(&p,&p,&p);
 
157
                }
 
158
                r->real = aa.real; r->imag = aa.imag;
 
159
                if (br < 0) nc_quot(&nc_1, r, r);
 
160
                return;
 
161
            }
 
162
        }
 
163
        /* complexobect.c uses an inline version of this formula
 
164
           investigate whether this had better performance or accuracy */
 
165
        nc_log(a, r);
 
166
        nc_prod(r, b, r);
 
167
        nc_exp(r, r);
 
168
        return;
 
169
}
 
170
 
 
171
 
 
172
static void
 
173
nc_prodi(cdouble *x, cdouble *r)
 
174
{
 
175
        r->real = -x->imag;
 
176
        r->imag = x->real;
 
177
        return;
 
178
}
 
179
 
 
180
 
 
181
static void
 
182
nc_acos(cdouble *x, cdouble *r)
 
183
{
 
184
        nc_prod(x,x,r);
 
185
        nc_diff(&nc_1, r, r);
 
186
        nc_sqrt(r, r);
 
187
        nc_prodi(r, r);
 
188
        nc_sum(x, r, r);
 
189
        nc_log(r, r);
 
190
        nc_prodi(r, r);
 
191
        nc_neg(r, r);
 
192
        return;
 
193
        /* return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i,
 
194
           nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))));
 
195
        */
 
196
}
 
197
 
 
198
static void
 
199
nc_acosh(cdouble *x, cdouble *r)
 
200
{
 
201
        nc_prod(x, x, r);
 
202
        nc_diff(&nc_1, r, r);
 
203
        nc_sqrt(r, r);
 
204
        nc_prodi(r, r);
 
205
        nc_sum(x, r, r);
 
206
        nc_log(r, r);
 
207
        return;
 
208
        /*
 
209
          return nc_log(nc_sum(x,nc_prod(nc_i,
 
210
          nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))));
 
211
        */
 
212
}
 
213
 
 
214
static void
 
215
nc_asin(cdouble *x, cdouble *r)
 
216
{
 
217
        cdouble a, *pa=&a;
 
218
        nc_prod(x, x, r);
 
219
        nc_diff(&nc_1, r, r);
 
220
        nc_sqrt(r, r);
 
221
        nc_prodi(x, pa);
 
222
        nc_sum(pa, r, r);
 
223
        nc_log(r, r);
 
224
        nc_prodi(r, r);
 
225
        nc_neg(r, r);
 
226
        return;
 
227
        /*
 
228
              return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x),
 
229
              nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))));
 
230
        */
 
231
}
 
232
 
 
233
 
 
234
static void
 
235
nc_asinh(cdouble *x, cdouble *r)
 
236
{
 
237
        nc_prod(x, x, r);
 
238
        nc_sum(&nc_1, r, r);
 
239
        nc_sqrt(r, r);
 
240
        nc_diff(r, x, r);
 
241
        nc_log(r, r);
 
242
        nc_neg(r, r);
 
243
        return;
 
244
        /*
 
245
          return nc_neg(nc_log(nc_diff(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)));
 
246
        */
 
247
}
 
248
 
 
249
static void
 
250
nc_atan(cdouble *x, cdouble *r)
 
251
{
 
252
        cdouble a, *pa=&a;
 
253
        nc_diff(&nc_i, x, pa);
 
254
        nc_sum(&nc_i, x, r);
 
255
        nc_quot(r, pa, r);
 
256
        nc_log(r,r);
 
257
        nc_prod(&nc_i2, r, r);
 
258
        return;
 
259
        /*
 
260
          return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
 
261
        */
 
262
}
 
263
 
 
264
static void
 
265
nc_atanh(cdouble *x, cdouble *r)
 
266
{
 
267
        cdouble a, *pa=&a;
 
268
        nc_diff(&nc_1, x, r);
 
269
        nc_sum(&nc_1, x, pa);
 
270
        nc_quot(pa, r, r);
 
271
        nc_log(r, r);
 
272
        nc_prod(&nc_half, r, r);
 
273
        return;
 
274
        /*
 
275
          return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
 
276
        */
 
277
}
 
278
 
 
279
static void
 
280
nc_cos(cdouble *x, cdouble *r)
 
281
{
 
282
        double xr=x->real, xi=x->imag;
 
283
        r->real = cos(xr)*cosh(xi);
 
284
        r->imag = -sin(xr)*sinh(xi);
 
285
        return;
 
286
}
 
287
 
 
288
static void
 
289
nc_cosh(cdouble *x, cdouble *r)
 
290
{
 
291
        double xr=x->real, xi=x->imag;
 
292
        r->real = cos(xi)*cosh(xr);
 
293
        r->imag = sin(xi)*sinh(xr);
 
294
        return;
 
295
}
 
296
 
 
297
 
 
298
#define M_LOG10_E 0.434294481903251827651128918916605082294397
 
299
 
 
300
static void
 
301
nc_log10(cdouble *x, cdouble *r)
 
302
{
 
303
        nc_log(x, r);
 
304
        r->real *= M_LOG10_E;
 
305
        r->imag *= M_LOG10_E;
 
306
        return;
 
307
}
 
308
 
 
309
static void
 
310
nc_sin(cdouble *x, cdouble *r)
 
311
{
 
312
        double xr=x->real, xi=x->imag;
 
313
        r->real = sin(xr)*cosh(xi);
 
314
        r->imag = cos(xr)*sinh(xi);
 
315
        return;
 
316
}
 
317
 
 
318
static void
 
319
nc_sinh(cdouble *x, cdouble *r)
 
320
{
 
321
        double xr=x->real, xi=x->imag;
 
322
        r->real = cos(xi)*sinh(xr);
 
323
        r->imag = sin(xi)*cosh(xr);
 
324
        return;
 
325
}
 
326
 
 
327
static void
 
328
nc_tan(cdouble *x, cdouble *r)
 
329
{
 
330
        double sr,cr,shi,chi;
 
331
        double rs,is,rc,ic;
 
332
        double d;
 
333
        double xr=x->real, xi=x->imag;
 
334
        sr = sin(xr);
 
335
        cr = cos(xr);
 
336
        shi = sinh(xi);
 
337
        chi = cosh(xi);
 
338
        rs = sr*chi;
 
339
        is = cr*shi;
 
340
        rc = cr*chi;
 
341
        ic = -sr*shi;
 
342
        d = rc*rc + ic*ic;
 
343
        r->real = (rs*rc+is*ic)/d;
 
344
        r->imag = (is*rc-rs*ic)/d;
 
345
        return;
 
346
}
 
347
 
 
348
static void
 
349
nc_tanh(cdouble *x, cdouble *r)
 
350
{
 
351
        double si,ci,shr,chr;
 
352
        double rs,is,rc,ic;
 
353
        double d;
 
354
        double xr=x->real, xi=x->imag;
 
355
        si = sin(xi);
 
356
        ci = cos(xi);
 
357
        shr = sinh(xr);
 
358
        chr = cosh(xr);
 
359
        rs = ci*shr;
 
360
        is = si*chr;
 
361
        rc = ci*chr;
 
362
        ic = si*shr;
 
363
        d = rc*rc + ic*ic;
 
364
        r->real = (rs*rc+is*ic)/d;
 
365
        r->imag = (is*rc-rs*ic)/d;
 
366
        return;
 
367
}