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 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 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 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 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 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 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 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 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 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 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 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 231 Math_getFPcontrol(void *fp) 232 { 233 F_Math_getFPcontrol *f; 234 235 f = fp; 236 237 *f->ret = getFPcontrol(); 238 } 239 240 void 241 Math_getFPstatus(void *fp) 242 { 243 F_Math_getFPstatus *f; 244 245 f = fp; 246 247 *f->ret = getFPstatus(); 248 } 249 250 void 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 * 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 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 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 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 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 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 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 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