summaryrefslogtreecommitdiff
path: root/src/regress/lib/libc/cephes/ieee.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/regress/lib/libc/cephes/ieee.c')
-rw-r--r--src/regress/lib/libc/cephes/ieee.c4153
1 files changed, 0 insertions, 4153 deletions
diff --git a/src/regress/lib/libc/cephes/ieee.c b/src/regress/lib/libc/cephes/ieee.c
deleted file mode 100644
index e2b8aa7b99..0000000000
--- a/src/regress/lib/libc/cephes/ieee.c
+++ /dev/null
@@ -1,4153 +0,0 @@
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}
generated by cgit v1.2.3-55-g6feb (git 2.39.1) at 2025-03-21 08:52:52 +0000