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