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