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