1 /**
2 * This file has no copyright assigned and is placed in the Public Domain.
3 * This file is part of the mingw-w64 runtime package.
4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
5 */
6 #ifndef _CEPHES_EMATH_H
7 #define _CEPHES_EMATH_H
8
9 /**
10 * This is a workaround for a gcc bug
11 */
12 #define __restrict__
13
14 /* This file is extracted from S L Moshier's ioldoubl.c,
15 * modified for use in MinGW
16 *
17 * Extended precision arithmetic functions for long double I/O.
18 * This program has been placed in the public domain.
19 */
20
21
22 /*
23 * Revision history:
24 *
25 * 5 Jan 84 PDP-11 assembly language version
26 * 6 Dec 86 C language version
27 * 30 Aug 88 100 digit version, improved rounding
28 * 15 May 92 80-bit long double support
29 *
30 * Author: S. L. Moshier.
31 *
32 * 6 Oct 02 Modified for MinGW by inlining utility routines,
33 * removing global variables, and splitting out strtold
34 * from _IO_ldtoa and _IO_ldtostr.
35 *
36 * Danny Smith <dannysmith@users.sourceforge.net>
37 *
38 */
39
40
41 /* ieee.c
42 *
43 * Extended precision IEEE binary floating point arithmetic routines
44 *
45 * Numbers are stored in C language as arrays of 16-bit unsigned
46 * short integers. The arguments of the routines are pointers to
47 * the arrays.
48 *
49 *
50 * External e type data structure, simulates Intel 8087 chip
51 * temporary real format but possibly with a larger significand:
52 *
53 * NE-1 significand words (least significant word first,
54 * most significant bit is normally set)
55 * exponent (value = EXONE for 1.0,
56 * top bit is the sign)
57 *
58 *
59 * Internal data structure of a number (a "word" is 16 bits):
60 *
61 * ei[0] sign word (0 for positive, 0xffff for negative)
62 * ei[1] biased __exponent (value = EXONE for the number 1.0)
63 * ei[2] high guard word (always zero after normalization)
64 * ei[3]
65 * to ei[NI-2] significand (NI-4 significand words,
66 * most significant word first,
67 * most significant bit is set)
68 * ei[NI-1] low guard word (0x8000 bit is rounding place)
69 *
70 *
71 *
72 * Routines for external format numbers
73 *
74 * __asctoe64( string, &d ) ASCII string to long double
75 * __asctoeg( string, e, prec ) ASCII string to specified precision
76 * __e64toe( &d, e ) IEEE long double precision to e type
77 * __eadd( a, b, c ) c = b + a
78 * __eclear(e) e = 0
79 * __ecmp (a, b) Returns 1 if a > b, 0 if a == b,
80 * -1 if a < b, -2 if either a or b is a NaN.
81 * __ediv( a, b, c ) c = b / a
82 * __efloor( a, b ) truncate to integer, toward -infinity
83 * __efrexp( a, exp, s ) extract exponent and significand
84 * __eifrac( e, &l, frac ) e to long integer and e type fraction
85 * __euifrac( e, &l, frac ) e to unsigned long integer and e type fraction
86 * __einfin( e ) set e to infinity, leaving its sign alone
87 * __eldexp( a, n, b ) multiply by 2**n
88 * __emov( a, b ) b = a
89 * __emul( a, b, c ) c = b * a
90 * __eneg(e) e = -e
91 * __eround( a, b ) b = nearest integer value to a
92 * __esub( a, b, c ) c = b - a
93 * __e24toasc( &f, str, n ) single to ASCII string, n digits after decimal
94 * __e53toasc( &d, str, n ) double to ASCII string, n digits after decimal
95 * __e64toasc( &d, str, n ) long double to ASCII string
96 * __etoasc( e, str, n ) e to ASCII string, n digits after decimal
97 * __etoe24( e, &f ) convert e type to IEEE single precision
98 * __etoe53( e, &d ) convert e type to IEEE double precision
99 * __etoe64( e, &d ) convert e type to IEEE long double precision
100 * __eisneg( e ) 1 if sign bit of e != 0, else 0
101 * __eisinf( e ) 1 if e has maximum exponent (non-IEEE)
102 * or is infinite (IEEE)
103 * __eisnan( e ) 1 if e is a NaN
104 * __esqrt( a, b ) b = square root of a
105 *
106 *
107 * Routines for internal format numbers
108 *
109 * __eaddm( ai, bi ) add significands, bi = bi + ai
110 * __ecleaz(ei) ei = 0
111 * __ecleazs(ei) set ei = 0 but leave its sign alone
112 * __ecmpm( ai, bi ) compare significands, return 1, 0, or -1
113 * __edivm( ai, bi ) divide significands, bi = bi / ai
114 * __emdnorm(ai,l,s,exp) normalize and round off
115 * __emovi( a, ai ) convert external a to internal ai
116 * __emovo( ai, a ) convert internal ai to external a
117 * __emovz( ai, bi ) bi = ai, low guard word of bi = 0
118 * __emulm( ai, bi ) multiply significands, bi = bi * ai
119 * __enormlz(ei) left-justify the significand
120 * __eshdn1( ai ) shift significand and guards down 1 bit
121 * __eshdn8( ai ) shift down 8 bits
122 * __eshdn6( ai ) shift down 16 bits
123 * __eshift( ai, n ) shift ai n bits up (or down if n < 0)
124 * __eshup1( ai ) shift significand and guards up 1 bit
125 * __eshup8( ai ) shift up 8 bits
126 * __eshup6( ai ) shift up 16 bits
127 * __esubm( ai, bi ) subtract significands, bi = bi - ai
128 *
129 *
130 * The result is always normalized and rounded to NI-4 word precision
131 * after each arithmetic operation.
132 *
133 * Exception flags are NOT fully supported.
134 *
135 * Define INFINITY in mconf.h for support of infinity; otherwise a
136 * saturation arithmetic is implemented.
137 *
138 * Define NANS for support of Not-a-Number items; otherwise the
139 * arithmetic will never produce a NaN output, and might be confused
140 * by a NaN input.
141 * If NaN's are supported, the output of ecmp(a,b) is -2 if
142 * either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
143 * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
144 * if in doubt.
145 * Signaling NaN's are NOT supported; they are treated the same
146 * as quiet NaN's.
147 *
148 * Denormals are always supported here where appropriate (e.g., not
149 * for conversion to DEC numbers).
150 */
151
152 #include <stdio.h>
153 #include <stdlib.h>
154 #include <string.h>
155 #include <errno.h>
156 #include <math.h>
157 #include <locale.h>
158 #include <ctype.h>
159
160 #undef alloca
161 #define alloca __builtin_alloca
162
163 /* Don't build non-ANSI _IO_ldtoa. It is not thread safe. */
164 #ifndef USE_LDTOA
165 #define USE_LDTOA 0
166 #endif
167
168
169 /* Number of 16 bit words in external x type format */
170 #define NE 6
171
172 /* Number of 16 bit words in internal format */
173 #define NI (NE+3)
174
175 /* Array offset to exponent */
176 #define E 1
177
178 /* Array offset to high guard word */
179 #define M 2
180
181 /* Number of bits of precision */
182 #define NBITS ((NI-4)*16)
183
184 /* Maximum number of decimal digits in ASCII conversion
185 * = NBITS*log10(2)
186 */
187 #define NDEC (NBITS*8/27)
188
189 /* The exponent of 1.0 */
190 #define EXONE (0x3fff)
191
192
193 #define mtherr(fname, code)
194
195
196 extern long double strtold (const char * __restrict__ s, char ** __restrict__ se);
197 extern int __asctoe64(const char * __restrict__ ss,
198 short unsigned int * __restrict__ y);
199 extern void __emul(const short unsigned int * a,
200 const short unsigned int * b,
201 short unsigned int * c);
202 extern int __ecmp(const short unsigned int * __restrict__ a,
203 const short unsigned int * __restrict__ b);
204 extern int __enormlz(short unsigned int *x);
205 extern int __eshift(short unsigned int *x, int sc);
206 extern void __eaddm(const short unsigned int * __restrict__ x,
207 short unsigned int * __restrict__ y);
208 extern void __esubm(const short unsigned int * __restrict__ x,
209 short unsigned int * __restrict__ y);
210 extern void __emdnorm(short unsigned int *s, int lost, int subflg,
211 int exp, int rcntrl, const int rndprc);
212 extern void __toe64(short unsigned int * __restrict__ a,
213 short unsigned int * __restrict__ b);
214 extern int __edivm(short unsigned int * __restrict__ den,
215 short unsigned int * __restrict__ num);
216 extern int __emulm(const short unsigned int * __restrict__ a,
217 short unsigned int * __restrict__ b);
218 extern void __emovi(const short unsigned int * __restrict__ a,
219 short unsigned int * __restrict__ b);
220 extern void __emovo(const short unsigned int * __restrict__ a,
221 short unsigned int * __restrict__ b);
222
223 #if USE_LDTOA
224
225 extern char * _IO_ldtoa(long double, int, int, int *, int *, char **);
226 extern void _IO_ldtostr(long double *x, char *string, int ndigs,
227 int flags, char fmt);
228
229 extern void __eiremain(short unsigned int * __restrict__ den,
230 short unsigned int *__restrict__ num,
231 short unsigned int *__restrict__ equot);
232 extern void __efloor(short unsigned int *x, short unsigned int *y);
233 extern void __eadd1(const short unsigned int * __restrict__ a,
234 const short unsigned int * __restrict__ b,
235 short unsigned int * __restrict__ c,
236 int subflg);
237 extern void __esub(const short unsigned int *a, const short unsigned int *b,
238 short unsigned int *c);
239 extern void __ediv(const short unsigned int *a, const short unsigned int *b,
240 short unsigned int *c);
241 extern void __e64toe(short unsigned int *pe, short unsigned int *y);
242
243
244 #endif
245
246 static __inline__ int __eisneg(const short unsigned int *x);
247 static __inline__ int __eisinf(const short unsigned int *x);
248 static __inline__ int __eisnan(const short unsigned int *x);
249 static __inline__ int __eiszero(const short unsigned int *a);
250 static __inline__ void __emovz(register const short unsigned int * __restrict__ a,
251 register short unsigned int * __restrict__ b);
252 static __inline__ void __eclear(register short unsigned int *x);
253 static __inline__ void __ecleaz(register short unsigned int *xi);
254 static __inline__ void __ecleazs(register short unsigned int *xi);
255 static __inline__ int __eiisinf(const short unsigned int *x);
256 static __inline__ int __eiisnan(const short unsigned int *x);
257 static __inline__ int __eiiszero(const short unsigned int *x);
258 static __inline__ void __enan_64(short unsigned int *nanptr);
259 static __inline__ void __enan_NBITS (short unsigned int *nanptr);
260 static __inline__ void __enan_NI16 (short unsigned int *nanptr);
261 static __inline__ void __einfin(register short unsigned int *x);
262 static __inline__ void __eneg(short unsigned int *x);
263 static __inline__ void __eshup1(register short unsigned int *x);
264 static __inline__ void __eshup8(register short unsigned int *x);
265 static __inline__ void __eshup6(register short unsigned int *x);
266 static __inline__ void __eshdn1(register short unsigned int *x);
267 static __inline__ void __eshdn8(register short unsigned int *x);
268 static __inline__ void __eshdn6(register short unsigned int *x);
269
270
271
272 /* Intel IEEE, low order words come first:
273 */
274 #define IBMPC 1
275
276 /* Define 1 for ANSI C atan2() function
277 * See atan.c and clog.c.
278 */
279 #define ANSIC 1
280
281 /*define VOLATILE volatile*/
282 #define VOLATILE
283
284 /* For 12-byte long doubles on an i386, pad a 16-bit short 0
285 * to the end of real constants initialized by integer arrays.
286 *
287 * #define XPD 0,
288 *
289 * Otherwise, the type is 10 bytes long and XPD should be
290 * defined blank.
291 *
292 * #define XPD
293 */
294 #define XPD 0,
295 /* #define XPD */
296 #define NANS 1
297
298 /* NaN's require infinity support. */
299 #ifdef NANS
300 #ifndef INFINITY
301 #define INFINITY
302 #endif
303 #endif
304
305 /* This handles 64-bit long ints. */
306 #define LONGBITS (8 * sizeof(long))
307
308
309 #define NTEN 12
310 #define MAXP 4096
311
312 /*
313 ; Clear out entire external format number.
314 ;
315 ; unsigned short x[];
316 ; eclear( x );
317 */
318
__eclear(register short unsigned int * x)319 static __inline__ void __eclear(register short unsigned int *x)
320 {
321 memset(x, 0, NE * sizeof(unsigned short));
322 }
323
324
325 /* Move external format number from a to b.
326 *
327 * emov( a, b );
328 */
329
__emov(register const short unsigned int * __restrict__ a,register short unsigned int * __restrict__ b)330 static __inline__ void __emov(register const short unsigned int * __restrict__ a,
331 register short unsigned int * __restrict__ b)
332 {
333 memcpy(b, a, NE * sizeof(unsigned short));
334 }
335
336
337 /*
338 ; Negate external format number
339 ;
340 ; unsigned short x[NE];
341 ; eneg( x );
342 */
343
__eneg(short unsigned int * x)344 static __inline__ void __eneg(short unsigned int *x)
345 {
346 #ifdef NANS
347 if (__eisnan(x))
348 return;
349 #endif
350 x[NE-1] ^= 0x8000; /* Toggle the sign bit */
351 }
352
353
354 /* Return 1 if external format number is negative,
355 * else return zero.
356 */
__eisneg(const short unsigned int * x)357 static __inline__ int __eisneg(const short unsigned int *x)
358 {
359 #ifdef NANS
360 if (__eisnan(x))
361 return (0);
362 #endif
363 if (x[NE-1] & 0x8000)
364 return (1);
365 else
366 return (0);
367 }
368
369
370 /* Return 1 if external format number has maximum possible exponent,
371 * else return zero.
372 */
__eisinf(const short unsigned int * x)373 static __inline__ int __eisinf(const short unsigned int *x)
374 {
375 if ((x[NE - 1] & 0x7fff) == 0x7fff)
376 {
377 #ifdef NANS
378 if (__eisnan(x))
379 return (0);
380 #endif
381 return (1);
382 }
383 else
384 return (0);
385 }
386
387 /* Check if e-type number is not a number.
388 */
__eisnan(const short unsigned int * x)389 static __inline__ int __eisnan(const short unsigned int *x)
390 {
391 #ifdef NANS
392 int i;
393 /* NaN has maximum __exponent */
394 if ((x[NE - 1] & 0x7fff) == 0x7fff)
395 /* ... and non-zero significand field. */
396 for (i = 0; i < NE - 1; i++)
397 {
398 if (*x++ != 0)
399 return (1);
400 }
401 #endif
402 return (0);
403 }
404
405 /*
406 ; Fill __entire number, including __exponent and significand, with
407 ; largest possible number. These programs implement a saturation
408 ; value that is an ordinary, legal number. A special value
409 ; "infinity" may also be implemented; this would require tests
410 ; for that value and implementation of special rules for arithmetic
411 ; operations involving inifinity.
412 */
413
__einfin(register short unsigned int * x)414 static __inline__ void __einfin(register short unsigned int *x)
415 {
416 register int i;
417 #ifdef INFINITY
418 for (i = 0; i < NE - 1; i++)
419 *x++ = 0;
420 *x |= 32767;
421 #else
422 for (i = 0; i < NE - 1; i++)
423 *x++ = 0xffff;
424 *x |= 32766;
425 *(x - 5) = 0;
426 #endif
427 }
428
429 /* Clear out internal format number.
430 */
431
__ecleaz(register short unsigned int * xi)432 static __inline__ void __ecleaz(register short unsigned int *xi)
433 {
434 memset(xi, 0, NI * sizeof(unsigned short));
435 }
436
437 /* same, but don't touch the sign. */
438
__ecleazs(register short unsigned int * xi)439 static __inline__ void __ecleazs(register short unsigned int *xi)
440 {
441 ++xi;
442 memset(xi, 0, (NI-1) * sizeof(unsigned short));
443 }
444
445 /* Move internal format number from a to b.
446 */
__emovz(register const short unsigned int * __restrict__ a,register short unsigned int * __restrict__ b)447 static __inline__ void __emovz(register const short unsigned int * __restrict__ a,
448 register short unsigned int * __restrict__ b)
449 {
450 memcpy(b, a, (NI-1) * sizeof(unsigned short));
451 b[NI - 1] = 0;
452 }
453
454 /* Return nonzero if internal format number is a NaN.
455 */
456
__eiisnan(const short unsigned int * x)457 static __inline__ int __eiisnan (const short unsigned int *x)
458 {
459 int i;
460
461 if ((x[E] & 0x7fff) == 0x7fff)
462 {
463 for (i = M + 1; i < NI; i++ )
464 {
465 if (x[i] != 0)
466 return (1);
467 }
468 }
469 return (0);
470 }
471
472 /* Return nonzero if external format number is zero. */
473
474 static __inline__ int
__eiszero(const short unsigned int * a)475 __eiszero(const short unsigned int * a)
476 {
477 union {
478 long double ld;
479 unsigned short sh[8];
480 } av;
481 av.ld = 0.0;
482 memcpy (av.sh, a, 12);
483 if (av.ld == 0.0)
484 return (1);
485 return (0);
486 }
487
488 /* Return nonzero if internal format number is zero. */
489
490 static __inline__ int
__eiiszero(const short unsigned int * ai)491 __eiiszero(const short unsigned int * ai)
492 {
493 int i;
494 /* skip the sign word */
495 for (i = 1; i < NI - 1; i++ )
496 {
497 if (ai[i] != 0)
498 return (0);
499 }
500 return (1);
501 }
502
503
504 /* Return nonzero if internal format number is infinite. */
505
506 static __inline__ int
__eiisinf(const unsigned short * x)507 __eiisinf (const unsigned short *x)
508 {
509 #ifdef NANS
510 if (__eiisnan (x))
511 return (0);
512 #endif
513 if ((x[E] & 0x7fff) == 0x7fff)
514 return (1);
515 return (0);
516 }
517
518 /*
519 ; Compare significands of numbers in internal format.
520 ; Guard words are included in the comparison.
521 ;
522 ; unsigned short a[NI], b[NI];
523 ; cmpm( a, b );
524 ;
525 ; for the significands:
526 ; returns +1 if a > b
527 ; 0 if a == b
528 ; -1 if a < b
529 */
__ecmpm(register const short unsigned int * __restrict__ a,register const short unsigned int * __restrict__ b)530 static __inline__ int __ecmpm(register const short unsigned int * __restrict__ a,
531 register const short unsigned int * __restrict__ b)
532 {
533 int i;
534
535 a += M; /* skip up to significand area */
536 b += M;
537 for (i = M; i < NI; i++)
538 {
539 if( *a++ != *b++ )
540 goto difrnt;
541 }
542 return(0);
543
544 difrnt:
545 if ( *(--a) > *(--b) )
546 return (1);
547 else
548 return (-1);
549 }
550
551
552 /*
553 ; Shift significand down by 1 bit
554 */
555
__eshdn1(register short unsigned int * x)556 static __inline__ void __eshdn1(register short unsigned int *x)
557 {
558 register unsigned short bits;
559 int i;
560
561 x += M; /* point to significand area */
562
563 bits = 0;
564 for (i = M; i < NI; i++ )
565 {
566 if (*x & 1)
567 bits |= 1;
568 *x >>= 1;
569 if (bits & 2)
570 *x |= 0x8000;
571 bits <<= 1;
572 ++x;
573 }
574 }
575
576 /*
577 ; Shift significand up by 1 bit
578 */
579
__eshup1(register short unsigned int * x)580 static __inline__ void __eshup1(register short unsigned int *x)
581 {
582 register unsigned short bits;
583 int i;
584
585 x += NI-1;
586 bits = 0;
587
588 for (i = M; i < NI; i++)
589 {
590 if (*x & 0x8000)
591 bits |= 1;
592 *x <<= 1;
593 if (bits & 2)
594 *x |= 1;
595 bits <<= 1;
596 --x;
597 }
598 }
599
600
601 /*
602 ; Shift significand down by 8 bits
603 */
604
__eshdn8(register short unsigned int * x)605 static __inline__ void __eshdn8(register short unsigned int *x)
606 {
607 register unsigned short newbyt, oldbyt;
608 int i;
609
610 x += M;
611 oldbyt = 0;
612 for (i = M; i < NI; i++)
613 {
614 newbyt = *x << 8;
615 *x >>= 8;
616 *x |= oldbyt;
617 oldbyt = newbyt;
618 ++x;
619 }
620 }
621
622 /*
623 ; Shift significand up by 8 bits
624 */
625
__eshup8(register short unsigned int * x)626 static __inline__ void __eshup8(register short unsigned int *x)
627 {
628 int i;
629 register unsigned short newbyt, oldbyt;
630
631 x += NI - 1;
632 oldbyt = 0;
633
634 for (i = M; i < NI; i++)
635 {
636 newbyt = *x >> 8;
637 *x <<= 8;
638 *x |= oldbyt;
639 oldbyt = newbyt;
640 --x;
641 }
642 }
643
644 /*
645 ; Shift significand up by 16 bits
646 */
647
__eshup6(register short unsigned int * x)648 static __inline__ void __eshup6(register short unsigned int *x)
649 {
650 int i;
651 register unsigned short *p;
652
653 p = x + M;
654 x += M + 1;
655
656 for (i = M; i < NI - 1; i++)
657 *p++ = *x++;
658
659 *p = 0;
660 }
661
662 /*
663 ; Shift significand down by 16 bits
664 */
665
__eshdn6(register short unsigned int * x)666 static __inline__ void __eshdn6(register short unsigned int *x)
667 {
668 int i;
669 register unsigned short *p;
670
671 x += NI - 1;
672 p = x + 1;
673
674 for (i = M; i < NI - 1; i++)
675 *(--p) = *(--x);
676
677 *(--p) = 0;
678 }
679
680 /*
681 ; Add significands
682 ; x + y replaces y
683 */
684
__enan_64(unsigned short * nanptr)685 static __inline__ void __enan_64(unsigned short* nanptr)
686 {
687 int i;
688 for (i = 0; i < 3; i++)
689 *nanptr++ = 0;
690 *nanptr++ = 0xc000;
691 *nanptr++ = 0x7fff;
692 *nanptr = 0;
693 return;
694 }
695
__enan_NBITS(unsigned short * nanptr)696 static __inline__ void __enan_NBITS(unsigned short* nanptr)
697 {
698 int i;
699 for (i = 0; i < NE - 2; i++)
700 *nanptr++ = 0;
701 *nanptr++ = 0xc000;
702 *nanptr = 0x7fff;
703 return;
704 }
705
__enan_NI16(unsigned short * nanptr)706 static __inline__ void __enan_NI16(unsigned short* nanptr)
707 {
708 int i;
709 *nanptr++ = 0;
710 *nanptr++ = 0x7fff;
711 *nanptr++ = 0;
712 *nanptr++ = 0xc000;
713 for (i = 4; i < NI; i++)
714 *nanptr++ = 0;
715 return;
716 }
717
718 #endif /* _CEPHES_EMATH_H */
719
720