1 /* $OpenBSD: ieee.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* ieee.c 20 * 21 * Extended precision IEEE binary floating point arithmetic routines 22 * 23 * Numbers are stored in C language as arrays of 16-bit unsigned 24 * short integers. The arguments of the routines are pointers to 25 * the arrays. 26 * 27 * 28 * External e type data structure, simulates Intel 8087 chip 29 * temporary real format but possibly with a larger significand: 30 * 31 * NE-1 significand words (least significant word first, 32 * most significant bit is normally set) 33 * exponent (value = EXONE for 1.0, 34 * top bit is the sign) 35 * 36 * 37 * Internal data structure of a number (a "word" is 16 bits): 38 * 39 * ei[0] sign word (0 for positive, 0xffff for negative) 40 * ei[1] biased exponent (value = EXONE for the number 1.0) 41 * ei[2] high guard word (always zero after normalization) 42 * ei[3] 43 * to ei[NI-2] significand (NI-4 significand words, 44 * most significant word first, 45 * most significant bit is set) 46 * ei[NI-1] low guard word (0x8000 bit is rounding place) 47 * 48 * 49 * 50 * Routines for external format numbers 51 * 52 * asctoe( string, e ) ASCII string to extended double e type 53 * asctoe64( string, &d ) ASCII string to long double 54 * asctoe53( string, &d ) ASCII string to double 55 * asctoe24( string, &f ) ASCII string to single 56 * asctoeg( string, e, prec ) ASCII string to specified precision 57 * e24toe( &f, e ) IEEE single precision to e type 58 * e53toe( &d, e ) IEEE double precision to e type 59 * e64toe( &d, e ) IEEE long double precision to e type 60 * eabs(e) absolute value 61 * eadd( a, b, c ) c = b + a 62 * eclear(e) e = 0 63 * ecmp (a, b) Returns 1 if a > b, 0 if a == b, 64 * -1 if a < b, -2 if either a or b is a NaN. 65 * ediv( a, b, c ) c = b / a 66 * efloor( a, b ) truncate to integer, toward -infinity 67 * efrexp( a, exp, s ) extract exponent and significand 68 * eifrac( e, &l, frac ) e to long integer and e type fraction 69 * euifrac( e, &l, frac ) e to unsigned long integer and e type fraction 70 * einfin( e ) set e to infinity, leaving its sign alone 71 * eldexp( a, n, b ) multiply by 2**n 72 * emov( a, b ) b = a 73 * emul( a, b, c ) c = b * a 74 * eneg(e) e = -e 75 * eround( a, b ) b = nearest integer value to a 76 * esub( a, b, c ) c = b - a 77 * e24toasc( &f, str, n ) single to ASCII string, n digits after decimal 78 * e53toasc( &d, str, n ) double to ASCII string, n digits after decimal 79 * e64toasc( &d, str, n ) long double to ASCII string 80 * etoasc( e, str, n ) e to ASCII string, n digits after decimal 81 * etoe24( e, &f ) convert e type to IEEE single precision 82 * etoe53( e, &d ) convert e type to IEEE double precision 83 * etoe64( e, &d ) convert e type to IEEE long double precision 84 * ltoe( &l, e ) long (32 bit) integer to e type 85 * ultoe( &l, e ) unsigned long (32 bit) integer to e type 86 * eisneg( e ) 1 if sign bit of e != 0, else 0 87 * eisinf( e ) 1 if e has maximum exponent (non-IEEE) 88 * or is infinite (IEEE) 89 * eisnan( e ) 1 if e is a NaN 90 * esqrt( a, b ) b = square root of a 91 * 92 * 93 * Routines for internal format numbers 94 * 95 * eaddm( ai, bi ) add significands, bi = bi + ai 96 * ecleaz(ei) ei = 0 97 * ecleazs(ei) set ei = 0 but leave its sign alone 98 * ecmpm( ai, bi ) compare significands, return 1, 0, or -1 99 * edivm( ai, bi ) divide significands, bi = bi / ai 100 * emdnorm(ai,l,s,exp) normalize and round off 101 * emovi( a, ai ) convert external a to internal ai 102 * emovo( ai, a ) convert internal ai to external a 103 * emovz( ai, bi ) bi = ai, low guard word of bi = 0 104 * emulm( ai, bi ) multiply significands, bi = bi * ai 105 * enormlz(ei) left-justify the significand 106 * eshdn1( ai ) shift significand and guards down 1 bit 107 * eshdn8( ai ) shift down 8 bits 108 * eshdn6( ai ) shift down 16 bits 109 * eshift( ai, n ) shift ai n bits up (or down if n < 0) 110 * eshup1( ai ) shift significand and guards up 1 bit 111 * eshup8( ai ) shift up 8 bits 112 * eshup6( ai ) shift up 16 bits 113 * esubm( ai, bi ) subtract significands, bi = bi - ai 114 * 115 * 116 * The result is always normalized and rounded to NI-4 word precision 117 * after each arithmetic operation. 118 * 119 * Exception flags are NOT fully supported. 120 * 121 * Define INFINITY in mconf.h for support of infinity; otherwise a 122 * saturation arithmetic is implemented. 123 * 124 * Define NANS for support of Not-a-Number items; otherwise the 125 * arithmetic will never produce a NaN output, and might be confused 126 * by a NaN input. 127 * If NaN's are supported, the output of ecmp(a,b) is -2 if 128 * either a or b is a NaN. This means asking if(ecmp(a,b) < 0) 129 * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than 130 * if in doubt. 131 * Signaling NaN's are NOT supported; they are treated the same 132 * as quiet NaN's. 133 * 134 * Denormals are always supported here where appropriate (e.g., not 135 * for conversion to DEC numbers). 136 */ 137 138 /* 139 * Revision history: 140 * 141 * 5 Jan 84 PDP-11 assembly language version 142 * 2 Mar 86 fixed bug in asctoq() 143 * 6 Dec 86 C language version 144 * 30 Aug 88 100 digit version, improved rounding 145 * 15 May 92 80-bit long double support 146 * 147 * Author: S. L. Moshier. 148 */ 149 150 #include <stdio.h> 151 #include "mconf.h" 152 #include "ehead.h" 153 154 /* Change UNK into something else. */ 155 #ifdef UNK 156 #undef UNK 157 #if BIGENDIAN 158 #define MIEEE 1 159 #else 160 #define IBMPC 1 161 #endif 162 #endif 163 164 /* NaN's require infinity support. */ 165 #ifdef NANS 166 #ifndef INFINITY 167 #define INFINITY 168 #endif 169 #endif 170 171 /* This handles 64-bit long ints. */ 172 #define LONGBITS (8 * sizeof(long)) 173 174 /* Control register for rounding precision. 175 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits. 176 */ 177 int rndprc = NBITS; 178 extern int rndprc; 179 180 void eaddm(), esubm(), emdnorm(), asctoeg(), enan(); 181 static void toe24(), toe53(), toe64(), toe113(); 182 void eremain(), einit(), eiremain(); 183 int ecmpm(), edivm(), emulm(), eisneg(), eisinf(); 184 void emovi(), emovo(), emovz(), ecleaz(), eadd1(); 185 void etodec(), todec(), dectoe(); 186 int eisnan(), eiisnan(); 187 188 189 190 void einit() 191 { 192 } 193 194 /* 195 ; Clear out entire external format number. 196 ; 197 ; unsigned short x[]; 198 ; eclear( x ); 199 */ 200 201 void eclear( x ) 202 register unsigned short *x; 203 { 204 register int i; 205 206 for( i=0; i<NE; i++ ) 207 *x++ = 0; 208 } 209 210 211 212 /* Move external format number from a to b. 213 * 214 * emov( a, b ); 215 */ 216 217 void emov( a, b ) 218 register unsigned short *a, *b; 219 { 220 register int i; 221 222 for( i=0; i<NE; i++ ) 223 *b++ = *a++; 224 } 225 226 227 /* 228 ; Absolute value of external format number 229 ; 230 ; short x[NE]; 231 ; eabs( x ); 232 */ 233 234 void eabs(x) 235 unsigned short x[]; /* x is the memory address of a short */ 236 { 237 238 x[NE-1] &= 0x7fff; /* sign is top bit of last word of external format */ 239 } 240 241 242 243 244 /* 245 ; Negate external format number 246 ; 247 ; unsigned short x[NE]; 248 ; eneg( x ); 249 */ 250 251 void eneg(x) 252 unsigned short x[]; 253 { 254 255 #ifdef NANS 256 if( eisnan(x) ) 257 return; 258 #endif 259 x[NE-1] ^= 0x8000; /* Toggle the sign bit */ 260 } 261 262 263 264 /* Return 1 if external format number is negative, 265 * else return zero. 266 */ 267 int eisneg(x) 268 unsigned short x[]; 269 { 270 271 #ifdef NANS 272 if( eisnan(x) ) 273 return( 0 ); 274 #endif 275 if( x[NE-1] & 0x8000 ) 276 return( 1 ); 277 else 278 return( 0 ); 279 } 280 281 282 /* Return 1 if external format number has maximum possible exponent, 283 * else return zero. 284 */ 285 int eisinf(x) 286 unsigned short x[]; 287 { 288 289 if( (x[NE-1] & 0x7fff) == 0x7fff ) 290 { 291 #ifdef NANS 292 if( eisnan(x) ) 293 return( 0 ); 294 #endif 295 return( 1 ); 296 } 297 else 298 return( 0 ); 299 } 300 301 /* Check if e-type number is not a number. 302 */ 303 int eisnan(x) 304 unsigned short x[]; 305 { 306 307 #ifdef NANS 308 int i; 309 /* NaN has maximum exponent */ 310 if( (x[NE-1] & 0x7fff) != 0x7fff ) 311 return (0); 312 /* ... and non-zero significand field. */ 313 for( i=0; i<NE-1; i++ ) 314 { 315 if( *x++ != 0 ) 316 return (1); 317 } 318 #endif 319 return (0); 320 } 321 322 /* 323 ; Fill entire number, including exponent and significand, with 324 ; largest possible number. These programs implement a saturation 325 ; value that is an ordinary, legal number. A special value 326 ; "infinity" may also be implemented; this would require tests 327 ; for that value and implementation of special rules for arithmetic 328 ; operations involving inifinity. 329 */ 330 331 void einfin(x) 332 register unsigned short *x; 333 { 334 register int i; 335 336 #ifdef INFINITY 337 for( i=0; i<NE-1; i++ ) 338 *x++ = 0; 339 *x |= 32767; 340 #else 341 for( i=0; i<NE-1; i++ ) 342 *x++ = 0xffff; 343 *x |= 32766; 344 if( rndprc < NBITS ) 345 { 346 if (rndprc == 113) 347 { 348 *(x - 9) = 0; 349 *(x - 8) = 0; 350 } 351 if( rndprc == 64 ) 352 { 353 *(x-5) = 0; 354 } 355 if( rndprc == 53 ) 356 { 357 *(x-4) = 0xf800; 358 } 359 else 360 { 361 *(x-4) = 0; 362 *(x-3) = 0; 363 *(x-2) = 0xff00; 364 } 365 } 366 #endif 367 } 368 369 370 371 /* Move in external format number, 372 * converting it to internal format. 373 */ 374 void emovi( a, b ) 375 unsigned short *a, *b; 376 { 377 register unsigned short *p, *q; 378 int i; 379 380 q = b; 381 p = a + (NE-1); /* point to last word of external number */ 382 /* get the sign bit */ 383 if( *p & 0x8000 ) 384 *q++ = 0xffff; 385 else 386 *q++ = 0; 387 /* get the exponent */ 388 *q = *p--; 389 *q++ &= 0x7fff; /* delete the sign bit */ 390 #ifdef INFINITY 391 if( (*(q-1) & 0x7fff) == 0x7fff ) 392 { 393 #ifdef NANS 394 if( eisnan(a) ) 395 { 396 *q++ = 0; 397 for( i=3; i<NI; i++ ) 398 *q++ = *p--; 399 return; 400 } 401 #endif 402 for( i=2; i<NI; i++ ) 403 *q++ = 0; 404 return; 405 } 406 #endif 407 /* clear high guard word */ 408 *q++ = 0; 409 /* move in the significand */ 410 for( i=0; i<NE-1; i++ ) 411 *q++ = *p--; 412 /* clear low guard word */ 413 *q = 0; 414 } 415 416 417 /* Move internal format number out, 418 * converting it to external format. 419 */ 420 void emovo( a, b ) 421 unsigned short *a, *b; 422 { 423 register unsigned short *p, *q; 424 unsigned short i; 425 426 p = a; 427 q = b + (NE-1); /* point to output exponent */ 428 /* combine sign and exponent */ 429 i = *p++; 430 if( i ) 431 *q-- = *p++ | 0x8000; 432 else 433 *q-- = *p++; 434 #ifdef INFINITY 435 if( *(p-1) == 0x7fff ) 436 { 437 #ifdef NANS 438 if( eiisnan(a) ) 439 { 440 enan( b, NBITS ); 441 return; 442 } 443 #endif 444 einfin(b); 445 return; 446 } 447 #endif 448 /* skip over guard word */ 449 ++p; 450 /* move the significand */ 451 for( i=0; i<NE-1; i++ ) 452 *q-- = *p++; 453 } 454 455 456 457 458 /* Clear out internal format number. 459 */ 460 461 void ecleaz( xi ) 462 register unsigned short *xi; 463 { 464 register int i; 465 466 for( i=0; i<NI; i++ ) 467 *xi++ = 0; 468 } 469 470 /* same, but don't touch the sign. */ 471 472 void ecleazs( xi ) 473 register unsigned short *xi; 474 { 475 register int i; 476 477 ++xi; 478 for(i=0; i<NI-1; i++) 479 *xi++ = 0; 480 } 481 482 483 484 485 /* Move internal format number from a to b. 486 */ 487 void emovz( a, b ) 488 register unsigned short *a, *b; 489 { 490 register int i; 491 492 for( i=0; i<NI-1; i++ ) 493 *b++ = *a++; 494 /* clear low guard word */ 495 *b = 0; 496 } 497 498 /* Return nonzero if internal format number is a NaN. 499 */ 500 501 int eiisnan (x) 502 unsigned short x[]; 503 { 504 int i; 505 506 if( (x[E] & 0x7fff) == 0x7fff ) 507 { 508 for( i=M+1; i<NI; i++ ) 509 { 510 if( x[i] != 0 ) 511 return(1); 512 } 513 } 514 return(0); 515 } 516 517 #ifdef INFINITY 518 /* Return nonzero if internal format number is infinite. */ 519 520 static int 521 eiisinf (x) 522 unsigned short x[]; 523 { 524 525 #ifdef NANS 526 if (eiisnan (x)) 527 return (0); 528 #endif 529 if ((x[E] & 0x7fff) == 0x7fff) 530 return (1); 531 return (0); 532 } 533 #endif 534 535 /* 536 ; Compare significands of numbers in internal format. 537 ; Guard words are included in the comparison. 538 ; 539 ; unsigned short a[NI], b[NI]; 540 ; cmpm( a, b ); 541 ; 542 ; for the significands: 543 ; returns +1 if a > b 544 ; 0 if a == b 545 ; -1 if a < b 546 */ 547 int ecmpm( a, b ) 548 register unsigned short *a, *b; 549 { 550 int i; 551 552 a += M; /* skip up to significand area */ 553 b += M; 554 for( i=M; i<NI; i++ ) 555 { 556 if( *a++ != *b++ ) 557 goto difrnt; 558 } 559 return(0); 560 561 difrnt: 562 if( *(--a) > *(--b) ) 563 return(1); 564 else 565 return(-1); 566 } 567 568 569 /* 570 ; Shift significand down by 1 bit 571 */ 572 573 void eshdn1(x) 574 register unsigned short *x; 575 { 576 register unsigned short bits; 577 int i; 578 579 x += M; /* point to significand area */ 580 581 bits = 0; 582 for( i=M; i<NI; i++ ) 583 { 584 if( *x & 1 ) 585 bits |= 1; 586 *x >>= 1; 587 if( bits & 2 ) 588 *x |= 0x8000; 589 bits <<= 1; 590 ++x; 591 } 592 } 593 594 595 596 /* 597 ; Shift significand up by 1 bit 598 */ 599 600 void eshup1(x) 601 register unsigned short *x; 602 { 603 register unsigned short bits; 604 int i; 605 606 x += NI-1; 607 bits = 0; 608 609 for( i=M; i<NI; i++ ) 610 { 611 if( *x & 0x8000 ) 612 bits |= 1; 613 *x <<= 1; 614 if( bits & 2 ) 615 *x |= 1; 616 bits <<= 1; 617 --x; 618 } 619 } 620 621 622 623 /* 624 ; Shift significand down by 8 bits 625 */ 626 627 void eshdn8(x) 628 register unsigned short *x; 629 { 630 register unsigned short newbyt, oldbyt; 631 int i; 632 633 x += M; 634 oldbyt = 0; 635 for( i=M; i<NI; i++ ) 636 { 637 newbyt = *x << 8; 638 *x >>= 8; 639 *x |= oldbyt; 640 oldbyt = newbyt; 641 ++x; 642 } 643 } 644 645 /* 646 ; Shift significand up by 8 bits 647 */ 648 649 void eshup8(x) 650 register unsigned short *x; 651 { 652 int i; 653 register unsigned short newbyt, oldbyt; 654 655 x += NI-1; 656 oldbyt = 0; 657 658 for( i=M; i<NI; i++ ) 659 { 660 newbyt = *x >> 8; 661 *x <<= 8; 662 *x |= oldbyt; 663 oldbyt = newbyt; 664 --x; 665 } 666 } 667 668 /* 669 ; Shift significand up by 16 bits 670 */ 671 672 void eshup6(x) 673 register unsigned short *x; 674 { 675 int i; 676 register unsigned short *p; 677 678 p = x + M; 679 x += M + 1; 680 681 for( i=M; i<NI-1; i++ ) 682 *p++ = *x++; 683 684 *p = 0; 685 } 686 687 /* 688 ; Shift significand down by 16 bits 689 */ 690 691 void eshdn6(x) 692 register unsigned short *x; 693 { 694 int i; 695 register unsigned short *p; 696 697 x += NI-1; 698 p = x + 1; 699 700 for( i=M; i<NI-1; i++ ) 701 *(--p) = *(--x); 702 703 *(--p) = 0; 704 } 705 706 /* 707 ; Add significands 708 ; x + y replaces y 709 */ 710 711 void eaddm( x, y ) 712 unsigned short *x, *y; 713 { 714 register unsigned long a; 715 int i; 716 unsigned int carry; 717 718 x += NI-1; 719 y += NI-1; 720 carry = 0; 721 for( i=M; i<NI; i++ ) 722 { 723 a = (unsigned long )(*x) + (unsigned long )(*y) + carry; 724 if( a & 0x10000 ) 725 carry = 1; 726 else 727 carry = 0; 728 *y = (unsigned short )a; 729 --x; 730 --y; 731 } 732 } 733 734 /* 735 ; Subtract significands 736 ; y - x replaces y 737 */ 738 739 void esubm( x, y ) 740 unsigned short *x, *y; 741 { 742 unsigned long a; 743 int i; 744 unsigned int carry; 745 746 x += NI-1; 747 y += NI-1; 748 carry = 0; 749 for( i=M; i<NI; i++ ) 750 { 751 a = (unsigned long )(*y) - (unsigned long )(*x) - carry; 752 if( a & 0x10000 ) 753 carry = 1; 754 else 755 carry = 0; 756 *y = (unsigned short )a; 757 --x; 758 --y; 759 } 760 } 761 762 763 /* Divide significands */ 764 765 static unsigned short equot[NI] = {0}; /* was static */ 766 767 #if 0 768 int edivm( den, num ) 769 unsigned short den[], num[]; 770 { 771 int i; 772 register unsigned short *p, *q; 773 unsigned short j; 774 775 p = &equot[0]; 776 *p++ = num[0]; 777 *p++ = num[1]; 778 779 for( i=M; i<NI; i++ ) 780 { 781 *p++ = 0; 782 } 783 784 /* Use faster compare and subtraction if denominator 785 * has only 15 bits of significance. 786 */ 787 p = &den[M+2]; 788 if( *p++ == 0 ) 789 { 790 for( i=M+3; i<NI; i++ ) 791 { 792 if( *p++ != 0 ) 793 goto fulldiv; 794 } 795 if( (den[M+1] & 1) != 0 ) 796 goto fulldiv; 797 eshdn1(num); 798 eshdn1(den); 799 800 p = &den[M+1]; 801 q = &num[M+1]; 802 803 for( i=0; i<NBITS+2; i++ ) 804 { 805 if( *p <= *q ) 806 { 807 *q -= *p; 808 j = 1; 809 } 810 else 811 { 812 j = 0; 813 } 814 eshup1(equot); 815 equot[NI-2] |= j; 816 eshup1(num); 817 } 818 goto divdon; 819 } 820 821 /* The number of quotient bits to calculate is 822 * NBITS + 1 scaling guard bit + 1 roundoff bit. 823 */ 824 fulldiv: 825 826 p = &equot[NI-2]; 827 for( i=0; i<NBITS+2; i++ ) 828 { 829 if( ecmpm(den,num) <= 0 ) 830 { 831 esubm(den, num); 832 j = 1; /* quotient bit = 1 */ 833 } 834 else 835 j = 0; 836 eshup1(equot); 837 *p |= j; 838 eshup1(num); 839 } 840 841 divdon: 842 843 eshdn1( equot ); 844 eshdn1( equot ); 845 846 /* test for nonzero remainder after roundoff bit */ 847 p = &num[M]; 848 j = 0; 849 for( i=M; i<NI; i++ ) 850 { 851 j |= *p++; 852 } 853 if( j ) 854 j = 1; 855 856 857 for( i=0; i<NI; i++ ) 858 num[i] = equot[i]; 859 return( (int )j ); 860 } 861 862 /* Multiply significands */ 863 int emulm( a, b ) 864 unsigned short a[], b[]; 865 { 866 unsigned short *p, *q; 867 int i, j, k; 868 869 equot[0] = b[0]; 870 equot[1] = b[1]; 871 for( i=M; i<NI; i++ ) 872 equot[i] = 0; 873 874 p = &a[NI-2]; 875 k = NBITS; 876 while( *p == 0 ) /* significand is not supposed to be all zero */ 877 { 878 eshdn6(a); 879 k -= 16; 880 } 881 if( (*p & 0xff) == 0 ) 882 { 883 eshdn8(a); 884 k -= 8; 885 } 886 887 q = &equot[NI-1]; 888 j = 0; 889 for( i=0; i<k; i++ ) 890 { 891 if( *p & 1 ) 892 eaddm(b, equot); 893 /* remember if there were any nonzero bits shifted out */ 894 if( *q & 1 ) 895 j |= 1; 896 eshdn1(a); 897 eshdn1(equot); 898 } 899 900 for( i=0; i<NI; i++ ) 901 b[i] = equot[i]; 902 903 /* return flag for lost nonzero bits */ 904 return(j); 905 } 906 907 #else 908 909 /* Multiply significand of e-type number b 910 by 16-bit quantity a, e-type result to c. */ 911 912 void m16m( a, b, c ) 913 unsigned short a; 914 unsigned short b[], c[]; 915 { 916 register unsigned short *pp; 917 register unsigned long carry; 918 unsigned short *ps; 919 unsigned short p[NI]; 920 unsigned long aa, m; 921 int i; 922 923 aa = a; 924 pp = &p[NI-2]; 925 *pp++ = 0; 926 *pp = 0; 927 ps = &b[NI-1]; 928 929 for( i=M+1; i<NI; i++ ) 930 { 931 if( *ps == 0 ) 932 { 933 --ps; 934 --pp; 935 *(pp-1) = 0; 936 } 937 else 938 { 939 m = (unsigned long) aa * *ps--; 940 carry = (m & 0xffff) + *pp; 941 *pp-- = (unsigned short )carry; 942 carry = (carry >> 16) + (m >> 16) + *pp; 943 *pp = (unsigned short )carry; 944 *(pp-1) = carry >> 16; 945 } 946 } 947 for( i=M; i<NI; i++ ) 948 c[i] = p[i]; 949 } 950 951 952 /* Divide significands. Neither the numerator nor the denominator 953 is permitted to have its high guard word nonzero. */ 954 955 956 int edivm( den, num ) 957 unsigned short den[], num[]; 958 { 959 int i; 960 register unsigned short *p; 961 unsigned long tnum; 962 unsigned short j, tdenm, tquot; 963 unsigned short tprod[NI+1]; 964 965 p = &equot[0]; 966 *p++ = num[0]; 967 *p++ = num[1]; 968 969 for( i=M; i<NI; i++ ) 970 { 971 *p++ = 0; 972 } 973 eshdn1( num ); 974 tdenm = den[M+1]; 975 for( i=M; i<NI; i++ ) 976 { 977 /* Find trial quotient digit (the radix is 65536). */ 978 tnum = (((unsigned long) num[M]) << 16) + num[M+1]; 979 980 /* Do not execute the divide instruction if it will overflow. */ 981 if( (tdenm * 0xffffL) < tnum ) 982 tquot = 0xffff; 983 else 984 tquot = tnum / tdenm; 985 986 /* Prove that the divide worked. */ 987 /* 988 tcheck = (unsigned long )tquot * tdenm; 989 if( tnum - tcheck > tdenm ) 990 tquot = 0xffff; 991 */ 992 /* Multiply denominator by trial quotient digit. */ 993 m16m( tquot, den, tprod ); 994 /* The quotient digit may have been overestimated. */ 995 if( ecmpm( tprod, num ) > 0 ) 996 { 997 tquot -= 1; 998 esubm( den, tprod ); 999 if( ecmpm( tprod, num ) > 0 ) 1000 { 1001 tquot -= 1; 1002 esubm( den, tprod ); 1003 } 1004 } 1005 /* 1006 if( ecmpm( tprod, num ) > 0 ) 1007 { 1008 eshow( "tprod", tprod ); 1009 eshow( "num ", num ); 1010 printf( "tnum = %08lx, tden = %04x, tquot = %04x\n", 1011 tnum, den[M+1], tquot ); 1012 } 1013 */ 1014 esubm( tprod, num ); 1015 /* 1016 if( ecmpm( num, den ) >= 0 ) 1017 { 1018 eshow( "num ", num ); 1019 eshow( "den ", den ); 1020 printf( "tnum = %08lx, tden = %04x, tquot = %04x\n", 1021 tnum, den[M+1], tquot ); 1022 } 1023 */ 1024 equot[i] = tquot; 1025 eshup6(num); 1026 } 1027 /* test for nonzero remainder after roundoff bit */ 1028 p = &num[M]; 1029 j = 0; 1030 for( i=M; i<NI; i++ ) 1031 { 1032 j |= *p++; 1033 } 1034 if( j ) 1035 j = 1; 1036 1037 for( i=0; i<NI; i++ ) 1038 num[i] = equot[i]; 1039 1040 return( (int )j ); 1041 } 1042 1043 1044 1045 /* Multiply significands */ 1046 int emulm( a, b ) 1047 unsigned short a[], b[]; 1048 { 1049 unsigned short *p, *q; 1050 unsigned short pprod[NI]; 1051 unsigned short j; 1052 int i; 1053 1054 equot[0] = b[0]; 1055 equot[1] = b[1]; 1056 for( i=M; i<NI; i++ ) 1057 equot[i] = 0; 1058 1059 j = 0; 1060 p = &a[NI-1]; 1061 q = &equot[NI-1]; 1062 for( i=M+1; i<NI; i++ ) 1063 { 1064 if( *p == 0 ) 1065 { 1066 --p; 1067 } 1068 else 1069 { 1070 m16m( *p--, b, pprod ); 1071 eaddm(pprod, equot); 1072 } 1073 j |= *q; 1074 eshdn6(equot); 1075 } 1076 1077 for( i=0; i<NI; i++ ) 1078 b[i] = equot[i]; 1079 1080 /* return flag for lost nonzero bits */ 1081 return( (int)j ); 1082 } 1083 1084 1085 /* 1086 eshow(str, x) 1087 char *str; 1088 unsigned short *x; 1089 { 1090 int i; 1091 1092 printf( "%s ", str ); 1093 for( i=0; i<NI; i++ ) 1094 printf( "%04x ", *x++ ); 1095 printf( "\n" ); 1096 } 1097 */ 1098 #endif 1099 1100 1101 1102 /* 1103 * Normalize and round off. 1104 * 1105 * The internal format number to be rounded is "s". 1106 * Input "lost" indicates whether the number is exact. 1107 * This is the so-called sticky bit. 1108 * 1109 * Input "subflg" indicates whether the number was obtained 1110 * by a subtraction operation. In that case if lost is nonzero 1111 * then the number is slightly smaller than indicated. 1112 * 1113 * Input "exp" is the biased exponent, which may be negative. 1114 * the exponent field of "s" is ignored but is replaced by 1115 * "exp" as adjusted by normalization and rounding. 1116 * 1117 * Input "rcntrl" is the rounding control. 1118 */ 1119 1120 static int rlast = -1; 1121 static int rw = 0; 1122 static unsigned short rmsk = 0; 1123 static unsigned short rmbit = 0; 1124 static unsigned short rebit = 0; 1125 static int re = 0; 1126 static unsigned short rbit[NI] = {0,0,0,0,0,0,0,0}; 1127 1128 void emdnorm( s, lost, subflg, exp, rcntrl ) 1129 unsigned short s[]; 1130 int lost; 1131 int subflg; 1132 long exp; 1133 int rcntrl; 1134 { 1135 int i, j; 1136 unsigned short r; 1137 1138 /* Normalize */ 1139 j = enormlz( s ); 1140 1141 /* a blank significand could mean either zero or infinity. */ 1142 #ifndef INFINITY 1143 if( j > NBITS ) 1144 { 1145 ecleazs( s ); 1146 return; 1147 } 1148 #endif 1149 exp -= j; 1150 #ifndef INFINITY 1151 if( exp >= 32767L ) 1152 goto overf; 1153 #else 1154 if( (j > NBITS) && (exp < 32767L) ) 1155 { 1156 ecleazs( s ); 1157 return; 1158 } 1159 #endif 1160 if( exp < 0L ) 1161 { 1162 if( exp > (long )(-NBITS-1) ) 1163 { 1164 j = (int )exp; 1165 i = eshift( s, j ); 1166 if( i ) 1167 lost = 1; 1168 } 1169 else 1170 { 1171 ecleazs( s ); 1172 return; 1173 } 1174 } 1175 /* Round off, unless told not to by rcntrl. */ 1176 if( rcntrl == 0 ) 1177 goto mdfin; 1178 /* Set up rounding parameters if the control register changed. */ 1179 if( rndprc != rlast ) 1180 { 1181 ecleaz( rbit ); 1182 switch( rndprc ) 1183 { 1184 default: 1185 case NBITS: 1186 rw = NI-1; /* low guard word */ 1187 rmsk = 0xffff; 1188 rmbit = 0x8000; 1189 rebit = 1; 1190 re = rw - 1; 1191 break; 1192 case 113: 1193 rw = 10; 1194 rmsk = 0x7fff; 1195 rmbit = 0x4000; 1196 rebit = 0x8000; 1197 re = rw; 1198 break; 1199 case 64: 1200 rw = 7; 1201 rmsk = 0xffff; 1202 rmbit = 0x8000; 1203 rebit = 1; 1204 re = rw-1; 1205 break; 1206 /* For DEC arithmetic */ 1207 case 56: 1208 rw = 6; 1209 rmsk = 0xff; 1210 rmbit = 0x80; 1211 rebit = 0x100; 1212 re = rw; 1213 break; 1214 case 53: 1215 rw = 6; 1216 rmsk = 0x7ff; 1217 rmbit = 0x0400; 1218 rebit = 0x800; 1219 re = rw; 1220 break; 1221 case 24: 1222 rw = 4; 1223 rmsk = 0xff; 1224 rmbit = 0x80; 1225 rebit = 0x100; 1226 re = rw; 1227 break; 1228 } 1229 rbit[re] = rebit; 1230 rlast = rndprc; 1231 } 1232 1233 /* Shift down 1 temporarily if the data structure has an implied 1234 * most significant bit and the number is denormal. 1235 * For rndprc = 64 or NBITS, there is no implied bit. 1236 * But Intel long double denormals lose one bit of significance even so. 1237 */ 1238 #ifdef IBMPC 1239 if( (exp <= 0) && (rndprc != NBITS) ) 1240 #else 1241 if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) ) 1242 #endif 1243 { 1244 lost |= s[NI-1] & 1; 1245 eshdn1(s); 1246 } 1247 /* Clear out all bits below the rounding bit, 1248 * remembering in r if any were nonzero. 1249 */ 1250 r = s[rw] & rmsk; 1251 if( rndprc < NBITS ) 1252 { 1253 i = rw + 1; 1254 while( i < NI ) 1255 { 1256 if( s[i] ) 1257 r |= 1; 1258 s[i] = 0; 1259 ++i; 1260 } 1261 } 1262 s[rw] &= ~rmsk; 1263 if( (r & rmbit) != 0 ) 1264 { 1265 if( r == rmbit ) 1266 { 1267 if( lost == 0 ) 1268 { /* round to even */ 1269 if( (s[re] & rebit) == 0 ) 1270 goto mddone; 1271 } 1272 else 1273 { 1274 if( subflg != 0 ) 1275 goto mddone; 1276 } 1277 } 1278 eaddm( rbit, s ); 1279 } 1280 mddone: 1281 #ifdef IBMPC 1282 if( (exp <= 0) && (rndprc != NBITS) ) 1283 #else 1284 if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) ) 1285 #endif 1286 { 1287 eshup1(s); 1288 } 1289 if( s[2] != 0 ) 1290 { /* overflow on roundoff */ 1291 eshdn1(s); 1292 exp += 1; 1293 } 1294 mdfin: 1295 s[NI-1] = 0; 1296 if( exp >= 32767L ) 1297 { 1298 #ifndef INFINITY 1299 overf: 1300 #endif 1301 #ifdef INFINITY 1302 s[1] = 32767; 1303 for( i=2; i<NI-1; i++ ) 1304 s[i] = 0; 1305 #else 1306 s[1] = 32766; 1307 s[2] = 0; 1308 for( i=M+1; i<NI-1; i++ ) 1309 s[i] = 0xffff; 1310 s[NI-1] = 0; 1311 if( (rndprc < 64) || (rndprc == 113) ) 1312 { 1313 s[rw] &= ~rmsk; 1314 if( rndprc == 24 ) 1315 { 1316 s[5] = 0; 1317 s[6] = 0; 1318 } 1319 } 1320 #endif 1321 return; 1322 } 1323 if( exp < 0 ) 1324 s[1] = 0; 1325 else 1326 s[1] = (unsigned short )exp; 1327 } 1328 1329 1330 1331 /* 1332 ; Subtract external format numbers. 1333 ; 1334 ; unsigned short a[NE], b[NE], c[NE]; 1335 ; esub( a, b, c ); c = b - a 1336 */ 1337 1338 static int subflg = 0; 1339 1340 void esub( a, b, c ) 1341 unsigned short *a, *b, *c; 1342 { 1343 1344 #ifdef NANS 1345 if( eisnan(a) ) 1346 { 1347 emov (a, c); 1348 return; 1349 } 1350 if( eisnan(b) ) 1351 { 1352 emov(b,c); 1353 return; 1354 } 1355 /* Infinity minus infinity is a NaN. 1356 * Test for subtracting infinities of the same sign. 1357 */ 1358 if( eisinf(a) && eisinf(b) && ((eisneg (a) ^ eisneg (b)) == 0)) 1359 { 1360 mtherr( "esub", DOMAIN ); 1361 enan( c, NBITS ); 1362 return; 1363 } 1364 #endif 1365 subflg = 1; 1366 eadd1( a, b, c ); 1367 } 1368 1369 1370 /* 1371 ; Add. 1372 ; 1373 ; unsigned short a[NE], b[NE], c[NE]; 1374 ; eadd( a, b, c ); c = b + a 1375 */ 1376 void eadd( a, b, c ) 1377 unsigned short *a, *b, *c; 1378 { 1379 1380 #ifdef NANS 1381 /* NaN plus anything is a NaN. */ 1382 if( eisnan(a) ) 1383 { 1384 emov(a,c); 1385 return; 1386 } 1387 if( eisnan(b) ) 1388 { 1389 emov(b,c); 1390 return; 1391 } 1392 /* Infinity minus infinity is a NaN. 1393 * Test for adding infinities of opposite signs. 1394 */ 1395 if( eisinf(a) && eisinf(b) 1396 && ((eisneg(a) ^ eisneg(b)) != 0) ) 1397 { 1398 mtherr( "eadd", DOMAIN ); 1399 enan( c, NBITS ); 1400 return; 1401 } 1402 #endif 1403 subflg = 0; 1404 eadd1( a, b, c ); 1405 } 1406 1407 void eadd1( a, b, c ) 1408 unsigned short *a, *b, *c; 1409 { 1410 unsigned short ai[NI], bi[NI], ci[NI]; 1411 int i, lost, j, k; 1412 long lt, lta, ltb; 1413 1414 #ifdef INFINITY 1415 if( eisinf(a) ) 1416 { 1417 emov(a,c); 1418 if( subflg ) 1419 eneg(c); 1420 return; 1421 } 1422 if( eisinf(b) ) 1423 { 1424 emov(b,c); 1425 return; 1426 } 1427 #endif 1428 emovi( a, ai ); 1429 emovi( b, bi ); 1430 if( subflg ) 1431 ai[0] = ~ai[0]; 1432 1433 /* compare exponents */ 1434 lta = ai[E]; 1435 ltb = bi[E]; 1436 lt = lta - ltb; 1437 if( lt > 0L ) 1438 { /* put the larger number in bi */ 1439 emovz( bi, ci ); 1440 emovz( ai, bi ); 1441 emovz( ci, ai ); 1442 ltb = bi[E]; 1443 lt = -lt; 1444 } 1445 lost = 0; 1446 if( lt != 0L ) 1447 { 1448 if( lt < (long )(-NBITS-1) ) 1449 goto done; /* answer same as larger addend */ 1450 k = (int )lt; 1451 lost = eshift( ai, k ); /* shift the smaller number down */ 1452 } 1453 else 1454 { 1455 /* exponents were the same, so must compare significands */ 1456 i = ecmpm( ai, bi ); 1457 if( i == 0 ) 1458 { /* the numbers are identical in magnitude */ 1459 /* if different signs, result is zero */ 1460 if( ai[0] != bi[0] ) 1461 { 1462 eclear(c); 1463 return; 1464 } 1465 /* if same sign, result is double */ 1466 /* double denomalized tiny number */ 1467 if( (bi[E] == 0) && ((bi[3] & 0x8000) == 0) ) 1468 { 1469 eshup1( bi ); 1470 goto done; 1471 } 1472 /* add 1 to exponent unless both are zero! */ 1473 for( j=1; j<NI-1; j++ ) 1474 { 1475 if( bi[j] != 0 ) 1476 { 1477 ltb += 1; 1478 if( ltb >= 0x7fff ) 1479 { 1480 eclear(c); 1481 einfin(c); 1482 if( ai[0] != 0 ) 1483 eneg(c); 1484 return; 1485 } 1486 break; 1487 } 1488 } 1489 bi[E] = (unsigned short )ltb; 1490 goto done; 1491 } 1492 if( i > 0 ) 1493 { /* put the larger number in bi */ 1494 emovz( bi, ci ); 1495 emovz( ai, bi ); 1496 emovz( ci, ai ); 1497 } 1498 } 1499 if( ai[0] == bi[0] ) 1500 { 1501 eaddm( ai, bi ); 1502 subflg = 0; 1503 } 1504 else 1505 { 1506 esubm( ai, bi ); 1507 subflg = 1; 1508 } 1509 emdnorm( bi, lost, subflg, ltb, 64 ); 1510 1511 done: 1512 emovo( bi, c ); 1513 } 1514 1515 1516 1517 /* 1518 ; Divide. 1519 ; 1520 ; unsigned short a[NE], b[NE], c[NE]; 1521 ; ediv( a, b, c ); c = b / a 1522 */ 1523 void ediv( a, b, c ) 1524 unsigned short *a, *b, *c; 1525 { 1526 unsigned short ai[NI], bi[NI]; 1527 int i, sign; 1528 long lt, lta, ltb; 1529 1530 /* IEEE says if result is not a NaN, the sign is "-" if and only if 1531 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */ 1532 sign = eisneg(a) ^ eisneg(b); 1533 1534 #ifdef NANS 1535 /* Return any NaN input. */ 1536 if( eisnan(a) ) 1537 { 1538 emov(a,c); 1539 return; 1540 } 1541 if( eisnan(b) ) 1542 { 1543 emov(b,c); 1544 return; 1545 } 1546 /* Zero over zero, or infinity over infinity, is a NaN. */ 1547 if( ((ecmp(a,ezero) == 0) && (ecmp(b,ezero) == 0)) 1548 || (eisinf (a) && eisinf (b)) ) 1549 { 1550 mtherr( "ediv", DOMAIN ); 1551 enan( c, NBITS ); 1552 return; 1553 } 1554 #endif 1555 /* Infinity over anything else is infinity. */ 1556 #ifdef INFINITY 1557 if( eisinf(b) ) 1558 { 1559 einfin(c); 1560 goto divsign; 1561 } 1562 if( eisinf(a) ) 1563 { 1564 eclear(c); 1565 goto divsign; 1566 } 1567 #endif 1568 emovi( a, ai ); 1569 emovi( b, bi ); 1570 lta = ai[E]; 1571 ltb = bi[E]; 1572 if( bi[E] == 0 ) 1573 { /* See if numerator is zero. */ 1574 for( i=1; i<NI-1; i++ ) 1575 { 1576 if( bi[i] != 0 ) 1577 { 1578 ltb -= enormlz( bi ); 1579 goto dnzro1; 1580 } 1581 } 1582 eclear(c); 1583 goto divsign; 1584 } 1585 dnzro1: 1586 1587 if( ai[E] == 0 ) 1588 { /* possible divide by zero */ 1589 for( i=1; i<NI-1; i++ ) 1590 { 1591 if( ai[i] != 0 ) 1592 { 1593 lta -= enormlz( ai ); 1594 goto dnzro2; 1595 } 1596 } 1597 einfin(c); 1598 mtherr( "ediv", SING ); 1599 goto divsign; 1600 } 1601 dnzro2: 1602 1603 i = edivm( ai, bi ); 1604 /* calculate exponent */ 1605 lt = ltb - lta + EXONE; 1606 emdnorm( bi, i, 0, lt, 64 ); 1607 emovo( bi, c ); 1608 1609 divsign: 1610 1611 if( sign ) 1612 *(c+(NE-1)) |= 0x8000; 1613 else 1614 *(c+(NE-1)) &= ~0x8000; 1615 } 1616 1617 1618 1619 /* 1620 ; Multiply. 1621 ; 1622 ; unsigned short a[NE], b[NE], c[NE]; 1623 ; emul( a, b, c ); c = b * a 1624 */ 1625 void emul( a, b, c ) 1626 unsigned short *a, *b, *c; 1627 { 1628 unsigned short ai[NI], bi[NI]; 1629 int i, j, sign; 1630 long lt, lta, ltb; 1631 1632 /* IEEE says if result is not a NaN, the sign is "-" if and only if 1633 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */ 1634 sign = eisneg(a) ^ eisneg(b); 1635 1636 #ifdef NANS 1637 /* NaN times anything is the same NaN. */ 1638 if( eisnan(a) ) 1639 { 1640 emov(a,c); 1641 return; 1642 } 1643 if( eisnan(b) ) 1644 { 1645 emov(b,c); 1646 return; 1647 } 1648 /* Zero times infinity is a NaN. */ 1649 if( (eisinf(a) && (ecmp(b,ezero) == 0)) 1650 || (eisinf(b) && (ecmp(a,ezero) == 0)) ) 1651 { 1652 mtherr( "emul", DOMAIN ); 1653 enan( c, NBITS ); 1654 return; 1655 } 1656 #endif 1657 /* Infinity times anything else is infinity. */ 1658 #ifdef INFINITY 1659 if( eisinf(a) || eisinf(b) ) 1660 { 1661 einfin(c); 1662 goto mulsign; 1663 } 1664 #endif 1665 emovi( a, ai ); 1666 emovi( b, bi ); 1667 lta = ai[E]; 1668 ltb = bi[E]; 1669 if( ai[E] == 0 ) 1670 { 1671 for( i=1; i<NI-1; i++ ) 1672 { 1673 if( ai[i] != 0 ) 1674 { 1675 lta -= enormlz( ai ); 1676 goto mnzer1; 1677 } 1678 } 1679 eclear(c); 1680 goto mulsign; 1681 } 1682 mnzer1: 1683 1684 if( bi[E] == 0 ) 1685 { 1686 for( i=1; i<NI-1; i++ ) 1687 { 1688 if( bi[i] != 0 ) 1689 { 1690 ltb -= enormlz( bi ); 1691 goto mnzer2; 1692 } 1693 } 1694 eclear(c); 1695 goto mulsign; 1696 } 1697 mnzer2: 1698 1699 /* Multiply significands */ 1700 j = emulm( ai, bi ); 1701 /* calculate exponent */ 1702 lt = lta + ltb - (EXONE - 1); 1703 emdnorm( bi, j, 0, lt, 64 ); 1704 emovo( bi, c ); 1705 /* IEEE says sign is "-" if and only if operands have opposite signs. */ 1706 mulsign: 1707 if( sign ) 1708 *(c+(NE-1)) |= 0x8000; 1709 else 1710 *(c+(NE-1)) &= ~0x8000; 1711 } 1712 1713 1714 1715 1716 /* 1717 ; Convert IEEE double precision to e type 1718 ; double d; 1719 ; unsigned short x[N+2]; 1720 ; e53toe( &d, x ); 1721 */ 1722 void e53toe( pe, y ) 1723 unsigned short *pe, *y; 1724 { 1725 #ifdef DEC 1726 1727 dectoe( pe, y ); /* see etodec.c */ 1728 1729 #else 1730 1731 register unsigned short r; 1732 register unsigned short *p, *e; 1733 unsigned short yy[NI]; 1734 int denorm, k; 1735 1736 e = pe; 1737 denorm = 0; /* flag if denormalized number */ 1738 ecleaz(yy); 1739 #ifdef IBMPC 1740 e += 3; 1741 #endif 1742 r = *e; 1743 yy[0] = 0; 1744 if( r & 0x8000 ) 1745 yy[0] = 0xffff; 1746 yy[M] = (r & 0x0f) | 0x10; 1747 r &= ~0x800f; /* strip sign and 4 significand bits */ 1748 #ifdef INFINITY 1749 if( r == 0x7ff0 ) 1750 { 1751 #ifdef NANS 1752 #ifdef IBMPC 1753 if( ((pe[3] & 0xf) != 0) || (pe[2] != 0) 1754 || (pe[1] != 0) || (pe[0] != 0) ) 1755 { 1756 enan( y, NBITS ); 1757 return; 1758 } 1759 #else 1760 if( ((pe[0] & 0xf) != 0) || (pe[1] != 0) 1761 || (pe[2] != 0) || (pe[3] != 0) ) 1762 { 1763 enan( y, NBITS ); 1764 return; 1765 } 1766 #endif 1767 #endif /* NANS */ 1768 eclear( y ); 1769 einfin( y ); 1770 if( yy[0] ) 1771 eneg(y); 1772 return; 1773 } 1774 #endif 1775 r >>= 4; 1776 /* If zero exponent, then the significand is denormalized. 1777 * So, take back the understood high significand bit. */ 1778 if( r == 0 ) 1779 { 1780 denorm = 1; 1781 yy[M] &= ~0x10; 1782 } 1783 r += EXONE - 01777; 1784 yy[E] = r; 1785 p = &yy[M+1]; 1786 #ifdef IBMPC 1787 *p++ = *(--e); 1788 *p++ = *(--e); 1789 *p++ = *(--e); 1790 #endif 1791 #ifdef MIEEE 1792 ++e; 1793 *p++ = *e++; 1794 *p++ = *e++; 1795 *p++ = *e++; 1796 #endif 1797 (void )eshift( yy, -5 ); 1798 if( denorm ) 1799 { /* if zero exponent, then normalize the significand */ 1800 if( (k = enormlz(yy)) > NBITS ) 1801 ecleazs(yy); 1802 else 1803 yy[E] -= (unsigned short )(k-1); 1804 } 1805 emovo( yy, y ); 1806 #endif /* not DEC */ 1807 } 1808 1809 void e64toe( pe, y ) 1810 unsigned short *pe, *y; 1811 { 1812 unsigned short yy[NI]; 1813 unsigned short *p, *q, *e; 1814 int i; 1815 1816 e = pe; 1817 p = yy; 1818 for( i=0; i<NE-5; i++ ) 1819 *p++ = 0; 1820 #ifdef IBMPC 1821 for( i=0; i<5; i++ ) 1822 *p++ = *e++; 1823 #endif 1824 #ifdef DEC 1825 for( i=0; i<5; i++ ) 1826 *p++ = *e++; 1827 #endif 1828 #ifdef MIEEE 1829 p = &yy[0] + (NE-1); 1830 *p-- = *e++; 1831 ++e; 1832 for( i=0; i<4; i++ ) 1833 *p-- = *e++; 1834 #endif 1835 1836 #ifdef IBMPC 1837 /* For Intel long double, shift denormal significand up 1 1838 -- but only if the top significand bit is zero. */ 1839 if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0) 1840 { 1841 unsigned short temp[NI+1]; 1842 emovi(yy, temp); 1843 eshup1(temp); 1844 emovo(temp,y); 1845 return; 1846 } 1847 #endif 1848 #ifdef INFINITY 1849 /* Point to the exponent field. */ 1850 p = &yy[NE-1]; 1851 if ((*p & 0x7fff) == 0x7fff) 1852 { 1853 #ifdef NANS 1854 #ifdef IBMPC 1855 for( i=0; i<4; i++ ) 1856 { 1857 if((i != 3 && pe[i] != 0) 1858 /* Check for Intel long double infinity pattern. */ 1859 || (i == 3 && pe[i] != 0x8000)) 1860 { 1861 enan( y, NBITS ); 1862 return; 1863 } 1864 } 1865 #else 1866 /* In Motorola extended precision format, the most significant 1867 bit of an infinity mantissa could be either 1 or 0. It is 1868 the lower order bits that tell whether the value is a NaN. */ 1869 if ((pe[2] & 0x7fff) != 0) 1870 goto bigend_nan; 1871 1872 for( i=3; i<=5; i++ ) 1873 { 1874 if( pe[i] != 0 ) 1875 { 1876 bigend_nan: 1877 enan( y, NBITS ); 1878 return; 1879 } 1880 } 1881 #endif 1882 #endif /* NANS */ 1883 eclear( y ); 1884 einfin( y ); 1885 if( *p & 0x8000 ) 1886 eneg(y); 1887 return; 1888 } 1889 #endif 1890 p = yy; 1891 q = y; 1892 for( i=0; i<NE; i++ ) 1893 *q++ = *p++; 1894 } 1895 1896 void e113toe(pe,y) 1897 unsigned short *pe, *y; 1898 { 1899 register unsigned short r; 1900 unsigned short *e, *p; 1901 unsigned short yy[NI]; 1902 int denorm, i; 1903 1904 e = pe; 1905 denorm = 0; 1906 ecleaz(yy); 1907 #ifdef IBMPC 1908 e += 7; 1909 #endif 1910 r = *e; 1911 yy[0] = 0; 1912 if( r & 0x8000 ) 1913 yy[0] = 0xffff; 1914 r &= 0x7fff; 1915 #ifdef INFINITY 1916 if( r == 0x7fff ) 1917 { 1918 #ifdef NANS 1919 #ifdef IBMPC 1920 for( i=0; i<7; i++ ) 1921 { 1922 if( pe[i] != 0 ) 1923 { 1924 enan( y, NBITS ); 1925 return; 1926 } 1927 } 1928 #else 1929 for( i=1; i<8; i++ ) 1930 { 1931 if( pe[i] != 0 ) 1932 { 1933 enan( y, NBITS ); 1934 return; 1935 } 1936 } 1937 #endif 1938 #endif /* NANS */ 1939 eclear( y ); 1940 einfin( y ); 1941 if( *e & 0x8000 ) 1942 eneg(y); 1943 return; 1944 } 1945 #endif /* INFINITY */ 1946 yy[E] = r; 1947 p = &yy[M + 1]; 1948 #ifdef IBMPC 1949 for( i=0; i<7; i++ ) 1950 *p++ = *(--e); 1951 #endif 1952 #ifdef MIEEE 1953 ++e; 1954 for( i=0; i<7; i++ ) 1955 *p++ = *e++; 1956 #endif 1957 /* If denormal, remove the implied bit; else shift down 1. */ 1958 if( r == 0 ) 1959 { 1960 yy[M] = 0; 1961 } 1962 else 1963 { 1964 yy[M] = 1; 1965 eshift( yy, -1 ); 1966 } 1967 emovo(yy,y); 1968 } 1969 1970 1971 /* 1972 ; Convert IEEE single precision to e type 1973 ; float d; 1974 ; unsigned short x[N+2]; 1975 ; dtox( &d, x ); 1976 */ 1977 void e24toe( pe, y ) 1978 unsigned short *pe, *y; 1979 { 1980 register unsigned short r; 1981 register unsigned short *p, *e; 1982 unsigned short yy[NI]; 1983 int denorm, k; 1984 1985 e = pe; 1986 denorm = 0; /* flag if denormalized number */ 1987 ecleaz(yy); 1988 #ifdef IBMPC 1989 e += 1; 1990 #endif 1991 #ifdef DEC 1992 e += 1; 1993 #endif 1994 r = *e; 1995 yy[0] = 0; 1996 if( r & 0x8000 ) 1997 yy[0] = 0xffff; 1998 yy[M] = (r & 0x7f) | 0200; 1999 r &= ~0x807f; /* strip sign and 7 significand bits */ 2000 #ifdef INFINITY 2001 if( r == 0x7f80 ) 2002 { 2003 #ifdef NANS 2004 #ifdef MIEEE 2005 if( ((pe[0] & 0x7f) != 0) || (pe[1] != 0) ) 2006 { 2007 enan( y, NBITS ); 2008 return; 2009 } 2010 #else 2011 if( ((pe[1] & 0x7f) != 0) || (pe[0] != 0) ) 2012 { 2013 enan( y, NBITS ); 2014 return; 2015 } 2016 #endif 2017 #endif /* NANS */ 2018 eclear( y ); 2019 einfin( y ); 2020 if( yy[0] ) 2021 eneg(y); 2022 return; 2023 } 2024 #endif 2025 r >>= 7; 2026 /* If zero exponent, then the significand is denormalized. 2027 * So, take back the understood high significand bit. */ 2028 if( r == 0 ) 2029 { 2030 denorm = 1; 2031 yy[M] &= ~0200; 2032 } 2033 r += EXONE - 0177; 2034 yy[E] = r; 2035 p = &yy[M+1]; 2036 #ifdef IBMPC 2037 *p++ = *(--e); 2038 #endif 2039 #ifdef DEC 2040 *p++ = *(--e); 2041 #endif 2042 #ifdef MIEEE 2043 ++e; 2044 *p++ = *e++; 2045 #endif 2046 (void )eshift( yy, -8 ); 2047 if( denorm ) 2048 { /* if zero exponent, then normalize the significand */ 2049 if( (k = enormlz(yy)) > NBITS ) 2050 ecleazs(yy); 2051 else 2052 yy[E] -= (unsigned short )(k-1); 2053 } 2054 emovo( yy, y ); 2055 } 2056 2057 void etoe113(x,e) 2058 unsigned short *x, *e; 2059 { 2060 unsigned short xi[NI]; 2061 long exp; 2062 int rndsav; 2063 2064 #ifdef NANS 2065 if( eisnan(x) ) 2066 { 2067 enan( e, 113 ); 2068 return; 2069 } 2070 #endif 2071 emovi( x, xi ); 2072 exp = (long )xi[E]; 2073 #ifdef INFINITY 2074 if( eisinf(x) ) 2075 goto nonorm; 2076 #endif 2077 /* round off to nearest or even */ 2078 rndsav = rndprc; 2079 rndprc = 113; 2080 emdnorm( xi, 0, 0, exp, 64 ); 2081 rndprc = rndsav; 2082 nonorm: 2083 toe113 (xi, e); 2084 } 2085 2086 /* move out internal format to ieee long double */ 2087 static void toe113(a,b) 2088 unsigned short *a, *b; 2089 { 2090 register unsigned short *p, *q; 2091 unsigned short i; 2092 2093 #ifdef NANS 2094 if( eiisnan(a) ) 2095 { 2096 enan( b, 113 ); 2097 return; 2098 } 2099 #endif 2100 p = a; 2101 #ifdef MIEEE 2102 q = b; 2103 #else 2104 q = b + 7; /* point to output exponent */ 2105 #endif 2106 2107 /* If not denormal, delete the implied bit. */ 2108 if( a[E] != 0 ) 2109 { 2110 eshup1 (a); 2111 } 2112 /* combine sign and exponent */ 2113 i = *p++; 2114 #ifdef MIEEE 2115 if( i ) 2116 *q++ = *p++ | 0x8000; 2117 else 2118 *q++ = *p++; 2119 #else 2120 if( i ) 2121 *q-- = *p++ | 0x8000; 2122 else 2123 *q-- = *p++; 2124 #endif 2125 /* skip over guard word */ 2126 ++p; 2127 /* move the significand */ 2128 #ifdef MIEEE 2129 for (i = 0; i < 7; i++) 2130 *q++ = *p++; 2131 #else 2132 for (i = 0; i < 7; i++) 2133 *q-- = *p++; 2134 #endif 2135 } 2136 2137 2138 void etoe64( x, e ) 2139 unsigned short *x, *e; 2140 { 2141 unsigned short xi[NI]; 2142 long exp; 2143 int rndsav; 2144 2145 #ifdef NANS 2146 if( eisnan(x) ) 2147 { 2148 enan( e, 64 ); 2149 return; 2150 } 2151 #endif 2152 emovi( x, xi ); 2153 exp = (long )xi[E]; /* adjust exponent for offset */ 2154 #ifdef INFINITY 2155 if( eisinf(x) ) 2156 goto nonorm; 2157 #endif 2158 /* round off to nearest or even */ 2159 rndsav = rndprc; 2160 rndprc = 64; 2161 emdnorm( xi, 0, 0, exp, 64 ); 2162 rndprc = rndsav; 2163 nonorm: 2164 toe64( xi, e ); 2165 } 2166 2167 /* move out internal format to ieee long double */ 2168 static void toe64( a, b ) 2169 unsigned short *a, *b; 2170 { 2171 register unsigned short *p, *q; 2172 unsigned short i; 2173 2174 #ifdef NANS 2175 if( eiisnan(a) ) 2176 { 2177 enan( b, 64 ); 2178 return; 2179 } 2180 #endif 2181 #ifdef IBMPC 2182 /* Shift Intel denormal significand down 1. */ 2183 if( a[E] == 0 ) 2184 eshdn1(a); 2185 #endif 2186 p = a; 2187 #ifdef MIEEE 2188 q = b; 2189 #else 2190 q = b + 4; /* point to output exponent */ 2191 #if 1 2192 /* NOTE: if data type is 96 bits wide, clear the last word here. */ 2193 *(q+1)= 0; 2194 #endif 2195 #endif 2196 2197 /* combine sign and exponent */ 2198 i = *p++; 2199 #ifdef MIEEE 2200 if( i ) 2201 *q++ = *p++ | 0x8000; 2202 else 2203 *q++ = *p++; 2204 *q++ = 0; 2205 #else 2206 if( i ) 2207 *q-- = *p++ | 0x8000; 2208 else 2209 *q-- = *p++; 2210 #endif 2211 /* skip over guard word */ 2212 ++p; 2213 /* move the significand */ 2214 #ifdef MIEEE 2215 for( i=0; i<4; i++ ) 2216 *q++ = *p++; 2217 #else 2218 #ifdef INFINITY 2219 if (eiisinf (a)) 2220 { 2221 /* Intel long double infinity. */ 2222 *q-- = 0x8000; 2223 *q-- = 0; 2224 *q-- = 0; 2225 *q = 0; 2226 return; 2227 } 2228 #endif 2229 for( i=0; i<4; i++ ) 2230 *q-- = *p++; 2231 #endif 2232 } 2233 2234 2235 /* 2236 ; e type to IEEE double precision 2237 ; double d; 2238 ; unsigned short x[NE]; 2239 ; etoe53( x, &d ); 2240 */ 2241 2242 #ifdef DEC 2243 2244 void etoe53( x, e ) 2245 unsigned short *x, *e; 2246 { 2247 etodec( x, e ); /* see etodec.c */ 2248 } 2249 2250 static void toe53( x, y ) 2251 unsigned short *x, *y; 2252 { 2253 todec( x, y ); 2254 } 2255 2256 #else 2257 2258 void etoe53( x, e ) 2259 unsigned short *x, *e; 2260 { 2261 unsigned short xi[NI]; 2262 long exp; 2263 int rndsav; 2264 2265 #ifdef NANS 2266 if( eisnan(x) ) 2267 { 2268 enan( e, 53 ); 2269 return; 2270 } 2271 #endif 2272 emovi( x, xi ); 2273 exp = (long )xi[E] - (EXONE - 0x3ff); /* adjust exponent for offsets */ 2274 #ifdef INFINITY 2275 if( eisinf(x) ) 2276 goto nonorm; 2277 #endif 2278 /* round off to nearest or even */ 2279 rndsav = rndprc; 2280 rndprc = 53; 2281 emdnorm( xi, 0, 0, exp, 64 ); 2282 rndprc = rndsav; 2283 nonorm: 2284 toe53( xi, e ); 2285 } 2286 2287 2288 static void toe53( x, y ) 2289 unsigned short *x, *y; 2290 { 2291 unsigned short i; 2292 unsigned short *p; 2293 2294 2295 #ifdef NANS 2296 if( eiisnan(x) ) 2297 { 2298 enan( y, 53 ); 2299 return; 2300 } 2301 #endif 2302 p = &x[0]; 2303 #ifdef IBMPC 2304 y += 3; 2305 #endif 2306 *y = 0; /* output high order */ 2307 if( *p++ ) 2308 *y = 0x8000; /* output sign bit */ 2309 2310 i = *p++; 2311 if( i >= (unsigned int )2047 ) 2312 { /* Saturate at largest number less than infinity. */ 2313 #ifdef INFINITY 2314 *y |= 0x7ff0; 2315 #ifdef IBMPC 2316 *(--y) = 0; 2317 *(--y) = 0; 2318 *(--y) = 0; 2319 #endif 2320 #ifdef MIEEE 2321 ++y; 2322 *y++ = 0; 2323 *y++ = 0; 2324 *y++ = 0; 2325 #endif 2326 #else 2327 *y |= (unsigned short )0x7fef; 2328 #ifdef IBMPC 2329 *(--y) = 0xffff; 2330 *(--y) = 0xffff; 2331 *(--y) = 0xffff; 2332 #endif 2333 #ifdef MIEEE 2334 ++y; 2335 *y++ = 0xffff; 2336 *y++ = 0xffff; 2337 *y++ = 0xffff; 2338 #endif 2339 #endif 2340 return; 2341 } 2342 if( i == 0 ) 2343 { 2344 (void )eshift( x, 4 ); 2345 } 2346 else 2347 { 2348 i <<= 4; 2349 (void )eshift( x, 5 ); 2350 } 2351 i |= *p++ & (unsigned short )0x0f; /* *p = xi[M] */ 2352 *y |= (unsigned short )i; /* high order output already has sign bit set */ 2353 #ifdef IBMPC 2354 *(--y) = *p++; 2355 *(--y) = *p++; 2356 *(--y) = *p; 2357 #endif 2358 #ifdef MIEEE 2359 ++y; 2360 *y++ = *p++; 2361 *y++ = *p++; 2362 *y++ = *p++; 2363 #endif 2364 } 2365 2366 #endif /* not DEC */ 2367 2368 2369 2370 /* 2371 ; e type to IEEE single precision 2372 ; float d; 2373 ; unsigned short x[N+2]; 2374 ; xtod( x, &d ); 2375 */ 2376 void etoe24( x, e ) 2377 unsigned short *x, *e; 2378 { 2379 long exp; 2380 unsigned short xi[NI]; 2381 int rndsav; 2382 2383 #ifdef NANS 2384 if( eisnan(x) ) 2385 { 2386 enan( e, 24 ); 2387 return; 2388 } 2389 #endif 2390 emovi( x, xi ); 2391 exp = (long )xi[E] - (EXONE - 0177); /* adjust exponent for offsets */ 2392 #ifdef INFINITY 2393 if( eisinf(x) ) 2394 goto nonorm; 2395 #endif 2396 /* round off to nearest or even */ 2397 rndsav = rndprc; 2398 rndprc = 24; 2399 emdnorm( xi, 0, 0, exp, 64 ); 2400 rndprc = rndsav; 2401 nonorm: 2402 toe24( xi, e ); 2403 } 2404 2405 static void toe24( x, y ) 2406 unsigned short *x, *y; 2407 { 2408 unsigned short i; 2409 unsigned short *p; 2410 2411 #ifdef NANS 2412 if( eiisnan(x) ) 2413 { 2414 enan( y, 24 ); 2415 return; 2416 } 2417 #endif 2418 p = &x[0]; 2419 #ifdef IBMPC 2420 y += 1; 2421 #endif 2422 #ifdef DEC 2423 y += 1; 2424 #endif 2425 *y = 0; /* output high order */ 2426 if( *p++ ) 2427 *y = 0x8000; /* output sign bit */ 2428 2429 i = *p++; 2430 if( i >= 255 ) 2431 { /* Saturate at largest number less than infinity. */ 2432 #ifdef INFINITY 2433 *y |= (unsigned short )0x7f80; 2434 #ifdef IBMPC 2435 *(--y) = 0; 2436 #endif 2437 #ifdef DEC 2438 *(--y) = 0; 2439 #endif 2440 #ifdef MIEEE 2441 ++y; 2442 *y = 0; 2443 #endif 2444 #else 2445 *y |= (unsigned short )0x7f7f; 2446 #ifdef IBMPC 2447 *(--y) = 0xffff; 2448 #endif 2449 #ifdef DEC 2450 *(--y) = 0xffff; 2451 #endif 2452 #ifdef MIEEE 2453 ++y; 2454 *y = 0xffff; 2455 #endif 2456 #endif 2457 return; 2458 } 2459 if( i == 0 ) 2460 { 2461 (void )eshift( x, 7 ); 2462 } 2463 else 2464 { 2465 i <<= 7; 2466 (void )eshift( x, 8 ); 2467 } 2468 i |= *p++ & (unsigned short )0x7f; /* *p = xi[M] */ 2469 *y |= i; /* high order output already has sign bit set */ 2470 #ifdef IBMPC 2471 *(--y) = *p; 2472 #endif 2473 #ifdef DEC 2474 *(--y) = *p; 2475 #endif 2476 #ifdef MIEEE 2477 ++y; 2478 *y = *p; 2479 #endif 2480 } 2481 2482 2483 /* Compare two e type numbers. 2484 * 2485 * unsigned short a[NE], b[NE]; 2486 * ecmp( a, b ); 2487 * 2488 * returns +1 if a > b 2489 * 0 if a == b 2490 * -1 if a < b 2491 * -2 if either a or b is a NaN. 2492 */ 2493 int ecmp( a, b ) 2494 unsigned short *a, *b; 2495 { 2496 unsigned short ai[NI], bi[NI]; 2497 register unsigned short *p, *q; 2498 register int i; 2499 int msign; 2500 2501 #ifdef NANS 2502 if (eisnan (a) || eisnan (b)) 2503 return( -2 ); 2504 #endif 2505 emovi( a, ai ); 2506 p = ai; 2507 emovi( b, bi ); 2508 q = bi; 2509 2510 if( *p != *q ) 2511 { /* the signs are different */ 2512 /* -0 equals + 0 */ 2513 for( i=1; i<NI-1; i++ ) 2514 { 2515 if( ai[i] != 0 ) 2516 goto nzro; 2517 if( bi[i] != 0 ) 2518 goto nzro; 2519 } 2520 return(0); 2521 nzro: 2522 if( *p == 0 ) 2523 return( 1 ); 2524 else 2525 return( -1 ); 2526 } 2527 /* both are the same sign */ 2528 if( *p == 0 ) 2529 msign = 1; 2530 else 2531 msign = -1; 2532 i = NI-1; 2533 do 2534 { 2535 if( *p++ != *q++ ) 2536 { 2537 goto diff; 2538 } 2539 } 2540 while( --i > 0 ); 2541 2542 return(0); /* equality */ 2543 2544 2545 2546 diff: 2547 2548 if( *(--p) > *(--q) ) 2549 return( msign ); /* p is bigger */ 2550 else 2551 return( -msign ); /* p is littler */ 2552 } 2553 2554 2555 2556 2557 /* Find nearest integer to x = floor( x + 0.5 ) 2558 * 2559 * unsigned short x[NE], y[NE] 2560 * eround( x, y ); 2561 */ 2562 void eround( x, y ) 2563 unsigned short *x, *y; 2564 { 2565 2566 eadd( ehalf, x, y ); 2567 efloor( y, y ); 2568 } 2569 2570 2571 2572 2573 /* 2574 ; convert long (32-bit) integer to e type 2575 ; 2576 ; long l; 2577 ; unsigned short x[NE]; 2578 ; ltoe( &l, x ); 2579 ; note &l is the memory address of l 2580 */ 2581 void ltoe( lp, y ) 2582 long *lp; /* lp is the memory address of a long integer */ 2583 unsigned short *y; /* y is the address of a short */ 2584 { 2585 unsigned short yi[NI]; 2586 unsigned long ll; 2587 int k; 2588 2589 ecleaz( yi ); 2590 if( *lp < 0 ) 2591 { 2592 ll = (unsigned long )( -(*lp) ); /* make it positive */ 2593 yi[0] = 0xffff; /* put correct sign in the e type number */ 2594 } 2595 else 2596 { 2597 ll = (unsigned long )( *lp ); 2598 } 2599 /* move the long integer to yi significand area */ 2600 if( sizeof(long) == 8 ) 2601 { 2602 yi[M] = (unsigned short) (ll >> (LONGBITS - 16)); 2603 yi[M + 1] = (unsigned short) (ll >> (LONGBITS - 32)); 2604 yi[M + 2] = (unsigned short) (ll >> 16); 2605 yi[M + 3] = (unsigned short) ll; 2606 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */ 2607 } 2608 else 2609 { 2610 yi[M] = (unsigned short )(ll >> 16); 2611 yi[M+1] = (unsigned short )ll; 2612 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */ 2613 } 2614 if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */ 2615 ecleaz( yi ); /* it was zero */ 2616 else 2617 yi[E] -= (unsigned short )k; /* subtract shift count from exponent */ 2618 emovo( yi, y ); /* output the answer */ 2619 } 2620 2621 /* 2622 ; convert unsigned long (32-bit) integer to e type 2623 ; 2624 ; unsigned long l; 2625 ; unsigned short x[NE]; 2626 ; ltox( &l, x ); 2627 ; note &l is the memory address of l 2628 */ 2629 void ultoe( lp, y ) 2630 unsigned long *lp; /* lp is the memory address of a long integer */ 2631 unsigned short *y; /* y is the address of a short */ 2632 { 2633 unsigned short yi[NI]; 2634 unsigned long ll; 2635 int k; 2636 2637 ecleaz( yi ); 2638 ll = *lp; 2639 2640 /* move the long integer to ayi significand area */ 2641 if( sizeof(long) == 8 ) 2642 { 2643 yi[M] = (unsigned short) (ll >> (LONGBITS - 16)); 2644 yi[M + 1] = (unsigned short) (ll >> (LONGBITS - 32)); 2645 yi[M + 2] = (unsigned short) (ll >> 16); 2646 yi[M + 3] = (unsigned short) ll; 2647 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */ 2648 } 2649 else 2650 { 2651 yi[M] = (unsigned short )(ll >> 16); 2652 yi[M+1] = (unsigned short )ll; 2653 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */ 2654 } 2655 if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */ 2656 ecleaz( yi ); /* it was zero */ 2657 else 2658 yi[E] -= (unsigned short )k; /* subtract shift count from exponent */ 2659 emovo( yi, y ); /* output the answer */ 2660 } 2661 2662 2663 /* 2664 ; Find long integer and fractional parts 2665 2666 ; long i; 2667 ; unsigned short x[NE], frac[NE]; 2668 ; xifrac( x, &i, frac ); 2669 2670 The integer output has the sign of the input. The fraction is 2671 the positive fractional part of abs(x). 2672 */ 2673 void eifrac( x, i, frac ) 2674 unsigned short *x; 2675 long *i; 2676 unsigned short *frac; 2677 { 2678 unsigned short xi[NI]; 2679 int j, k; 2680 unsigned long ll; 2681 2682 emovi( x, xi ); 2683 k = (int )xi[E] - (EXONE - 1); 2684 if( k <= 0 ) 2685 { 2686 /* if exponent <= 0, integer = 0 and real output is fraction */ 2687 *i = 0L; 2688 emovo( xi, frac ); 2689 return; 2690 } 2691 if( k > (8 * sizeof(long) - 1) ) 2692 { 2693 /* 2694 ; long integer overflow: output large integer 2695 ; and correct fraction 2696 */ 2697 j = 8 * sizeof(long) - 1; 2698 if( xi[0] ) 2699 *i = (long) ((unsigned long) 1) << j; 2700 else 2701 *i = (long) (((unsigned long) (~(0L))) >> 1); 2702 (void )eshift( xi, k ); 2703 } 2704 if( k > 16 ) 2705 { 2706 /* 2707 Shift more than 16 bits: shift up k-16 mod 16 2708 then shift by 16's. 2709 */ 2710 j = k - ((k >> 4) << 4); 2711 eshift (xi, j); 2712 ll = xi[M]; 2713 k -= j; 2714 do 2715 { 2716 eshup6 (xi); 2717 ll = (ll << 16) | xi[M]; 2718 } 2719 while ((k -= 16) > 0); 2720 *i = ll; 2721 if (xi[0]) 2722 *i = -(*i); 2723 } 2724 else 2725 { 2726 /* shift not more than 16 bits */ 2727 eshift( xi, k ); 2728 *i = (long )xi[M] & 0xffff; 2729 if( xi[0] ) 2730 *i = -(*i); 2731 } 2732 xi[0] = 0; 2733 xi[E] = EXONE - 1; 2734 xi[M] = 0; 2735 if( (k = enormlz( xi )) > NBITS ) 2736 ecleaz( xi ); 2737 else 2738 xi[E] -= (unsigned short )k; 2739 2740 emovo( xi, frac ); 2741 } 2742 2743 2744 /* 2745 ; Find unsigned long integer and fractional parts 2746 2747 ; unsigned long i; 2748 ; unsigned short x[NE], frac[NE]; 2749 ; xifrac( x, &i, frac ); 2750 2751 A negative e type input yields integer output = 0 2752 but correct fraction. 2753 */ 2754 void euifrac( x, i, frac ) 2755 unsigned short *x; 2756 unsigned long *i; 2757 unsigned short *frac; 2758 { 2759 unsigned short xi[NI]; 2760 int j, k; 2761 unsigned long ll; 2762 2763 emovi( x, xi ); 2764 k = (int )xi[E] - (EXONE - 1); 2765 if( k <= 0 ) 2766 { 2767 /* if exponent <= 0, integer = 0 and argument is fraction */ 2768 *i = 0L; 2769 emovo( xi, frac ); 2770 return; 2771 } 2772 if( k > (8 * sizeof(long)) ) 2773 { 2774 /* 2775 ; long integer overflow: output large integer 2776 ; and correct fraction 2777 */ 2778 *i = ~(0L); 2779 (void )eshift( xi, k ); 2780 } 2781 else if( k > 16 ) 2782 { 2783 /* 2784 Shift more than 16 bits: shift up k-16 mod 16 2785 then shift up by 16's. 2786 */ 2787 j = k - ((k >> 4) << 4); 2788 eshift (xi, j); 2789 ll = xi[M]; 2790 k -= j; 2791 do 2792 { 2793 eshup6 (xi); 2794 ll = (ll << 16) | xi[M]; 2795 } 2796 while ((k -= 16) > 0); 2797 *i = ll; 2798 } 2799 else 2800 { 2801 /* shift not more than 16 bits */ 2802 eshift( xi, k ); 2803 *i = (long )xi[M] & 0xffff; 2804 } 2805 2806 if( xi[0] ) /* A negative value yields unsigned integer 0. */ 2807 *i = 0L; 2808 2809 xi[0] = 0; 2810 xi[E] = EXONE - 1; 2811 xi[M] = 0; 2812 if( (k = enormlz( xi )) > NBITS ) 2813 ecleaz( xi ); 2814 else 2815 xi[E] -= (unsigned short )k; 2816 2817 emovo( xi, frac ); 2818 } 2819 2820 2821 2822 /* 2823 ; Shift significand 2824 ; 2825 ; Shifts significand area up or down by the number of bits 2826 ; given by the variable sc. 2827 */ 2828 int eshift( x, sc ) 2829 unsigned short *x; 2830 int sc; 2831 { 2832 unsigned short lost; 2833 unsigned short *p; 2834 2835 if( sc == 0 ) 2836 return( 0 ); 2837 2838 lost = 0; 2839 p = x + NI-1; 2840 2841 if( sc < 0 ) 2842 { 2843 sc = -sc; 2844 while( sc >= 16 ) 2845 { 2846 lost |= *p; /* remember lost bits */ 2847 eshdn6(x); 2848 sc -= 16; 2849 } 2850 2851 while( sc >= 8 ) 2852 { 2853 lost |= *p & 0xff; 2854 eshdn8(x); 2855 sc -= 8; 2856 } 2857 2858 while( sc > 0 ) 2859 { 2860 lost |= *p & 1; 2861 eshdn1(x); 2862 sc -= 1; 2863 } 2864 } 2865 else 2866 { 2867 while( sc >= 16 ) 2868 { 2869 eshup6(x); 2870 sc -= 16; 2871 } 2872 2873 while( sc >= 8 ) 2874 { 2875 eshup8(x); 2876 sc -= 8; 2877 } 2878 2879 while( sc > 0 ) 2880 { 2881 eshup1(x); 2882 sc -= 1; 2883 } 2884 } 2885 if( lost ) 2886 lost = 1; 2887 return( (int )lost ); 2888 } 2889 2890 2891 2892 /* 2893 ; normalize 2894 ; 2895 ; Shift normalizes the significand area pointed to by argument 2896 ; shift count (up = positive) is returned. 2897 */ 2898 int enormlz(x) 2899 unsigned short x[]; 2900 { 2901 register unsigned short *p; 2902 int sc; 2903 2904 sc = 0; 2905 p = &x[M]; 2906 if( *p != 0 ) 2907 goto normdn; 2908 ++p; 2909 if( *p & 0x8000 ) 2910 return( 0 ); /* already normalized */ 2911 while( *p == 0 ) 2912 { 2913 eshup6(x); 2914 sc += 16; 2915 /* With guard word, there are NBITS+16 bits available. 2916 * return true if all are zero. 2917 */ 2918 if( sc > NBITS ) 2919 return( sc ); 2920 } 2921 /* see if high byte is zero */ 2922 while( (*p & 0xff00) == 0 ) 2923 { 2924 eshup8(x); 2925 sc += 8; 2926 } 2927 /* now shift 1 bit at a time */ 2928 while( (*p & 0x8000) == 0) 2929 { 2930 eshup1(x); 2931 sc += 1; 2932 if( sc > (NBITS+16) ) 2933 { 2934 mtherr( "enormlz", UNDERFLOW ); 2935 return( sc ); 2936 } 2937 } 2938 return( sc ); 2939 2940 /* Normalize by shifting down out of the high guard word 2941 of the significand */ 2942 normdn: 2943 2944 if( *p & 0xff00 ) 2945 { 2946 eshdn8(x); 2947 sc -= 8; 2948 } 2949 while( *p != 0 ) 2950 { 2951 eshdn1(x); 2952 sc -= 1; 2953 2954 if( sc < -NBITS ) 2955 { 2956 mtherr( "enormlz", OVERFLOW ); 2957 return( sc ); 2958 } 2959 } 2960 return( sc ); 2961 } 2962 2963 2964 2965 2966 /* Convert e type number to decimal format ASCII string. 2967 * The constants are for 64 bit precision. 2968 */ 2969 2970 #define NTEN 12 2971 #define MAXP 4096 2972 2973 #if NE == 10 2974 static unsigned short etens[NTEN + 1][NE] = 2975 { 2976 {0x6576, 0x4a92, 0x804a, 0x153f, 2977 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */ 2978 {0x6a32, 0xce52, 0x329a, 0x28ce, 2979 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */ 2980 {0x526c, 0x50ce, 0xf18b, 0x3d28, 2981 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,}, 2982 {0x9c66, 0x58f8, 0xbc50, 0x5c54, 2983 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,}, 2984 {0x851e, 0xeab7, 0x98fe, 0x901b, 2985 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,}, 2986 {0x0235, 0x0137, 0x36b1, 0x336c, 2987 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,}, 2988 {0x50f8, 0x25fb, 0xc76b, 0x6b71, 2989 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,}, 2990 {0x0000, 0x0000, 0x0000, 0x0000, 2991 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,}, 2992 {0x0000, 0x0000, 0x0000, 0x0000, 2993 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,}, 2994 {0x0000, 0x0000, 0x0000, 0x0000, 2995 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,}, 2996 {0x0000, 0x0000, 0x0000, 0x0000, 2997 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,}, 2998 {0x0000, 0x0000, 0x0000, 0x0000, 2999 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,}, 3000 {0x0000, 0x0000, 0x0000, 0x0000, 3001 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */ 3002 }; 3003 3004 static unsigned short emtens[NTEN + 1][NE] = 3005 { 3006 {0x2030, 0xcffc, 0xa1c3, 0x8123, 3007 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */ 3008 {0x8264, 0xd2cb, 0xf2ea, 0x12d4, 3009 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */ 3010 {0xf53f, 0xf698, 0x6bd3, 0x0158, 3011 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,}, 3012 {0xe731, 0x04d4, 0xe3f2, 0xd332, 3013 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,}, 3014 {0xa23e, 0x5308, 0xfefb, 0x1155, 3015 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,}, 3016 {0xe26d, 0xdbde, 0xd05d, 0xb3f6, 3017 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,}, 3018 {0x2a20, 0x6224, 0x47b3, 0x98d7, 3019 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,}, 3020 {0x0b5b, 0x4af2, 0xa581, 0x18ed, 3021 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,}, 3022 {0xbf71, 0xa9b3, 0x7989, 0xbe68, 3023 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,}, 3024 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b, 3025 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,}, 3026 {0xc155, 0xa4a8, 0x404e, 0x6113, 3027 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,}, 3028 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 3029 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,}, 3030 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 3031 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */ 3032 }; 3033 #else 3034 static unsigned short etens[NTEN+1][NE] = { 3035 {0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */ 3036 {0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */ 3037 {0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,}, 3038 {0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,}, 3039 {0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,}, 3040 {0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,}, 3041 {0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,}, 3042 {0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,}, 3043 {0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,}, 3044 {0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,}, 3045 {0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,}, 3046 {0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,}, 3047 {0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */ 3048 }; 3049 3050 static unsigned short emtens[NTEN+1][NE] = { 3051 {0x2de4,0x9fde,0xd2ce,0x04c8,0xa6dd,0x0ad8,}, /* 10**-4096 */ 3052 {0x4925,0x2de4,0x3436,0x534f,0xceae,0x256b,}, /* 10**-2048 */ 3053 {0x87a6,0xc0bd,0xda57,0x82a5,0xa2a6,0x32b5,}, 3054 {0x7133,0xd21c,0xdb23,0xee32,0x9049,0x395a,}, 3055 {0xfa91,0x1939,0x637a,0x4325,0xc031,0x3cac,}, 3056 {0xac7d,0xe4a0,0x64bc,0x467c,0xddd0,0x3e55,}, 3057 {0x3f24,0xe9a5,0xa539,0xea27,0xa87f,0x3f2a,}, 3058 {0x67de,0x94ba,0x4539,0x1ead,0xcfb1,0x3f94,}, 3059 {0x4c2f,0xe15b,0xc44d,0x94be,0xe695,0x3fc9,}, 3060 {0xfdc2,0xcefc,0x8461,0x7711,0xabcc,0x3fe4,}, 3061 {0xd3c3,0x652b,0xe219,0x1758,0xd1b7,0x3ff1,}, 3062 {0x3d71,0xd70a,0x70a3,0x0a3d,0xa3d7,0x3ff8,}, 3063 {0xcccd,0xcccc,0xcccc,0xcccc,0xcccc,0x3ffb,}, /* 10**-1 */ 3064 }; 3065 #endif 3066 3067 void e24toasc( x, string, ndigs ) 3068 unsigned short x[]; 3069 char *string; 3070 int ndigs; 3071 { 3072 unsigned short w[NI]; 3073 3074 e24toe( x, w ); 3075 etoasc( w, string, ndigs ); 3076 } 3077 3078 3079 void e53toasc( x, string, ndigs ) 3080 unsigned short x[]; 3081 char *string; 3082 int ndigs; 3083 { 3084 unsigned short w[NI]; 3085 3086 e53toe( x, w ); 3087 etoasc( w, string, ndigs ); 3088 } 3089 3090 3091 void e64toasc( x, string, ndigs ) 3092 unsigned short x[]; 3093 char *string; 3094 int ndigs; 3095 { 3096 unsigned short w[NI]; 3097 3098 e64toe( x, w ); 3099 etoasc( w, string, ndigs ); 3100 } 3101 3102 void e113toasc (x, string, ndigs) 3103 unsigned short x[]; 3104 char *string; 3105 int ndigs; 3106 { 3107 unsigned short w[NI]; 3108 3109 e113toe (x, w); 3110 etoasc (w, string, ndigs); 3111 } 3112 3113 3114 void etoasc( x, string, ndigs ) 3115 unsigned short x[]; 3116 char *string; 3117 int ndigs; 3118 { 3119 long digit; 3120 unsigned short y[NI], t[NI], u[NI], w[NI]; 3121 unsigned short *p, *r, *ten; 3122 unsigned short sign; 3123 int i, j, k, expon, rndsav; 3124 char *s, *ss; 3125 unsigned short m; 3126 3127 rndsav = rndprc; 3128 #ifdef NANS 3129 if( eisnan(x) ) 3130 { 3131 sprintf( string, " NaN " ); 3132 goto bxit; 3133 } 3134 #endif 3135 rndprc = NBITS; /* set to full precision */ 3136 emov( x, y ); /* retain external format */ 3137 if( y[NE-1] & 0x8000 ) 3138 { 3139 sign = 0xffff; 3140 y[NE-1] &= 0x7fff; 3141 } 3142 else 3143 { 3144 sign = 0; 3145 } 3146 expon = 0; 3147 ten = &etens[NTEN][0]; 3148 emov( eone, t ); 3149 /* Test for zero exponent */ 3150 if( y[NE-1] == 0 ) 3151 { 3152 for( k=0; k<NE-1; k++ ) 3153 { 3154 if( y[k] != 0 ) 3155 goto tnzro; /* denormalized number */ 3156 } 3157 goto isone; /* legal all zeros */ 3158 } 3159 tnzro: 3160 3161 /* Test for infinity. 3162 */ 3163 if( y[NE-1] == 0x7fff ) 3164 { 3165 if( sign ) 3166 sprintf( string, " -Infinity " ); 3167 else 3168 sprintf( string, " Infinity " ); 3169 goto bxit; 3170 } 3171 3172 /* Test for exponent nonzero but significand denormalized. 3173 * This is an error condition. 3174 */ 3175 if( (y[NE-1] != 0) && ((y[NE-2] & 0x8000) == 0) ) 3176 { 3177 mtherr( "etoasc", DOMAIN ); 3178 sprintf( string, "NaN" ); 3179 goto bxit; 3180 } 3181 3182 /* Compare to 1.0 */ 3183 i = ecmp( eone, y ); 3184 if( i == 0 ) 3185 goto isone; 3186 3187 if( i < 0 ) 3188 { /* Number is greater than 1 */ 3189 /* Convert significand to an integer and strip trailing decimal zeros. */ 3190 emov( y, u ); 3191 u[NE-1] = EXONE + NBITS - 1; 3192 3193 p = &etens[NTEN-4][0]; 3194 m = 16; 3195 do 3196 { 3197 ediv( p, u, t ); 3198 efloor( t, w ); 3199 for( j=0; j<NE-1; j++ ) 3200 { 3201 if( t[j] != w[j] ) 3202 goto noint; 3203 } 3204 emov( t, u ); 3205 expon += (int )m; 3206 noint: 3207 p += NE; 3208 m >>= 1; 3209 } 3210 while( m != 0 ); 3211 3212 /* Rescale from integer significand */ 3213 u[NE-1] += y[NE-1] - (unsigned int )(EXONE + NBITS - 1); 3214 emov( u, y ); 3215 /* Find power of 10 */ 3216 emov( eone, t ); 3217 m = MAXP; 3218 p = &etens[0][0]; 3219 while( ecmp( ten, u ) <= 0 ) 3220 { 3221 if( ecmp( p, u ) <= 0 ) 3222 { 3223 ediv( p, u, u ); 3224 emul( p, t, t ); 3225 expon += (int )m; 3226 } 3227 m >>= 1; 3228 if( m == 0 ) 3229 break; 3230 p += NE; 3231 } 3232 } 3233 else 3234 { /* Number is less than 1.0 */ 3235 /* Pad significand with trailing decimal zeros. */ 3236 if( y[NE-1] == 0 ) 3237 { 3238 while( (y[NE-2] & 0x8000) == 0 ) 3239 { 3240 emul( ten, y, y ); 3241 expon -= 1; 3242 } 3243 } 3244 else 3245 { 3246 emovi( y, w ); 3247 for( i=0; i<NDEC+1; i++ ) 3248 { 3249 if( (w[NI-1] & 0x7) != 0 ) 3250 break; 3251 /* multiply by 10 */ 3252 emovz( w, u ); 3253 eshdn1( u ); 3254 eshdn1( u ); 3255 eaddm( w, u ); 3256 u[1] += 3; 3257 while( u[2] != 0 ) 3258 { 3259 eshdn1(u); 3260 u[1] += 1; 3261 } 3262 if( u[NI-1] != 0 ) 3263 break; 3264 if( eone[NE-1] <= u[1] ) 3265 break; 3266 emovz( u, w ); 3267 expon -= 1; 3268 } 3269 emovo( w, y ); 3270 } 3271 k = -MAXP; 3272 p = &emtens[0][0]; 3273 r = &etens[0][0]; 3274 emov( y, w ); 3275 emov( eone, t ); 3276 while( ecmp( eone, w ) > 0 ) 3277 { 3278 if( ecmp( p, w ) >= 0 ) 3279 { 3280 emul( r, w, w ); 3281 emul( r, t, t ); 3282 expon += k; 3283 } 3284 k /= 2; 3285 if( k == 0 ) 3286 break; 3287 p += NE; 3288 r += NE; 3289 } 3290 ediv( t, eone, t ); 3291 } 3292 isone: 3293 /* Find the first (leading) digit. */ 3294 emovi( t, w ); 3295 emovz( w, t ); 3296 emovi( y, w ); 3297 emovz( w, y ); 3298 eiremain( t, y ); 3299 digit = equot[NI-1]; 3300 while( (digit == 0) && (ecmp(y,ezero) != 0) ) 3301 { 3302 eshup1( y ); 3303 emovz( y, u ); 3304 eshup1( u ); 3305 eshup1( u ); 3306 eaddm( u, y ); 3307 eiremain( t, y ); 3308 digit = equot[NI-1]; 3309 expon -= 1; 3310 } 3311 s = string; 3312 if( sign ) 3313 *s++ = '-'; 3314 else 3315 *s++ = ' '; 3316 /* Examine number of digits requested by caller. */ 3317 if( ndigs < 0 ) 3318 ndigs = 0; 3319 if( ndigs > NDEC ) 3320 ndigs = NDEC; 3321 if( digit == 10 ) 3322 { 3323 *s++ = '1'; 3324 *s++ = '.'; 3325 if( ndigs > 0 ) 3326 { 3327 *s++ = '0'; 3328 ndigs -= 1; 3329 } 3330 expon += 1; 3331 } 3332 else 3333 { 3334 *s++ = (char )digit + '0'; 3335 *s++ = '.'; 3336 } 3337 /* Generate digits after the decimal point. */ 3338 for( k=0; k<=ndigs; k++ ) 3339 { 3340 /* multiply current number by 10, without normalizing */ 3341 eshup1( y ); 3342 emovz( y, u ); 3343 eshup1( u ); 3344 eshup1( u ); 3345 eaddm( u, y ); 3346 eiremain( t, y ); 3347 *s++ = (char )equot[NI-1] + '0'; 3348 } 3349 digit = equot[NI-1]; 3350 --s; 3351 ss = s; 3352 /* round off the ASCII string */ 3353 if( digit > 4 ) 3354 { 3355 /* Test for critical rounding case in ASCII output. */ 3356 if( digit == 5 ) 3357 { 3358 emovo( y, t ); 3359 if( ecmp(t,ezero) != 0 ) 3360 goto roun; /* round to nearest */ 3361 if( (*(s-1) & 1) == 0 ) 3362 goto doexp; /* round to even */ 3363 } 3364 /* Round up and propagate carry-outs */ 3365 roun: 3366 --s; 3367 k = *s & 0x7f; 3368 /* Carry out to most significant digit? */ 3369 if( k == '.' ) 3370 { 3371 --s; 3372 k = *s; 3373 k += 1; 3374 *s = (char )k; 3375 /* Most significant digit carries to 10? */ 3376 if( k > '9' ) 3377 { 3378 expon += 1; 3379 *s = '1'; 3380 } 3381 goto doexp; 3382 } 3383 /* Round up and carry out from less significant digits */ 3384 k += 1; 3385 *s = (char )k; 3386 if( k > '9' ) 3387 { 3388 *s = '0'; 3389 goto roun; 3390 } 3391 } 3392 doexp: 3393 /* 3394 if( expon >= 0 ) 3395 sprintf( ss, "e+%d", expon ); 3396 else 3397 sprintf( ss, "e%d", expon ); 3398 */ 3399 sprintf( ss, "E%d", expon ); 3400 bxit: 3401 rndprc = rndsav; 3402 } 3403 3404 3405 3406 3407 /* 3408 ; ASCTOQ 3409 ; ASCTOQ.MAC LATEST REV: 11 JAN 84 3410 ; SLM, 3 JAN 78 3411 ; 3412 ; Convert ASCII string to quadruple precision floating point 3413 ; 3414 ; Numeric input is free field decimal number 3415 ; with max of 15 digits with or without 3416 ; decimal point entered as ASCII from teletype. 3417 ; Entering E after the number followed by a second 3418 ; number causes the second number to be interpreted 3419 ; as a power of 10 to be multiplied by the first number 3420 ; (i.e., "scientific" notation). 3421 ; 3422 ; Usage: 3423 ; asctoq( string, q ); 3424 */ 3425 3426 /* ASCII to single */ 3427 void asctoe24( s, y ) 3428 char *s; 3429 unsigned short *y; 3430 { 3431 asctoeg( s, y, 24 ); 3432 } 3433 3434 3435 /* ASCII to double */ 3436 void asctoe53( s, y ) 3437 char *s; 3438 unsigned short *y; 3439 { 3440 #ifdef DEC 3441 asctoeg( s, y, 56 ); 3442 #else 3443 asctoeg( s, y, 53 ); 3444 #endif 3445 } 3446 3447 3448 /* ASCII to long double */ 3449 void asctoe64( s, y ) 3450 char *s; 3451 unsigned short *y; 3452 { 3453 asctoeg( s, y, 64 ); 3454 } 3455 3456 /* ASCII to 128-bit long double */ 3457 void asctoe113 (s, y) 3458 char *s; 3459 unsigned short *y; 3460 { 3461 asctoeg( s, y, 113 ); 3462 } 3463 3464 /* ASCII to super double */ 3465 void asctoe( s, y ) 3466 char *s; 3467 unsigned short *y; 3468 { 3469 asctoeg( s, y, NBITS ); 3470 } 3471 3472 /* Space to make a copy of the input string: */ 3473 static char lstr[82] = {0}; 3474 3475 void asctoeg( ss, y, oprec ) 3476 char *ss; 3477 unsigned short *y; 3478 int oprec; 3479 { 3480 unsigned short yy[NI], xt[NI], tt[NI]; 3481 int esign, decflg, sgnflg, nexp, exp, prec, lost; 3482 int k, trail, c, rndsav; 3483 long lexp; 3484 unsigned short nsign, *p; 3485 char *sp, *s; 3486 3487 /* Copy the input string. */ 3488 s = ss; 3489 while( *s == ' ' ) /* skip leading spaces */ 3490 ++s; 3491 sp = lstr; 3492 for( k=0; k<79; k++ ) 3493 { 3494 if( (*sp++ = *s++) == '\0' ) 3495 break; 3496 } 3497 *sp = '\0'; 3498 s = lstr; 3499 3500 rndsav = rndprc; 3501 rndprc = NBITS; /* Set to full precision */ 3502 lost = 0; 3503 nsign = 0; 3504 decflg = 0; 3505 sgnflg = 0; 3506 nexp = 0; 3507 exp = 0; 3508 prec = 0; 3509 ecleaz( yy ); 3510 trail = 0; 3511 3512 nxtcom: 3513 k = *s - '0'; 3514 if( (k >= 0) && (k <= 9) ) 3515 { 3516 /* Ignore leading zeros */ 3517 if( (prec == 0) && (decflg == 0) && (k == 0) ) 3518 goto donchr; 3519 /* Identify and strip trailing zeros after the decimal point. */ 3520 if( (trail == 0) && (decflg != 0) ) 3521 { 3522 sp = s; 3523 while( (*sp >= '0') && (*sp <= '9') ) 3524 ++sp; 3525 /* Check for syntax error */ 3526 c = *sp & 0x7f; 3527 if( (c != 'e') && (c != 'E') && (c != '\0') 3528 && (c != '\n') && (c != '\r') && (c != ' ') 3529 && (c != ',') ) 3530 goto error; 3531 --sp; 3532 while( *sp == '0' ) 3533 *sp-- = 'z'; 3534 trail = 1; 3535 if( *s == 'z' ) 3536 goto donchr; 3537 } 3538 /* If enough digits were given to more than fill up the yy register, 3539 * continuing until overflow into the high guard word yy[2] 3540 * guarantees that there will be a roundoff bit at the top 3541 * of the low guard word after normalization. 3542 */ 3543 if( yy[2] == 0 ) 3544 { 3545 if( decflg ) 3546 nexp += 1; /* count digits after decimal point */ 3547 eshup1( yy ); /* multiply current number by 10 */ 3548 emovz( yy, xt ); 3549 eshup1( xt ); 3550 eshup1( xt ); 3551 eaddm( xt, yy ); 3552 ecleaz( xt ); 3553 xt[NI-2] = (unsigned short )k; 3554 eaddm( xt, yy ); 3555 } 3556 else 3557 { 3558 /* Mark any lost non-zero digit. */ 3559 lost |= k; 3560 /* Count lost digits before the decimal point. */ 3561 if (decflg == 0) 3562 nexp -= 1; 3563 } 3564 prec += 1; 3565 goto donchr; 3566 } 3567 3568 switch( *s ) 3569 { 3570 case 'z': 3571 break; 3572 case 'E': 3573 case 'e': 3574 goto expnt; 3575 case '.': /* decimal point */ 3576 if( decflg ) 3577 goto error; 3578 ++decflg; 3579 break; 3580 case '-': 3581 nsign = 0xffff; 3582 if( sgnflg ) 3583 goto error; 3584 ++sgnflg; 3585 break; 3586 case '+': 3587 if( sgnflg ) 3588 goto error; 3589 ++sgnflg; 3590 break; 3591 case ',': 3592 case ' ': 3593 case '\0': 3594 case '\n': 3595 case '\r': 3596 goto daldone; 3597 case 'i': 3598 case 'I': 3599 goto infinite; 3600 default: 3601 error: 3602 #ifdef NANS 3603 enan( yy, NI*16 ); 3604 #else 3605 mtherr( "asctoe", DOMAIN ); 3606 ecleaz(yy); 3607 #endif 3608 goto aexit; 3609 } 3610 donchr: 3611 ++s; 3612 goto nxtcom; 3613 3614 /* Exponent interpretation */ 3615 expnt: 3616 3617 esign = 1; 3618 exp = 0; 3619 ++s; 3620 /* check for + or - */ 3621 if( *s == '-' ) 3622 { 3623 esign = -1; 3624 ++s; 3625 } 3626 if( *s == '+' ) 3627 ++s; 3628 while( (*s >= '0') && (*s <= '9') ) 3629 { 3630 exp *= 10; 3631 exp += *s++ - '0'; 3632 if (exp > 4977) 3633 { 3634 if (esign < 0) 3635 goto zero; 3636 else 3637 goto infinite; 3638 } 3639 } 3640 if( esign < 0 ) 3641 exp = -exp; 3642 if( exp > 4932 ) 3643 { 3644 infinite: 3645 ecleaz(yy); 3646 yy[E] = 0x7fff; /* infinity */ 3647 goto aexit; 3648 } 3649 if( exp < -4977 ) 3650 { 3651 zero: 3652 ecleaz(yy); 3653 goto aexit; 3654 } 3655 3656 daldone: 3657 nexp = exp - nexp; 3658 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */ 3659 while( (nexp > 0) && (yy[2] == 0) ) 3660 { 3661 emovz( yy, xt ); 3662 eshup1( xt ); 3663 eshup1( xt ); 3664 eaddm( yy, xt ); 3665 eshup1( xt ); 3666 if( xt[2] != 0 ) 3667 break; 3668 nexp -= 1; 3669 emovz( xt, yy ); 3670 } 3671 if( (k = enormlz(yy)) > NBITS ) 3672 { 3673 ecleaz(yy); 3674 goto aexit; 3675 } 3676 lexp = (EXONE - 1 + NBITS) - k; 3677 emdnorm( yy, lost, 0, lexp, 64 ); 3678 /* convert to external format */ 3679 3680 3681 /* Multiply by 10**nexp. If precision is 64 bits, 3682 * the maximum relative error incurred in forming 10**n 3683 * for 0 <= n <= 324 is 8.2e-20, at 10**180. 3684 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947. 3685 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435. 3686 */ 3687 lexp = yy[E]; 3688 if( nexp == 0 ) 3689 { 3690 k = 0; 3691 goto expdon; 3692 } 3693 esign = 1; 3694 if( nexp < 0 ) 3695 { 3696 nexp = -nexp; 3697 esign = -1; 3698 if( nexp > 4096 ) 3699 { /* Punt. Can't handle this without 2 divides. */ 3700 emovi( etens[0], tt ); 3701 lexp -= tt[E]; 3702 k = edivm( tt, yy ); 3703 lexp += EXONE; 3704 nexp -= 4096; 3705 } 3706 } 3707 p = &etens[NTEN][0]; 3708 emov( eone, xt ); 3709 exp = 1; 3710 do 3711 { 3712 if( exp & nexp ) 3713 emul( p, xt, xt ); 3714 p -= NE; 3715 exp = exp + exp; 3716 } 3717 while( exp <= MAXP ); 3718 3719 emovi( xt, tt ); 3720 if( esign < 0 ) 3721 { 3722 lexp -= tt[E]; 3723 k = edivm( tt, yy ); 3724 lexp += EXONE; 3725 } 3726 else 3727 { 3728 lexp += tt[E]; 3729 k = emulm( tt, yy ); 3730 lexp -= EXONE - 1; 3731 } 3732 3733 expdon: 3734 3735 /* Round and convert directly to the destination type */ 3736 if( oprec == 53 ) 3737 lexp -= EXONE - 0x3ff; 3738 else if( oprec == 24 ) 3739 lexp -= EXONE - 0177; 3740 #ifdef DEC 3741 else if( oprec == 56 ) 3742 lexp -= EXONE - 0201; 3743 #endif 3744 rndprc = oprec; 3745 emdnorm( yy, k, 0, lexp, 64 ); 3746 3747 aexit: 3748 3749 rndprc = rndsav; 3750 yy[0] = nsign; 3751 switch( oprec ) 3752 { 3753 #ifdef DEC 3754 case 56: 3755 todec( yy, y ); /* see etodec.c */ 3756 break; 3757 #endif 3758 case 53: 3759 toe53( yy, y ); 3760 break; 3761 case 24: 3762 toe24( yy, y ); 3763 break; 3764 case 64: 3765 toe64( yy, y ); 3766 break; 3767 case 113: 3768 toe113( yy, y ); 3769 break; 3770 case NBITS: 3771 emovo( yy, y ); 3772 break; 3773 } 3774 } 3775 3776 3777 3778 /* y = largest integer not greater than x 3779 * (truncated toward minus infinity) 3780 * 3781 * unsigned short x[NE], y[NE] 3782 * 3783 * efloor( x, y ); 3784 */ 3785 static unsigned short bmask[] = { 3786 0xffff, 3787 0xfffe, 3788 0xfffc, 3789 0xfff8, 3790 0xfff0, 3791 0xffe0, 3792 0xffc0, 3793 0xff80, 3794 0xff00, 3795 0xfe00, 3796 0xfc00, 3797 0xf800, 3798 0xf000, 3799 0xe000, 3800 0xc000, 3801 0x8000, 3802 0x0000, 3803 }; 3804 3805 void efloor( x, y ) 3806 unsigned short x[], y[]; 3807 { 3808 register unsigned short *p; 3809 int e, expon, i; 3810 unsigned short f[NE]; 3811 3812 emov( x, f ); /* leave in external format */ 3813 expon = (int )f[NE-1]; 3814 e = (expon & 0x7fff) - (EXONE - 1); 3815 if( e <= 0 ) 3816 { 3817 eclear(y); 3818 goto isitneg; 3819 } 3820 /* number of bits to clear out */ 3821 e = NBITS - e; 3822 emov( f, y ); 3823 if( e <= 0 ) 3824 return; 3825 3826 p = &y[0]; 3827 while( e >= 16 ) 3828 { 3829 *p++ = 0; 3830 e -= 16; 3831 } 3832 /* clear the remaining bits */ 3833 *p &= bmask[e]; 3834 /* truncate negatives toward minus infinity */ 3835 isitneg: 3836 3837 if( (unsigned short )expon & (unsigned short )0x8000 ) 3838 { 3839 for( i=0; i<NE-1; i++ ) 3840 { 3841 if( f[i] != y[i] ) 3842 { 3843 esub( eone, y, y ); 3844 break; 3845 } 3846 } 3847 } 3848 } 3849 3850 3851 /* unsigned short x[], s[]; 3852 * long *exp; 3853 * 3854 * efrexp( x, exp, s ); 3855 * 3856 * Returns s and exp such that s * 2**exp = x and .5 <= s < 1. 3857 * For example, 1.1 = 0.55 * 2**1 3858 * Handles denormalized numbers properly using long integer exp. 3859 */ 3860 void efrexp( x, exp, s ) 3861 unsigned short x[]; 3862 long *exp; 3863 unsigned short s[]; 3864 { 3865 unsigned short xi[NI]; 3866 long li; 3867 3868 emovi( x, xi ); 3869 li = (long )((short )xi[1]); 3870 3871 if( li == 0 ) 3872 { 3873 li -= enormlz( xi ); 3874 } 3875 xi[1] = 0x3ffe; 3876 emovo( xi, s ); 3877 *exp = li - 0x3ffe; 3878 } 3879 3880 3881 3882 /* unsigned short x[], y[]; 3883 * long pwr2; 3884 * 3885 * eldexp( x, pwr2, y ); 3886 * 3887 * Returns y = x * 2**pwr2. 3888 */ 3889 void eldexp( x, pwr2, y ) 3890 unsigned short x[]; 3891 long pwr2; 3892 unsigned short y[]; 3893 { 3894 unsigned short xi[NI]; 3895 long li; 3896 int i; 3897 3898 emovi( x, xi ); 3899 li = xi[1]; 3900 li += pwr2; 3901 i = 0; 3902 emdnorm( xi, i, i, li, 64 ); 3903 emovo( xi, y ); 3904 } 3905 3906 3907 /* c = remainder after dividing b by a 3908 * Least significant integer quotient bits left in equot[]. 3909 */ 3910 void eremain( a, b, c ) 3911 unsigned short a[], b[], c[]; 3912 { 3913 unsigned short den[NI], num[NI]; 3914 3915 #ifdef NANS 3916 if( eisinf(b) || (ecmp(a,ezero) == 0) || eisnan(a) || eisnan(b)) 3917 { 3918 enan( c, NBITS ); 3919 return; 3920 } 3921 #endif 3922 if( ecmp(a,ezero) == 0 ) 3923 { 3924 mtherr( "eremain", SING ); 3925 eclear( c ); 3926 return; 3927 } 3928 emovi( a, den ); 3929 emovi( b, num ); 3930 eiremain( den, num ); 3931 /* Sign of remainder = sign of quotient */ 3932 if( a[0] == b[0] ) 3933 num[0] = 0; 3934 else 3935 num[0] = 0xffff; 3936 emovo( num, c ); 3937 } 3938 3939 3940 void eiremain( den, num ) 3941 unsigned short den[], num[]; 3942 { 3943 long ld, ln; 3944 unsigned short j; 3945 3946 ld = den[E]; 3947 ld -= enormlz( den ); 3948 ln = num[E]; 3949 ln -= enormlz( num ); 3950 ecleaz( equot ); 3951 while( ln >= ld ) 3952 { 3953 if( ecmpm(den,num) <= 0 ) 3954 { 3955 esubm(den, num); 3956 j = 1; 3957 } 3958 else 3959 { 3960 j = 0; 3961 } 3962 eshup1(equot); 3963 equot[NI-1] |= j; 3964 eshup1(num); 3965 ln -= 1; 3966 } 3967 emdnorm( num, 0, 0, ln, 0 ); 3968 } 3969 3970 /* NaN bit patterns 3971 */ 3972 #ifdef MIEEE 3973 unsigned short nan113[8] = { 3974 0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}; 3975 unsigned short nan64[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}; 3976 unsigned short nan53[4] = {0x7fff, 0xffff, 0xffff, 0xffff}; 3977 unsigned short nan24[2] = {0x7fff, 0xffff}; 3978 #endif 3979 3980 #ifdef IBMPC 3981 unsigned short nan113[8] = {0, 0, 0, 0, 0, 0, 0xc000, 0xffff}; 3982 unsigned short nan64[6] = {0, 0, 0, 0xc000, 0xffff, 0}; 3983 unsigned short nan53[4] = {0, 0, 0, 0xfff8}; 3984 unsigned short nan24[2] = {0, 0xffc0}; 3985 #endif 3986 3987 3988 void enan (nan, size) 3989 unsigned short *nan; 3990 int size; 3991 { 3992 int i, n; 3993 unsigned short *p; 3994 3995 switch( size ) 3996 { 3997 #ifndef DEC 3998 case 113: 3999 n = 8; 4000 p = nan113; 4001 break; 4002 4003 case 64: 4004 n = 6; 4005 p = nan64; 4006 break; 4007 4008 case 53: 4009 n = 4; 4010 p = nan53; 4011 break; 4012 4013 case 24: 4014 n = 2; 4015 p = nan24; 4016 break; 4017 4018 case NBITS: 4019 for( i=0; i<NE-2; i++ ) 4020 *nan++ = 0; 4021 *nan++ = 0xc000; 4022 *nan++ = 0x7fff; 4023 return; 4024 4025 case NI*16: 4026 *nan++ = 0; 4027 *nan++ = 0x7fff; 4028 *nan++ = 0; 4029 *nan++ = 0xc000; 4030 for( i=4; i<NI; i++ ) 4031 *nan++ = 0; 4032 return; 4033 #endif 4034 default: 4035 mtherr( "enan", DOMAIN ); 4036 return; 4037 } 4038 for (i=0; i < n; i++) 4039 *nan++ = *p++; 4040 } 4041 4042 4043 4044 /* Longhand square root. */ 4045 4046 static int esqinited = 0; 4047 static unsigned short sqrndbit[NI]; 4048 4049 void esqrt( x, y ) 4050 short *x, *y; 4051 { 4052 unsigned short temp[NI], num[NI], sq[NI], xx[NI]; 4053 int i, j, k, n, nlups; 4054 long m, exp; 4055 4056 if( esqinited == 0 ) 4057 { 4058 ecleaz( sqrndbit ); 4059 sqrndbit[NI-2] = 1; 4060 esqinited = 1; 4061 } 4062 /* Check for arg <= 0 */ 4063 i = ecmp( x, ezero ); 4064 if( i <= 0 ) 4065 { 4066 #ifdef NANS 4067 if (i == -2) 4068 { 4069 enan (y, NBITS); 4070 return; 4071 } 4072 #endif 4073 eclear(y); 4074 if( i < 0 ) 4075 mtherr( "esqrt", DOMAIN ); 4076 return; 4077 } 4078 4079 #ifdef INFINITY 4080 if( eisinf(x) ) 4081 { 4082 eclear(y); 4083 einfin(y); 4084 return; 4085 } 4086 #endif 4087 /* Bring in the arg and renormalize if it is denormal. */ 4088 emovi( x, xx ); 4089 m = (long )xx[1]; /* local long word exponent */ 4090 if( m == 0 ) 4091 m -= enormlz( xx ); 4092 4093 /* Divide exponent by 2 */ 4094 m -= 0x3ffe; 4095 exp = (unsigned short )( (m / 2) + 0x3ffe ); 4096 4097 /* Adjust if exponent odd */ 4098 if( (m & 1) != 0 ) 4099 { 4100 if( m > 0 ) 4101 exp += 1; 4102 eshdn1( xx ); 4103 } 4104 4105 ecleaz( sq ); 4106 ecleaz( num ); 4107 n = 8; /* get 8 bits of result per inner loop */ 4108 nlups = rndprc; 4109 j = 0; 4110 4111 while( nlups > 0 ) 4112 { 4113 /* bring in next word of arg */ 4114 if( j < NE ) 4115 num[NI-1] = xx[j+3]; 4116 /* Do additional bit on last outer loop, for roundoff. */ 4117 if( nlups <= 8 ) 4118 n = nlups + 1; 4119 for( i=0; i<n; i++ ) 4120 { 4121 /* Next 2 bits of arg */ 4122 eshup1( num ); 4123 eshup1( num ); 4124 /* Shift up answer */ 4125 eshup1( sq ); 4126 /* Make trial divisor */ 4127 for( k=0; k<NI; k++ ) 4128 temp[k] = sq[k]; 4129 eshup1( temp ); 4130 eaddm( sqrndbit, temp ); 4131 /* Subtract and insert answer bit if it goes in */ 4132 if( ecmpm( temp, num ) <= 0 ) 4133 { 4134 esubm( temp, num ); 4135 sq[NI-2] |= 1; 4136 } 4137 } 4138 nlups -= n; 4139 j += 1; 4140 } 4141 4142 /* Adjust for extra, roundoff loop done. */ 4143 exp += (NBITS - 1) - rndprc; 4144 4145 /* Sticky bit = 1 if the remainder is nonzero. */ 4146 k = 0; 4147 for( i=3; i<NI; i++ ) 4148 k |= (int )num[i]; 4149 4150 /* Renormalize and round off. */ 4151 emdnorm( sq, k, 0, exp, 64 ); 4152 emovo( sq, y ); 4153 } 4154