xref: /plan9/sys/src/ape/lib/ap/stdio/fconv.h (revision 7dd7cddf99dd7472612f1413b4da293630e6b1bc)
1 /****************************************************************
2 
3   The author of this software (_dtoa, strtod) is David M. Gay.
4   Please send bug reports to
5 	David M. Gay
6 	Bell Laboratories, Room 2C-463
7 	600 Mountain Avenue
8 	Murray Hill, NJ 07974-2070
9 	U.S.A.
10 	dmg@research.bell-labs.com
11  */
12 #include <stdlib.h>
13 #include <string.h>
14 #define _RESEARCH_SOURCE
15 #include <errno.h>
16 #include <float.h>
17 #include <math.h>
18 #define CONST const
19 
20 /*
21  * #define IEEE_8087 for IEEE-arithmetic machines where the least
22  *	significant byte has the lowest address.
23  * #define IEEE_MC68k for IEEE-arithmetic machines where the most
24  *	significant byte has the lowest address.
25  * #define Sudden_Underflow for IEEE-format machines without gradual
26  *	underflow (i.e., that flush to zero on underflow).
27  * #define IBM for IBM mainframe-style floating-point arithmetic.
28  * #define VAX for VAX-style floating-point arithmetic.
29  * #define Unsigned_Shifts if >> does treats its left operand as unsigned.
30  * #define No_leftright to omit left-right logic in fast floating-point
31  *	computation of dtoa.
32  * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
33  * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
34  *	that use extended-precision instructions to compute rounded
35  *	products and quotients) with IBM.
36  * #define ROUND_BIASED for IEEE-format with biased rounding.
37  * #define Inaccurate_Divide for IEEE-format with correctly rounded
38  *	products but inaccurate quotients, e.g., for Intel i860.
39  * #define Just_16 to store 16 bits per 32-bit long when doing high-precision
40  *	integer arithmetic.  Whether this speeds things up or slows things
41  *	down depends on the machine and the number being converted.
42  */
43 
44 #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
45 Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
46 #endif
47 
48 typedef union {
49 	double d;
50 	unsigned long ul[2];
51 } Dul;
52 
53 #ifdef IEEE_8087
54 #define word0(x) ((x).ul[1])
55 #define word1(x) ((x).ul[0])
56 #else
57 #define word0(x) ((x).ul[0])
58 #define word1(x) ((x).ul[1])
59 #endif
60 
61 #ifdef Unsigned_Shifts
62 #define Sign_Extend(a,b) if (b < 0) a |= 0xffff0000;
63 #else
64 #define Sign_Extend(a,b) /*no-op*/
65 #endif
66 
67 /* The following definition of Storeinc is appropriate for MIPS processors.
68  * An alternative that might be better on some machines is
69  * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
70  */
71 #if defined(IEEE_8087) + defined(VAX)
72 #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
73 ((unsigned short *)a)[0] = (unsigned short)c, a++)
74 #else
75 #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
76 ((unsigned short *)a)[1] = (unsigned short)c, a++)
77 #endif
78 
79 /* #define P DBL_MANT_DIG */
80 /* Ten_pmax = floor(P*log(2)/log(5)) */
81 /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
82 /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
83 /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
84 
85 #if defined(IEEE_8087) + defined(IEEE_MC68k)
86 #define Exp_shift  20
87 #define Exp_shift1 20
88 #define Exp_msk1    0x100000
89 #define Exp_msk11   0x100000
90 #define Exp_mask  0x7ff00000
91 #define P 53
92 #define Bias 1023
93 #define IEEE_Arith
94 #define Emin (-1022)
95 #define Exp_1  0x3ff00000
96 #define Exp_11 0x3ff00000
97 #define Ebits 11
98 #define Frac_mask  0xfffff
99 #define Frac_mask1 0xfffff
100 #define Ten_pmax 22
101 #define Bletch 0x10
102 #define Bndry_mask  0xfffff
103 #define Bndry_mask1 0xfffff
104 #define LSB 1
105 #define Sign_bit 0x80000000
106 #define Log2P 1
107 #define Tiny0 0
108 #define Tiny1 1
109 #define Quick_max 14
110 #define Int_max 14
111 #define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
112 #else
113 #undef  Sudden_Underflow
114 #define Sudden_Underflow
115 #ifdef IBM
116 #define Exp_shift  24
117 #define Exp_shift1 24
118 #define Exp_msk1   0x1000000
119 #define Exp_msk11  0x1000000
120 #define Exp_mask  0x7f000000
121 #define P 14
122 #define Bias 65
123 #define Exp_1  0x41000000
124 #define Exp_11 0x41000000
125 #define Ebits 8	/* exponent has 7 bits, but 8 is the right value in b2d */
126 #define Frac_mask  0xffffff
127 #define Frac_mask1 0xffffff
128 #define Bletch 4
129 #define Ten_pmax 22
130 #define Bndry_mask  0xefffff
131 #define Bndry_mask1 0xffffff
132 #define LSB 1
133 #define Sign_bit 0x80000000
134 #define Log2P 4
135 #define Tiny0 0x100000
136 #define Tiny1 0
137 #define Quick_max 14
138 #define Int_max 15
139 #else /* VAX */
140 #define Exp_shift  23
141 #define Exp_shift1 7
142 #define Exp_msk1    0x80
143 #define Exp_msk11   0x800000
144 #define Exp_mask  0x7f80
145 #define P 56
146 #define Bias 129
147 #define Exp_1  0x40800000
148 #define Exp_11 0x4080
149 #define Ebits 8
150 #define Frac_mask  0x7fffff
151 #define Frac_mask1 0xffff007f
152 #define Ten_pmax 24
153 #define Bletch 2
154 #define Bndry_mask  0xffff007f
155 #define Bndry_mask1 0xffff007f
156 #define LSB 0x10000
157 #define Sign_bit 0x8000
158 #define Log2P 1
159 #define Tiny0 0x80
160 #define Tiny1 0
161 #define Quick_max 15
162 #define Int_max 15
163 #endif
164 #endif
165 
166 #ifndef IEEE_Arith
167 #define ROUND_BIASED
168 #endif
169 
170 #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
171 #define Big1 0xffffffff
172 
173 #ifndef Just_16
174 /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
175  * This makes some inner loops simpler and sometimes saves work
176  * during multiplications, but it often seems to make things slightly
177  * slower.  Hence the default is now to store 32 bits per long.
178  */
179 #ifndef Pack_32
180 #define Pack_32
181 #endif
182 #endif
183 
184 #define Kmax 15
185 
186  struct
187 Bigint {
188 	struct Bigint *next;
189 	int k, maxwds, sign, wds;
190 	unsigned long x[1];
191 	};
192 
193  typedef struct Bigint Bigint;
194 
195 /* This routines shouldn't be visible externally */
196 extern Bigint	*_Balloc(int);
197 extern void	_Bfree(Bigint *);
198 extern Bigint	*_multadd(Bigint *, int, int);
199 extern int	_hi0bits(unsigned long);
200 extern Bigint	*_mult(Bigint *, Bigint *);
201 extern Bigint	*_pow5mult(Bigint *, int);
202 extern Bigint	*_lshift(Bigint *, int);
203 extern int	_cmp(Bigint *, Bigint *);
204 extern Bigint	*_diff(Bigint *, Bigint *);
205 extern Bigint	*_d2b(double, int *, int *);
206 extern Bigint	*_i2b(int);
207 
208 extern double	_tens[], _bigtens[], _tinytens[];
209 
210 #ifdef IEEE_Arith
211 #define n_bigtens 5
212 #else
213 #ifdef IBM
214 #define n_bigtens 3
215 #else
216 #define n_bigtens 2
217 #endif
218 #endif
219 
220 #define Balloc(x) _Balloc(x)
221 #define Bfree(x) _Bfree(x)
222 #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
223 y->wds*sizeof(long) + 2*sizeof(int))
224 #define multadd(x,y,z) _multadd(x,y,z)
225 #define hi0bits(x) _hi0bits(x)
226 #define i2b(x) _i2b(x)
227 #define mult(x,y) _mult(x,y)
228 #define pow5mult(x,y) _pow5mult(x,y)
229 #define lshift(x,y) _lshift(x,y)
230 #define cmp(x,y) _cmp(x,y)
231 #define diff(x,y) _diff(x,y)
232 #define d2b(x,y,z) _d2b(x,y,z)
233 
234 #define tens _tens
235 #define bigtens _bigtens
236 #define tinytens _tinytens
237