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 */ 177int rndprc = NBITS; 178extern int rndprc; 179 180void eaddm(), esubm(), emdnorm(), asctoeg(), enan(); 181static void toe24(), toe53(), toe64(), toe113(); 182void eremain(), einit(), eiremain(); 183int ecmpm(), edivm(), emulm(), eisneg(), eisinf(); 184void emovi(), emovo(), emovz(), ecleaz(), eadd1(); 185void etodec(), todec(), dectoe(); 186int eisnan(), eiisnan(); 187 188 189 190void einit() 191{ 192} 193 194/* 195; Clear out entire external format number. 196; 197; unsigned short x[]; 198; eclear( x ); 199*/ 200 201void eclear( x ) 202register unsigned short *x; 203{ 204register int i; 205 206for( 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 217void emov( a, b ) 218register unsigned short *a, *b; 219{ 220register int i; 221 222for( 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 234void eabs(x) 235unsigned short x[]; /* x is the memory address of a short */ 236{ 237 238x[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 251void eneg(x) 252unsigned short x[]; 253{ 254 255#ifdef NANS 256if( eisnan(x) ) 257 return; 258#endif 259x[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 */ 267int eisneg(x) 268unsigned short x[]; 269{ 270 271#ifdef NANS 272if( eisnan(x) ) 273 return( 0 ); 274#endif 275if( x[NE-1] & 0x8000 ) 276 return( 1 ); 277else 278 return( 0 ); 279} 280 281 282/* Return 1 if external format number has maximum possible exponent, 283 * else return zero. 284 */ 285int eisinf(x) 286unsigned short x[]; 287{ 288 289if( (x[NE-1] & 0x7fff) == 0x7fff ) 290 { 291#ifdef NANS 292 if( eisnan(x) ) 293 return( 0 ); 294#endif 295 return( 1 ); 296 } 297else 298 return( 0 ); 299} 300 301/* Check if e-type number is not a number. 302 */ 303int eisnan(x) 304unsigned short x[]; 305{ 306 307#ifdef NANS 308int i; 309/* NaN has maximum exponent */ 310if( (x[NE-1] & 0x7fff) != 0x7fff ) 311 return (0); 312/* ... and non-zero significand field. */ 313for( i=0; i<NE-1; i++ ) 314 { 315 if( *x++ != 0 ) 316 return (1); 317 } 318#endif 319return (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 331void einfin(x) 332register unsigned short *x; 333{ 334register int i; 335 336#ifdef INFINITY 337for( i=0; i<NE-1; i++ ) 338 *x++ = 0; 339*x |= 32767; 340#else 341for( i=0; i<NE-1; i++ ) 342 *x++ = 0xffff; 343*x |= 32766; 344if( 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 */ 374void emovi( a, b ) 375unsigned short *a, *b; 376{ 377register unsigned short *p, *q; 378int i; 379 380q = b; 381p = a + (NE-1); /* point to last word of external number */ 382/* get the sign bit */ 383if( *p & 0x8000 ) 384 *q++ = 0xffff; 385else 386 *q++ = 0; 387/* get the exponent */ 388*q = *p--; 389*q++ &= 0x7fff; /* delete the sign bit */ 390#ifdef INFINITY 391if( (*(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 */ 410for( 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 */ 420void emovo( a, b ) 421unsigned short *a, *b; 422{ 423register unsigned short *p, *q; 424unsigned short i; 425 426p = a; 427q = b + (NE-1); /* point to output exponent */ 428/* combine sign and exponent */ 429i = *p++; 430if( i ) 431 *q-- = *p++ | 0x8000; 432else 433 *q-- = *p++; 434#ifdef INFINITY 435if( *(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 */ 451for( i=0; i<NE-1; i++ ) 452 *q-- = *p++; 453} 454 455 456 457 458/* Clear out internal format number. 459 */ 460 461void ecleaz( xi ) 462register unsigned short *xi; 463{ 464register int i; 465 466for( i=0; i<NI; i++ ) 467 *xi++ = 0; 468} 469 470/* same, but don't touch the sign. */ 471 472void ecleazs( xi ) 473register unsigned short *xi; 474{ 475register int i; 476 477++xi; 478for(i=0; i<NI-1; i++) 479 *xi++ = 0; 480} 481 482 483 484 485/* Move internal format number from a to b. 486 */ 487void emovz( a, b ) 488register unsigned short *a, *b; 489{ 490register int i; 491 492for( 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 501int eiisnan (x) 502unsigned short x[]; 503{ 504int i; 505 506if( (x[E] & 0x7fff) == 0x7fff ) 507 { 508 for( i=M+1; i<NI; i++ ) 509 { 510 if( x[i] != 0 ) 511 return(1); 512 } 513 } 514return(0); 515} 516 517#ifdef INFINITY 518/* Return nonzero if internal format number is infinite. */ 519 520static int 521eiisinf (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*/ 547int ecmpm( a, b ) 548register unsigned short *a, *b; 549{ 550int i; 551 552a += M; /* skip up to significand area */ 553b += M; 554for( i=M; i<NI; i++ ) 555 { 556 if( *a++ != *b++ ) 557 goto difrnt; 558 } 559return(0); 560 561difrnt: 562if( *(--a) > *(--b) ) 563 return(1); 564else 565 return(-1); 566} 567 568 569/* 570; Shift significand down by 1 bit 571*/ 572 573void eshdn1(x) 574register unsigned short *x; 575{ 576register unsigned short bits; 577int i; 578 579x += M; /* point to significand area */ 580 581bits = 0; 582for( 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 600void eshup1(x) 601register unsigned short *x; 602{ 603register unsigned short bits; 604int i; 605 606x += NI-1; 607bits = 0; 608 609for( 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 627void eshdn8(x) 628register unsigned short *x; 629{ 630register unsigned short newbyt, oldbyt; 631int i; 632 633x += M; 634oldbyt = 0; 635for( 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 649void eshup8(x) 650register unsigned short *x; 651{ 652int i; 653register unsigned short newbyt, oldbyt; 654 655x += NI-1; 656oldbyt = 0; 657 658for( 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 672void eshup6(x) 673register unsigned short *x; 674{ 675int i; 676register unsigned short *p; 677 678p = x + M; 679x += M + 1; 680 681for( 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 691void eshdn6(x) 692register unsigned short *x; 693{ 694int i; 695register unsigned short *p; 696 697x += NI-1; 698p = x + 1; 699 700for( 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 711void eaddm( x, y ) 712unsigned short *x, *y; 713{ 714register unsigned long a; 715int i; 716unsigned int carry; 717 718x += NI-1; 719y += NI-1; 720carry = 0; 721for( 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 739void esubm( x, y ) 740unsigned short *x, *y; 741{ 742unsigned long a; 743int i; 744unsigned int carry; 745 746x += NI-1; 747y += NI-1; 748carry = 0; 749for( 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 765static unsigned short equot[NI] = {0}; /* was static */ 766 767#if 0 768int edivm( den, num ) 769unsigned short den[], num[]; 770{ 771int i; 772register unsigned short *p, *q; 773unsigned short j; 774 775p = &equot[0]; 776*p++ = num[0]; 777*p++ = num[1]; 778 779for( 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 */ 787p = &den[M+2]; 788if( *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 */ 824fulldiv: 825 826p = &equot[NI-2]; 827for( 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 841divdon: 842 843eshdn1( equot ); 844eshdn1( equot ); 845 846/* test for nonzero remainder after roundoff bit */ 847p = &num[M]; 848j = 0; 849for( i=M; i<NI; i++ ) 850 { 851 j |= *p++; 852 } 853if( j ) 854 j = 1; 855 856 857for( i=0; i<NI; i++ ) 858 num[i] = equot[i]; 859return( (int )j ); 860} 861 862/* Multiply significands */ 863int emulm( a, b ) 864unsigned short a[], b[]; 865{ 866unsigned short *p, *q; 867int i, j, k; 868 869equot[0] = b[0]; 870equot[1] = b[1]; 871for( i=M; i<NI; i++ ) 872 equot[i] = 0; 873 874p = &a[NI-2]; 875k = NBITS; 876while( *p == 0 ) /* significand is not supposed to be all zero */ 877 { 878 eshdn6(a); 879 k -= 16; 880 } 881if( (*p & 0xff) == 0 ) 882 { 883 eshdn8(a); 884 k -= 8; 885 } 886 887q = &equot[NI-1]; 888j = 0; 889for( 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 900for( i=0; i<NI; i++ ) 901 b[i] = equot[i]; 902 903/* return flag for lost nonzero bits */ 904return(j); 905} 906 907#else 908 909/* Multiply significand of e-type number b 910by 16-bit quantity a, e-type result to c. */ 911 912void m16m( a, b, c ) 913unsigned short a; 914unsigned short b[], c[]; 915{ 916register unsigned short *pp; 917register unsigned long carry; 918unsigned short *ps; 919unsigned short p[NI]; 920unsigned long aa, m; 921int i; 922 923aa = a; 924pp = &p[NI-2]; 925*pp++ = 0; 926*pp = 0; 927ps = &b[NI-1]; 928 929for( 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 } 947for( i=M; i<NI; i++ ) 948 c[i] = p[i]; 949} 950 951 952/* Divide significands. Neither the numerator nor the denominator 953is permitted to have its high guard word nonzero. */ 954 955 956int edivm( den, num ) 957unsigned short den[], num[]; 958{ 959int i; 960register unsigned short *p; 961unsigned long tnum; 962unsigned short j, tdenm, tquot; 963unsigned short tprod[NI+1]; 964 965p = &equot[0]; 966*p++ = num[0]; 967*p++ = num[1]; 968 969for( i=M; i<NI; i++ ) 970 { 971 *p++ = 0; 972 } 973eshdn1( num ); 974tdenm = den[M+1]; 975for( 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 */ 1028p = &num[M]; 1029j = 0; 1030for( i=M; i<NI; i++ ) 1031 { 1032 j |= *p++; 1033 } 1034if( j ) 1035 j = 1; 1036 1037for( i=0; i<NI; i++ ) 1038 num[i] = equot[i]; 1039 1040return( (int )j ); 1041} 1042 1043 1044 1045/* Multiply significands */ 1046int emulm( a, b ) 1047unsigned short a[], b[]; 1048{ 1049unsigned short *p, *q; 1050unsigned short pprod[NI]; 1051unsigned short j; 1052int i; 1053 1054equot[0] = b[0]; 1055equot[1] = b[1]; 1056for( i=M; i<NI; i++ ) 1057 equot[i] = 0; 1058 1059j = 0; 1060p = &a[NI-1]; 1061q = &equot[NI-1]; 1062for( 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 1077for( i=0; i<NI; i++ ) 1078 b[i] = equot[i]; 1079 1080/* return flag for lost nonzero bits */ 1081return( (int)j ); 1082} 1083 1084 1085/* 1086eshow(str, x) 1087char *str; 1088unsigned short *x; 1089{ 1090int i; 1091 1092printf( "%s ", str ); 1093for( i=0; i<NI; i++ ) 1094 printf( "%04x ", *x++ ); 1095printf( "\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 1120static int rlast = -1; 1121static int rw = 0; 1122static unsigned short rmsk = 0; 1123static unsigned short rmbit = 0; 1124static unsigned short rebit = 0; 1125static int re = 0; 1126static unsigned short rbit[NI] = {0,0,0,0,0,0,0,0}; 1127 1128void emdnorm( s, lost, subflg, exp, rcntrl ) 1129unsigned short s[]; 1130int lost; 1131int subflg; 1132long exp; 1133int rcntrl; 1134{ 1135int i, j; 1136unsigned short r; 1137 1138/* Normalize */ 1139j = enormlz( s ); 1140 1141/* a blank significand could mean either zero or infinity. */ 1142#ifndef INFINITY 1143if( j > NBITS ) 1144 { 1145 ecleazs( s ); 1146 return; 1147 } 1148#endif 1149exp -= j; 1150#ifndef INFINITY 1151if( exp >= 32767L ) 1152 goto overf; 1153#else 1154if( (j > NBITS) && (exp < 32767L) ) 1155 { 1156 ecleazs( s ); 1157 return; 1158 } 1159#endif 1160if( 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. */ 1176if( rcntrl == 0 ) 1177 goto mdfin; 1178/* Set up rounding parameters if the control register changed. */ 1179if( 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 1239if( (exp <= 0) && (rndprc != NBITS) ) 1240#else 1241if( (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 */ 1250r = s[rw] & rmsk; 1251if( 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 } 1262s[rw] &= ~rmsk; 1263if( (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 } 1280mddone: 1281#ifdef IBMPC 1282if( (exp <= 0) && (rndprc != NBITS) ) 1283#else 1284if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) ) 1285#endif 1286 { 1287 eshup1(s); 1288 } 1289if( s[2] != 0 ) 1290 { /* overflow on roundoff */ 1291 eshdn1(s); 1292 exp += 1; 1293 } 1294mdfin: 1295s[NI-1] = 0; 1296if( exp >= 32767L ) 1297 { 1298#ifndef INFINITY 1299overf: 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 } 1323if( exp < 0 ) 1324 s[1] = 0; 1325else 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 1338static int subflg = 0; 1339 1340void esub( a, b, c ) 1341unsigned short *a, *b, *c; 1342{ 1343 1344#ifdef NANS 1345if( eisnan(a) ) 1346 { 1347 emov (a, c); 1348 return; 1349 } 1350if( 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 */ 1358if( eisinf(a) && eisinf(b) && ((eisneg (a) ^ eisneg (b)) == 0)) 1359 { 1360 mtherr( "esub", DOMAIN ); 1361 enan( c, NBITS ); 1362 return; 1363 } 1364#endif 1365subflg = 1; 1366eadd1( 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*/ 1376void eadd( a, b, c ) 1377unsigned short *a, *b, *c; 1378{ 1379 1380#ifdef NANS 1381/* NaN plus anything is a NaN. */ 1382if( eisnan(a) ) 1383 { 1384 emov(a,c); 1385 return; 1386 } 1387if( 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 */ 1395if( eisinf(a) && eisinf(b) 1396 && ((eisneg(a) ^ eisneg(b)) != 0) ) 1397 { 1398 mtherr( "eadd", DOMAIN ); 1399 enan( c, NBITS ); 1400 return; 1401 } 1402#endif 1403subflg = 0; 1404eadd1( a, b, c ); 1405} 1406 1407void eadd1( a, b, c ) 1408unsigned short *a, *b, *c; 1409{ 1410unsigned short ai[NI], bi[NI], ci[NI]; 1411int i, lost, j, k; 1412long lt, lta, ltb; 1413 1414#ifdef INFINITY 1415if( eisinf(a) ) 1416 { 1417 emov(a,c); 1418 if( subflg ) 1419 eneg(c); 1420 return; 1421 } 1422if( eisinf(b) ) 1423 { 1424 emov(b,c); 1425 return; 1426 } 1427#endif 1428emovi( a, ai ); 1429emovi( b, bi ); 1430if( subflg ) 1431 ai[0] = ~ai[0]; 1432 1433/* compare exponents */ 1434lta = ai[E]; 1435ltb = bi[E]; 1436lt = lta - ltb; 1437if( 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 } 1445lost = 0; 1446if( 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 } 1453else 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 } 1499if( ai[0] == bi[0] ) 1500 { 1501 eaddm( ai, bi ); 1502 subflg = 0; 1503 } 1504else 1505 { 1506 esubm( ai, bi ); 1507 subflg = 1; 1508 } 1509emdnorm( bi, lost, subflg, ltb, 64 ); 1510 1511done: 1512emovo( 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*/ 1523void ediv( a, b, c ) 1524unsigned short *a, *b, *c; 1525{ 1526unsigned short ai[NI], bi[NI]; 1527int i, sign; 1528long 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. */ 1532sign = eisneg(a) ^ eisneg(b); 1533 1534#ifdef NANS 1535/* Return any NaN input. */ 1536if( eisnan(a) ) 1537 { 1538 emov(a,c); 1539 return; 1540 } 1541if( eisnan(b) ) 1542 { 1543 emov(b,c); 1544 return; 1545 } 1546/* Zero over zero, or infinity over infinity, is a NaN. */ 1547if( ((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 1557if( eisinf(b) ) 1558 { 1559 einfin(c); 1560 goto divsign; 1561 } 1562if( eisinf(a) ) 1563 { 1564 eclear(c); 1565 goto divsign; 1566 } 1567#endif 1568emovi( a, ai ); 1569emovi( b, bi ); 1570lta = ai[E]; 1571ltb = bi[E]; 1572if( 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 } 1585dnzro1: 1586 1587if( 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 } 1601dnzro2: 1602 1603i = edivm( ai, bi ); 1604/* calculate exponent */ 1605lt = ltb - lta + EXONE; 1606emdnorm( bi, i, 0, lt, 64 ); 1607emovo( bi, c ); 1608 1609divsign: 1610 1611if( sign ) 1612 *(c+(NE-1)) |= 0x8000; 1613else 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*/ 1625void emul( a, b, c ) 1626unsigned short *a, *b, *c; 1627{ 1628unsigned short ai[NI], bi[NI]; 1629int i, j, sign; 1630long 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. */ 1634sign = eisneg(a) ^ eisneg(b); 1635 1636#ifdef NANS 1637/* NaN times anything is the same NaN. */ 1638if( eisnan(a) ) 1639 { 1640 emov(a,c); 1641 return; 1642 } 1643if( eisnan(b) ) 1644 { 1645 emov(b,c); 1646 return; 1647 } 1648/* Zero times infinity is a NaN. */ 1649if( (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 1659if( eisinf(a) || eisinf(b) ) 1660 { 1661 einfin(c); 1662 goto mulsign; 1663 } 1664#endif 1665emovi( a, ai ); 1666emovi( b, bi ); 1667lta = ai[E]; 1668ltb = bi[E]; 1669if( 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 } 1682mnzer1: 1683 1684if( 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 } 1697mnzer2: 1698 1699/* Multiply significands */ 1700j = emulm( ai, bi ); 1701/* calculate exponent */ 1702lt = lta + ltb - (EXONE - 1); 1703emdnorm( bi, j, 0, lt, 64 ); 1704emovo( bi, c ); 1705/* IEEE says sign is "-" if and only if operands have opposite signs. */ 1706mulsign: 1707if( sign ) 1708 *(c+(NE-1)) |= 0x8000; 1709else 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*/ 1722void e53toe( pe, y ) 1723unsigned short *pe, *y; 1724{ 1725#ifdef DEC 1726 1727dectoe( pe, y ); /* see etodec.c */ 1728 1729#else 1730 1731register unsigned short r; 1732register unsigned short *p, *e; 1733unsigned short yy[NI]; 1734int denorm, k; 1735 1736e = pe; 1737denorm = 0; /* flag if denormalized number */ 1738ecleaz(yy); 1739#ifdef IBMPC 1740e += 3; 1741#endif 1742r = *e; 1743yy[0] = 0; 1744if( r & 0x8000 ) 1745 yy[0] = 0xffff; 1746yy[M] = (r & 0x0f) | 0x10; 1747r &= ~0x800f; /* strip sign and 4 significand bits */ 1748#ifdef INFINITY 1749if( 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 1775r >>= 4; 1776/* If zero exponent, then the significand is denormalized. 1777 * So, take back the understood high significand bit. */ 1778if( r == 0 ) 1779 { 1780 denorm = 1; 1781 yy[M] &= ~0x10; 1782 } 1783r += EXONE - 01777; 1784yy[E] = r; 1785p = &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 ); 1798if( 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 } 1805emovo( yy, y ); 1806#endif /* not DEC */ 1807} 1808 1809void e64toe( pe, y ) 1810unsigned short *pe, *y; 1811{ 1812unsigned short yy[NI]; 1813unsigned short *p, *q, *e; 1814int i; 1815 1816e = pe; 1817p = yy; 1818for( i=0; i<NE-5; i++ ) 1819 *p++ = 0; 1820#ifdef IBMPC 1821for( i=0; i<5; i++ ) 1822 *p++ = *e++; 1823#endif 1824#ifdef DEC 1825for( i=0; i<5; i++ ) 1826 *p++ = *e++; 1827#endif 1828#ifdef MIEEE 1829p = &yy[0] + (NE-1); 1830*p-- = *e++; 1831++e; 1832for( 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. */ 1839if((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. */ 1850p = &yy[NE-1]; 1851if ((*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 { 1876bigend_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 1890p = yy; 1891q = y; 1892for( i=0; i<NE; i++ ) 1893 *q++ = *p++; 1894} 1895 1896void e113toe(pe,y) 1897unsigned short *pe, *y; 1898{ 1899register unsigned short r; 1900unsigned short *e, *p; 1901unsigned short yy[NI]; 1902int denorm, i; 1903 1904e = pe; 1905denorm = 0; 1906ecleaz(yy); 1907#ifdef IBMPC 1908e += 7; 1909#endif 1910r = *e; 1911yy[0] = 0; 1912if( r & 0x8000 ) 1913 yy[0] = 0xffff; 1914r &= 0x7fff; 1915#ifdef INFINITY 1916if( 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 */ 1946yy[E] = r; 1947p = &yy[M + 1]; 1948#ifdef IBMPC 1949for( i=0; i<7; i++ ) 1950 *p++ = *(--e); 1951#endif 1952#ifdef MIEEE 1953++e; 1954for( i=0; i<7; i++ ) 1955 *p++ = *e++; 1956#endif 1957/* If denormal, remove the implied bit; else shift down 1. */ 1958if( r == 0 ) 1959 { 1960 yy[M] = 0; 1961 } 1962else 1963 { 1964 yy[M] = 1; 1965 eshift( yy, -1 ); 1966 } 1967emovo(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*/ 1977void e24toe( pe, y ) 1978unsigned short *pe, *y; 1979{ 1980register unsigned short r; 1981register unsigned short *p, *e; 1982unsigned short yy[NI]; 1983int denorm, k; 1984 1985e = pe; 1986denorm = 0; /* flag if denormalized number */ 1987ecleaz(yy); 1988#ifdef IBMPC 1989e += 1; 1990#endif 1991#ifdef DEC 1992e += 1; 1993#endif 1994r = *e; 1995yy[0] = 0; 1996if( r & 0x8000 ) 1997 yy[0] = 0xffff; 1998yy[M] = (r & 0x7f) | 0200; 1999r &= ~0x807f; /* strip sign and 7 significand bits */ 2000#ifdef INFINITY 2001if( 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 2025r >>= 7; 2026/* If zero exponent, then the significand is denormalized. 2027 * So, take back the understood high significand bit. */ 2028if( r == 0 ) 2029 { 2030 denorm = 1; 2031 yy[M] &= ~0200; 2032 } 2033r += EXONE - 0177; 2034yy[E] = r; 2035p = &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 ); 2047if( 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 } 2054emovo( yy, y ); 2055} 2056 2057void etoe113(x,e) 2058unsigned short *x, *e; 2059{ 2060unsigned short xi[NI]; 2061long exp; 2062int rndsav; 2063 2064#ifdef NANS 2065if( eisnan(x) ) 2066 { 2067 enan( e, 113 ); 2068 return; 2069 } 2070#endif 2071emovi( x, xi ); 2072exp = (long )xi[E]; 2073#ifdef INFINITY 2074if( eisinf(x) ) 2075 goto nonorm; 2076#endif 2077/* round off to nearest or even */ 2078rndsav = rndprc; 2079rndprc = 113; 2080emdnorm( xi, 0, 0, exp, 64 ); 2081rndprc = rndsav; 2082nonorm: 2083toe113 (xi, e); 2084} 2085 2086/* move out internal format to ieee long double */ 2087static void toe113(a,b) 2088unsigned short *a, *b; 2089{ 2090register unsigned short *p, *q; 2091unsigned short i; 2092 2093#ifdef NANS 2094if( eiisnan(a) ) 2095 { 2096 enan( b, 113 ); 2097 return; 2098 } 2099#endif 2100p = a; 2101#ifdef MIEEE 2102q = b; 2103#else 2104q = b + 7; /* point to output exponent */ 2105#endif 2106 2107/* If not denormal, delete the implied bit. */ 2108if( a[E] != 0 ) 2109 { 2110 eshup1 (a); 2111 } 2112/* combine sign and exponent */ 2113i = *p++; 2114#ifdef MIEEE 2115if( i ) 2116 *q++ = *p++ | 0x8000; 2117else 2118 *q++ = *p++; 2119#else 2120if( i ) 2121 *q-- = *p++ | 0x8000; 2122else 2123 *q-- = *p++; 2124#endif 2125/* skip over guard word */ 2126++p; 2127/* move the significand */ 2128#ifdef MIEEE 2129for (i = 0; i < 7; i++) 2130 *q++ = *p++; 2131#else 2132for (i = 0; i < 7; i++) 2133 *q-- = *p++; 2134#endif 2135} 2136 2137 2138void etoe64( x, e ) 2139unsigned short *x, *e; 2140{ 2141unsigned short xi[NI]; 2142long exp; 2143int rndsav; 2144 2145#ifdef NANS 2146if( eisnan(x) ) 2147 { 2148 enan( e, 64 ); 2149 return; 2150 } 2151#endif 2152emovi( x, xi ); 2153exp = (long )xi[E]; /* adjust exponent for offset */ 2154#ifdef INFINITY 2155if( eisinf(x) ) 2156 goto nonorm; 2157#endif 2158/* round off to nearest or even */ 2159rndsav = rndprc; 2160rndprc = 64; 2161emdnorm( xi, 0, 0, exp, 64 ); 2162rndprc = rndsav; 2163nonorm: 2164toe64( xi, e ); 2165} 2166 2167/* move out internal format to ieee long double */ 2168static void toe64( a, b ) 2169unsigned short *a, *b; 2170{ 2171register unsigned short *p, *q; 2172unsigned short i; 2173 2174#ifdef NANS 2175if( eiisnan(a) ) 2176 { 2177 enan( b, 64 ); 2178 return; 2179 } 2180#endif 2181#ifdef IBMPC 2182/* Shift Intel denormal significand down 1. */ 2183if( a[E] == 0 ) 2184 eshdn1(a); 2185#endif 2186p = a; 2187#ifdef MIEEE 2188q = b; 2189#else 2190q = 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 */ 2198i = *p++; 2199#ifdef MIEEE 2200if( i ) 2201 *q++ = *p++ | 0x8000; 2202else 2203 *q++ = *p++; 2204*q++ = 0; 2205#else 2206if( i ) 2207 *q-- = *p++ | 0x8000; 2208else 2209 *q-- = *p++; 2210#endif 2211/* skip over guard word */ 2212++p; 2213/* move the significand */ 2214#ifdef MIEEE 2215for( i=0; i<4; i++ ) 2216 *q++ = *p++; 2217#else 2218#ifdef INFINITY 2219if (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 2229for( 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 2244void etoe53( x, e ) 2245unsigned short *x, *e; 2246{ 2247etodec( x, e ); /* see etodec.c */ 2248} 2249 2250static void toe53( x, y ) 2251unsigned short *x, *y; 2252{ 2253todec( x, y ); 2254} 2255 2256#else 2257 2258void etoe53( x, e ) 2259unsigned short *x, *e; 2260{ 2261unsigned short xi[NI]; 2262long exp; 2263int rndsav; 2264 2265#ifdef NANS 2266if( eisnan(x) ) 2267 { 2268 enan( e, 53 ); 2269 return; 2270 } 2271#endif 2272emovi( x, xi ); 2273exp = (long )xi[E] - (EXONE - 0x3ff); /* adjust exponent for offsets */ 2274#ifdef INFINITY 2275if( eisinf(x) ) 2276 goto nonorm; 2277#endif 2278/* round off to nearest or even */ 2279rndsav = rndprc; 2280rndprc = 53; 2281emdnorm( xi, 0, 0, exp, 64 ); 2282rndprc = rndsav; 2283nonorm: 2284toe53( xi, e ); 2285} 2286 2287 2288static void toe53( x, y ) 2289unsigned short *x, *y; 2290{ 2291unsigned short i; 2292unsigned short *p; 2293 2294 2295#ifdef NANS 2296if( eiisnan(x) ) 2297 { 2298 enan( y, 53 ); 2299 return; 2300 } 2301#endif 2302p = &x[0]; 2303#ifdef IBMPC 2304y += 3; 2305#endif 2306*y = 0; /* output high order */ 2307if( *p++ ) 2308 *y = 0x8000; /* output sign bit */ 2309 2310i = *p++; 2311if( 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 } 2342if( i == 0 ) 2343 { 2344 (void )eshift( x, 4 ); 2345 } 2346else 2347 { 2348 i <<= 4; 2349 (void )eshift( x, 5 ); 2350 } 2351i |= *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*/ 2376void etoe24( x, e ) 2377unsigned short *x, *e; 2378{ 2379long exp; 2380unsigned short xi[NI]; 2381int rndsav; 2382 2383#ifdef NANS 2384if( eisnan(x) ) 2385 { 2386 enan( e, 24 ); 2387 return; 2388 } 2389#endif 2390emovi( x, xi ); 2391exp = (long )xi[E] - (EXONE - 0177); /* adjust exponent for offsets */ 2392#ifdef INFINITY 2393if( eisinf(x) ) 2394 goto nonorm; 2395#endif 2396/* round off to nearest or even */ 2397rndsav = rndprc; 2398rndprc = 24; 2399emdnorm( xi, 0, 0, exp, 64 ); 2400rndprc = rndsav; 2401nonorm: 2402toe24( xi, e ); 2403} 2404 2405static void toe24( x, y ) 2406unsigned short *x, *y; 2407{ 2408unsigned short i; 2409unsigned short *p; 2410 2411#ifdef NANS 2412if( eiisnan(x) ) 2413 { 2414 enan( y, 24 ); 2415 return; 2416 } 2417#endif 2418p = &x[0]; 2419#ifdef IBMPC 2420y += 1; 2421#endif 2422#ifdef DEC 2423y += 1; 2424#endif 2425*y = 0; /* output high order */ 2426if( *p++ ) 2427 *y = 0x8000; /* output sign bit */ 2428 2429i = *p++; 2430if( 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 } 2459if( i == 0 ) 2460 { 2461 (void )eshift( x, 7 ); 2462 } 2463else 2464 { 2465 i <<= 7; 2466 (void )eshift( x, 8 ); 2467 } 2468i |= *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 */ 2493int ecmp( a, b ) 2494unsigned short *a, *b; 2495{ 2496unsigned short ai[NI], bi[NI]; 2497register unsigned short *p, *q; 2498register int i; 2499int msign; 2500 2501#ifdef NANS 2502if (eisnan (a) || eisnan (b)) 2503 return( -2 ); 2504#endif 2505emovi( a, ai ); 2506p = ai; 2507emovi( b, bi ); 2508q = bi; 2509 2510if( *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); 2521nzro: 2522 if( *p == 0 ) 2523 return( 1 ); 2524 else 2525 return( -1 ); 2526 } 2527/* both are the same sign */ 2528if( *p == 0 ) 2529 msign = 1; 2530else 2531 msign = -1; 2532i = NI-1; 2533do 2534 { 2535 if( *p++ != *q++ ) 2536 { 2537 goto diff; 2538 } 2539 } 2540while( --i > 0 ); 2541 2542return(0); /* equality */ 2543 2544 2545 2546diff: 2547 2548if( *(--p) > *(--q) ) 2549 return( msign ); /* p is bigger */ 2550else 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 */ 2562void eround( x, y ) 2563unsigned short *x, *y; 2564{ 2565 2566eadd( ehalf, x, y ); 2567efloor( 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*/ 2581void ltoe( lp, y ) 2582long *lp; /* lp is the memory address of a long integer */ 2583unsigned short *y; /* y is the address of a short */ 2584{ 2585unsigned short yi[NI]; 2586unsigned long ll; 2587int k; 2588 2589ecleaz( yi ); 2590if( *lp < 0 ) 2591 { 2592 ll = (unsigned long )( -(*lp) ); /* make it positive */ 2593 yi[0] = 0xffff; /* put correct sign in the e type number */ 2594 } 2595else 2596 { 2597 ll = (unsigned long )( *lp ); 2598 } 2599/* move the long integer to yi significand area */ 2600if( 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 } 2608else 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 } 2614if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */ 2615 ecleaz( yi ); /* it was zero */ 2616else 2617 yi[E] -= (unsigned short )k; /* subtract shift count from exponent */ 2618emovo( 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*/ 2629void ultoe( lp, y ) 2630unsigned long *lp; /* lp is the memory address of a long integer */ 2631unsigned short *y; /* y is the address of a short */ 2632{ 2633unsigned short yi[NI]; 2634unsigned long ll; 2635int k; 2636 2637ecleaz( yi ); 2638ll = *lp; 2639 2640/* move the long integer to ayi significand area */ 2641if( 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 } 2649else 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 } 2655if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */ 2656 ecleaz( yi ); /* it was zero */ 2657else 2658 yi[E] -= (unsigned short )k; /* subtract shift count from exponent */ 2659emovo( 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*/ 2673void eifrac( x, i, frac ) 2674unsigned short *x; 2675long *i; 2676unsigned short *frac; 2677{ 2678unsigned short xi[NI]; 2679int j, k; 2680unsigned long ll; 2681 2682emovi( x, xi ); 2683k = (int )xi[E] - (EXONE - 1); 2684if( k <= 0 ) 2685 { 2686/* if exponent <= 0, integer = 0 and real output is fraction */ 2687 *i = 0L; 2688 emovo( xi, frac ); 2689 return; 2690 } 2691if( 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 } 2704if( 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 } 2724else 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 } 2732xi[0] = 0; 2733xi[E] = EXONE - 1; 2734xi[M] = 0; 2735if( (k = enormlz( xi )) > NBITS ) 2736 ecleaz( xi ); 2737else 2738 xi[E] -= (unsigned short )k; 2739 2740emovo( 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*/ 2754void euifrac( x, i, frac ) 2755unsigned short *x; 2756unsigned long *i; 2757unsigned short *frac; 2758{ 2759unsigned short xi[NI]; 2760int j, k; 2761unsigned long ll; 2762 2763emovi( x, xi ); 2764k = (int )xi[E] - (EXONE - 1); 2765if( k <= 0 ) 2766 { 2767/* if exponent <= 0, integer = 0 and argument is fraction */ 2768 *i = 0L; 2769 emovo( xi, frac ); 2770 return; 2771 } 2772if( 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 } 2781else 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 } 2799else 2800 { 2801/* shift not more than 16 bits */ 2802 eshift( xi, k ); 2803 *i = (long )xi[M] & 0xffff; 2804 } 2805 2806if( xi[0] ) /* A negative value yields unsigned integer 0. */ 2807 *i = 0L; 2808 2809xi[0] = 0; 2810xi[E] = EXONE - 1; 2811xi[M] = 0; 2812if( (k = enormlz( xi )) > NBITS ) 2813 ecleaz( xi ); 2814else 2815 xi[E] -= (unsigned short )k; 2816 2817emovo( 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*/ 2828int eshift( x, sc ) 2829unsigned short *x; 2830int sc; 2831{ 2832unsigned short lost; 2833unsigned short *p; 2834 2835if( sc == 0 ) 2836 return( 0 ); 2837 2838lost = 0; 2839p = x + NI-1; 2840 2841if( 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 } 2865else 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 } 2885if( lost ) 2886 lost = 1; 2887return( (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*/ 2898int enormlz(x) 2899unsigned short x[]; 2900{ 2901register unsigned short *p; 2902int sc; 2903 2904sc = 0; 2905p = &x[M]; 2906if( *p != 0 ) 2907 goto normdn; 2908++p; 2909if( *p & 0x8000 ) 2910 return( 0 ); /* already normalized */ 2911while( *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 */ 2922while( (*p & 0xff00) == 0 ) 2923 { 2924 eshup8(x); 2925 sc += 8; 2926 } 2927/* now shift 1 bit at a time */ 2928while( (*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 } 2938return( sc ); 2939 2940/* Normalize by shifting down out of the high guard word 2941 of the significand */ 2942normdn: 2943 2944if( *p & 0xff00 ) 2945 { 2946 eshdn8(x); 2947 sc -= 8; 2948 } 2949while( *p != 0 ) 2950 { 2951 eshdn1(x); 2952 sc -= 1; 2953 2954 if( sc < -NBITS ) 2955 { 2956 mtherr( "enormlz", OVERFLOW ); 2957 return( sc ); 2958 } 2959 } 2960return( 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 2974static 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 3004static 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 3034static 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 3050static 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 3067void e24toasc( x, string, ndigs ) 3068unsigned short x[]; 3069char *string; 3070int ndigs; 3071{ 3072unsigned short w[NI]; 3073 3074e24toe( x, w ); 3075etoasc( w, string, ndigs ); 3076} 3077 3078 3079void e53toasc( x, string, ndigs ) 3080unsigned short x[]; 3081char *string; 3082int ndigs; 3083{ 3084unsigned short w[NI]; 3085 3086e53toe( x, w ); 3087etoasc( w, string, ndigs ); 3088} 3089 3090 3091void e64toasc( x, string, ndigs ) 3092unsigned short x[]; 3093char *string; 3094int ndigs; 3095{ 3096unsigned short w[NI]; 3097 3098e64toe( x, w ); 3099etoasc( w, string, ndigs ); 3100} 3101 3102void e113toasc (x, string, ndigs) 3103unsigned short x[]; 3104char *string; 3105int ndigs; 3106{ 3107unsigned short w[NI]; 3108 3109e113toe (x, w); 3110etoasc (w, string, ndigs); 3111} 3112 3113 3114void etoasc( x, string, ndigs ) 3115unsigned short x[]; 3116char *string; 3117int ndigs; 3118{ 3119long digit; 3120unsigned short y[NI], t[NI], u[NI], w[NI]; 3121unsigned short *p, *r, *ten; 3122unsigned short sign; 3123int i, j, k, expon, rndsav; 3124char *s, *ss; 3125unsigned short m; 3126 3127rndsav = rndprc; 3128#ifdef NANS 3129if( eisnan(x) ) 3130 { 3131 sprintf( string, " NaN " ); 3132 goto bxit; 3133 } 3134#endif 3135rndprc = NBITS; /* set to full precision */ 3136emov( x, y ); /* retain external format */ 3137if( y[NE-1] & 0x8000 ) 3138 { 3139 sign = 0xffff; 3140 y[NE-1] &= 0x7fff; 3141 } 3142else 3143 { 3144 sign = 0; 3145 } 3146expon = 0; 3147ten = &etens[NTEN][0]; 3148emov( eone, t ); 3149/* Test for zero exponent */ 3150if( 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 } 3159tnzro: 3160 3161/* Test for infinity. 3162 */ 3163if( 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 */ 3175if( (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 */ 3183i = ecmp( eone, y ); 3184if( i == 0 ) 3185 goto isone; 3186 3187if( 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; 3195do 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; 3206noint: 3207 p += NE; 3208 m >>= 1; 3209 } 3210while( 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 } 3233else 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 } 3292isone: 3293/* Find the first (leading) digit. */ 3294emovi( t, w ); 3295emovz( w, t ); 3296emovi( y, w ); 3297emovz( w, y ); 3298eiremain( t, y ); 3299digit = equot[NI-1]; 3300while( (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 } 3311s = string; 3312if( sign ) 3313 *s++ = '-'; 3314else 3315 *s++ = ' '; 3316/* Examine number of digits requested by caller. */ 3317if( ndigs < 0 ) 3318 ndigs = 0; 3319if( ndigs > NDEC ) 3320 ndigs = NDEC; 3321if( digit == 10 ) 3322 { 3323 *s++ = '1'; 3324 *s++ = '.'; 3325 if( ndigs > 0 ) 3326 { 3327 *s++ = '0'; 3328 ndigs -= 1; 3329 } 3330 expon += 1; 3331 } 3332else 3333 { 3334 *s++ = (char )digit + '0'; 3335 *s++ = '.'; 3336 } 3337/* Generate digits after the decimal point. */ 3338for( 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 } 3349digit = equot[NI-1]; 3350--s; 3351ss = s; 3352/* round off the ASCII string */ 3353if( 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 */ 3365roun: 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 } 3392doexp: 3393/* 3394if( expon >= 0 ) 3395 sprintf( ss, "e+%d", expon ); 3396else 3397 sprintf( ss, "e%d", expon ); 3398*/ 3399 sprintf( ss, "E%d", expon ); 3400bxit: 3401rndprc = 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 */ 3427void asctoe24( s, y ) 3428char *s; 3429unsigned short *y; 3430{ 3431asctoeg( s, y, 24 ); 3432} 3433 3434 3435/* ASCII to double */ 3436void asctoe53( s, y ) 3437char *s; 3438unsigned short *y; 3439{ 3440#ifdef DEC 3441asctoeg( s, y, 56 ); 3442#else 3443asctoeg( s, y, 53 ); 3444#endif 3445} 3446 3447 3448/* ASCII to long double */ 3449void asctoe64( s, y ) 3450char *s; 3451unsigned short *y; 3452{ 3453asctoeg( s, y, 64 ); 3454} 3455 3456/* ASCII to 128-bit long double */ 3457void asctoe113 (s, y) 3458char *s; 3459unsigned short *y; 3460{ 3461asctoeg( s, y, 113 ); 3462} 3463 3464/* ASCII to super double */ 3465void asctoe( s, y ) 3466char *s; 3467unsigned short *y; 3468{ 3469asctoeg( s, y, NBITS ); 3470} 3471 3472/* Space to make a copy of the input string: */ 3473static char lstr[82] = {0}; 3474 3475void asctoeg( ss, y, oprec ) 3476char *ss; 3477unsigned short *y; 3478int oprec; 3479{ 3480unsigned short yy[NI], xt[NI], tt[NI]; 3481int esign, decflg, sgnflg, nexp, exp, prec, lost; 3482int k, trail, c, rndsav; 3483long lexp; 3484unsigned short nsign, *p; 3485char *sp, *s; 3486 3487/* Copy the input string. */ 3488s = ss; 3489while( *s == ' ' ) /* skip leading spaces */ 3490 ++s; 3491sp = lstr; 3492for( k=0; k<79; k++ ) 3493 { 3494 if( (*sp++ = *s++) == '\0' ) 3495 break; 3496 } 3497*sp = '\0'; 3498s = lstr; 3499 3500rndsav = rndprc; 3501rndprc = NBITS; /* Set to full precision */ 3502lost = 0; 3503nsign = 0; 3504decflg = 0; 3505sgnflg = 0; 3506nexp = 0; 3507exp = 0; 3508prec = 0; 3509ecleaz( yy ); 3510trail = 0; 3511 3512nxtcom: 3513k = *s - '0'; 3514if( (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 3568switch( *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 } 3610donchr: 3611++s; 3612goto nxtcom; 3613 3614/* Exponent interpretation */ 3615expnt: 3616 3617esign = 1; 3618exp = 0; 3619++s; 3620/* check for + or - */ 3621if( *s == '-' ) 3622 { 3623 esign = -1; 3624 ++s; 3625 } 3626if( *s == '+' ) 3627 ++s; 3628while( (*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 } 3640if( esign < 0 ) 3641 exp = -exp; 3642if( exp > 4932 ) 3643 { 3644infinite: 3645 ecleaz(yy); 3646 yy[E] = 0x7fff; /* infinity */ 3647 goto aexit; 3648 } 3649if( exp < -4977 ) 3650 { 3651zero: 3652 ecleaz(yy); 3653 goto aexit; 3654 } 3655 3656daldone: 3657nexp = exp - nexp; 3658/* Pad trailing zeros to minimize power of 10, per IEEE spec. */ 3659while( (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 } 3671if( (k = enormlz(yy)) > NBITS ) 3672 { 3673 ecleaz(yy); 3674 goto aexit; 3675 } 3676lexp = (EXONE - 1 + NBITS) - k; 3677emdnorm( 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 */ 3687lexp = yy[E]; 3688if( nexp == 0 ) 3689 { 3690 k = 0; 3691 goto expdon; 3692 } 3693esign = 1; 3694if( 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 } 3707p = &etens[NTEN][0]; 3708emov( eone, xt ); 3709exp = 1; 3710do 3711 { 3712 if( exp & nexp ) 3713 emul( p, xt, xt ); 3714 p -= NE; 3715 exp = exp + exp; 3716 } 3717while( exp <= MAXP ); 3718 3719emovi( xt, tt ); 3720if( esign < 0 ) 3721 { 3722 lexp -= tt[E]; 3723 k = edivm( tt, yy ); 3724 lexp += EXONE; 3725 } 3726else 3727 { 3728 lexp += tt[E]; 3729 k = emulm( tt, yy ); 3730 lexp -= EXONE - 1; 3731 } 3732 3733expdon: 3734 3735/* Round and convert directly to the destination type */ 3736if( oprec == 53 ) 3737 lexp -= EXONE - 0x3ff; 3738else if( oprec == 24 ) 3739 lexp -= EXONE - 0177; 3740#ifdef DEC 3741else if( oprec == 56 ) 3742 lexp -= EXONE - 0201; 3743#endif 3744rndprc = oprec; 3745emdnorm( yy, k, 0, lexp, 64 ); 3746 3747aexit: 3748 3749rndprc = rndsav; 3750yy[0] = nsign; 3751switch( 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 */ 3785static unsigned short bmask[] = { 37860xffff, 37870xfffe, 37880xfffc, 37890xfff8, 37900xfff0, 37910xffe0, 37920xffc0, 37930xff80, 37940xff00, 37950xfe00, 37960xfc00, 37970xf800, 37980xf000, 37990xe000, 38000xc000, 38010x8000, 38020x0000, 3803}; 3804 3805void efloor( x, y ) 3806unsigned short x[], y[]; 3807{ 3808register unsigned short *p; 3809int e, expon, i; 3810unsigned short f[NE]; 3811 3812emov( x, f ); /* leave in external format */ 3813expon = (int )f[NE-1]; 3814e = (expon & 0x7fff) - (EXONE - 1); 3815if( e <= 0 ) 3816 { 3817 eclear(y); 3818 goto isitneg; 3819 } 3820/* number of bits to clear out */ 3821e = NBITS - e; 3822emov( f, y ); 3823if( e <= 0 ) 3824 return; 3825 3826p = &y[0]; 3827while( 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 */ 3835isitneg: 3836 3837if( (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 */ 3860void efrexp( x, exp, s ) 3861unsigned short x[]; 3862long *exp; 3863unsigned short s[]; 3864{ 3865unsigned short xi[NI]; 3866long li; 3867 3868emovi( x, xi ); 3869li = (long )((short )xi[1]); 3870 3871if( li == 0 ) 3872 { 3873 li -= enormlz( xi ); 3874 } 3875xi[1] = 0x3ffe; 3876emovo( 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 */ 3889void eldexp( x, pwr2, y ) 3890unsigned short x[]; 3891long pwr2; 3892unsigned short y[]; 3893{ 3894unsigned short xi[NI]; 3895long li; 3896int i; 3897 3898emovi( x, xi ); 3899li = xi[1]; 3900li += pwr2; 3901i = 0; 3902emdnorm( xi, i, i, li, 64 ); 3903emovo( xi, y ); 3904} 3905 3906 3907/* c = remainder after dividing b by a 3908 * Least significant integer quotient bits left in equot[]. 3909 */ 3910void eremain( a, b, c ) 3911unsigned short a[], b[], c[]; 3912{ 3913unsigned short den[NI], num[NI]; 3914 3915#ifdef NANS 3916if( eisinf(b) || (ecmp(a,ezero) == 0) || eisnan(a) || eisnan(b)) 3917 { 3918 enan( c, NBITS ); 3919 return; 3920 } 3921#endif 3922if( ecmp(a,ezero) == 0 ) 3923 { 3924 mtherr( "eremain", SING ); 3925 eclear( c ); 3926 return; 3927 } 3928emovi( a, den ); 3929emovi( b, num ); 3930eiremain( den, num ); 3931/* Sign of remainder = sign of quotient */ 3932if( a[0] == b[0] ) 3933 num[0] = 0; 3934else 3935 num[0] = 0xffff; 3936emovo( num, c ); 3937} 3938 3939 3940void eiremain( den, num ) 3941unsigned short den[], num[]; 3942{ 3943long ld, ln; 3944unsigned short j; 3945 3946ld = den[E]; 3947ld -= enormlz( den ); 3948ln = num[E]; 3949ln -= enormlz( num ); 3950ecleaz( equot ); 3951while( 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 } 3967emdnorm( num, 0, 0, ln, 0 ); 3968} 3969 3970/* NaN bit patterns 3971 */ 3972#ifdef MIEEE 3973unsigned short nan113[8] = { 3974 0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}; 3975unsigned short nan64[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff}; 3976unsigned short nan53[4] = {0x7fff, 0xffff, 0xffff, 0xffff}; 3977unsigned short nan24[2] = {0x7fff, 0xffff}; 3978#endif 3979 3980#ifdef IBMPC 3981unsigned short nan113[8] = {0, 0, 0, 0, 0, 0, 0xc000, 0xffff}; 3982unsigned short nan64[6] = {0, 0, 0, 0xc000, 0xffff, 0}; 3983unsigned short nan53[4] = {0, 0, 0, 0xfff8}; 3984unsigned short nan24[2] = {0, 0xffc0}; 3985#endif 3986 3987 3988void enan (nan, size) 3989unsigned short *nan; 3990int size; 3991{ 3992int i, n; 3993unsigned short *p; 3994 3995switch( 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 } 4038for (i=0; i < n; i++) 4039 *nan++ = *p++; 4040} 4041 4042 4043 4044/* Longhand square root. */ 4045 4046static int esqinited = 0; 4047static unsigned short sqrndbit[NI]; 4048 4049void esqrt( x, y ) 4050short *x, *y; 4051{ 4052unsigned short temp[NI], num[NI], sq[NI], xx[NI]; 4053int i, j, k, n, nlups; 4054long m, exp; 4055 4056if( esqinited == 0 ) 4057 { 4058 ecleaz( sqrndbit ); 4059 sqrndbit[NI-2] = 1; 4060 esqinited = 1; 4061 } 4062/* Check for arg <= 0 */ 4063i = ecmp( x, ezero ); 4064if( 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 4080if( 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. */ 4088emovi( x, xx ); 4089m = (long )xx[1]; /* local long word exponent */ 4090if( m == 0 ) 4091 m -= enormlz( xx ); 4092 4093/* Divide exponent by 2 */ 4094m -= 0x3ffe; 4095exp = (unsigned short )( (m / 2) + 0x3ffe ); 4096 4097/* Adjust if exponent odd */ 4098if( (m & 1) != 0 ) 4099 { 4100 if( m > 0 ) 4101 exp += 1; 4102 eshdn1( xx ); 4103 } 4104 4105ecleaz( sq ); 4106ecleaz( num ); 4107n = 8; /* get 8 bits of result per inner loop */ 4108nlups = rndprc; 4109j = 0; 4110 4111while( 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. */ 4143exp += (NBITS - 1) - rndprc; 4144 4145/* Sticky bit = 1 if the remainder is nonzero. */ 4146k = 0; 4147for( i=3; i<NI; i++ ) 4148 k |= (int )num[i]; 4149 4150/* Renormalize and round off. */ 4151emdnorm( sq, k, 0, exp, 64 ); 4152emovo( sq, y ); 4153} 4154