1 /* Complex Data Type definition for S-Lang */
2 /*
3 Copyright (C) 2004-2017,2018 John E. Davis
4
5 This file is part of the S-Lang Library.
6
7 The S-Lang Library is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 2 of the
10 License, or (at your option) any later version.
11
12 The S-Lang Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this library; if not, write to the Free Software
19 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
20 USA.
21 */
22
23 #include "slinclud.h"
24
25 #include "slang.h"
26 #include "_slang.h"
27
28 /* The rest of the file is enclosed in this #if */
29 #if SLANG_HAS_COMPLEX
30
31 #if SLANG_HAS_FLOAT
32 # include <math.h>
33 #endif
34
35 #ifdef PI
36 # undef PI
37 #endif
38 #define PI 3.14159265358979323846
39
SLang_pop_complex(double * r,double * i)40 int SLang_pop_complex (double *r, double *i)
41 {
42 double *c;
43
44 switch (SLang_peek_at_stack ())
45 {
46 case SLANG_COMPLEX_TYPE:
47 if (-1 == SLclass_pop_ptr_obj (SLANG_COMPLEX_TYPE, VOID_STAR_STAR(&c)))
48 return -1;
49 *r = c[0];
50 *i = c[1];
51 SLfree ((char *) c);
52 break;
53
54 default:
55 *i = 0.0;
56 if (-1 == SLang_pop_double (r))
57 return -1;
58 break;
59
60 case -1:
61 return -1;
62 }
63 return 0;
64 }
65
SLang_push_complex(double r,double i)66 int SLang_push_complex (double r, double i)
67 {
68 double *c;
69
70 c = (double *) SLmalloc (2 * sizeof (double));
71 if (c == NULL)
72 return -1;
73
74 c[0] = r;
75 c[1] = i;
76
77 if (-1 == SLclass_push_ptr_obj (SLANG_COMPLEX_TYPE, (VOID_STAR) c))
78 {
79 SLfree ((char *) c);
80 return -1;
81 }
82 return 0;
83 }
84
SLcomplex_times(double * c,double * a,double * b)85 double *SLcomplex_times (double *c, double *a, double *b)
86 {
87 double a_real, b_real, a_imag, b_imag;
88
89 a_real = a[0];
90 b_real = b[0];
91 a_imag = a[1];
92 b_imag = b[1];
93
94 c[0] = a_real * b_real - a_imag * b_imag;
95 c[1] = a_imag * b_real + a_real * b_imag;
96
97 return c;
98 }
99
SLcomplex_divide(double * c,double * a,double * b)100 double *SLcomplex_divide (double *c, double *a, double *b)
101 {
102 double a_real, b_real, a_imag, b_imag;
103 double ratio, invden;
104
105 a_real = a[0];
106 b_real = b[0];
107 a_imag = a[1];
108 b_imag = b[1];
109
110 /* Do it this way to avoid overflow in the denom */
111 if (fabs(b_real) > fabs(b_imag))
112 {
113 ratio = b_imag / b_real;
114 invden = 1.0 / (b_real + b_imag * ratio);
115 c[0] = (a_real + ratio * a_imag) * invden;
116 c[1] = (a_imag - a_real * ratio) * invden;
117 }
118 else
119 {
120 ratio = b_real / b_imag;
121 invden = 1.0 / (b_real * ratio + b_imag);
122 c[0] = (a_real * ratio + a_imag) * invden;
123 c[1] = (a_imag * ratio - a_real) * invden;
124 }
125 return c;
126 }
127
128 /* a^b = exp (b log a); */
SLcomplex_pow(double * c,double * a,double * b)129 double *SLcomplex_pow (double *c, double *a, double *b)
130 {
131 if ((a[0] == 0.0) && (b[0] == 0.0) && (a[1] == 0.0) && (b[1] == 0.0))
132 {
133 c[0] = 1.0;
134 c[1] = 0.0;
135 return c;
136 }
137 return SLcomplex_exp (c, SLcomplex_times (c, b, SLcomplex_log (c, a)));
138 }
139
complex_dpow(double * c,double * a,double b)140 static double *complex_dpow (double *c, double *a, double b)
141 {
142 if ((b == 0.0) && (a[0] == 0.0) && (a[1] == 0.0))
143 {
144 c[0] = 1.0;
145 c[1] = 0.0;
146 return c;
147 }
148 SLcomplex_log (c, a);
149 c[0] *= b;
150 c[1] *= b;
151 return SLcomplex_exp (c, c);
152 }
153
dcomplex_pow(double * c,double a,double * b)154 static double *dcomplex_pow (double *c, double a, double *b)
155 {
156 if ((a == 0.0) && (b[0] == 0.0) && (b[1] == 0.0))
157 {
158 c[0] = 1.0;
159 c[1] = 0.0;
160 return c;
161 }
162 a = log (a);
163 c[0] = a * b[0];
164 c[1] = a * b[1];
165 return SLcomplex_exp (c, c);
166 }
167
SLcomplex_abs(double * z)168 double SLcomplex_abs (double *z)
169 {
170 return SLmath_hypot (z[0], z[1]);
171 }
172
173 /* It appears that FORTRAN assumes that the branch cut for the log function
174 * is along the -x axis. So, use this for atan2:
175 */
my_atan2(double y,double x)176 static double my_atan2 (double y, double x)
177 {
178 double val;
179
180 val = atan (y/x);
181
182 if (x >= 0)
183 return val; /* I, IV */
184
185 if (y <= 0) /* III */
186 return val - PI;
187
188 return PI + val; /* II */
189 }
190
polar_form(double * r,double * theta,double * z)191 static void polar_form (double *r, double *theta, double *z)
192 {
193 double x, y;
194
195 *r = SLcomplex_abs (z);
196
197 x = z[0];
198 y = z[1];
199
200 if (x == 0.0)
201 {
202 if (y >= 0)
203 *theta = 0.5 * PI;
204 else
205 *theta = 1.5 * PI;
206 }
207 else *theta = my_atan2 (y, x);
208 }
209
SLcomplex_sin(double * sinz,double * z)210 double *SLcomplex_sin (double *sinz, double *z)
211 {
212 double x, y;
213
214 x = z[0]; y = z[1];
215 sinz[0] = sin (x) * cosh (y);
216 sinz[1] = cos (x) * sinh (y);
217 return sinz;
218 }
219
SLcomplex_cos(double * cosz,double * z)220 double *SLcomplex_cos (double *cosz, double *z)
221 {
222 double x, y;
223
224 x = z[0]; y = z[1];
225 cosz[0] = cos (x) * cosh (y);
226 cosz[1] = -sin (x) * sinh (y);
227 return cosz;
228 }
229
_pSLcomplex_expm1(double * f,double * z)230 double *_pSLcomplex_expm1 (double *f, double *z)
231 {
232 /* exp(z)-1 = (e^x-1) + (cos(y)-1)*e^x + ie^x*sin(y) */
233 double ex;
234 double x = z[0], y = z[1];
235 double siny2 = sin(0.5*y);
236
237 ex = exp(x);
238 f[0] = EXPM1_FUNC(x) - ex*(2.0*siny2*siny2); /* cos(y)-1 = 2*sin^2(y/2) */
239 f[1] = ex*sin(y);
240 return f;
241 }
242
SLcomplex_exp(double * expz,double * z)243 double *SLcomplex_exp (double *expz, double *z)
244 {
245 double r, i;
246
247 r = exp (z[0]);
248 i = z[1];
249 expz[0] = r * cos (i);
250 expz[1] = r * sin (i);
251 return expz;
252 }
253
SLcomplex_log(double * logz,double * z)254 double *SLcomplex_log (double *logz, double *z)
255 {
256 double r, theta;
257
258 polar_form (&r, &theta, z); /* log R.e^(ix) = log R + ix */
259 logz[0] = log(r);
260 logz[1] = theta;
261 return logz;
262 }
263
_pSLcomplex_log1p(double * f,double * z)264 double *_pSLcomplex_log1p (double *f, double *z)
265 {
266 /* f(z) = log(1+z)
267 * = log(1+x+iy)
268 * = log(sqrt((1+x)^2 + y^2)) + i*theta
269 * = 0.5*log((1+x)^2 + y^2) + i*theta
270 * = 0.5*log(1 + 2*x + x^2 + y^2) + i*theta
271 * = 0.5*log(1 + r^2 + 2*x) + i*theta
272 */
273 double r, theta;
274
275 polar_form (&r, &theta, z);
276 if (r < 1)
277 f[0] = 0.5*LOG1P_FUNC(r*r + 2.0*z[0]);
278 else
279 f[0] = log(SLmath_hypot(1.0+z[0], z[1]));
280
281 f[1] = theta;
282 return f;
283 }
284
SLcomplex_log10(double * log10z,double * z)285 double *SLcomplex_log10 (double *log10z, double *z)
286 {
287 double l10 = log (10.0);
288 (void) SLcomplex_log (log10z, z);
289 log10z[0] = log10z[0] / l10;
290 log10z[1] = log10z[1] / l10;
291 return log10z;
292 }
293
SLcomplex_sqrt(double * sqrtz,double * z)294 double *SLcomplex_sqrt (double *sqrtz, double *z)
295 {
296 double r, x, y;
297
298 x = z[0];
299 y = z[1];
300
301 r = SLmath_hypot (x, y);
302
303 if (r == 0.0)
304 {
305 sqrtz [0] = sqrtz [1] = 0.0;
306 return sqrtz;
307 }
308
309 if (x >= 0.0)
310 {
311 x = sqrt (0.5 * (r + x));
312 y = 0.5 * y / x;
313 }
314 else
315 {
316 r = sqrt (0.5 * (r - x));
317 x = 0.5 * y / r;
318 y = r;
319
320 if (x < 0.0)
321 {
322 x = -x;
323 y = -y;
324 }
325 }
326
327 sqrtz[0] = x;
328 sqrtz[1] = y;
329
330 return sqrtz;
331 }
332
SLcomplex_tan(double * tanz,double * z)333 double *SLcomplex_tan (double *tanz, double *z)
334 {
335 double x, y, invden;
336
337 x = 2 * z[0];
338 y = 2 * z[1];
339 invden = 1.0 / (cos (x) + cosh (y));
340 tanz[0] = invden * sin (x);
341 tanz[1] = invden * sinh (y);
342 return tanz;
343 }
344
345 /* Utility Function */
compute_alpha_beta(double * z,double * alpha,double * beta)346 static void compute_alpha_beta (double *z, double *alpha, double *beta)
347 {
348 double x, y, a, b;
349
350 x = z[0];
351 y = z[1];
352 a = 0.5 * SLmath_hypot (x + 1, y);
353 b = 0.5 * SLmath_hypot (x - 1, y);
354
355 *alpha = a + b;
356 *beta = a - b;
357 }
358
SLcomplex_asin(double * asinz,double * z)359 double *SLcomplex_asin (double *asinz, double *z)
360 {
361 double alpha, beta;
362
363 compute_alpha_beta (z, &alpha, &beta);
364 asinz[0] = asin (beta);
365 asinz[1] = log (alpha + sqrt (alpha * alpha - 1));
366 return asinz;
367 }
368
SLcomplex_acos(double * acosz,double * z)369 double *SLcomplex_acos (double *acosz, double *z)
370 {
371 double alpha, beta;
372
373 compute_alpha_beta (z, &alpha, &beta);
374 acosz[0] = acos (beta);
375 acosz[1] = -log (alpha + sqrt (alpha * alpha - 1));
376 return acosz;
377 }
378
SLcomplex_atan(double * atanz,double * z)379 double *SLcomplex_atan (double *atanz, double *z)
380 {
381 double x, y;
382 double z1[2], z2[2];
383
384 x = z[0]; y = z[1];
385 z1[0] = x;
386 z1[1] = 1 + y;
387 z2[0] = -x;
388 z2[1] = 1 - y;
389
390 SLcomplex_log (z1, SLcomplex_divide (z2, z1, z2));
391 atanz[0] = -0.5 * z1[1];
392 atanz[1] = 0.5 * z1[0];
393
394 return atanz;
395 }
396
SLcomplex_sinh(double * sinhz,double * z)397 double *SLcomplex_sinh (double *sinhz, double *z)
398 {
399 double x, y;
400 x = z[0]; y = z[1];
401 sinhz[0] = sinh (x) * cos (y);
402 sinhz[1] = cosh (x) * sin (y);
403 return sinhz;
404 }
405
SLcomplex_cosh(double * coshz,double * z)406 double *SLcomplex_cosh (double *coshz, double *z)
407 {
408 double x, y;
409 x = z[0]; y = z[1];
410 coshz[0] = cosh (x) * cos (y);
411 coshz[1] = sinh (x) * sin (y);
412 return coshz;
413 }
414
SLcomplex_tanh(double * tanhz,double * z)415 double *SLcomplex_tanh (double *tanhz, double *z)
416 {
417 double x, y, invden;
418 x = 2 * z[0];
419 y = 2 * z[1];
420 invden = 1.0 / (cosh (x) + cos (y));
421 tanhz[0] = invden * sinh (x);
422 tanhz[1] = invden * sin (y);
423 return tanhz;
424 }
425 #if 0
426 static double *not_implemented (char *fun, double *p)
427 {
428 _pSLang_verror (SL_NOT_IMPLEMENTED, "%s for complex numbers has not been implemented",
429 fun);
430 *p = -1.0;
431 return p;
432 }
433 #endif
434 /* Use: asinh(z) = -i asin(iz) */
SLcomplex_asinh(double * asinhz,double * z)435 double *SLcomplex_asinh (double *asinhz, double *z)
436 {
437 double iz[2];
438
439 iz[0] = -z[1];
440 iz[1] = z[0];
441
442 (void) SLcomplex_asin (iz, iz);
443 asinhz[0] = iz[1];
444 asinhz[1] = -iz[0];
445
446 return asinhz;
447 }
448
449 /* Use: acosh (z) = i acos(z) */
SLcomplex_acosh(double * acoshz,double * z)450 double *SLcomplex_acosh (double *acoshz, double *z)
451 {
452 double iz[2];
453
454 (void) SLcomplex_acos (iz, z);
455 acoshz[0] = -iz[1];
456 acoshz[1] = iz[0];
457
458 return acoshz;
459 }
460
461 /* Use: atanh(z) = -i atan(iz) */
SLcomplex_atanh(double * atanhz,double * z)462 double *SLcomplex_atanh (double *atanhz, double *z)
463 {
464 double iz[2];
465
466 iz[0] = -z[1];
467 iz[1] = z[0];
468
469 (void) SLcomplex_atan (iz, iz);
470 atanhz[0] = iz[1];
471 atanhz[1] = -iz[0];
472
473 return atanhz;
474 }
475
complex_binary_result(int op,SLtype a,SLtype b,SLtype * c)476 static int complex_binary_result (int op, SLtype a, SLtype b,
477 SLtype *c)
478 {
479 (void) a; (void) b;
480
481 switch (op)
482 {
483 default:
484 case SLANG_POW:
485 case SLANG_PLUS:
486 case SLANG_MINUS:
487 case SLANG_TIMES:
488 case SLANG_DIVIDE:
489 *c = SLANG_COMPLEX_TYPE;
490 break;
491
492 case SLANG_EQ:
493 case SLANG_NE:
494 *c = SLANG_CHAR_TYPE;
495 break;
496 }
497 return 1;
498 }
499
complex_complex_binary(int op,SLtype a_type,VOID_STAR ap,SLuindex_Type na,SLtype b_type,VOID_STAR bp,SLuindex_Type nb,VOID_STAR cp)500 static int complex_complex_binary (int op,
501 SLtype a_type, VOID_STAR ap, SLuindex_Type na,
502 SLtype b_type, VOID_STAR bp, SLuindex_Type nb,
503 VOID_STAR cp)
504 {
505 char *ic;
506 double *a, *b, *c;
507 SLuindex_Type n, n_max;
508 SLuindex_Type da, db;
509
510 (void) a_type;
511 (void) b_type;
512
513 a = (double *) ap;
514 b = (double *) bp;
515 c = (double *) cp;
516 ic = (char *) cp;
517
518 if (na == 1) da = 0; else da = 2;
519 if (nb == 1) db = 0; else db = 2;
520
521 if (na > nb) n_max = na; else n_max = nb;
522 n_max = 2 * n_max;
523
524 switch (op)
525 {
526 default:
527 return 0;
528
529 case SLANG_PLUS:
530 for (n = 0; n < n_max; n += 2)
531 {
532 c[n] = a[0] + b[0];
533 c[n + 1] = a[1] + b[1];
534 a += da; b += db;
535 }
536 break;
537
538 case SLANG_MINUS:
539 for (n = 0; n < n_max; n += 2)
540 {
541 c[n] = a[0] - b[0];
542 c[n + 1] = a[1] - b[1];
543 a += da; b += db;
544 }
545 break;
546
547 case SLANG_TIMES:
548 for (n = 0; n < n_max; n += 2)
549 {
550 SLcomplex_times (c + n, a, b);
551 a += da; b += db;
552 }
553 break;
554
555 case SLANG_DIVIDE: /* / */
556 for (n = 0; n < n_max; n += 2)
557 {
558 #if 0
559 if ((b[0] == 0.0) && (b[1] == 0.0))
560 {
561 SLang_set_error (SL_DIVIDE_ERROR);
562 return -1;
563 }
564 #endif
565 SLcomplex_divide (c + n, a, b);
566 a += da; b += db;
567 }
568 break;
569
570 case SLANG_EQ: /* == */
571 for (n = 0; n < n_max; n += 2)
572 {
573 ic[n/2] = ((a[0] == b[0]) && (a[1] == b[1]));
574 a += da; b += db;
575 }
576 break;
577
578 case SLANG_NE: /* != */
579 for (n = 0; n < n_max; n += 2)
580 {
581 ic[n/2] = ((a[0] != b[0]) || (a[1] != b[1]));
582 a += da; b += db;
583 }
584 break;
585
586 case SLANG_POW:
587 for (n = 0; n < n_max; n += 2)
588 {
589 SLcomplex_pow (c + n, a, b);
590 a += da; b += db;
591 }
592 break;
593
594 }
595
596 return 1;
597 }
598
complex_double_binary(int op,SLtype a_type,VOID_STAR ap,SLuindex_Type na,SLtype b_type,VOID_STAR bp,SLuindex_Type nb,VOID_STAR cp)599 static int complex_double_binary (int op,
600 SLtype a_type, VOID_STAR ap, SLuindex_Type na,
601 SLtype b_type, VOID_STAR bp, SLuindex_Type nb,
602 VOID_STAR cp)
603 {
604 char *ic;
605 double *a, *b, *c;
606 SLuindex_Type n, n_max;
607 SLuindex_Type da, db;
608
609 (void) a_type;
610 (void) b_type;
611
612 a = (double *) ap;
613 b = (double *) bp;
614 c = (double *) cp;
615 ic = (char *) cp;
616
617 if (na == 1) da = 0; else da = 2;
618 if (nb == 1) db = 0; else db = 1;
619
620 if (na > nb) n_max = na; else n_max = nb;
621 n_max = 2 * n_max;
622
623 switch (op)
624 {
625 default:
626 return 0;
627
628 case SLANG_PLUS:
629 for (n = 0; n < n_max; n += 2)
630 {
631 c[n] = a[0] + b[0];
632 c[n + 1] = a[1];
633 a += da; b += db;
634 }
635 break;
636
637 case SLANG_MINUS:
638 for (n = 0; n < n_max; n += 2)
639 {
640 c[n] = a[0] - b[0];
641 c[n + 1] = a[1];
642 a += da; b += db;
643 }
644 break;
645
646 case SLANG_TIMES:
647 for (n = 0; n < n_max; n += 2)
648 {
649 double b0 = b[0];
650 c[n] = a[0] * b0;
651 c[n + 1] = a[1] * b0;
652 a += da; b += db;
653 }
654 break;
655
656 case SLANG_DIVIDE: /* / */
657 for (n = 0; n < n_max; n += 2)
658 {
659 double b0 = b[0];
660 #if 0
661 if (b0 == 0.0)
662 {
663 SLang_set_error (SL_DIVIDE_ERROR);
664 return -1;
665 }
666 #endif
667 c[n] = a[0] / b0;
668 c[n + 1] = a[1] / b0;
669 a += da; b += db;
670 }
671 break;
672
673 case SLANG_EQ: /* == */
674 for (n = 0; n < n_max; n += 2)
675 {
676 ic[n/2] = ((a[0] == b[0]) && (a[1] == 0.0));
677 a += da; b += db;
678 }
679 break;
680
681 case SLANG_NE: /* != */
682 for (n = 0; n < n_max; n += 2)
683 {
684 ic[n/2] = ((a[0] != b[0]) || (a[1] != 0.0));
685 a += da; b += db;
686 }
687 break;
688
689 case SLANG_POW:
690 for (n = 0; n < n_max; n += 2)
691 {
692 complex_dpow (c + n, a, b[0]);
693 a += da; b += db;
694 }
695 break;
696 }
697
698 return 1;
699 }
700
double_complex_binary(int op,SLtype a_type,VOID_STAR ap,SLuindex_Type na,SLtype b_type,VOID_STAR bp,SLuindex_Type nb,VOID_STAR cp)701 static int double_complex_binary (int op,
702 SLtype a_type, VOID_STAR ap, SLuindex_Type na,
703 SLtype b_type, VOID_STAR bp, SLuindex_Type nb,
704 VOID_STAR cp)
705 {
706 char *ic;
707 double *a, *b, *c;
708 SLuindex_Type n, n_max;
709 SLuindex_Type da, db;
710
711 (void) a_type;
712 (void) b_type;
713
714 a = (double *) ap;
715 b = (double *) bp;
716 c = (double *) cp;
717 ic = (char *) cp;
718
719 if (na == 1) da = 0; else da = 1;
720 if (nb == 1) db = 0; else db = 2;
721
722 if (na > nb) n_max = na; else n_max = nb;
723 n_max = 2 * n_max;
724
725 switch (op)
726 {
727 default:
728 return 0;
729
730 case SLANG_PLUS:
731 for (n = 0; n < n_max; n += 2)
732 {
733 c[n] = a[0] + b[0];
734 c[n + 1] = b[1];
735 a += da; b += db;
736 }
737 break;
738
739 case SLANG_MINUS:
740 for (n = 0; n < n_max; n += 2)
741 {
742 c[n] = a[0] - b[0];
743 c[n + 1] = -b[1];
744 a += da; b += db;
745 }
746 break;
747
748 case SLANG_TIMES:
749 for (n = 0; n < n_max; n += 2)
750 {
751 double a0 = a[0];
752 c[n] = a0 * b[0];
753 c[n + 1] = a0 * b[1];
754 a += da; b += db;
755 }
756 break;
757
758 case SLANG_DIVIDE: /* / */
759 for (n = 0; n < n_max; n += 2)
760 {
761 double z[2];
762 #if 0
763 if ((b[0] == 0.0) && (b[1] == 0.0))
764 {
765 SLang_set_error (SL_DIVIDE_ERROR);
766 return -1;
767 }
768 #endif
769 z[0] = a[0];
770 z[1] = 0.0;
771 SLcomplex_divide (c + n, z, b);
772 a += da; b += db;
773 }
774 break;
775
776 case SLANG_EQ: /* == */
777 for (n = 0; n < n_max; n += 2)
778 {
779 ic[n/2] = ((a[0] == b[0]) && (0.0 == b[1]));
780 a += da; b += db;
781 }
782 break;
783
784 case SLANG_NE: /* != */
785 for (n = 0; n < n_max; n += 2)
786 {
787 ic[n/2] = ((a[0] != b[0]) || (0.0 != b[1]));
788 a += da; b += db;
789 }
790 break;
791
792 case SLANG_POW:
793 for (n = 0; n < n_max; n += 2)
794 {
795 dcomplex_pow (c + n, a[0], b);
796 a += da; b += db;
797 }
798 break;
799 }
800
801 return 1;
802 }
803
complex_generic_binary(int op,SLtype a_type,VOID_STAR ap,SLuindex_Type na,SLtype b_type,VOID_STAR bp,SLuindex_Type nb,VOID_STAR cp)804 static int complex_generic_binary (int op,
805 SLtype a_type, VOID_STAR ap, SLuindex_Type na,
806 SLtype b_type, VOID_STAR bp, SLuindex_Type nb,
807 VOID_STAR cp)
808 {
809 char *ic;
810 char *b;
811 double *a, *c;
812 SLuindex_Type n, n_max;
813 SLuindex_Type da, db;
814 unsigned int sizeof_b;
815 SLang_To_Double_Fun_Type to_double;
816
817 if (NULL == (to_double = SLarith_get_to_double_fun (b_type, &sizeof_b)))
818 return 0;
819
820 (void) a_type;
821
822 a = (double *) ap;
823 b = (char *) bp;
824 c = (double *) cp;
825 ic = (char *) cp;
826
827 if (na == 1) da = 0; else da = 2;
828 if (nb == 1) db = 0; else db = sizeof_b;
829
830 if (na > nb) n_max = na; else n_max = nb;
831 n_max = 2 * n_max;
832
833 switch (op)
834 {
835 default:
836 return 0;
837
838 case SLANG_POW:
839 for (n = 0; n < n_max; n += 2)
840 {
841 complex_dpow (c + n, a, to_double((VOID_STAR)b));
842 a += da; b += db;
843 }
844 break;
845
846 case SLANG_PLUS:
847 for (n = 0; n < n_max; n += 2)
848 {
849 c[n] = a[0] + to_double((VOID_STAR)b);
850 c[n + 1] = a[1];
851 a += da; b += db;
852 }
853 break;
854
855 case SLANG_MINUS:
856 for (n = 0; n < n_max; n += 2)
857 {
858 c[n] = a[0] - to_double((VOID_STAR)b);
859 c[n + 1] = a[1];
860 a += da; b += db;
861 }
862 break;
863
864 case SLANG_TIMES:
865 for (n = 0; n < n_max; n += 2)
866 {
867 double b0 = to_double((VOID_STAR)b);
868 c[n] = a[0] * b0;
869 c[n + 1] = a[1] * b0;
870 a += da; b += db;
871 }
872 break;
873
874 case SLANG_DIVIDE: /* / */
875 for (n = 0; n < n_max; n += 2)
876 {
877 double b0 = to_double((VOID_STAR)b);
878 #if 0
879 if (b0 == 0)
880 {
881 SLang_set_error (SL_DIVIDE_ERROR);
882 return -1;
883 }
884 #endif
885 c[n] = a[0] / b0;
886 c[n + 1] = a[1] / b0;
887 a += da; b += db;
888 }
889 break;
890
891 case SLANG_EQ: /* == */
892 for (n = 0; n < n_max; n += 2)
893 {
894 ic[n/2] = ((a[0] == to_double((VOID_STAR)b)) && (a[1] == 0.0));
895 a += da; b += db;
896 }
897 break;
898
899 case SLANG_NE: /* != */
900 for (n = 0; n < n_max; n += 2)
901 {
902 ic[n/2] = ((a[0] != to_double((VOID_STAR)b)) || (a[1] != 0.0));
903 a += da; b += db;
904 }
905 break;
906 }
907
908 return 1;
909 }
910
generic_complex_binary(int op,SLtype a_type,VOID_STAR ap,SLuindex_Type na,SLtype b_type,VOID_STAR bp,SLuindex_Type nb,VOID_STAR cp)911 static int generic_complex_binary (int op,
912 SLtype a_type, VOID_STAR ap, SLuindex_Type na,
913 SLtype b_type, VOID_STAR bp, SLuindex_Type nb,
914 VOID_STAR cp)
915 {
916 double *b, *c;
917 char *a, *ic;
918 SLuindex_Type n, n_max;
919 SLuindex_Type da, db;
920 unsigned int sizeof_a;
921 SLang_To_Double_Fun_Type to_double;
922
923 if (NULL == (to_double = SLarith_get_to_double_fun (a_type, &sizeof_a)))
924 return 0;
925
926 (void) b_type;
927
928 a = (char *) ap;
929 b = (double *) bp;
930 c = (double *) cp;
931 ic = (char *) cp;
932
933 if (na == 1) da = 0; else da = sizeof_a;
934 if (nb == 1) db = 0; else db = 2;
935
936 if (na > nb) n_max = na; else n_max = nb;
937 n_max = 2 * n_max;
938
939 switch (op)
940 {
941 default:
942 return 0;
943 case SLANG_POW:
944 for (n = 0; n < n_max; n += 2)
945 {
946 dcomplex_pow (c + n, to_double((VOID_STAR)a), b);
947 a += da; b += db;
948 }
949 break;
950
951 case SLANG_PLUS:
952 for (n = 0; n < n_max; n += 2)
953 {
954 c[n] = to_double((VOID_STAR)a) + b[0];
955 c[n + 1] = b[1];
956 a += da; b += db;
957 }
958 break;
959
960 case SLANG_MINUS:
961 for (n = 0; n < n_max; n += 2)
962 {
963 c[n] = to_double((VOID_STAR)a) - b[0];
964 c[n + 1] = -b[1];
965 a += da; b += db;
966 }
967 break;
968
969 case SLANG_TIMES:
970 for (n = 0; n < n_max; n += 2)
971 {
972 double a0 = to_double((VOID_STAR)a);
973 c[n] = a0 * b[0];
974 c[n + 1] = a0 * b[1];
975 a += da; b += db;
976 }
977 break;
978
979 case SLANG_DIVIDE: /* / */
980 for (n = 0; n < n_max; n += 2)
981 {
982 double z[2];
983 #if 0
984 if ((b[0] == 0.0) && (b[1] == 0.0))
985 {
986 SLang_set_error (SL_DIVIDE_ERROR);
987 return -1;
988 }
989 #endif
990 z[0] = to_double((VOID_STAR)a);
991 z[1] = 0.0;
992 SLcomplex_divide (c + n, z, b);
993 a += da; b += db;
994 }
995 break;
996
997 case SLANG_EQ: /* == */
998 for (n = 0; n < n_max; n += 2)
999 {
1000 ic[n/2] = ((to_double((VOID_STAR)a) == b[0]) && (0.0 == b[1]));
1001 a += da; b += db;
1002 }
1003 break;
1004
1005 case SLANG_NE: /* != */
1006 for (n = 0; n < n_max; n += 2)
1007 {
1008 ic[n/2] = ((to_double((VOID_STAR)a) != b[0]) || (0.0 != b[1]));
1009 a += da; b += db;
1010 }
1011 break;
1012 }
1013
1014 return 1;
1015 }
1016
complex_unary_result(int op,SLtype a,SLtype * b)1017 static int complex_unary_result (int op, SLtype a, SLtype *b)
1018 {
1019 (void) a;
1020
1021 switch (op)
1022 {
1023 default:
1024 return 0;
1025
1026 case SLANG_PLUSPLUS:
1027 case SLANG_MINUSMINUS:
1028 case SLANG_CHS:
1029 case SLANG_MUL2:
1030 *b = SLANG_COMPLEX_TYPE;
1031 break;
1032
1033 case SLANG_SQR: /* |Real|^2 + |Imag|^2 ==> double */
1034 case SLANG_ABS: /* |z| ==> double */
1035 *b = SLANG_DOUBLE_TYPE;
1036 break;
1037
1038 case SLANG_SIGN:
1039 *b = SLANG_INT_TYPE;
1040 break;
1041 }
1042 return 1;
1043 }
1044
complex_unary(int op,SLtype a_type,VOID_STAR ap,SLuindex_Type na,VOID_STAR bp)1045 static int complex_unary (int op,
1046 SLtype a_type, VOID_STAR ap, SLuindex_Type na,
1047 VOID_STAR bp)
1048 {
1049 SLuindex_Type n;
1050 double *a, *b;
1051 int *ic;
1052
1053 (void) a_type;
1054
1055 a = (double *) ap;
1056 b = (double *) bp;
1057 ic = (int *) bp;
1058
1059 na = 2 * na;
1060
1061 switch (op)
1062 {
1063 default:
1064 return 0;
1065
1066 case SLANG_PLUSPLUS:
1067 for (n = 0; n < na; n += 2) b[n] = (a[n] + 1);
1068 break;
1069 case SLANG_MINUSMINUS:
1070 for (n = 0; n < na; n += 2) b[n] = (a[n] - 1);
1071 break;
1072 case SLANG_CHS:
1073 for (n = 0; n < na; n += 2)
1074 {
1075 b[n] = -(a[n]);
1076 b[n + 1] = -(a[n + 1]);
1077 }
1078 break;
1079 case SLANG_SQR: /* |Real|^2 + |Imag|^2 ==> double */
1080 for (n = 0; n < na; n += 2)
1081 b[n/2] = (a[n] * a[n] + a[n + 1] * a[n + 1]);
1082 break;
1083
1084 case SLANG_MUL2:
1085 for (n = 0; n < na; n += 2)
1086 {
1087 b[n] = (2 * a[n]);
1088 b[n + 1] = (2 * a[n + 1]);
1089 }
1090 break;
1091
1092 case SLANG_ABS: /* |z| ==> double */
1093 for (n = 0; n < na; n += 2)
1094 b[n/2] = SLcomplex_abs (a + n);
1095 break;
1096
1097 case SLANG_SIGN:
1098 /* Another creative extension. Lets return an integer which indicates
1099 * whether the complex number is in the upperhalf plane or not.
1100 */
1101 for (n = 0; n < na; n += 2)
1102 {
1103 if (a[n + 1] < 0.0) ic[n/2] = -1;
1104 else if (a[n + 1] > 0.0) ic[n/2] = 1;
1105 else ic[n/2] = 0;
1106 }
1107 break;
1108 }
1109
1110 return 1;
1111 }
1112
1113 static int
complex_typecast(SLtype from_type,VOID_STAR from,SLuindex_Type num,SLtype to_type,VOID_STAR to)1114 complex_typecast (SLtype from_type, VOID_STAR from, SLuindex_Type num,
1115 SLtype to_type, VOID_STAR to)
1116 {
1117 double *z;
1118 double *d;
1119 char *i;
1120 SLuindex_Type n;
1121 unsigned int sizeof_i;
1122 SLang_To_Double_Fun_Type to_double;
1123
1124 (void) to_type;
1125
1126 z = (double *) to;
1127
1128 switch (from_type)
1129 {
1130 default:
1131 if (NULL == (to_double = SLarith_get_to_double_fun (from_type, &sizeof_i)))
1132 return 0;
1133 i = (char *) from;
1134 for (n = 0; n < num; n++)
1135 {
1136 *z++ = to_double ((VOID_STAR) i);
1137 *z++ = 0.0;
1138
1139 i += sizeof_i;
1140 }
1141 break;
1142
1143 case SLANG_DOUBLE_TYPE:
1144 d = (double *) from;
1145 for (n = 0; n < num; n++)
1146 {
1147 *z++ = d[n];
1148 *z++ = 0.0;
1149 }
1150 break;
1151 }
1152
1153 return 1;
1154 }
1155
complex_destroy(SLtype type,VOID_STAR ptr)1156 static void complex_destroy (SLtype type, VOID_STAR ptr)
1157 {
1158 (void) type;
1159 SLfree ((char *)*(double **) ptr);
1160 }
1161
complex_push(SLtype type,VOID_STAR ptr)1162 static int complex_push (SLtype type, VOID_STAR ptr)
1163 {
1164 double *z;
1165
1166 (void) type;
1167 z = *(double **) ptr;
1168 return SLang_push_complex (z[0], z[1]);
1169 }
1170
complex_pop(SLtype type,VOID_STAR ptr)1171 static int complex_pop (SLtype type, VOID_STAR ptr)
1172 {
1173 double *z;
1174
1175 (void) type;
1176 z = *(double **) ptr;
1177 return SLang_pop_complex (&z[0], &z[1]);
1178 }
1179
_pSLinit_slcomplex(void)1180 int _pSLinit_slcomplex (void)
1181 {
1182 SLang_Class_Type *cl;
1183 SLtype *types;
1184
1185 if (NULL == (cl = SLclass_allocate_class ("Complex_Type")))
1186 return -1;
1187
1188 (void) SLclass_set_destroy_function (cl, complex_destroy);
1189 (void) SLclass_set_push_function (cl, complex_push);
1190 (void) SLclass_set_pop_function (cl, complex_pop);
1191
1192 if (-1 == SLclass_register_class (cl, SLANG_COMPLEX_TYPE, 2 * sizeof (double),
1193 SLANG_CLASS_TYPE_VECTOR))
1194 return -1;
1195
1196 types = _pSLarith_Arith_Types;
1197 while (*types != SLANG_DOUBLE_TYPE)
1198 {
1199 SLtype t = *types++;
1200
1201 if ((-1 == SLclass_add_binary_op (t, SLANG_COMPLEX_TYPE, generic_complex_binary, complex_binary_result))
1202 || (-1 == SLclass_add_binary_op (SLANG_COMPLEX_TYPE, t, complex_generic_binary, complex_binary_result))
1203 || (-1 == (SLclass_add_typecast (t, SLANG_COMPLEX_TYPE, complex_typecast, 1))))
1204 return -1;
1205 }
1206
1207 if ((-1 == (SLclass_add_binary_op (SLANG_COMPLEX_TYPE, SLANG_COMPLEX_TYPE, complex_complex_binary, complex_binary_result)))
1208 || (-1 == (SLclass_add_binary_op (SLANG_COMPLEX_TYPE, SLANG_DOUBLE_TYPE, complex_double_binary, complex_binary_result)))
1209 || (-1 == (SLclass_add_binary_op (SLANG_DOUBLE_TYPE, SLANG_COMPLEX_TYPE, double_complex_binary, complex_binary_result)))
1210 || (-1 == (SLclass_add_unary_op (SLANG_COMPLEX_TYPE, complex_unary, complex_unary_result)))
1211 || (-1 == (SLclass_add_typecast (SLANG_DOUBLE_TYPE, SLANG_COMPLEX_TYPE, complex_typecast, 1))))
1212 return -1;
1213
1214 return 0;
1215 }
1216
1217 #endif /* if SLANG_HAS_COMPLEX */
1218
1219