1 /* sin, cos, etc, 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 #ifndef _GNU_SOURCE
24 # define _GNU_SOURCE /* sincos */
25 #endif
26 #include "slinclud.h"
27 #include <float.h>
28 #include <math.h>
29
30 #if SLANG_HAS_FLOAT
31
32 #ifdef HAVE_FLOATINGPOINT_H
33 # include <floatingpoint.h>
34 #endif
35
36 #ifdef HAVE_FENV_H
37 # include <fenv.h>
38 #endif
39
40 #ifdef HAVE_IEEEFP_H
41 # include <ieeefp.h>
42 #endif
43
44 #ifdef HAVE_NAN_H
45 # include <nan.h>
46 #endif
47
48 #include "slang.h"
49 #include "_slang.h"
50
51 #ifdef PI
52 # undef PI
53 #endif
54 #define PI 3.14159265358979323846264338327950288
55
56 #if defined(__unix__)
57 #include <signal.h>
58 #include <errno.h>
59
math_floating_point_exception(int sig)60 static void math_floating_point_exception (int sig)
61 {
62 sig = errno;
63 (void) SLsignal (SIGFPE, math_floating_point_exception);
64 errno = sig;
65 }
66 #endif
67
SLmath_hypot(double x,double y)68 double SLmath_hypot (double x, double y)
69 {
70 #ifdef HAVE_HYPOT
71 return hypot (x,y);
72 #else
73 double fr, fi, ratio;
74
75 fr = fabs(x);
76 fi = fabs(y);
77
78 if (fr > fi)
79 {
80 ratio = y / x;
81 x = fr * sqrt (1.0 + ratio * ratio);
82 }
83 else if (fi == 0.0) x = 0.0;
84 else
85 {
86 ratio = x / y;
87 x = fi * sqrt (1.0 + ratio * ratio);
88 }
89
90 return x;
91 #endif
92 }
93
double_math_op_result(int op,SLtype a,SLtype * b)94 static int double_math_op_result (int op, SLtype a, SLtype *b)
95 {
96 switch (op)
97 {
98 case SLMATH_ISINF:
99 case SLMATH_ISNAN:
100 *b = SLANG_CHAR_TYPE;
101 break;
102
103 default:
104 if (a != SLANG_FLOAT_TYPE)
105 *b = SLANG_DOUBLE_TYPE;
106 else
107 *b = a;
108 break;
109 }
110 return 1;
111 }
112
113 #ifdef HAVE_ISNAN
114 # define ISNAN_FUN isnan
115 #else
116 # define ISNAN_FUN my_isnan
my_isnan(double x)117 static int my_isnan (double x)
118 {
119 return (volatile double) x != (volatile double) x;
120 }
121 #endif
122
_pSLmath_isnan(double x)123 int _pSLmath_isnan (double x)
124 {
125 return ISNAN_FUN (x);
126 }
127
128 #ifdef HAVE_ISINF
129 # define ISINF_FUN isinf
130 #else
131 # define ISINF_FUN my_isinf
my_isinf(double x)132 static int my_isinf (double x)
133 {
134 #ifdef HAVE_FINITE
135 return (0 == finite (x)) && (0 == ISNAN_FUN (x));
136 #else
137 double y;
138 if (x == 0.0)
139 return 0;
140 y = x * 0.5;
141 return (x == y);
142 #endif
143 }
144 #endif
145
_pSLmath_isinf(double x)146 int _pSLmath_isinf (double x)
147 {
148 return ISINF_FUN (x);
149 }
150
151 #ifdef HAVE_ROUND
152 # define ROUND_FUN round
153 #else
154 # define ROUND_FUN my_round
my_round(double x)155 static double my_round (double x)
156 {
157 double xf, xi;
158
159 xf = modf (x, &xi); /* x = xi + xf */
160 if (xi > 0)
161 {
162 if (xf >= 0.5)
163 return xi + 1.0;
164 }
165 else if (xi < 0)
166 {
167 if (xf <= -0.5)
168 return xi - 1.0;
169 }
170 else if (xf < 0) /* xi=0 */
171 {
172 if (xf <= -0.5)
173 return -1.0;
174 }
175 else if (xf >= 0.5) /* xi=0 */
176 return 1.0;
177
178 return xi;
179 }
180 #endif
181
182 #ifndef HAVE_LOG1P
_pSLmath_log1p(double x)183 double _pSLmath_log1p (double x)
184 {
185 /*
186 * log(1+x) ~= x for small x
187 * Suppose the error in 1+x is e. That is, 1+x produces 1+x+e.
188 * Then log(1+x) will produce log(1+x+e). The error is
189 * df = log(1+x+e)-log(1+x)
190 * = log (1 + e/(1+x))
191 * ~= e/(1+x)
192 * Here e = (1+x) - (1+x)
193 * = ((1+x)-1) - x
194 * So compute:
195 * log1p(x) = log(1+x) - (((1+x)-1)-x)/(1+x)
196 */
197 volatile double xp1, xx;
198 if (ISINF_FUN(x))
199 {
200 if (x < 0) return _pSLang_NaN;
201 return _pSLang_Inf;
202 }
203 xp1 = 1.0 + x;
204 if (xp1 == 0.0)
205 return -_pSLang_Inf;
206 xx = xp1-1.0;
207 return log(xp1) - (xx-x)/xp1;
208 }
209 #endif
210
211 #ifndef HAVE_EXPM1
_pSLmath_expm1(double x)212 double _pSLmath_expm1 (double x)
213 {
214 /* f = exp(x)-1
215 * = (exp(x)-1) + x - x
216 * = (exp(x)-(1+x)) + x
217 *
218 * While the above works, it is not as accurate as the algroithm below,
219 * which is due to Kahan (yes, that Kahan).
220 */
221 volatile double u;
222 if (ISINF_FUN(x))
223 {
224 if (x < 0) return -1.0;
225 return _pSLang_Inf;
226 }
227
228 u = exp(x);
229 if (u == 1.0)
230 return x;
231 if (u - 1.0 == -1.0)
232 return -1.0;
233 return (u-1.0)*x/log(u);
234 }
235 #endif
236
237 #ifdef HAVE_ASINH
238 # define ASINH_FUN asinh
239 #else
240 # define ASINH_FUN my_asinh
my_asinh(double x)241 static double my_asinh (double x)
242 {
243 /* f(x) = log(x + sqrt(x^2+1))
244 * Note: f(x)+f(-x)
245 * = log(x+sqrt(x^2+1))+log(-x+sqrt(x^2+1))
246 * = log(x^2+1 - x^2)
247 * = 0
248 * Thus, f(-x)=-f(x) as required.
249 *
250 * So consider x>=0.
251 * For large x
252 * = log(2x + sqrt(x^2+1)-x)
253 * = log(2x + ((x^2+1)-x^2)/(x+sqrt(x^2+1)))
254 * = log(2x + 1/(x+sqrt(x^2+1)))
255 * = log(2x + 1/x/(1+sqrt(1+1/x^2)))
256 * For small x
257 * = log(1 + x + sqrt(1+x^2)-1)
258 * = log(1 + x + (1+x^2-1)/(sqrt(1+x^2)+1))
259 * = log(1 + x + x^2/(1+sqrt(1+x^2)))
260 */
261 double s = 1.0;
262 double x2 = x*x;
263
264 if (x < 0.0)
265 {
266 s = -1.0;
267 x = -x;
268 }
269
270 if (x > 1.0)
271 return s*log (2.0*x + 1.0/x/(1.0+sqrt(1.0+1.0/x2)));
272
273 return LOG1P_FUNC(x+x2/(1.0+sqrt(1.0+x2)));
274 }
275 #endif
276
277 #ifdef HAVE_ACOSH
278 # define ACOSH_FUN acosh
279 #else
280 # define ACOSH_FUN my_acosh
my_acosh(double x)281 static double my_acosh (double x)
282 {
283 /* log (x + sqrt(x^2 - 1); x >= 1
284 * = log(2*x + sqrt(x^2-1)-x)
285 * = log(2*x + (x^2-1-x^2)/(x+sqrt(x^2-1)))
286 * = log(2*x - 1/x/(1+sqrt(1-1/x^2)))
287 */
288 return log(2.0*x - 1.0/x/(1.0+sqrt(1.0-1/x/x)));
289 }
290 #endif
291 #ifdef HAVE_ATANH
292 # define ATANH_FUN atanh
293 #else
294 # define ATANH_FUN my_atanh
295
my_atanh(double x)296 static double my_atanh (double x)
297 {
298 /* 0.5 * log ((1.0 + x)/(1.0 - x)); (0 <= x^2 < 1)
299 * = 0.5*log((1-x+2x)/(1-x))
300 * = 0.5*log(1 + 2x/(1-x))
301 */
302 return 0.5 * LOG1P_FUNC(2.0*x/(1-x));
303 }
304 #endif
305
double_math_op(int op,SLtype type,VOID_STAR ap,SLuindex_Type na,VOID_STAR bp)306 static int double_math_op (int op,
307 SLtype type, VOID_STAR ap, SLuindex_Type na,
308 VOID_STAR bp)
309 {
310 double *a, *b;
311 unsigned int i;
312 char *c;
313
314 (void) type;
315 a = (double *) ap;
316 b = (double *) bp;
317
318 switch (op)
319 {
320 default:
321 return 0;
322
323 case SLMATH_SINH:
324 for (i = 0; i < na; i++) b[i] = sinh(a[i]);
325 break;
326 case SLMATH_COSH:
327 for (i = 0; i < na; i++) b[i] = cosh(a[i]);
328 break;
329 case SLMATH_TANH:
330 for (i = 0; i < na; i++) b[i] = tanh(a[i]);
331 break;
332 case SLMATH_TAN:
333 for (i = 0; i < na; i++) b[i] = tan(a[i]);
334 break;
335 case SLMATH_ASIN:
336 for (i = 0; i < na; i++) b[i] = asin(a[i]);
337 break;
338 case SLMATH_ACOS:
339 for (i = 0; i < na; i++) b[i] = acos(a[i]);
340 break;
341 case SLMATH_ATAN:
342 for (i = 0; i < na; i++) b[i] = atan(a[i]);
343 break;
344 case SLMATH_EXP:
345 for (i = 0; i < na; i++) b[i] = exp(a[i]);
346 break;
347 case SLMATH_LOG:
348 for (i = 0; i < na; i++) b[i] = log(a[i]);
349 break;
350 case SLMATH_LOG10:
351 for (i = 0; i < na; i++) b[i] = log10(a[i]);
352 break;
353 case SLMATH_SQRT:
354 for (i = 0; i < na; i++) b[i] = sqrt (a[i]);
355 return 1;
356 case SLMATH_SIN:
357 for (i = 0; i < na; i++) b[i] = sin(a[i]);
358 break;
359 case SLMATH_COS:
360 for (i = 0; i < na; i++) b[i] = cos(a[i]);
361 break;
362 case SLMATH_ASINH:
363 for (i = 0; i < na; i++) b[i] = ASINH_FUN(a[i]);
364 break;
365 case SLMATH_ATANH:
366 for (i = 0; i < na; i++) b[i] = ATANH_FUN(a[i]);
367 break;
368 case SLMATH_ACOSH:
369 for (i = 0; i < na; i++) b[i] = ACOSH_FUN(a[i]);
370 break;
371
372 case SLMATH_CONJ:
373 case SLMATH_REAL:
374 for (i = 0; i < na; i++)
375 b[i] = a[i];
376 return 1;
377 case SLMATH_IMAG:
378 for (i = 0; i < na; i++)
379 b[i] = 0.0;
380 return 1;
381
382 case SLMATH_ISINF:
383 c = (char *) bp;
384 for (i = 0; i < na; i++)
385 c[i] = (char) ISINF_FUN(a[i]);
386 return 1;
387
388 case SLMATH_ISNAN:
389 c = (char *) bp;
390 for (i = 0; i < na; i++)
391 c[i] = (char) ISNAN_FUN(a[i]);
392 return 1;
393
394 case SLMATH_FLOOR:
395 for (i = 0; i < na; i++) b[i] = floor(a[i]);
396 break;
397
398 case SLMATH_CEIL:
399 for (i = 0; i < na; i++) b[i] = ceil(a[i]);
400 break;
401
402 case SLMATH_ROUND:
403 for (i = 0; i < na; i++) b[i] = ROUND_FUN(a[i]);
404 break;
405
406 case SLMATH_EXPM1:
407 for (i = 0; i < na; i++) b[i] = EXPM1_FUNC(a[i]);
408 break;
409
410 case SLMATH_LOG1P:
411 for (i = 0; i < na; i++) b[i] = LOG1P_FUNC(a[i]);
412 break;
413 }
414
415 return 1;
416 }
417
float_math_op(int op,SLtype type,VOID_STAR ap,SLuindex_Type na,VOID_STAR bp)418 static int float_math_op (int op,
419 SLtype type, VOID_STAR ap, SLuindex_Type na,
420 VOID_STAR bp)
421 {
422 float *a, *b;
423 unsigned int i;
424 char *c;
425
426 (void) type;
427 a = (float *) ap;
428 b = (float *) bp;
429
430 switch (op)
431 {
432 default:
433 return 0;
434
435 case SLMATH_SINH:
436 for (i = 0; i < na; i++) b[i] = sinh(a[i]);
437 break;
438 case SLMATH_COSH:
439 for (i = 0; i < na; i++) b[i] = cosh(a[i]);
440 break;
441 case SLMATH_TANH:
442 for (i = 0; i < na; i++) b[i] = tanh(a[i]);
443 break;
444 case SLMATH_TAN:
445 for (i = 0; i < na; i++) b[i] = tan(a[i]);
446 break;
447 case SLMATH_ASIN:
448 for (i = 0; i < na; i++) b[i] = asin(a[i]);
449 break;
450 case SLMATH_ACOS:
451 for (i = 0; i < na; i++) b[i] = acos(a[i]);
452 break;
453 case SLMATH_ATAN:
454 for (i = 0; i < na; i++) b[i] = atan(a[i]);
455 break;
456 case SLMATH_EXP:
457 for (i = 0; i < na; i++) b[i] = exp(a[i]);
458 break;
459 case SLMATH_LOG:
460 for (i = 0; i < na; i++) b[i] = log(a[i]);
461 break;
462 case SLMATH_LOG10:
463 for (i = 0; i < na; i++) b[i] = log10(a[i]);
464 break;
465 case SLMATH_SQRT:
466 for (i = 0; i < na; i++)
467 b[i] = (float) sqrt ((double) a[i]);
468 return 1;
469
470 case SLMATH_SIN:
471 for (i = 0; i < na; i++) b[i] = sin(a[i]);
472 break;
473 case SLMATH_COS:
474 for (i = 0; i < na; i++) b[i] = cos(a[i]);
475 break;
476
477 case SLMATH_ASINH:
478 for (i = 0; i < na; i++) b[i] = ASINH_FUN(a[i]);
479 break;
480 case SLMATH_ATANH:
481 for (i = 0; i < na; i++) b[i] = ATANH_FUN(a[i]);
482 break;
483 case SLMATH_ACOSH:
484 for (i = 0; i < na; i++) b[i] = ACOSH_FUN(a[i]);
485 break;
486
487 case SLMATH_CONJ:
488 case SLMATH_REAL:
489 for (i = 0; i < na; i++)
490 b[i] = a[i];
491 return 1;
492 case SLMATH_IMAG:
493 for (i = 0; i < na; i++)
494 b[i] = 0.0;
495 return 1;
496 case SLMATH_ISINF:
497 c = (char *) bp;
498 for (i = 0; i < na; i++)
499 c[i] = (char) ISINF_FUN((double) a[i]);
500 return 1;
501
502 case SLMATH_ISNAN:
503 c = (char *) bp;
504 for (i = 0; i < na; i++)
505 c[i] = (char) ISNAN_FUN((double) a[i]);
506 return 1;
507
508 case SLMATH_FLOOR:
509 for (i = 0; i < na; i++) b[i] = floor(a[i]);
510 break;
511 case SLMATH_CEIL:
512 for (i = 0; i < na; i++) b[i] = ceil(a[i]);
513 break;
514 case SLMATH_ROUND:
515 for (i = 0; i < na; i++) b[i] = ROUND_FUN(a[i]);
516 break;
517
518 case SLMATH_EXPM1:
519 for (i = 0; i < na; i++) b[i] = EXPM1_FUNC(a[i]);
520 break;
521
522 case SLMATH_LOG1P:
523 for (i = 0; i < na; i++) b[i] = LOG1P_FUNC(a[i]);
524 break;
525 }
526
527 return 1;
528 }
529
generic_math_op(int op,SLtype type,VOID_STAR ap,SLuindex_Type na,VOID_STAR bp)530 static int generic_math_op (int op,
531 SLtype type, VOID_STAR ap, SLuindex_Type na,
532 VOID_STAR bp)
533 {
534 double *b;
535 unsigned int i;
536 SLang_To_Double_Fun_Type to_double;
537 unsigned int da;
538 char *a, *c;
539
540 if (NULL == (to_double = SLarith_get_to_double_fun (type, &da)))
541 return 0;
542
543 b = (double *) bp;
544 a = (char *) ap;
545
546 switch (op)
547 {
548 default:
549 return 0;
550
551 case SLMATH_SINH:
552 for (i = 0; i < na; i++)
553 {
554 b[i] = sinh (to_double ((VOID_STAR) a));
555 a += da;
556 }
557 break;
558 case SLMATH_COSH:
559 for (i = 0; i < na; i++)
560 {
561 b[i] = cosh (to_double ((VOID_STAR) a));
562 a += da;
563 }
564 break;
565 case SLMATH_TANH:
566 for (i = 0; i < na; i++)
567 {
568 b[i] = tanh (to_double ((VOID_STAR) a));
569 a += da;
570 }
571 break;
572 case SLMATH_TAN:
573 for (i = 0; i < na; i++)
574 {
575 b[i] = tan (to_double ((VOID_STAR) a));
576 a += da;
577 }
578 break;
579 case SLMATH_ASIN:
580 for (i = 0; i < na; i++)
581 {
582 b[i] = asin (to_double ((VOID_STAR) a));
583 a += da;
584 }
585 break;
586 case SLMATH_ACOS:
587 for (i = 0; i < na; i++)
588 {
589 b[i] = acos (to_double ((VOID_STAR) a));
590 a += da;
591 }
592 break;
593 case SLMATH_ATAN:
594 for (i = 0; i < na; i++)
595 {
596 b[i] = atan (to_double ((VOID_STAR) a));
597 a += da;
598 }
599 break;
600 case SLMATH_EXP:
601 for (i = 0; i < na; i++)
602 {
603 b[i] = exp (to_double ((VOID_STAR) a));
604 a += da;
605 }
606 break;
607 case SLMATH_LOG:
608 for (i = 0; i < na; i++)
609 {
610 b[i] = log (to_double ((VOID_STAR) a));
611 a += da;
612 }
613 break;
614 case SLMATH_LOG10:
615 for (i = 0; i < na; i++)
616 {
617 b[i] = log10 (to_double ((VOID_STAR) a));
618 a += da;
619 }
620 break;
621 case SLMATH_SQRT:
622 for (i = 0; i < na; i++)
623 {
624 b[i] = sqrt (to_double ((VOID_STAR) a));
625 a += da;
626 }
627 break;
628 case SLMATH_SIN:
629 for (i = 0; i < na; i++)
630 {
631 b[i] = sin (to_double ((VOID_STAR) a));
632 a += da;
633 }
634 break;
635 case SLMATH_COS:
636 for (i = 0; i < na; i++)
637 {
638 b[i] = cos (to_double ((VOID_STAR) a));
639 a += da;
640 }
641 break;
642
643 case SLMATH_ASINH:
644 for (i = 0; i < na; i++)
645 {
646 b[i] = ASINH_FUN (to_double ((VOID_STAR) a));
647 a += da;
648 }
649 break;
650 case SLMATH_ATANH:
651 for (i = 0; i < na; i++)
652 {
653 b[i] = ATANH_FUN (to_double ((VOID_STAR) a));
654 a += da;
655 }
656 break;
657 case SLMATH_ACOSH:
658 for (i = 0; i < na; i++)
659 {
660 b[i] = ACOSH_FUN (to_double ((VOID_STAR) a));
661 a += da;
662 }
663 break;
664
665 case SLMATH_CONJ:
666 case SLMATH_REAL:
667 for (i = 0; i < na; i++)
668 {
669 b[i] = to_double((VOID_STAR) a);
670 a += da;
671 }
672 return 1;
673
674 case SLMATH_IMAG:
675 for (i = 0; i < na; i++)
676 b[i] = 0.0;
677 return 1;
678
679 case SLMATH_ISINF:
680 c = (char *) bp;
681 for (i = 0; i < na; i++)
682 {
683 c[i] = (char) ISINF_FUN(to_double((VOID_STAR) a));
684 a += da;
685 }
686 return 1;
687 case SLMATH_ISNAN:
688 c = (char *) bp;
689 for (i = 0; i < na; i++)
690 {
691 c[i] = (char) ISNAN_FUN(to_double((VOID_STAR) a));
692 a += da;
693 }
694 return 1;
695 case SLMATH_FLOOR:
696 for (i = 0; i < na; i++)
697 {
698 b[i] = floor (to_double ((VOID_STAR) a));
699 a += da;
700 }
701 break;
702 case SLMATH_CEIL:
703 for (i = 0; i < na; i++)
704 {
705 b[i] = ceil (to_double ((VOID_STAR) a));
706 a += da;
707 }
708 break;
709 case SLMATH_ROUND:
710 for (i = 0; i < na; i++)
711 {
712 b[i] = ROUND_FUN (to_double ((VOID_STAR) a));
713 a += da;
714 }
715 break;
716
717 case SLMATH_EXPM1:
718 for (i = 0; i < na; i++)
719 {
720 b[i] = EXPM1_FUNC (to_double ((VOID_STAR) a));
721 a += da;
722 }
723 break;
724
725 case SLMATH_LOG1P:
726 for (i = 0; i < na; i++)
727 {
728 b[i] = LOG1P_FUNC (to_double ((VOID_STAR) a));
729 a += da;
730 }
731 break;
732 }
733
734 return 1;
735 }
736
737 #if SLANG_HAS_COMPLEX
complex_math_op_result(int op,SLtype a,SLtype * b)738 static int complex_math_op_result (int op, SLtype a, SLtype *b)
739 {
740 (void) a;
741 switch (op)
742 {
743 default:
744 *b = SLANG_COMPLEX_TYPE;
745 break;
746
747 case SLMATH_ISINF:
748 case SLMATH_ISNAN:
749 *b = SLANG_CHAR_TYPE;
750 break;
751
752 case SLMATH_REAL:
753 case SLMATH_IMAG:
754 *b = SLANG_DOUBLE_TYPE;
755 break;
756 }
757 return 1;
758 }
759
complex_math_op(int op,SLtype type,VOID_STAR ap,SLuindex_Type na,VOID_STAR bp)760 static int complex_math_op (int op,
761 SLtype type, VOID_STAR ap, SLuindex_Type na,
762 VOID_STAR bp)
763 {
764 double *a, *b;
765 SLuindex_Type i;
766 SLuindex_Type na2 = na * 2;
767 double *(*fun) (double *, double *);
768 char *c;
769
770 (void) type;
771 a = (double *) ap;
772 b = (double *) bp;
773
774 switch (op)
775 {
776 default:
777 return 0;
778
779 case SLMATH_REAL:
780 for (i = 0; i < na; i++)
781 b[i] = a[2 * i];
782 return 1;
783
784 case SLMATH_IMAG:
785 for (i = 0; i < na; i++)
786 b[i] = a[2 * i + 1];
787 return 1;
788
789 case SLMATH_CONJ:
790 for (i = 0; i < na2; i += 2)
791 {
792 b[i] = a[i];
793 b[i+1] = -a[i+1];
794 }
795 return 1;
796
797 case SLMATH_ATANH:
798 fun = SLcomplex_atanh;
799 break;
800 case SLMATH_ACOSH:
801 fun = SLcomplex_acosh;
802 break;
803 case SLMATH_ASINH:
804 fun = SLcomplex_asinh;
805 break;
806 case SLMATH_EXP:
807 fun = SLcomplex_exp;
808 break;
809 case SLMATH_LOG:
810 fun = SLcomplex_log;
811 break;
812 case SLMATH_LOG10:
813 fun = SLcomplex_log10;
814 break;
815 case SLMATH_SQRT:
816 fun = SLcomplex_sqrt;
817 break;
818 case SLMATH_SIN:
819 fun = SLcomplex_sin;
820 break;
821 case SLMATH_COS:
822 fun = SLcomplex_cos;
823 break;
824 case SLMATH_SINH:
825 fun = SLcomplex_sinh;
826 break;
827 case SLMATH_COSH:
828 fun = SLcomplex_cosh;
829 break;
830 case SLMATH_TANH:
831 fun = SLcomplex_tanh;
832 break;
833 case SLMATH_TAN:
834 fun = SLcomplex_tan;
835 break;
836 case SLMATH_ASIN:
837 fun = SLcomplex_asin;
838 break;
839 case SLMATH_ACOS:
840 fun = SLcomplex_acos;
841 break;
842 case SLMATH_ATAN:
843 fun = SLcomplex_atan;
844 break;
845 case SLMATH_ISINF:
846 c = (char *) bp;
847 for (i = 0; i < na; i++)
848 {
849 SLuindex_Type j = 2*i;
850 c[i] = (char) (ISINF_FUN(a[j]) || ISINF_FUN(a[j+1]));
851 }
852 return 1;
853 case SLMATH_ISNAN:
854 c = (char *) bp;
855 for (i = 0; i < na; i++)
856 {
857 SLuindex_Type j = 2*i;
858 c[i] = (char) (ISNAN_FUN(a[j]) || ISNAN_FUN(a[j+1]));
859 }
860 return 1;
861
862 case SLMATH_FLOOR:
863 case SLMATH_ROUND:
864 case SLMATH_CEIL:
865 return double_math_op (op, type, ap, na2, bp);
866 }
867
868 for (i = 0; i < na2; i += 2)
869 (void) (*fun) (b + i, a + i);
870
871 return 1;
872 }
873 #endif
874
875 typedef struct
876 {
877 SLang_Array_Type *at;
878 int is_float; /* if non-zero, use fptr */
879 float f;
880 double d;
881 char c;
882 float *fptr; /* points at at->data or &f */
883 double *dptr;
884 char *cptr;
885 unsigned int inc; /* inc = 0, if at==NULL */
886 SLuindex_Type num;
887 }
888 Array_Or_Scalar_Type;
889
free_array_or_scalar(Array_Or_Scalar_Type * ast)890 static void free_array_or_scalar (Array_Or_Scalar_Type *ast)
891 {
892 if (ast->at != NULL)
893 SLang_free_array (ast->at);
894 }
895
pop_array_or_scalar(Array_Or_Scalar_Type * ast)896 static int pop_array_or_scalar (Array_Or_Scalar_Type *ast)
897 {
898 SLang_Array_Type *at;
899 SLtype dtype;
900
901 ast->at = NULL;
902 ast->inc = 0;
903 ast->num = 1;
904 switch (_pSLang_peek_at_stack2 (&dtype))
905 {
906 case -1:
907 return -1;
908
909 case SLANG_ARRAY_TYPE:
910 if (dtype == SLANG_FLOAT_TYPE)
911 {
912 if (-1 == SLang_pop_array_of_type (&at, SLANG_FLOAT_TYPE))
913 return -1;
914 ast->is_float = 1;
915 ast->fptr = (float *) at->data;
916 }
917 else
918 {
919 if (-1 == SLang_pop_array_of_type (&at, SLANG_DOUBLE_TYPE))
920 return -1;
921 ast->is_float = 0;
922 ast->dptr = (double *) at->data;
923 }
924 ast->inc = 1;
925 ast->num = at->num_elements;
926 ast->at = at;
927 return 0;
928
929 case SLANG_FLOAT_TYPE:
930 ast->is_float = 1;
931 ast->fptr = &ast->f;
932 if (-1 == SLang_pop_float (ast->fptr))
933 return -1;
934 return 0;
935
936 default:
937 ast->is_float = 0;
938 ast->dptr = &ast->d;
939 if (-1 == SLang_pop_double (ast->dptr))
940 return -1;
941 return 0;
942 }
943 }
944
pop_2_arrays_or_scalar(Array_Or_Scalar_Type * a,Array_Or_Scalar_Type * b)945 static int pop_2_arrays_or_scalar (Array_Or_Scalar_Type *a, Array_Or_Scalar_Type *b)
946 {
947 if (-1 == pop_array_or_scalar (b))
948 return -1;
949
950 if (-1 == pop_array_or_scalar (a))
951 {
952 free_array_or_scalar (b);
953 return -1;
954 }
955
956 if ((a->at != NULL) && (b->at != NULL))
957 {
958 if (a->num != b->num)
959 {
960 _pSLang_verror (SL_TypeMismatch_Error, "Array sizes do not match");
961 free_array_or_scalar (a);
962 free_array_or_scalar (b);
963 return -1;
964 }
965 }
966 return 0;
967 }
968
create_from_tmp_array(SLang_Array_Type * a,SLang_Array_Type * b,SLtype type)969 static SLang_Array_Type *create_from_tmp_array (SLang_Array_Type *a, SLang_Array_Type *b, SLtype type)
970 {
971 SLang_Array_Type *c;
972
973 if ((a != NULL) && (a->data_type == type) && (a->num_refs == 1))
974 {
975 a->num_refs += 1;
976 return a;
977 }
978 if ((b != NULL) && (b->data_type == type) && (b->num_refs == 1))
979 {
980 b->num_refs += 1;
981 return b;
982 }
983
984 if (a != NULL)
985 c = a;
986 else
987 c = b;
988
989 return SLang_create_array1 (type, 0, NULL, c->dims, c->num_dims, 1);
990 }
991
do_dd_fun(double (* f)(double,double),Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)992 static int do_dd_fun (double (*f)(double, double),
993 Array_Or_Scalar_Type *a_ast,
994 Array_Or_Scalar_Type *b_ast,
995 Array_Or_Scalar_Type *c_ast)
996 {
997 double *a = a_ast->dptr;
998 double *b = b_ast->dptr;
999 double *c = c_ast->dptr;
1000 unsigned int a_inc = a_ast->inc;
1001 unsigned int b_inc = b_ast->inc;
1002 SLuindex_Type i, n = c_ast->num;
1003
1004 for (i = 0; i < n; i++)
1005 {
1006 c[i] = (*f)(*a, *b);
1007 a += a_inc;
1008 b += b_inc;
1009 }
1010 return 0;
1011 }
1012
do_ff_fun(double (* f)(double,double),Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1013 static int do_ff_fun (double (*f)(double, double),
1014 Array_Or_Scalar_Type *a_ast,
1015 Array_Or_Scalar_Type *b_ast,
1016 Array_Or_Scalar_Type *c_ast)
1017 {
1018 float *a = a_ast->fptr;
1019 float *b = b_ast->fptr;
1020 float *c = c_ast->fptr;
1021 unsigned int a_inc = a_ast->inc;
1022 unsigned int b_inc = b_ast->inc;
1023 SLuindex_Type i, n = c_ast->num;
1024
1025 for (i = 0; i < n; i++)
1026 {
1027 c[i] = (float) (*f)(*a, *b);
1028 a += a_inc;
1029 b += b_inc;
1030 }
1031 return 0;
1032 }
1033
do_fd_fun(double (* f)(double,double),Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1034 static int do_fd_fun (double (*f)(double, double),
1035 Array_Or_Scalar_Type *a_ast,
1036 Array_Or_Scalar_Type *b_ast,
1037 Array_Or_Scalar_Type *c_ast)
1038 {
1039 float *a = a_ast->fptr;
1040 double *b = b_ast->dptr;
1041 double *c = c_ast->dptr;
1042 unsigned int a_inc = a_ast->inc;
1043 unsigned int b_inc = b_ast->inc;
1044 SLuindex_Type i, n = c_ast->num;
1045
1046 for (i = 0; i < n; i++)
1047 {
1048 c[i] = (*f)(*a, *b);
1049 a += a_inc;
1050 b += b_inc;
1051 }
1052 return 0;
1053 }
1054
do_df_fun(double (* f)(double,double),Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1055 static int do_df_fun (double (*f)(double, double),
1056 Array_Or_Scalar_Type *a_ast,
1057 Array_Or_Scalar_Type *b_ast,
1058 Array_Or_Scalar_Type *c_ast)
1059 {
1060 double *a = a_ast->dptr;
1061 float *b = b_ast->fptr;
1062 double *c = c_ast->dptr;
1063 unsigned int a_inc = a_ast->inc;
1064 unsigned int b_inc = b_ast->inc;
1065 SLuindex_Type i, n = c_ast->num;
1066
1067 for (i = 0; i < n; i++)
1068 {
1069 c[i] = (*f)(*a, *b);
1070 a += a_inc;
1071 b += b_inc;
1072 }
1073 return 0;
1074 }
1075
do_binary_function(double (* f)(double,double))1076 static int do_binary_function (double (*f)(double, double))
1077 {
1078 SLtype type;
1079 Array_Or_Scalar_Type a_ast, b_ast, c_ast;
1080
1081 if (-1 == pop_2_arrays_or_scalar (&a_ast, &b_ast))
1082 return -1;
1083
1084 c_ast.is_float = (a_ast.is_float && b_ast.is_float);
1085 c_ast.at = NULL;
1086 c_ast.num = 1;
1087 c_ast.inc = 0;
1088 if (c_ast.is_float)
1089 {
1090 type = SLANG_FLOAT_TYPE;
1091 c_ast.fptr = &c_ast.f;
1092 }
1093 else
1094 {
1095 type = SLANG_DOUBLE_TYPE;
1096 c_ast.dptr = &c_ast.d;
1097 }
1098
1099 if ((a_ast.at != NULL) || (b_ast.at != NULL))
1100 {
1101 if (NULL == (c_ast.at = create_from_tmp_array (a_ast.at, b_ast.at, type)))
1102 {
1103 free_array_or_scalar (&a_ast);
1104 free_array_or_scalar (&b_ast);
1105 return -1;
1106 }
1107 c_ast.fptr = (float *) c_ast.at->data;
1108 c_ast.dptr = (double *) c_ast.at->data;
1109 c_ast.num = c_ast.at->num_elements;
1110 c_ast.inc = 1;
1111 }
1112
1113 if (a_ast.is_float)
1114 {
1115 if (b_ast.is_float)
1116 (void) do_ff_fun (f, &a_ast, &b_ast, &c_ast);
1117 else
1118 (void) do_fd_fun (f, &a_ast, &b_ast, &c_ast);
1119 }
1120 else if (b_ast.is_float)
1121 (void) do_df_fun (f, &a_ast, &b_ast, &c_ast);
1122 else
1123 (void) do_dd_fun (f, &a_ast, &b_ast, &c_ast);
1124
1125 free_array_or_scalar (&a_ast);
1126 free_array_or_scalar (&b_ast);
1127
1128 if (c_ast.at != NULL)
1129 return SLang_push_array (c_ast.at, 1);
1130
1131 if (c_ast.is_float)
1132 return SLang_push_float (c_ast.f);
1133
1134 return SLang_push_double (c_ast.d);
1135 }
1136
1137 /* num is assumed to be at least 2 here */
do_binary_function_on_nargs(double (* f)(double,double),int num)1138 static int do_binary_function_on_nargs (double (*f)(double, double), int num)
1139 {
1140 int depth = SLstack_depth () - num;
1141
1142 num--; /* first time consumes 2 args */
1143 while (num > 0)
1144 {
1145 if (-1 == do_binary_function (f))
1146 {
1147 num = SLstack_depth () - depth;
1148 if (num > 0)
1149 (void) SLdo_pop_n ((unsigned int) num);
1150 return -1;
1151 }
1152 num--;
1153 }
1154 return 0;
1155 }
1156
hypot_fun(void)1157 static void hypot_fun (void)
1158 {
1159 Array_Or_Scalar_Type ast;
1160 SLuindex_Type num;
1161
1162 if (SLang_Num_Function_Args >= 2)
1163 {
1164 (void) do_binary_function_on_nargs (SLmath_hypot, SLang_Num_Function_Args);
1165 return;
1166 }
1167
1168 if (-1 == pop_array_or_scalar (&ast))
1169 return;
1170
1171 if (0 == (num = ast.num))
1172 {
1173 SLang_verror (SL_InvalidParm_Error, "An empty array was passed to hypot");
1174 free_array_or_scalar (&ast);
1175 return;
1176 }
1177
1178 if (ast.is_float)
1179 {
1180 float *f = ast.fptr;
1181 double sum, esum, max;
1182 SLuindex_Type i, imax;
1183
1184 max = fabs((double)*f);
1185 imax = 0;
1186 for (i = 1; i < num; i++)
1187 {
1188 double max1 = fabs((double)f[i]);
1189 if (max1 > max)
1190 {
1191 imax = i;
1192 max = max1;
1193 }
1194 }
1195 sum = 0.0; esum = 0.0;
1196 if (max > 0.0)
1197 {
1198 for (i = 0; i < imax; i++)
1199 {
1200 double term = f[i]/max;
1201 double new_sum;
1202
1203 term = term*term - esum;
1204 new_sum = sum + term;
1205 esum = (new_sum - sum) - term;
1206 sum = new_sum;
1207 }
1208 for (i = imax+1; i < num; i++)
1209 {
1210 double term = f[i]/max;
1211 double new_sum;
1212
1213 term = term*term - esum;
1214 new_sum = sum + term;
1215 esum = (new_sum - sum) - term;
1216 sum = new_sum;
1217 }
1218 }
1219 (void) SLang_push_float ((float) max*sqrt(1.0+sum));
1220 }
1221 else
1222 {
1223 unsigned int i, imax;
1224 double *d = ast.dptr;
1225 double sum, esum, max;
1226 max = fabs(*d);
1227 imax = 0;
1228 for (i = 1; i < num; i++)
1229 {
1230 double max1 = fabs(d[i]);
1231 if (max1 > max)
1232 {
1233 max = max1;
1234 imax = i;
1235 }
1236 }
1237 sum = 0.0; esum = 0.0;
1238 if (max > 0.0)
1239 {
1240 for (i = 0; i < imax; i++)
1241 {
1242 double term = d[i]/max;
1243 double new_sum;
1244
1245 term = term*term - esum;
1246 new_sum = sum + term;
1247 esum = (new_sum - sum) - term;
1248 sum = new_sum;
1249 }
1250 for (i = imax+1; i < num; i++)
1251 {
1252 double term = d[i]/max;
1253 double new_sum;
1254
1255 term = term*term - esum;
1256 new_sum = sum + term;
1257 esum = (new_sum - sum) - term;
1258 sum = new_sum;
1259 }
1260 }
1261 (void) SLang_push_double (max*sqrt(1.0+sum));
1262 }
1263 free_array_or_scalar (&ast);
1264 }
1265
math_poly(void)1266 static void math_poly (void)
1267 {
1268 Array_Or_Scalar_Type ast;
1269 SLang_Array_Type *coeff_at, *y_at;
1270 SLuindex_Type i, k, n;
1271 SLuindex_Type num;
1272 double *a;
1273 double x, y;
1274 int use_factorial = 0;
1275
1276 switch (SLang_Num_Function_Args)
1277 {
1278 case 3:
1279 if (-1 == SLang_pop_int (&use_factorial))
1280 return;
1281 /* drop */
1282 case 2:
1283 break;
1284
1285 default:
1286 SLang_verror (SL_Usage_Error, "\
1287 Usage: y = polynom([a0,a1,...], x [,use_factorial_form])");
1288 return;
1289 }
1290
1291 if (-1 == pop_array_or_scalar (&ast))
1292 return;
1293
1294 if (-1 == SLang_pop_array_of_type (&coeff_at, SLANG_DOUBLE_TYPE))
1295 {
1296 free_array_or_scalar (&ast);
1297 return;
1298 }
1299 a = (double *) coeff_at->data;
1300 n = coeff_at->num_elements;
1301
1302 if (ast.inc == 0)
1303 {
1304 if (ast.is_float)
1305 x = (double) ast.f;
1306 else
1307 x = ast.d;
1308
1309 y = 0.0;
1310 k = n;
1311
1312 if (use_factorial)
1313 {
1314 while (k != 0)
1315 {
1316 y = a[k-1] + (x/k)*y;
1317 k--;
1318 }
1319 }
1320 else while (k != 0)
1321 {
1322 k--;
1323 y = a[k] + x*y;
1324 }
1325
1326 if (ast.is_float)
1327 (void) SLang_push_float ((float) y);
1328 else
1329 (void) SLang_push_double (y);
1330
1331 goto free_and_return;
1332 }
1333
1334 if (NULL == (y_at = create_from_tmp_array (ast.at, NULL, ast.at->data_type)))
1335 goto free_and_return;
1336
1337 num = ast.num;
1338
1339 if (ast.is_float)
1340 {
1341 float *f = ast.fptr;
1342 float *yf = (float *)y_at->data;
1343
1344 for (i = 0; i < num; i++)
1345 {
1346 x = (double) f[i];
1347 y = 0.0;
1348 k = n;
1349 if (use_factorial)
1350 {
1351 while (k != 0)
1352 {
1353 y = a[k-1] + (x/k)*y;
1354 k--;
1355 }
1356 }
1357 else while (k != 0)
1358 {
1359 k--;
1360 y = a[k] + x*y;
1361 }
1362 yf[i] = (float) y;
1363 }
1364 }
1365 else
1366 {
1367 double *d = ast.dptr;
1368 double *yd = (double *)y_at->data;
1369
1370 for (i = 0; i < num; i++)
1371 {
1372 x = d[i];
1373 y = 0.0;
1374 k = n;
1375 if (use_factorial)
1376 {
1377 while (k != 0)
1378 {
1379 y = a[k-1] + (x/k)*y;
1380 k--;
1381 }
1382 }
1383 else while (k != 0)
1384 {
1385 k--;
1386 y = a[k] + x*y;
1387 }
1388 yd[i] = y;
1389 }
1390 }
1391
1392 (void) SLang_push_array (y_at, 1);
1393 /* drop */
1394 free_and_return:
1395 free_array_or_scalar (&ast);
1396 SLang_free_array (coeff_at);
1397 }
1398
atan2_fun(void)1399 static void atan2_fun (void)
1400 {
1401 (void) do_binary_function (atan2);
1402 }
1403
do_min(double a,double b)1404 static double do_min (double a, double b)
1405 {
1406 if ((a >= b) || ISNAN_FUN(a))
1407 return b;
1408 return a;
1409 }
do_max(double a,double b)1410 static double do_max (double a, double b)
1411 {
1412 if ((a <= b) || ISNAN_FUN(a))
1413 return b;
1414 return a;
1415 }
do_diff(double a,double b)1416 static double do_diff (double a, double b)
1417 {
1418 return fabs(a-b);
1419 }
1420
min_fun(void)1421 static void min_fun (void)
1422 {
1423 int num = SLang_Num_Function_Args;
1424
1425 if (num <= 0)
1426 {
1427 SLang_verror (SL_Usage_Error, "_min: Expecting at least one argument");
1428 return;
1429 }
1430
1431 if (num > 1)
1432 (void) do_binary_function_on_nargs (do_min, num);
1433 }
1434
max_fun(void)1435 static void max_fun (void)
1436 {
1437 int num = SLang_Num_Function_Args;
1438
1439 if (num <= 0)
1440 {
1441 SLang_verror (SL_Usage_Error, "_max: Expecting at least one argument");
1442 return;
1443 }
1444
1445 if (num > 1)
1446 (void) do_binary_function_on_nargs (do_max, num);
1447 }
1448
diff_fun(void)1449 static void diff_fun (void)
1450 {
1451 (void) do_binary_function (do_diff);
1452 }
1453
1454 #if 1
1455
do_c_dd_fun(int (* f)(double,double,VOID_STAR),VOID_STAR cd,Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1456 static int do_c_dd_fun (int (*f)(double, double, VOID_STAR),
1457 VOID_STAR cd,
1458 Array_Or_Scalar_Type *a_ast,
1459 Array_Or_Scalar_Type *b_ast,
1460 Array_Or_Scalar_Type *c_ast)
1461 {
1462 double *a = a_ast->dptr;
1463 double *b = b_ast->dptr;
1464 char *c = c_ast->cptr;
1465 unsigned int a_inc = a_ast->inc;
1466 unsigned int b_inc = b_ast->inc;
1467 SLuindex_Type i, n = c_ast->num;
1468
1469 for (i = 0; i < n; i++)
1470 {
1471 c[i] = (char) (*f)(*a, *b, cd);
1472 a += a_inc;
1473 b += b_inc;
1474 }
1475 return 0;
1476 }
1477
do_c_ff_fun(int (* f)(double,double,VOID_STAR),VOID_STAR cd,Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1478 static int do_c_ff_fun (int (*f)(double, double, VOID_STAR),
1479 VOID_STAR cd,
1480 Array_Or_Scalar_Type *a_ast,
1481 Array_Or_Scalar_Type *b_ast,
1482 Array_Or_Scalar_Type *c_ast)
1483 {
1484 float *a = a_ast->fptr;
1485 float *b = b_ast->fptr;
1486 char *c = c_ast->cptr;
1487 unsigned int a_inc = a_ast->inc;
1488 unsigned int b_inc = b_ast->inc;
1489 SLuindex_Type i, n = c_ast->num;
1490
1491 for (i = 0; i < n; i++)
1492 {
1493 c[i] = (char) (*f)(*a, *b, cd);
1494 a += a_inc;
1495 b += b_inc;
1496 }
1497 return 0;
1498 }
1499
do_c_fd_fun(int (* f)(double,double,VOID_STAR),VOID_STAR cd,Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1500 static int do_c_fd_fun (int (*f)(double, double, VOID_STAR),
1501 VOID_STAR cd,
1502 Array_Or_Scalar_Type *a_ast,
1503 Array_Or_Scalar_Type *b_ast,
1504 Array_Or_Scalar_Type *c_ast)
1505 {
1506 float *a = a_ast->fptr;
1507 double *b = b_ast->dptr;
1508 char *c = c_ast->cptr;
1509 unsigned int a_inc = a_ast->inc;
1510 unsigned int b_inc = b_ast->inc;
1511 SLuindex_Type i, n = c_ast->num;
1512
1513 for (i = 0; i < n; i++)
1514 {
1515 c[i] = (char) (*f)(*a, *b, cd);
1516 a += a_inc;
1517 b += b_inc;
1518 }
1519 return 0;
1520 }
1521
do_c_df_fun(int (* f)(double,double,VOID_STAR),VOID_STAR cd,Array_Or_Scalar_Type * a_ast,Array_Or_Scalar_Type * b_ast,Array_Or_Scalar_Type * c_ast)1522 static int do_c_df_fun (int (*f)(double, double, VOID_STAR),
1523 VOID_STAR cd,
1524 Array_Or_Scalar_Type *a_ast,
1525 Array_Or_Scalar_Type *b_ast,
1526 Array_Or_Scalar_Type *c_ast)
1527 {
1528 double *a = a_ast->dptr;
1529 float *b = b_ast->fptr;
1530 char *c = c_ast->cptr;
1531 unsigned int a_inc = a_ast->inc;
1532 unsigned int b_inc = b_ast->inc;
1533 SLuindex_Type i, n = c_ast->num;
1534
1535 for (i = 0; i < n; i++)
1536 {
1537 c[i] = (char) (*f)(*a, *b, cd);
1538 a += a_inc;
1539 b += b_inc;
1540 }
1541 return 0;
1542 }
1543
do_binary_function_c(int (* f)(double,double,VOID_STAR),VOID_STAR cd)1544 static int do_binary_function_c (int (*f)(double, double,VOID_STAR), VOID_STAR cd)
1545 {
1546 Array_Or_Scalar_Type a_ast, b_ast, c_ast;
1547
1548 if (-1 == pop_2_arrays_or_scalar (&a_ast, &b_ast))
1549 return -1;
1550
1551 c_ast.at = NULL;
1552 c_ast.num = 1;
1553 c_ast.inc = 0;
1554 c_ast.cptr = &c_ast.c;
1555
1556 if ((a_ast.at != NULL) || (b_ast.at != NULL))
1557 {
1558 if (a_ast.at != NULL)
1559 c_ast.at = SLang_create_array1 (SLANG_CHAR_TYPE, 0, NULL, a_ast.at->dims, a_ast.at->num_dims, 1);
1560 else
1561 c_ast.at = SLang_create_array1 (SLANG_CHAR_TYPE, 0, NULL, b_ast.at->dims, b_ast.at->num_dims, 1);
1562
1563 if (c_ast.at == NULL)
1564 {
1565 free_array_or_scalar (&a_ast);
1566 free_array_or_scalar (&b_ast);
1567 return -1;
1568 }
1569 c_ast.cptr = (char *) c_ast.at->data;
1570 c_ast.num = c_ast.at->num_elements;
1571 c_ast.inc = 1;
1572 }
1573
1574 if (a_ast.is_float)
1575 {
1576 if (b_ast.is_float)
1577 (void) do_c_ff_fun (f, cd, &a_ast, &b_ast, &c_ast);
1578 else
1579 (void) do_c_fd_fun (f, cd, &a_ast, &b_ast, &c_ast);
1580 }
1581 else if (b_ast.is_float)
1582 (void) do_c_df_fun (f, cd, &a_ast, &b_ast, &c_ast);
1583 else
1584 (void) do_c_dd_fun (f, cd, &a_ast, &b_ast, &c_ast);
1585
1586 free_array_or_scalar (&a_ast);
1587 free_array_or_scalar (&b_ast);
1588
1589 if (c_ast.at != NULL)
1590 return SLang_push_array (c_ast.at, 1);
1591
1592 return SLang_push_char (c_ast.c);
1593 }
1594
1595 typedef struct
1596 {
1597 double relerr;
1598 double abserr;
1599 }
1600 Feqs_Err_Type;
1601
get_tolerances(int nargs,Feqs_Err_Type * ep)1602 static int get_tolerances (int nargs, Feqs_Err_Type *ep)
1603 {
1604 switch (nargs)
1605 {
1606 case 2:
1607 if ((-1 == SLang_pop_double (&ep->abserr))
1608 || (-1 == SLang_pop_double (&ep->relerr)))
1609 return -1;
1610 break;
1611
1612 case 1:
1613 if (-1 == SLang_pop_double (&ep->relerr))
1614 return -1;
1615 ep->abserr = 0.0;
1616 break;
1617
1618 default:
1619 ep->relerr = 0.01;
1620 ep->abserr = 1e-6;
1621 break;
1622 }
1623 return 0;
1624 }
1625
do_feqs(double a,double b,VOID_STAR cd)1626 static int do_feqs (double a, double b, VOID_STAR cd)
1627 {
1628 Feqs_Err_Type *e = (Feqs_Err_Type *)cd;
1629
1630 if (fabs (a-b) <= e->abserr)
1631 return 1;
1632
1633 if (fabs(a) > fabs(b))
1634 return fabs ((a-b)/a) <= e->relerr;
1635
1636 return fabs ((b-a)/b) <= e->relerr;
1637 }
1638
do_fneqs(double a,double b,VOID_STAR cd)1639 static int do_fneqs (double a, double b, VOID_STAR cd)
1640 {
1641 return !do_feqs (a, b, cd);
1642 }
1643
do_fleqs(double a,double b,VOID_STAR cd)1644 static int do_fleqs (double a, double b, VOID_STAR cd)
1645 {
1646 return (a < b) || do_feqs (a, b, cd);
1647 }
1648
do_fgeqs(double a,double b,VOID_STAR cd)1649 static int do_fgeqs (double a, double b, VOID_STAR cd)
1650 {
1651 return (a > b) || do_feqs (a, b, cd);
1652 }
1653
do_an_feqs_fun(int (* f)(double,double,VOID_STAR))1654 static void do_an_feqs_fun (int (*f)(double, double, VOID_STAR))
1655 {
1656 Feqs_Err_Type e;
1657 if (-1 == get_tolerances (SLang_Num_Function_Args-2, &e))
1658 return;
1659
1660 do_binary_function_c (f, (VOID_STAR)&e);
1661 }
1662
feqs_fun(void)1663 static void feqs_fun (void)
1664 {
1665 do_an_feqs_fun (do_feqs);
1666 }
1667
fgeqs_fun(void)1668 static void fgeqs_fun (void)
1669 {
1670 do_an_feqs_fun (do_fgeqs);
1671 }
1672
fleqs_fun(void)1673 static void fleqs_fun (void)
1674 {
1675 do_an_feqs_fun (do_fleqs);
1676 }
1677
fneqs_fun(void)1678 static void fneqs_fun (void)
1679 {
1680 do_an_feqs_fun (do_fneqs);
1681 }
1682 #endif
1683
do_nint(double x)1684 static int do_nint (double x)
1685 {
1686 double xf, xi;
1687
1688 xf = modf (x, &xi); /* x = xi + xf */
1689 if (x >= 0)
1690 {
1691 if (xf >= 0.5)
1692 return xi + 1;
1693 }
1694 else
1695 {
1696 if (xf <= -0.5)
1697 return xi - 1;
1698 }
1699 return xi;
1700 }
1701
float_to_nint(SLang_Array_Type * at,SLang_Array_Type * bt)1702 static int float_to_nint (SLang_Array_Type *at, SLang_Array_Type *bt)
1703 {
1704 SLuindex_Type n, i;
1705 int *ip;
1706 float *fp;
1707
1708 fp = (float *) at->data;
1709 ip = (int *) bt->data;
1710 n = at->num_elements;
1711
1712 for (i = 0; i < n; i++)
1713 ip[i] = do_nint ((double) fp[i]);
1714
1715 return 0;
1716 }
1717
double_to_nint(SLang_Array_Type * at,SLang_Array_Type * bt)1718 static int double_to_nint (SLang_Array_Type *at, SLang_Array_Type *bt)
1719 {
1720 SLuindex_Type n, i;
1721 int *ip;
1722 double *dp;
1723
1724 dp = (double *) at->data;
1725 ip = (int *) bt->data;
1726 n = at->num_elements;
1727
1728 for (i = 0; i < n; i++)
1729 ip[i] = do_nint (dp[i]);
1730
1731 return 0;
1732 }
1733
nint_intrin(void)1734 static void nint_intrin (void)
1735 {
1736 double x;
1737 SLang_Array_Type *at, *bt;
1738 int (*at_to_int_fun)(SLang_Array_Type *, SLang_Array_Type *);
1739 SLtype dtype;
1740
1741 if (SLANG_ARRAY_TYPE != _pSLang_peek_at_stack2 (&dtype))
1742 {
1743 if (-1 == SLang_pop_double (&x))
1744 return;
1745 (void) SLang_push_int (do_nint (x));
1746 return;
1747 }
1748 switch (dtype)
1749 {
1750 case SLANG_INT_TYPE:
1751 return;
1752
1753 case SLANG_FLOAT_TYPE:
1754 if (-1 == SLang_pop_array_of_type (&at, SLANG_FLOAT_TYPE))
1755 return;
1756 at_to_int_fun = float_to_nint;
1757 break;
1758
1759 case SLANG_DOUBLE_TYPE:
1760 default:
1761 if (-1 == SLang_pop_array_of_type (&at, SLANG_DOUBLE_TYPE))
1762 return;
1763 at_to_int_fun = double_to_nint;
1764 break;
1765 }
1766
1767 if (NULL == (bt = SLang_create_array1 (SLANG_INT_TYPE, 0, NULL, at->dims, at->num_dims, 1)))
1768 {
1769 SLang_free_array (at);
1770 return;
1771 }
1772 if (0 == (*at_to_int_fun) (at, bt))
1773 (void) SLang_push_array (bt, 0);
1774
1775 SLang_free_array (bt);
1776 SLang_free_array (at);
1777 }
1778
1779 #ifdef HAVE_FREXP
1780 # define FREXP_FUN(x,e) frexp(x,e)
1781 # ifdef HAVE_FREXPF
1782 # define FREXPF_FUN(x,e) frexpf(x,e)
1783 # else
1784 # define FREXPF_FUN(x,e) (float)FREXP_FUN(x,e)
1785 # endif
1786
frexp_intrin(void)1787 static void frexp_intrin (void)
1788 {
1789 double d;
1790 float f;
1791 int e, *ep;
1792 SLuindex_Type i, imax;
1793 SLang_Array_Type *at, *bt, *et;
1794 SLtype dtype;
1795
1796 switch (_pSLang_peek_at_stack2 (&dtype))
1797 {
1798 case SLANG_ARRAY_TYPE:
1799 switch (dtype)
1800 {
1801 case SLANG_FLOAT_TYPE:
1802 if (-1 == SLang_pop_array_of_type (&at, SLANG_FLOAT_TYPE))
1803 return;
1804 break;
1805
1806 default:
1807 case SLANG_DOUBLE_TYPE:
1808 if (-1 == SLang_pop_array_of_type (&at, SLANG_DOUBLE_TYPE))
1809 return;
1810 break;
1811 }
1812 break;
1813
1814 case SLANG_FLOAT_TYPE:
1815 if (-1 == SLang_pop_float (&f))
1816 return;
1817 f = FREXPF_FUN(f, &e);
1818 (void) SLang_push_float (f);
1819 (void) SLang_push_int (e);
1820 return;
1821
1822 default:
1823 case SLANG_DOUBLE_TYPE:
1824 if (-1 == SLang_pop_double (&d))
1825 return;
1826 d = FREXP_FUN(d, &e);
1827 (void) SLang_push_double (d);
1828 (void) SLang_push_int (e);
1829 return;
1830 }
1831
1832 if (NULL == (bt = SLang_create_array1 (at->data_type, 0, NULL, at->dims, at->num_dims, 1)))
1833 {
1834 SLang_free_array (at);
1835 return;
1836 }
1837 if (NULL == (et = SLang_create_array1 (SLANG_INT_TYPE, 0, NULL, at->dims, at->num_dims, 1)))
1838 {
1839 SLang_free_array (at);
1840 SLang_free_array (bt);
1841 return;
1842 }
1843
1844 imax = at->num_elements;
1845 ep = (int *)et->data;
1846
1847 if (at->data_type == SLANG_DOUBLE_TYPE)
1848 {
1849 double *a = (double *)at->data;
1850 double *b = (double *)bt->data;
1851 for (i = 0; i < imax; i++)
1852 {
1853 b[i] = FREXP_FUN(a[i], ep+i);
1854 }
1855 }
1856 else
1857 {
1858 float *a = (float *)at->data;
1859 float *b = (float *)bt->data;
1860 for (i = 0; i < imax; i++)
1861 {
1862 b[i] = FREXPF_FUN(a[i], ep+i);
1863 }
1864 }
1865
1866 (void) SLang_push_array (bt, 0);
1867 (void) SLang_push_array (et, 0);
1868 SLang_free_array (et);
1869 SLang_free_array (bt);
1870 SLang_free_array (at);
1871 }
1872 #endif /* HAVE_FREXP */
1873
1874 #ifdef HAVE_LDEXP
1875 # define LDEXP_FUN(a,b) ldexp(a,b)
1876 # ifdef HAVE_LDEXPF
1877 # define LDEXPF_FUN(a,b) ldexpf(a,b)
1878 # else
1879 # define LDEXPF_FUN(a,b) (float)LDEXP_FUN(a,b)
1880 # endif
ldexp_intrin(void)1881 static void ldexp_intrin (void)
1882 {
1883 Array_Or_Scalar_Type ast;
1884 SLang_Array_Type *e_at = NULL, *c_at;
1885 int *e_ptr, e_data;
1886 SLuindex_Type i, imax;
1887
1888 if (SLang_peek_at_stack () == SLANG_ARRAY_TYPE)
1889 {
1890 if (-1 == SLang_pop_array_of_type (&e_at, SLANG_INT_TYPE))
1891 return;
1892 e_ptr = (int *)e_at->data;
1893 }
1894 else
1895 {
1896 if (-1 == SLang_pop_int (&e_data))
1897 return;
1898 e_ptr = &e_data;
1899 }
1900
1901 if (-1 == pop_array_or_scalar (&ast))
1902 {
1903 if (e_at != NULL)
1904 SLang_free_array (e_at);
1905 return;
1906 }
1907
1908 if ((e_at == NULL) && (ast.at == NULL))
1909 {
1910 /* Scalar case */
1911 if (ast.is_float)
1912 (void) SLang_push_float (LDEXPF_FUN(ast.f, *e_ptr));
1913 else
1914 (void) SLang_push_double (LDEXP_FUN(ast.d, *e_ptr));
1915 /* free_array_or_scalar (&ast); Nothing to free */
1916 return;
1917 }
1918
1919 if (NULL == (c_at = create_from_tmp_array (ast.at, e_at, ast.is_float ? SLANG_FLOAT_TYPE : SLANG_DOUBLE_TYPE)))
1920 {
1921 free_array_or_scalar (&ast);
1922 SLang_free_array (e_at);
1923 return;
1924 }
1925
1926 if (e_at == NULL)
1927 {
1928 int e = *e_ptr;
1929 imax = ast.num;
1930
1931 if (ast.is_float)
1932 {
1933 float *c, *a;
1934 a = ast.fptr;
1935 c = (float *)c_at->data;
1936 for (i = 0; i < imax; i++)
1937 c[i] = LDEXPF_FUN(a[i], e);
1938 }
1939 else
1940 {
1941 double *c, *a;
1942 a = ast.dptr;
1943 c = (double *)c_at->data;
1944 for (i = 0; i < imax; i++)
1945 c[i] = LDEXP_FUN(a[i], e);
1946 }
1947 goto push_free_and_return;
1948 }
1949
1950 if (ast.at == NULL)
1951 {
1952 imax = e_at->num_elements;
1953
1954 if (ast.is_float)
1955 {
1956 float *c, a;
1957 a = ast.f;
1958 c = (float *)c_at->data;
1959 for (i = 0; i < imax; i++)
1960 c[i] = LDEXPF_FUN(a, e_ptr[i]);
1961 }
1962 else
1963 {
1964 double *c, a;
1965 a = ast.d;
1966 c = (double *)c_at->data;
1967 for (i = 0; i < imax; i++)
1968 c[i] = LDEXP_FUN(a, e_ptr[i]);
1969 }
1970 goto push_free_and_return;
1971 }
1972
1973 if (e_at->num_elements != ast.num)
1974 {
1975 SLang_verror (SL_TypeMismatch_Error, "ldexp: Array sizes do not match");
1976 goto free_and_return;
1977 }
1978
1979 imax = ast.num;
1980 if (ast.is_float)
1981 {
1982 float *c = (float *)c_at->data;
1983 float *a = ast.fptr;
1984 for (i = 0; i < imax; i++)
1985 c[i] = LDEXPF_FUN(a[i], e_ptr[i]);
1986 }
1987 else
1988 {
1989 double *c = (double *)c_at->data;
1990 double *a = ast.dptr;
1991 for (i = 0; i < imax; i++)
1992 c[i] = LDEXP_FUN(a[i], e_ptr[i]);
1993 }
1994 /* drop */
1995 push_free_and_return:
1996 (void) SLang_push_array (c_at, 0);
1997 /* drop */
1998 free_and_return:
1999 if (e_at != NULL) SLang_free_array (e_at);
2000 SLang_free_array (c_at);
2001 free_array_or_scalar (&ast);
2002 }
2003 #endif /* HAVE_LDEXPF */
2004
2005
2006 #ifdef HAVE_SINCOSF
2007 # define SINCOSF(x, sp, cp) sincosf (x, sp, cp)
2008 #else
2009 # define SINCOSF(x, sp, cp) \
2010 *(sp) = (float)sin(x); *(cp) = (float)cos(x)
2011 #endif
2012 #ifdef HAVE_SINCOS
2013 # define SINCOS(x, sp, cp) sincos (x, sp, cp)
2014 #else
2015 # define SINCOS(x, sp, cp) \
2016 *(sp) = sin(x); *(cp) = cos(x)
2017 #endif
2018
sincos_intrin(void)2019 static void sincos_intrin (void)
2020 {
2021 Array_Or_Scalar_Type ast;
2022 SLtype type;
2023 SLang_Array_Type *s_at, *c_at;
2024 SLuindex_Type num;
2025
2026 if (-1 == pop_array_or_scalar (&ast))
2027 return;
2028
2029 if (ast.inc == 0)
2030 {
2031 if (ast.is_float)
2032 {
2033 float s, c;
2034 SINCOSF(ast.f, &s, &c);
2035 (void) SLang_push_float (s);
2036 (void) SLang_push_float (c);
2037 }
2038 else
2039 {
2040 double s, c;
2041 SINCOS(ast.d, &s, &c);
2042 (void) SLang_push_double (s);
2043 (void) SLang_push_double (c);
2044 }
2045 free_array_or_scalar (&ast);
2046 return;
2047 }
2048
2049 c_at = s_at = NULL;
2050 num = ast.num;
2051 type = ast.is_float ? SLANG_FLOAT_TYPE : SLANG_DOUBLE_TYPE;
2052
2053 s_at = SLang_create_array1(type, 0, NULL, ast.at->dims, ast.at->num_dims, 1);
2054 if (s_at == NULL)
2055 goto free_and_return;
2056 c_at = SLang_create_array1(type, 0, NULL, ast.at->dims, ast.at->num_dims, 1);
2057 if (c_at == NULL)
2058 goto free_and_return;
2059
2060 if (ast.is_float)
2061 {
2062 SLuindex_Type i;
2063 float *x = ast.fptr;
2064 float *s = (float*) s_at->data;
2065 float *c = (float*) c_at->data;
2066 for (i=0; i<num; i++)
2067 {
2068 SINCOSF(x[i], s+i, c+i);
2069 }
2070 }
2071 else
2072 {
2073 SLuindex_Type i;
2074 double *x = ast.dptr;
2075 double *s = (double*) s_at->data;
2076 double *c = (double*) c_at->data;
2077 for (i=0; i<num; i++)
2078 {
2079 SINCOS (x[i], s+i, c+i);
2080 }
2081 }
2082
2083 if (0 == SLang_push_array (s_at, 0))
2084 (void) SLang_push_array (c_at, 0);
2085
2086 /* drop */
2087 free_and_return:
2088 if (c_at != NULL) SLang_free_array (c_at);
2089 if (s_at != NULL) SLang_free_array (s_at);
2090 free_array_or_scalar (&ast);
2091 }
2092
fpu_clear_except_bits(void)2093 static void fpu_clear_except_bits (void)
2094 {
2095 SLfpu_clear_except_bits ();
2096 }
2097
fpu_test_except_bits(int * bits)2098 static int fpu_test_except_bits (int *bits)
2099 {
2100 return (int) SLfpu_test_except_bits ((unsigned int) *bits);
2101 }
2102
2103 static SLang_DConstant_Type DConst_Table [] =
2104 {
2105 MAKE_DCONSTANT("E", 2.718281828459045),
2106 MAKE_DCONSTANT("PI", PI),
2107 SLANG_END_DCONST_TABLE
2108 };
2109
2110 static SLang_Math_Unary_Type SLmath_Table [] =
2111 {
2112 MAKE_MATH_UNARY("sinh", SLMATH_SINH),
2113 MAKE_MATH_UNARY("asinh", SLMATH_ASINH),
2114 MAKE_MATH_UNARY("cosh", SLMATH_COSH),
2115 MAKE_MATH_UNARY("acosh", SLMATH_ACOSH),
2116 MAKE_MATH_UNARY("tanh", SLMATH_TANH),
2117 MAKE_MATH_UNARY("atanh", SLMATH_ATANH),
2118 MAKE_MATH_UNARY("sin", SLMATH_SIN),
2119 MAKE_MATH_UNARY("cos", SLMATH_COS),
2120 MAKE_MATH_UNARY("tan", SLMATH_TAN),
2121 MAKE_MATH_UNARY("atan", SLMATH_ATAN),
2122 MAKE_MATH_UNARY("acos", SLMATH_ACOS),
2123 MAKE_MATH_UNARY("asin", SLMATH_ASIN),
2124 MAKE_MATH_UNARY("exp", SLMATH_EXP),
2125 MAKE_MATH_UNARY("log", SLMATH_LOG),
2126 MAKE_MATH_UNARY("sqrt", SLMATH_SQRT),
2127 MAKE_MATH_UNARY("log10", SLMATH_LOG10),
2128 MAKE_MATH_UNARY("isinf", SLMATH_ISINF),
2129 MAKE_MATH_UNARY("isnan", SLMATH_ISNAN),
2130 MAKE_MATH_UNARY("floor", SLMATH_FLOOR),
2131 MAKE_MATH_UNARY("ceil", SLMATH_CEIL),
2132 MAKE_MATH_UNARY("round", SLMATH_ROUND),
2133 MAKE_MATH_UNARY("expm1", SLMATH_EXPM1),
2134 MAKE_MATH_UNARY("log1p", SLMATH_LOG1P),
2135
2136 #if SLANG_HAS_COMPLEX
2137 MAKE_MATH_UNARY("Real", SLMATH_REAL),
2138 MAKE_MATH_UNARY("Imag", SLMATH_IMAG),
2139 MAKE_MATH_UNARY("Conj", SLMATH_CONJ),
2140 #endif
2141 SLANG_END_MATH_UNARY_TABLE
2142 };
2143
2144 static SLang_Intrin_Fun_Type SLang_Math_Table [] =
2145 {
2146 MAKE_INTRINSIC_0("nint", nint_intrin, SLANG_VOID_TYPE),
2147 MAKE_INTRINSIC_0("polynom", math_poly, SLANG_VOID_TYPE),
2148 MAKE_INTRINSIC_0("hypot", hypot_fun, SLANG_VOID_TYPE),
2149 MAKE_INTRINSIC_0("atan2", atan2_fun, SLANG_VOID_TYPE),
2150 MAKE_INTRINSIC_0("_min", min_fun, SLANG_VOID_TYPE),
2151 MAKE_INTRINSIC_0("_max", max_fun, SLANG_VOID_TYPE),
2152 MAKE_INTRINSIC_0("_diff", diff_fun, SLANG_VOID_TYPE),
2153 MAKE_INTRINSIC_0("feqs", feqs_fun, SLANG_VOID_TYPE),
2154 MAKE_INTRINSIC_0("fneqs", fneqs_fun, SLANG_VOID_TYPE),
2155 MAKE_INTRINSIC_0("flteqs", fleqs_fun, SLANG_VOID_TYPE),
2156 MAKE_INTRINSIC_0("fgteqs", fgeqs_fun, SLANG_VOID_TYPE),
2157 MAKE_INTRINSIC_0("fpu_clear_except_bits", fpu_clear_except_bits, SLANG_VOID_TYPE),
2158 MAKE_INTRINSIC_1("fpu_test_except_bits", fpu_test_except_bits, _pSLANG_LONG_TYPE, _pSLANG_LONG_TYPE),
2159 #ifdef HAVE_FREXP
2160 MAKE_INTRINSIC_0("frexp", frexp_intrin, SLANG_VOID_TYPE),
2161 #endif
2162 #ifdef HAVE_LDEXP
2163 MAKE_INTRINSIC_0("ldexp", ldexp_intrin, SLANG_VOID_TYPE),
2164 #endif
2165 MAKE_INTRINSIC_0("sincos", sincos_intrin, SLANG_VOID_TYPE),
2166 SLANG_END_INTRIN_FUN_TABLE
2167 };
2168
2169 static SLang_IConstant_Type IConsts [] =
2170 {
2171 MAKE_ICONSTANT("FE_DIVBYZERO", SL_FE_DIVBYZERO),
2172 MAKE_ICONSTANT("FE_INVALID", SL_FE_INVALID),
2173 MAKE_ICONSTANT("FE_OVERFLOW", SL_FE_OVERFLOW),
2174 MAKE_ICONSTANT("FE_UNDERFLOW", SL_FE_UNDERFLOW),
2175 MAKE_ICONSTANT("FE_INEXACT", SL_FE_INEXACT),
2176 MAKE_ICONSTANT("FE_ALL_EXCEPT", SL_FE_ALLEXCEPT),
2177 SLANG_END_ICONST_TABLE
2178 };
2179
add_nan_and_inf(void)2180 static int add_nan_and_inf (void)
2181 {
2182 if ((-1 == SLns_add_dconstant (NULL, "_NaN", _pSLang_NaN))
2183 || (-1 == SLns_add_dconstant (NULL, "_Inf", _pSLang_Inf)))
2184 return -1;
2185
2186 SLfpu_clear_except_bits ();
2187
2188 return 0;
2189 }
2190
SLang_init_slmath(void)2191 int SLang_init_slmath (void)
2192 {
2193 SLtype *int_types;
2194
2195 #if SLANG_HAS_COMPLEX
2196 if (-1 == _pSLinit_slcomplex ())
2197 return -1;
2198 #endif
2199 int_types = _pSLarith_Arith_Types;
2200
2201 while (*int_types != SLANG_FLOAT_TYPE)
2202 {
2203 if (-1 == SLclass_add_math_op (*int_types, generic_math_op, double_math_op_result))
2204 return -1;
2205 int_types++;
2206 }
2207
2208 if ((-1 == SLclass_add_math_op (SLANG_FLOAT_TYPE, float_math_op, double_math_op_result))
2209 || (-1 == SLclass_add_math_op (SLANG_DOUBLE_TYPE, double_math_op, double_math_op_result))
2210 #if SLANG_HAS_COMPLEX
2211 || (-1 == SLclass_add_math_op (SLANG_COMPLEX_TYPE, complex_math_op, complex_math_op_result))
2212 #endif
2213 )
2214 return -1;
2215
2216 if ((-1 == SLadd_math_unary_table (SLmath_Table, "__SLMATH__"))
2217 || (-1 == SLadd_intrin_fun_table (SLang_Math_Table, NULL))
2218 || (-1 == SLadd_dconstant_table (DConst_Table, NULL))
2219 || (-1 == SLadd_iconstant_table (IConsts, NULL))
2220 || (-1 == add_nan_and_inf ()))
2221 return -1;
2222
2223 #if defined(__unix__)
2224 (void) SLsignal (SIGFPE, math_floating_point_exception);
2225 #endif
2226
2227 return 0;
2228 }
2229 #endif /* SLANG_HAS_FLOAT */
2230