1 #include "schpriv.h"
2 #include <ctype.h>
3 #include <math.h>
4
5 #define zero scheme_exact_zero
6
make_complex(const Scheme_Object * r,const Scheme_Object * i,int normalize)7 static Scheme_Object *make_complex(const Scheme_Object *r, const Scheme_Object *i,
8 int normalize)
9 {
10 Scheme_Complex *c;
11
12 c = (Scheme_Complex *)scheme_malloc_small_dirty_tagged(sizeof(Scheme_Complex));
13 CLEAR_KEY_FIELD(&c->so);
14 c->so.type = scheme_complex_type;
15 c->r = (Scheme_Object *)r;
16 c->i = (Scheme_Object *)i;
17
18 if (normalize)
19 return scheme_complex_normalize((Scheme_Object *)c);
20 else
21 return (Scheme_Object *)c;
22 }
23
scheme_make_complex(const Scheme_Object * r,const Scheme_Object * i)24 Scheme_Object *scheme_make_complex(const Scheme_Object *r, const Scheme_Object *i)
25 {
26 return make_complex(r, i, 1);
27 }
28
scheme_real_to_complex(const Scheme_Object * n)29 Scheme_Object *scheme_real_to_complex(const Scheme_Object *n)
30 {
31 return make_complex(n, zero, 0);
32 }
33
scheme_make_small_complex(const Scheme_Object * n,Small_Complex * s)34 Scheme_Object *scheme_make_small_complex(const Scheme_Object *n, Small_Complex *s)
35 XFORM_SKIP_PROC
36 {
37 s->so.type = scheme_complex_type;
38 s->r = (Scheme_Object *)n;
39 s->i = zero;
40
41 return (Scheme_Object *)s;
42 }
43
scheme_is_complex_exact(const Scheme_Object * o)44 int scheme_is_complex_exact(const Scheme_Object *o)
45 {
46 Scheme_Complex *c = (Scheme_Complex *)o;
47
48 return !SCHEME_FLOATP(c->r) && !SCHEME_FLOATP(c->i);
49 }
50
scheme_complex_normalize(const Scheme_Object * o)51 Scheme_Object *scheme_complex_normalize(const Scheme_Object *o)
52 {
53 Scheme_Complex *c = (Scheme_Complex *)o;
54
55 if (c->i == zero)
56 return c->r;
57 if (c->r == zero) {
58 /* No coercions */
59 return (Scheme_Object *)c;
60 }
61
62 /* Coercions: Exact -> float -> double
63 If the complex contains a float and an exact, we coerce the exact
64 to a float, etc. */
65
66 #ifdef MZ_USE_SINGLE_FLOATS
67 if (SCHEME_FLTP(c->i)) {
68 if (!SCHEME_FLTP(c->r)) {
69 Scheme_Object *v;
70 if (SCHEME_DBLP(c->r)) {
71 v = scheme_make_double(SCHEME_FLT_VAL(c->i));
72 c->i = v;
73 } else {
74 v = scheme_make_float(scheme_get_val_as_float(c->r));
75 c->r = v;
76 }
77 }
78 } else if (SCHEME_FLTP(c->r)) {
79 Scheme_Object *v;
80 /* Imag part can't be a float, or we'd be in the previous case */
81 if (SCHEME_DBLP(c->i)) {
82 v = scheme_make_double(SCHEME_FLT_VAL(c->r));
83 c->r = v;
84 } else {
85 v = scheme_make_float(scheme_get_val_as_float(c->i));
86 c->i = v;
87 }
88 } else
89 #endif
90
91 if (SCHEME_DBLP(c->i)) {
92 if (!SCHEME_DBLP(c->r)) {
93 Scheme_Object *r;
94 r = scheme_make_double(scheme_get_val_as_double(c->r));
95 c->r = r;
96 }
97 } else if (SCHEME_DBLP(c->r)) {
98 Scheme_Object *i;
99 i = scheme_make_double(scheme_get_val_as_double(c->i));
100 c->i = i;
101 }
102
103 return (Scheme_Object *)c;
104 }
105
scheme_complex_real_part(const Scheme_Object * n)106 Scheme_Object *scheme_complex_real_part(const Scheme_Object *n)
107 {
108 return ((Scheme_Complex *)n)->r;
109 }
110
scheme_complex_imaginary_part(const Scheme_Object * n)111 Scheme_Object *scheme_complex_imaginary_part(const Scheme_Object *n)
112 {
113 return ((Scheme_Complex *)n)->i;
114 }
115
scheme_complex_eq(const Scheme_Object * a,const Scheme_Object * b)116 int scheme_complex_eq(const Scheme_Object *a, const Scheme_Object *b)
117 {
118 Scheme_Complex *ca = (Scheme_Complex *)a;
119 Scheme_Complex *cb = (Scheme_Complex *)b;
120 return scheme_bin_eq(ca->r, cb->r) && scheme_bin_eq(ca->i, cb->i);
121 }
122
scheme_complex_negate(const Scheme_Object * o)123 Scheme_Object *scheme_complex_negate(const Scheme_Object *o)
124 {
125 Scheme_Complex *c = (Scheme_Complex *)o;
126
127 return make_complex(scheme_bin_minus(scheme_make_integer(0),
128 c->r),
129 scheme_bin_minus(scheme_make_integer(0),
130 c->i),
131 0);
132 }
133
scheme_complex_add(const Scheme_Object * a,const Scheme_Object * b)134 Scheme_Object *scheme_complex_add(const Scheme_Object *a, const Scheme_Object *b)
135 {
136 Scheme_Complex *ca = (Scheme_Complex *)a;
137 Scheme_Complex *cb = (Scheme_Complex *)b;
138
139 return scheme_make_complex(scheme_bin_plus(ca->r, cb->r),
140 scheme_bin_plus(ca->i, cb->i));
141 }
142
scheme_complex_subtract(const Scheme_Object * a,const Scheme_Object * b)143 Scheme_Object *scheme_complex_subtract(const Scheme_Object *a, const Scheme_Object *b)
144 {
145 Scheme_Complex *ca = (Scheme_Complex *)a;
146 Scheme_Complex *cb = (Scheme_Complex *)b;
147
148 return scheme_make_complex(scheme_bin_minus(ca->r, cb->r),
149 scheme_bin_minus(ca->i, cb->i));
150 }
151
scheme_complex_add1(const Scheme_Object * n)152 Scheme_Object *scheme_complex_add1(const Scheme_Object *n)
153 {
154 Small_Complex s;
155
156 return scheme_complex_add(scheme_make_small_complex(scheme_make_integer(1), &s),
157 n);
158 }
159
scheme_complex_sub1(const Scheme_Object * n)160 Scheme_Object *scheme_complex_sub1(const Scheme_Object *n)
161 {
162 Small_Complex s;
163
164 return scheme_complex_add(n, scheme_make_small_complex(scheme_make_integer(-1),
165 &s));
166 }
167
scheme_complex_multiply(const Scheme_Object * a,const Scheme_Object * b)168 Scheme_Object *scheme_complex_multiply(const Scheme_Object *a, const Scheme_Object *b)
169 {
170 Scheme_Complex *ca = (Scheme_Complex *)a;
171 Scheme_Complex *cb = (Scheme_Complex *)b;
172
173 return scheme_make_complex(scheme_bin_minus(scheme_bin_mult(ca->r, cb->r),
174 scheme_bin_mult(ca->i, cb->i)),
175 scheme_bin_plus(scheme_bin_mult(ca->r, cb->i),
176 scheme_bin_mult(ca->i, cb->r)));
177
178 }
179
simple_complex_divide(Scheme_Object * a,Scheme_Object * b,Scheme_Object * c,Scheme_Object * d,int swap)180 static Scheme_Object *simple_complex_divide(Scheme_Object *a, Scheme_Object *b,
181 Scheme_Object *c, Scheme_Object *d,
182 int swap)
183 {
184 Scheme_Object *r, *i, *cm, *cb, *da, *ci;
185
186 cm = scheme_bin_plus(scheme_bin_mult(c, c),
187 scheme_bin_mult(d, d));
188
189 r = scheme_bin_div(scheme_bin_plus(scheme_bin_mult(c, a),
190 scheme_bin_mult(d, b)),
191 cm);
192
193 cb = scheme_bin_mult(c, b);
194 da = scheme_bin_mult(d, a);
195 if (swap)
196 ci = scheme_bin_minus(da, cb);
197 else
198 ci = scheme_bin_minus(cb, da);
199 i = scheme_bin_div(ci, cm);
200
201 return scheme_make_complex(r, i);
202 }
203
scheme_complex_divide(const Scheme_Object * _n,const Scheme_Object * _d)204 Scheme_Object *scheme_complex_divide(const Scheme_Object *_n, const Scheme_Object *_d)
205 {
206 Scheme_Complex *cn = (Scheme_Complex *)_n;
207 Scheme_Complex *cd = (Scheme_Complex *)_d;
208 Scheme_Object *den, *r, *i, *a, *b, *c, *d, *cm, *dm, *aa[1];
209 int swap;
210
211 if ((cn->r == zero) && (cn->i == zero))
212 return zero;
213
214 a = cn->r;
215 b = cn->i;
216 c = cd->r;
217 d = cd->i;
218
219 /* Check for exact-zero simplifications in d: */
220 if (c == zero) {
221 i = scheme_bin_minus(zero, scheme_bin_div(a, d));
222 r = scheme_bin_div(b, d);
223 return scheme_make_complex(r, i);
224 } else if (d == zero) {
225 r = scheme_bin_div(a, c);
226 i = scheme_bin_div(b, c);
227 return scheme_make_complex(r, i);
228 }
229
230 if (b == zero) {
231 /* As in Chez Scheme: a / c+di => c(a/(cc+dd)) + (-d(a/cc+dd))i */
232 cm = scheme_bin_div(a, scheme_bin_plus(scheme_bin_mult(c, c), scheme_bin_mult(d, d)));
233 return scheme_make_complex(scheme_bin_mult(c, cm),
234 scheme_bin_minus(zero, scheme_bin_mult(d, cm)));
235 }
236
237 if (!SCHEME_FLOATP(a) && !SCHEME_FLOATP(b) && !SCHEME_FLOATP(c) && !SCHEME_FLOATP(d))
238 return simple_complex_divide(a, b, c, d, 0);
239
240 aa[0] = c;
241 cm = scheme_abs(1, aa);
242 aa[0] = d;
243 dm = scheme_abs(1, aa);
244
245 if (scheme_bin_lt(cm, dm)) {
246 cm = a;
247 a = b;
248 b = cm;
249 cm = c;
250 c = d;
251 d = cm;
252 swap = 1;
253 } else
254 swap = 0;
255
256 r = scheme_bin_div(c, d);
257
258 if (!SCHEME_FLOATP(r) && (SCHEME_FLOATP(a) || SCHEME_FLOATP(b))) {
259 aa[0] = r;
260 r = scheme_exact_to_inexact(1, aa);
261 }
262
263 /* If r goes to infinity, try computing a different way to avoid overflow: */
264 if (SCHEME_FLOATP(r)) {
265 double v = SCHEME_FLOAT_VAL(r);
266 if (MZ_IS_POS_INFINITY(v) || MZ_IS_NEG_INFINITY(v)) {
267 /* This calculuation does not work as well for complex numbers with
268 large parts, such as `(/ 1e+300+1e+300i 4e+300+4e+300i)`, but it
269 works better for small parts, as in `(/ 0.0+0.0i 1+1e-320i)`. */
270 return simple_complex_divide(a, b, c, d, swap);
271 }
272 }
273
274 den = scheme_bin_plus(d, scheme_bin_mult(c, r));
275
276 if (swap)
277 i = scheme_bin_div(scheme_bin_minus(a, scheme_bin_mult(b, r)), den);
278 else
279 i = scheme_bin_div(scheme_bin_minus(scheme_bin_mult(b, r), a), den);
280
281 r = scheme_bin_div(scheme_bin_plus(b, scheme_bin_mult(a, r)), den);
282
283 return scheme_make_complex(r, i);
284 }
285
scheme_complex_power(const Scheme_Object * base,const Scheme_Object * exponent)286 Scheme_Object *scheme_complex_power(const Scheme_Object *base, const Scheme_Object *exponent)
287 {
288 Scheme_Complex *cb = (Scheme_Complex *)base;
289 Scheme_Complex *ce = (Scheme_Complex *)exponent;
290 double a, b, c, d, bm, ba, nm, na, r1, r2;
291 int d_is_zero;
292
293 if ((ce->i == zero) && !SCHEME_FLOATP(ce->r)) {
294 if (SCHEME_INTP(ce->r) || SCHEME_BIGNUMP(ce->r))
295 return scheme_generic_integer_power(base, ce->r);
296 }
297
298 a = scheme_get_val_as_double(cb->r);
299 b = scheme_get_val_as_double(cb->i);
300 c = scheme_get_val_as_double(ce->r);
301 d = scheme_get_val_as_double(ce->i);
302 d_is_zero = (ce->i == zero);
303
304 bm = sqrt(a * a + b * b);
305 ba = atan2(b, a);
306
307 /* New mag & angle */
308 nm = scheme_double_expt(bm, c) * exp(-(ba * d));
309 if (d_is_zero) /* precision here can avoid NaNs */
310 na = ba * c;
311 else
312 na = log(bm) * d + ba * c;
313
314 r1 = nm * cos(na);
315 r2 = nm * sin(na);
316
317 #ifdef MZ_USE_SINGLE_FLOATS
318 /* Coerce to double or float? */
319 if (!SCHEME_DBLP(cb->r) && !SCHEME_DBLP(cb->i)
320 && !SCHEME_DBLP(ce->r) && !SCHEME_DBLP(ce->i))
321 #ifndef USE_SINGLE_FLOATS_AS_DEFAULT
322 if (SCHEME_FLTP(cb->r) || SCHEME_FLTP(cb->i)
323 || SCHEME_FLTP(ce->r) || SCHEME_FLTP(ce->i))
324 #endif
325 return scheme_make_complex(scheme_make_float((float)r1),
326 scheme_make_float((float)r2));
327 #endif
328
329 return scheme_make_complex(scheme_make_double(r1),
330 scheme_make_double(r2));
331 }
332
scheme_complex_sqrt(const Scheme_Object * o)333 Scheme_Object *scheme_complex_sqrt(const Scheme_Object *o)
334 {
335 Scheme_Complex *c = (Scheme_Complex *)o;
336 Scheme_Object *r, *i, *ssq, *srssq, *nrsq, *prsq, *nr, *ni;
337
338 r = c->r;
339 i = c->i;
340
341 if (scheme_is_zero(i)) {
342 /* Special case for x+0.0i and x-0.0i: */
343 r = scheme_sqrt(1, &r);
344 if (!SCHEME_COMPLEXP(r))
345 return scheme_make_complex(r, i);
346 else {
347 c = (Scheme_Complex *)r;
348 if (SAME_OBJ(c->r, zero)) {
349 /* need an inexact-zero real part: */
350 #ifdef MZ_USE_SINGLE_FLOATS
351 if (SCHEME_FLTP(c->i))
352 r = scheme_make_float(0.0);
353 else
354 #endif
355 r = scheme_make_double(0.0);
356 if (scheme_minus_zero_p(scheme_real_to_double(i))) {
357 /* we started with x-0.0i */
358 return scheme_make_complex(r, scheme_bin_minus(scheme_make_integer(0), c->i));
359 }
360 else
361 return scheme_make_complex(r, c->i);
362 } else
363 return r;
364 }
365 }
366
367 ssq = scheme_bin_plus(scheme_bin_mult(r, r),
368 scheme_bin_mult(i, i));
369
370 srssq = scheme_sqrt(1, &ssq);
371
372 if (SCHEME_FLOATP(srssq)) {
373 /* We may have lost too much precision, if i << r. The result is
374 going to be inexact, anyway, so switch to using expt. */
375 Scheme_Object *a[2], *p;
376 a[0] = (Scheme_Object *)o;
377 #ifdef MZ_USE_SINGLE_FLOATS
378 if (SCHEME_FLTP(c->i))
379 p = scheme_make_float(0.5);
380 else
381 #endif
382 p = scheme_make_double(0.5);
383 a[1] = p;
384 return scheme_expt(2, a);
385 }
386
387 nrsq = scheme_bin_div(scheme_bin_minus(srssq, r),
388 scheme_make_integer(2));
389
390 nr = scheme_sqrt(1, &nrsq);
391 if (scheme_is_negative(i))
392 nr = scheme_bin_minus(zero, nr);
393
394 prsq = scheme_bin_div(scheme_bin_plus(srssq, r),
395 scheme_make_integer(2));
396
397 ni = scheme_sqrt(1, &prsq);
398
399 return scheme_make_complex(ni, nr);
400 }
401
scheme_complex_atan(const Scheme_Object * c)402 Scheme_Object *scheme_complex_atan(const Scheme_Object *c)
403 {
404 /* From Chez Scheme, which implements the principal expression from Kahan:
405 (log(1+z) - log(1-z))/2
406 */
407 # define OMEGA 1.7976931348623157e308
408 # define THETA 3.351951982485649e+153 /* = sqrt(OMEGA)/4 */
409 # define RHO 2.9833362924800834e-154 /* = 1/THETA */
410 # define HALF_PI 1.5707963267948966
411 # define IS_NEG(x) (((x) < 0.0) || ((x == 0.0) && scheme_minus_zero_p(x)))
412 double x, y, ay, r, i;
413 int negate;
414
415 /* Get components after multiplication by i */
416 x = -scheme_real_to_double(_scheme_complex_imaginary_part(c));
417 y = scheme_real_to_double(_scheme_complex_real_part(c));
418
419 /* Compute atanh */
420
421 if (x < 0.0) {
422 x = -x;
423 negate = 1;
424 } else {
425 y = -y;
426 negate = 0;
427 }
428
429 ay = fabs(y);
430
431 if ((x > THETA) || (y > THETA)) {
432 /* RP(1/z) +/- (pi/2)i */
433 if (x > ay)
434 r = 1/(x + ((y/x) * y));
435 else if (x < ay) {
436 r = y/x;
437 r = r/((x * r)+ y);
438 } else
439 r = 1/(x+ay);
440
441 i = (IS_NEG(y) ? HALF_PI : (-HALF_PI));
442 } else if (x == 1.0) {
443 double k = ay + RHO;
444 r = scheme_double_log(sqrt(sqrt((y * y) + 4.0)) / sqrt(k));
445 i = (HALF_PI + scheme_double_atan(k/2.0)) / (IS_NEG(y) ? 2.0 : -2.0);
446 } else {
447 double mx = 1.0 - x;
448 double k = ay + RHO;
449 k = k * k;
450
451 r = scheme_double_log(((4.0 * x) / ((mx * mx) + k)) + 1.0) / 4.0;
452 i = scheme_double_atan2(2.0 * y, (mx * (1.0 + x)) - k) / -2.0;
453 }
454
455 if (negate) {
456 i = -i;
457 r = -r;
458 }
459
460 /* Multiply by -i to get atan */
461 x = i;
462 y = -r;
463
464 #ifdef MZ_USE_SINGLE_FLOATS
465 if (SCHEME_FLTP(_scheme_complex_real_part(c))
466 || SCHEME_FLTP(_scheme_complex_imaginary_part(c))) {
467 return scheme_make_complex(scheme_make_float(x), scheme_make_float(y));
468 }
469 #endif
470
471 return scheme_make_complex(scheme_make_double(x), scheme_make_double(y));
472 }
473
scheme_complex_asin_or_acos(const Scheme_Object * z,int get_asin)474 Scheme_Object *scheme_complex_asin_or_acos(const Scheme_Object *z, int get_asin)
475 {
476 /* From Chez Scheme, which implements the principal expression from Kahan */
477 Scheme_Object *zp, *zm, *aa[1];
478 double a, b, c, d, r, i;
479
480 aa[0] = scheme_bin_minus(scheme_make_integer(1), z);
481 zm = scheme_sqrt(1, aa);
482
483 aa[0] = scheme_bin_plus(scheme_make_integer(1), z);
484 zp = scheme_sqrt(1, aa);
485
486 if (SCHEME_COMPLEXP(zm)) {
487 a = scheme_real_to_double(_scheme_complex_real_part(zm));
488 b = scheme_real_to_double(_scheme_complex_imaginary_part(zm));
489 } else {
490 a = scheme_real_to_double(zm);
491 b = 0.0;
492 }
493
494 if (SCHEME_COMPLEXP(zp)) {
495 c = scheme_real_to_double(_scheme_complex_real_part(zp));
496 d = scheme_real_to_double(_scheme_complex_imaginary_part(zp));
497 } else {
498 c = scheme_real_to_double(zp);
499 d = 0.0;
500 }
501
502 if (get_asin) {
503 if (SCHEME_COMPLEXP(z)) {
504 r = scheme_real_to_double(_scheme_complex_real_part(z));
505 r = scheme_double_atan2(r, (a*c)-(b*d));
506 } else {
507 r = scheme_real_to_double((Scheme_Object *)z);
508 r = scheme_double_atan2(r, 0.0); /* void +nan.0 from (a*c)-(b*d) */
509 }
510
511 i = asinh((a*d)-(b*c));
512 } else {
513 r = 2.0 * scheme_double_atan2(a, c);
514 i = asinh((b*c) - (a*d));
515 }
516
517 #ifdef MZ_USE_SINGLE_FLOATS
518 if (SCHEME_FLTP(_scheme_complex_real_part(z))
519 || SCHEME_FLTP(_scheme_complex_imaginary_part(z))) {
520 return scheme_make_complex(scheme_make_float(r), scheme_make_float(i));
521 }
522 #endif
523
524 return scheme_make_complex(scheme_make_double(r), scheme_make_double(i));
525 }
526
scheme_complex_asin(const Scheme_Object * c)527 Scheme_Object *scheme_complex_asin(const Scheme_Object *c)
528 {
529 return scheme_complex_asin_or_acos(c, 1);
530 }
531
scheme_complex_acos(const Scheme_Object * c)532 Scheme_Object *scheme_complex_acos(const Scheme_Object *c)
533 {
534 return scheme_complex_asin_or_acos(c, 0);
535 }
536