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
einit()190 void einit()
191 {
192 }
193
194 /*
195 ; Clear out entire external format number.
196 ;
197 ; unsigned short x[];
198 ; eclear( x );
199 */
200
eclear(x)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
emov(a,b)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
eabs(x)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
eneg(x)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 */
eisneg(x)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 */
eisinf(x)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 */
eisnan(x)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
einfin(x)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 */
emovi(a,b)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 */
emovo(a,b)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
ecleaz(xi)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
ecleazs(xi)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 */
emovz(a,b)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
eiisnan(x)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
eiisinf(x)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 */
ecmpm(a,b)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
eshdn1(x)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
eshup1(x)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
eshdn8(x)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
eshup8(x)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
eshup6(x)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
eshdn6(x)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
eaddm(x,y)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
esubm(x,y)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
m16m(a,b,c)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
edivm(den,num)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 */
emulm(a,b)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
emdnorm(s,lost,subflg,exp,rcntrl)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
esub(a,b,c)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 */
eadd(a,b,c)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
eadd1(a,b,c)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 */
ediv(a,b,c)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 */
emul(a,b,c)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 */
e53toe(pe,y)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
e64toe(pe,y)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
e113toe(pe,y)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 */
e24toe(pe,y)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
etoe113(x,e)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 */
toe113(a,b)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
etoe64(x,e)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 */
toe64(a,b)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
etoe53(x,e)2244 void etoe53( x, e )
2245 unsigned short *x, *e;
2246 {
2247 etodec( x, e ); /* see etodec.c */
2248 }
2249
toe53(x,y)2250 static void toe53( x, y )
2251 unsigned short *x, *y;
2252 {
2253 todec( x, y );
2254 }
2255
2256 #else
2257
etoe53(x,e)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
toe53(x,y)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 */
etoe24(x,e)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
toe24(x,y)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 */
ecmp(a,b)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 */
eround(x,y)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 */
ltoe(lp,y)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 */
ultoe(lp,y)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 */
eifrac(x,i,frac)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 */
euifrac(x,i,frac)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 */
eshift(x,sc)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 */
enormlz(x)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
e24toasc(x,string,ndigs)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
e53toasc(x,string,ndigs)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
e64toasc(x,string,ndigs)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
e113toasc(x,string,ndigs)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
etoasc(x,string,ndigs)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 */
asctoe24(s,y)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 */
asctoe53(s,y)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 */
asctoe64(s,y)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 */
asctoe113(s,y)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 */
asctoe(s,y)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
asctoeg(ss,y,oprec)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
efloor(x,y)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 */
efrexp(x,exp,s)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 */
eldexp(x,pwr2,y)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 */
eremain(a,b,c)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
eiremain(den,num)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
enan(nan,size)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
esqrt(x,y)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