1 #include "lib9.h"
2 #include "interp.h"
3 #include "isa.h"
4 #include "runt.h"
5 #include "raise.h"
6 #include "mathi.h"
7 #include "mathmod.h"
8
9 static union
10 {
11 double x;
12 uvlong u;
13 } bits64;
14
15 static union{
16 float x;
17 unsigned int u;
18 } bits32;
19
20 void
mathmodinit(void)21 mathmodinit(void)
22 {
23 builtinmod("$Math", Mathmodtab, Mathmodlen);
24 fmtinstall('g', gfltconv);
25 fmtinstall('G', gfltconv);
26 fmtinstall('e', gfltconv);
27 /* fmtinstall('E', gfltconv); */ /* avoid clash with ether address */
28 fmtinstall(0x00c9, gfltconv); /* L'É' */
29 fmtinstall('f', gfltconv);
30 }
31
32 void
Math_import_int(void * fp)33 Math_import_int(void *fp)
34 {
35 F_Math_import_int *f;
36 int i, n;
37 unsigned int u;
38 unsigned char *bp;
39 int *x;
40
41 f = fp;
42 n = f->x->len;
43 if(f->b->len!=4*n)
44 error(exMathia);
45 bp = (unsigned char *)(f->b->data);
46 x = (int*)(f->x->data);
47 for(i=0; i<n; i++){
48 u = *bp++;
49 u = (u<<8) | *bp++;
50 u = (u<<8) | *bp++;
51 u = (u<<8) | *bp++;
52 x[i] = u;
53 }
54 }
55
56 void
Math_import_real32(void * fp)57 Math_import_real32(void *fp)
58 {
59 F_Math_import_int *f;
60 int i, n;
61 unsigned int u;
62 unsigned char *bp;
63 double *x;
64
65 f = fp;
66 n = f->x->len;
67 if(f->b->len!=4*n)
68 error(exMathia);
69 bp = (unsigned char *)(f->b->data);
70 x = (double*)(f->x->data);
71 for(i=0; i<n; i++){
72 u = *bp++;
73 u = (u<<8) | *bp++;
74 u = (u<<8) | *bp++;
75 u = (u<<8) | *bp++;
76 bits32.u = u;
77 x[i] = bits32.x;
78 }
79 }
80
81 void
Math_import_real(void * fp)82 Math_import_real(void *fp)
83 {
84 F_Math_import_int *f;
85 int i, n;
86 uvlong u;
87 unsigned char *bp;
88 double *x;
89
90 f = fp;
91 n = f->x->len;
92 if(f->b->len!=8*n)
93 error(exMathia);
94 bp = f->b->data;
95 x = (double*)(f->x->data);
96 for(i=0; i<n; i++){
97 u = *bp++;
98 u = (u<<8) | *bp++;
99 u = (u<<8) | *bp++;
100 u = (u<<8) | *bp++;
101 u = (u<<8) | *bp++;
102 u = (u<<8) | *bp++;
103 u = (u<<8) | *bp++;
104 u = (u<<8) | *bp++;
105 bits64.u = u;
106 x[i] = bits64.x;
107 }
108 }
109
110 void
Math_export_int(void * fp)111 Math_export_int(void *fp)
112 {
113 F_Math_export_int *f;
114 int i, n;
115 unsigned int u;
116 unsigned char *bp;
117 int *x;
118
119 f = fp;
120 n = f->x->len;
121 if(f->b->len!=4*n)
122 error(exMathia);
123 bp = (unsigned char *)(f->b->data);
124 x = (int*)(f->x->data);
125 for(i=0; i<n; i++){
126 u = x[i];
127 *bp++ = u>>24;
128 *bp++ = u>>16;
129 *bp++ = u>>8;
130 *bp++ = u;
131 }
132 }
133
134 void
Math_export_real32(void * fp)135 Math_export_real32(void *fp)
136 {
137 F_Math_export_int *f;
138 int i, n;
139 unsigned int u;
140 unsigned char *bp;
141 double *x;
142
143 f = fp;
144 n = f->x->len;
145 if(f->b->len!=4*n)
146 error(exMathia);
147 bp = (unsigned char *)(f->b->data);
148 x = (double*)(f->x->data);
149 for(i=0; i<n; i++){
150 bits32.x = x[i];
151 u = bits32.u;
152 *bp++ = u>>24;
153 *bp++ = u>>16;
154 *bp++ = u>>8;
155 *bp++ = u;
156 }
157 }
158
159 void
Math_export_real(void * fp)160 Math_export_real(void *fp)
161 {
162 F_Math_export_int *f;
163 int i, n;
164 uvlong u;
165 unsigned char *bp;
166 double *x;
167
168 f = fp;
169 n = f->x->len;
170 if(f->b->len!=8*n)
171 error(exMathia);
172 bp = (unsigned char *)(f->b->data);
173 x = (double*)(f->x->data);
174 for(i=0; i<n; i++){
175 bits64.x = x[i];
176 u = bits64.u;
177 *bp++ = u>>56;
178 *bp++ = u>>48;
179 *bp++ = u>>40;
180 *bp++ = u>>32;
181 *bp++ = u>>24;
182 *bp++ = u>>16;
183 *bp++ = u>>8;
184 *bp++ = u;
185 }
186 }
187
188
189 void
Math_bits32real(void * fp)190 Math_bits32real(void *fp)
191 {
192 F_Math_bits32real *f;
193
194 f = fp;
195 bits32.u = f->b;
196 *f->ret = bits32.x;
197 }
198
199 void
Math_bits64real(void * fp)200 Math_bits64real(void *fp)
201 {
202 F_Math_bits64real *f;
203
204 f = fp;
205 bits64.u = f->b;
206 *f->ret = bits64.x;
207 }
208
209 void
Math_realbits32(void * fp)210 Math_realbits32(void *fp)
211 {
212 F_Math_realbits32 *f;
213
214 f = fp;
215 bits32.x = f->x;
216 *f->ret = bits32.u;
217 }
218
219 void
Math_realbits64(void * fp)220 Math_realbits64(void *fp)
221 {
222 F_Math_realbits64 *f;
223
224 f = fp;
225 bits64.x = f->x;
226 *f->ret = bits64.u;
227 }
228
229
230 void
Math_getFPcontrol(void * fp)231 Math_getFPcontrol(void *fp)
232 {
233 F_Math_getFPcontrol *f;
234
235 f = fp;
236
237 *f->ret = getFPcontrol();
238 }
239
240 void
Math_getFPstatus(void * fp)241 Math_getFPstatus(void *fp)
242 {
243 F_Math_getFPstatus *f;
244
245 f = fp;
246
247 *f->ret = getFPstatus();
248 }
249
250 void
Math_finite(void * fp)251 Math_finite(void *fp)
252 {
253 F_Math_finite *f;
254
255 f = fp;
256
257 *f->ret = finite(f->x);
258 }
259
260 void
Math_ilogb(void * fp)261 Math_ilogb(void *fp)
262 {
263 F_Math_ilogb *f;
264
265 f = fp;
266
267 *f->ret = ilogb(f->x);
268 }
269
270 void
Math_isnan(void * fp)271 Math_isnan(void *fp)
272 {
273 F_Math_isnan *f;
274
275 f = fp;
276
277 *f->ret = isNaN(f->x);
278 }
279
280 void
Math_acos(void * fp)281 Math_acos(void *fp)
282 {
283 F_Math_acos *f;
284
285 f = fp;
286
287 *f->ret = __ieee754_acos(f->x);
288 }
289
290 void
Math_acosh(void * fp)291 Math_acosh(void *fp)
292 {
293 F_Math_acosh *f;
294
295 f = fp;
296
297 *f->ret = __ieee754_acosh(f->x);
298 }
299
300 void
Math_asin(void * fp)301 Math_asin(void *fp)
302 {
303 F_Math_asin *f;
304
305 f = fp;
306
307 *f->ret = __ieee754_asin(f->x);
308 }
309
310 void
Math_asinh(void * fp)311 Math_asinh(void *fp)
312 {
313 F_Math_asinh *f;
314
315 f = fp;
316
317 *f->ret = asinh(f->x);
318 }
319
320 void
Math_atan(void * fp)321 Math_atan(void *fp)
322 {
323 F_Math_atan *f;
324
325 f = fp;
326
327 *f->ret = atan(f->x);
328 }
329
330 void
Math_atanh(void * fp)331 Math_atanh(void *fp)
332 {
333 F_Math_atanh *f;
334
335 f = fp;
336
337 *f->ret = __ieee754_atanh(f->x);
338 }
339
340 void
Math_cbrt(void * fp)341 Math_cbrt(void *fp)
342 {
343 F_Math_cbrt *f;
344
345 f = fp;
346
347 *f->ret = cbrt(f->x);
348 }
349
350 void
Math_ceil(void * fp)351 Math_ceil(void *fp)
352 {
353 F_Math_ceil *f;
354
355 f = fp;
356
357 *f->ret = ceil(f->x);
358 }
359
360 void
Math_cos(void * fp)361 Math_cos(void *fp)
362 {
363 F_Math_cos *f;
364
365 f = fp;
366
367 *f->ret = cos(f->x);
368 }
369
370 void
Math_cosh(void * fp)371 Math_cosh(void *fp)
372 {
373 F_Math_cosh *f;
374
375 f = fp;
376
377 *f->ret = __ieee754_cosh(f->x);
378 }
379
380 void
Math_erf(void * fp)381 Math_erf(void *fp)
382 {
383 F_Math_erf *f;
384
385 f = fp;
386
387 *f->ret = erf(f->x);
388 }
389
390 void
Math_erfc(void * fp)391 Math_erfc(void *fp)
392 {
393 F_Math_erfc *f;
394
395 f = fp;
396
397 *f->ret = erfc(f->x);
398 }
399
400 void
Math_exp(void * fp)401 Math_exp(void *fp)
402 {
403 F_Math_exp *f;
404
405 f = fp;
406
407 *f->ret = __ieee754_exp(f->x);
408 }
409
410 void
Math_expm1(void * fp)411 Math_expm1(void *fp)
412 {
413 F_Math_expm1 *f;
414
415 f = fp;
416
417 *f->ret = expm1(f->x);
418 }
419
420 void
Math_fabs(void * fp)421 Math_fabs(void *fp)
422 {
423 F_Math_fabs *f;
424
425 f = fp;
426
427 *f->ret = fabs(f->x);
428 }
429
430 void
Math_floor(void * fp)431 Math_floor(void *fp)
432 {
433 F_Math_floor *f;
434
435 f = fp;
436
437 *f->ret = floor(f->x);
438 }
439
440 void
Math_j0(void * fp)441 Math_j0(void *fp)
442 {
443 F_Math_j0 *f;
444
445 f = fp;
446
447 *f->ret = __ieee754_j0(f->x);
448 }
449
450 void
Math_j1(void * fp)451 Math_j1(void *fp)
452 {
453 F_Math_j1 *f;
454
455 f = fp;
456
457 *f->ret = __ieee754_j1(f->x);
458 }
459
460 void
Math_log(void * fp)461 Math_log(void *fp)
462 {
463 F_Math_log *f;
464
465 f = fp;
466
467 *f->ret = __ieee754_log(f->x);
468 }
469
470 void
Math_log10(void * fp)471 Math_log10(void *fp)
472 {
473 F_Math_log10 *f;
474
475 f = fp;
476
477 *f->ret = __ieee754_log10(f->x);
478 }
479
480 void
Math_log1p(void * fp)481 Math_log1p(void *fp)
482 {
483 F_Math_log1p *f;
484
485 f = fp;
486
487 *f->ret = log1p(f->x);
488 }
489
490 void
Math_rint(void * fp)491 Math_rint(void *fp)
492 {
493 F_Math_rint *f;
494
495 f = fp;
496
497 *f->ret = rint(f->x);
498 }
499
500 void
Math_sin(void * fp)501 Math_sin(void *fp)
502 {
503 F_Math_sin *f;
504
505 f = fp;
506
507 *f->ret = sin(f->x);
508 }
509
510 void
Math_sinh(void * fp)511 Math_sinh(void *fp)
512 {
513 F_Math_sinh *f;
514
515 f = fp;
516
517 *f->ret = __ieee754_sinh(f->x);
518 }
519
520 void
Math_sqrt(void * fp)521 Math_sqrt(void *fp)
522 {
523 F_Math_sqrt *f;
524
525 f = fp;
526
527 *f->ret = __ieee754_sqrt(f->x);
528 }
529
530 void
Math_tan(void * fp)531 Math_tan(void *fp)
532 {
533 F_Math_tan *f;
534
535 f = fp;
536
537 *f->ret = tan(f->x);
538 }
539
540 void
Math_tanh(void * fp)541 Math_tanh(void *fp)
542 {
543 F_Math_tanh *f;
544
545 f = fp;
546
547 *f->ret = tanh(f->x);
548 }
549
550 void
Math_y0(void * fp)551 Math_y0(void *fp)
552 {
553 F_Math_y0 *f;
554
555 f = fp;
556
557 *f->ret = __ieee754_y0(f->x);
558 }
559
560 void
Math_y1(void * fp)561 Math_y1(void *fp)
562 {
563 F_Math_y1 *f;
564
565 f = fp;
566
567 *f->ret = __ieee754_y1(f->x);
568 }
569
570 void
Math_fdim(void * fp)571 Math_fdim(void *fp)
572 {
573 F_Math_fdim *f;
574
575 f = fp;
576
577 *f->ret = fdim(f->x, f->y);
578 }
579
580 void
Math_fmax(void * fp)581 Math_fmax(void *fp)
582 {
583 F_Math_fmax *f;
584
585 f = fp;
586
587 *f->ret = fmax(f->x, f->y);
588 }
589
590 void
Math_fmin(void * fp)591 Math_fmin(void *fp)
592 {
593 F_Math_fmin *f;
594
595 f = fp;
596
597 *f->ret = fmin(f->x, f->y);
598 }
599
600 void
Math_fmod(void * fp)601 Math_fmod(void *fp)
602 {
603 F_Math_fmod *f;
604
605 f = fp;
606
607 *f->ret = __ieee754_fmod(f->x, f->y);
608 }
609
610 void
Math_hypot(void * fp)611 Math_hypot(void *fp)
612 {
613 F_Math_hypot *f;
614
615 f = fp;
616
617 *f->ret = __ieee754_hypot(f->x, f->y);
618 }
619
620 void
Math_nextafter(void * fp)621 Math_nextafter(void *fp)
622 {
623 F_Math_nextafter *f;
624
625 f = fp;
626
627 *f->ret = nextafter(f->x, f->y);
628 }
629
630 void
Math_pow(void * fp)631 Math_pow(void *fp)
632 {
633 F_Math_pow *f;
634
635 f = fp;
636
637 *f->ret = __ieee754_pow(f->x, f->y);
638 }
639
640
641
642 void
Math_FPcontrol(void * fp)643 Math_FPcontrol(void *fp)
644 {
645 F_Math_FPcontrol *f;
646
647 f = fp;
648
649 *f->ret = FPcontrol(f->r, f->mask);
650 }
651
652 void
Math_FPstatus(void * fp)653 Math_FPstatus(void *fp)
654 {
655 F_Math_FPstatus *f;
656
657 f = fp;
658
659 *f->ret = FPstatus(f->r, f->mask);
660 }
661
662 void
Math_atan2(void * fp)663 Math_atan2(void *fp)
664 {
665 F_Math_atan2 *f;
666
667 f = fp;
668
669 *f->ret = __ieee754_atan2(f->y, f->x);
670 }
671
672 void
Math_copysign(void * fp)673 Math_copysign(void *fp)
674 {
675 F_Math_copysign *f;
676
677 f = fp;
678
679 *f->ret = copysign(f->x, f->s);
680 }
681
682 void
Math_jn(void * fp)683 Math_jn(void *fp)
684 {
685 F_Math_jn *f;
686
687 f = fp;
688
689 *f->ret = __ieee754_jn(f->n, f->x);
690 }
691
692 void
Math_lgamma(void * fp)693 Math_lgamma(void *fp)
694 {
695 F_Math_lgamma *f;
696
697 f = fp;
698
699 f->ret->t1 = __ieee754_lgamma_r(f->x, &f->ret->t0);
700 }
701
702 void
Math_modf(void * fp)703 Math_modf(void *fp)
704 {
705 F_Math_modf *f;
706 double ipart;
707
708 f = fp;
709
710 f->ret->t1 = modf(f->x, &ipart);
711 f->ret->t0 = ipart;
712 }
713
714 void
Math_pow10(void * fp)715 Math_pow10(void *fp)
716 {
717 F_Math_pow10 *f;
718
719 f = fp;
720
721 *f->ret = ipow10(f->p);
722 }
723
724 void
Math_remainder(void * fp)725 Math_remainder(void *fp)
726 {
727 F_Math_remainder *f;
728
729 f = fp;
730
731 *f->ret = __ieee754_remainder(f->x, f->p);
732 }
733
734 void
Math_scalbn(void * fp)735 Math_scalbn(void *fp)
736 {
737 F_Math_scalbn *f;
738
739 f = fp;
740
741 *f->ret = scalbn(f->x, f->n);
742 }
743
744 void
Math_yn(void * fp)745 Math_yn(void *fp)
746 {
747 F_Math_yn *f;
748
749 f = fp;
750
751 *f->ret = __ieee754_yn(f->n, f->x);
752 }
753
754
755 /**** sorting real vectors through permutation vector ****/
756 /* qsort from coma:/usr/jlb/qsort/qsort.dir/qsort.c on 28 Sep '92
757 char* has been changed to uchar*, static internal functions.
758 specialized to swapping ints (which are 32-bit anyway in limbo).
759 converted uchar* to int* (and substituted 1 for es).
760 */
761
762 static int
cmp(int * u,int * v,double * x)763 cmp(int *u, int *v, double *x)
764 {
765 return ((x[*u]==x[*v])? 0 : ((x[*u]<x[*v])? -1 : 1));
766 }
767
768 #define swap(u, v) {int t = *(u); *(u) = *(v); *(v) = t;}
769
770 #define vecswap(u, v, n) if(n>0){ \
771 int i = n; \
772 register int *pi = u; \
773 register int *pj = v; \
774 do { \
775 register int t = *pi; \
776 *pi++ = *pj; \
777 *pj++ = t; \
778 } while (--i > 0); \
779 }
780
781 #define minimum(x, y) ((x)<=(y) ? (x) : (y))
782
783 static int *
med3(int * a,int * b,int * c,double * x)784 med3(int *a, int *b, int *c, double *x)
785 { return cmp(a, b, x) < 0 ?
786 (cmp(b, c, x) < 0 ? b : (cmp(a, c, x) < 0 ? c : a ) )
787 : (cmp(b, c, x) > 0 ? b : (cmp(a, c, x) < 0 ? a : c ) );
788 }
789
790 void
rqsort(int * a,int n,double * x)791 rqsort(int *a, int n, double *x)
792 {
793 int *pa, *pb, *pc, *pd, *pl, *pm, *pn;
794 int d, r;
795
796 if (n < 7) { /* Insertion sort on small arrays */
797 for (pm = a + 1; pm < a + n; pm++)
798 for (pl = pm; pl > a && cmp(pl-1, pl, x) > 0; pl--)
799 swap(pl, pl-1);
800 return;
801 }
802 pm = a + (n/2);
803 if (n > 7) {
804 pl = a;
805 pn = a + (n-1);
806 if (n > 40) { /* On big arrays, pseudomedian of 9 */
807 d = (n/8);
808 pl = med3(pl, pl+d, pl+2*d, x);
809 pm = med3(pm-d, pm, pm+d, x);
810 pn = med3(pn-2*d, pn-d, pn, x);
811 }
812 pm = med3(pl, pm, pn, x); /* On mid arrays, med of 3 */
813 }
814 swap(a, pm); /* On tiny arrays, partition around middle */
815 pa = pb = a + 1;
816 pc = pd = a + (n-1);
817 for (;;) {
818 while (pb <= pc && (r = cmp(pb, a, x)) <= 0) {
819 if (r == 0) { swap(pa, pb); pa++; }
820 pb++;
821 }
822 while (pb <= pc && (r = cmp(pc, a, x)) >= 0) {
823 if (r == 0) { swap(pc, pd); pd--; }
824 pc--;
825 }
826 if (pb > pc) break;
827 swap(pb, pc);
828 pb++;
829 pc--;
830 }
831 pn = a + n;
832 r = minimum(pa-a, pb-pa); vecswap(a, pb-r, r);
833 r = minimum(pd-pc, pn-pd-1); vecswap(pb, pn-r, r);
834 if ((r = pb-pa) > 1) rqsort(a, r, x);
835 if ((r = pd-pc) > 1) rqsort(pn-r, r, x);
836 }
837
838 void
Math_sort(void * fp)839 Math_sort(void*fp)
840 {
841 F_Math_sort *f;
842 int i, pilen, xlen, *p;
843
844 f = fp;
845
846 /* check that permutation contents are in [0,n-1] !!! */
847 p = (int*) (f->pi->data);
848 pilen = f->pi->len;
849 xlen = f->x->len - 1;
850
851 for(i = 0; i < pilen; i++) {
852 if((*p < 0) || (xlen < *p))
853 error(exMathia);
854 p++;
855 }
856
857 rqsort( (int*)(f->pi->data), f->pi->len, (double*)(f->x->data));
858 }
859
860
861 /************ BLAS ***************/
862
863 void
Math_dot(void * fp)864 Math_dot(void *fp)
865 {
866 F_Math_dot *f;
867
868 f = fp;
869 if(f->x->len!=f->y->len)
870 error(exMathia); /* incompatible lengths */
871 *f->ret = dot(f->x->len, (double*)(f->x->data), (double*)(f->y->data));
872 }
873
874 void
Math_iamax(void * fp)875 Math_iamax(void *fp)
876 {
877 F_Math_iamax *f;
878
879 f = fp;
880
881 *f->ret = iamax(f->x->len, (double*)(f->x->data));
882 }
883
884 void
Math_norm2(void * fp)885 Math_norm2(void *fp)
886 {
887 F_Math_norm2 *f;
888
889 f = fp;
890
891 *f->ret = norm2(f->x->len, (double*)(f->x->data));
892 }
893
894 void
Math_norm1(void * fp)895 Math_norm1(void *fp)
896 {
897 F_Math_norm1 *f;
898
899 f = fp;
900
901 *f->ret = norm1(f->x->len, (double*)(f->x->data));
902 }
903
904 void
Math_gemm(void * fp)905 Math_gemm(void *fp)
906 {
907 F_Math_gemm *f = fp;
908 int nrowa, ncola, nrowb, ncolb, mn, ld, m, n;
909 double *adata = 0, *bdata = 0, *cdata;
910 int nota = f->transa=='N';
911 int notb = f->transb=='N';
912 if(nota){
913 nrowa = f->m;
914 ncola = f->k;
915 }else{
916 nrowa = f->k;
917 ncola = f->m;
918 }
919 if(notb){
920 nrowb = f->k;
921 ncolb = f->n;
922 }else{
923 nrowb = f->n;
924 ncolb = f->k;
925 }
926 if( (!nota && f->transa!='C' && f->transa!='T') ||
927 (!notb && f->transb!='C' && f->transb!='T') ||
928 (f->m < 0 || f->n < 0 || f->k < 0) ){
929 error(exMathia);
930 }
931 if(f->a != H){
932 mn = f->a->len;
933 adata = (double*)(f->a->data);
934 ld = f->lda;
935 if(ld<nrowa || ld*(ncola-1)>mn)
936 error(exBounds);
937 }
938 if(f->b != H){
939 mn = f->b->len;
940 ld = f->ldb;
941 bdata = (double*)(f->b->data);
942 if(ld<nrowb || ld*(ncolb-1)>mn)
943 error(exBounds);
944 }
945 m = f->m;
946 n = f->n;
947 mn = f->c->len;
948 cdata = (double*)(f->c->data);
949 ld = f->ldc;
950 if(ld<m || ld*(n-1)>mn)
951 error(exBounds);
952
953 gemm(f->transa, f->transb, f->m, f->n, f->k, f->alpha,
954 adata, f->lda, bdata, f->ldb, f->beta, cdata, f->ldc);
955 }
956