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