Blame


1 b2cfc4e2 2003-09-30 devnull /*
2 b2cfc4e2 2003-09-30 devnull * The authors of this software are Rob Pike and Ken Thompson.
3 b2cfc4e2 2003-09-30 devnull * Copyright (c) 2002 by Lucent Technologies.
4 b2cfc4e2 2003-09-30 devnull * Permission to use, copy, modify, and distribute this software for any
5 b2cfc4e2 2003-09-30 devnull * purpose without fee is hereby granted, provided that this entire notice
6 b2cfc4e2 2003-09-30 devnull * is included in all copies of any software which is or includes a copy
7 b2cfc4e2 2003-09-30 devnull * or modification of this software and in all copies of the supporting
8 b2cfc4e2 2003-09-30 devnull * documentation for such software.
9 b2cfc4e2 2003-09-30 devnull * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
10 b2cfc4e2 2003-09-30 devnull * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR LUCENT TECHNOLOGIES MAKE ANY
11 b2cfc4e2 2003-09-30 devnull * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
12 b2cfc4e2 2003-09-30 devnull * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
13 b2cfc4e2 2003-09-30 devnull */
14 b2cfc4e2 2003-09-30 devnull #include <stdlib.h>
15 b2cfc4e2 2003-09-30 devnull #include <math.h>
16 b2cfc4e2 2003-09-30 devnull #include <ctype.h>
17 b2cfc4e2 2003-09-30 devnull #include <stdlib.h>
18 b2cfc4e2 2003-09-30 devnull #include <string.h>
19 b2cfc4e2 2003-09-30 devnull #include <errno.h>
20 b2cfc4e2 2003-09-30 devnull #include "fmt.h"
21 b2cfc4e2 2003-09-30 devnull #include "nan.h"
22 b2cfc4e2 2003-09-30 devnull
23 b2cfc4e2 2003-09-30 devnull #ifndef nelem
24 b2cfc4e2 2003-09-30 devnull #define nelem(x) (sizeof(x)/sizeof *(x))
25 b2cfc4e2 2003-09-30 devnull #endif
26 b2cfc4e2 2003-09-30 devnull #define nil ((void*)0)
27 b2cfc4e2 2003-09-30 devnull #define ulong _fmtulong
28 b2cfc4e2 2003-09-30 devnull typedef unsigned long ulong;
29 b2cfc4e2 2003-09-30 devnull
30 b2cfc4e2 2003-09-30 devnull static ulong
31 b2cfc4e2 2003-09-30 devnull umuldiv(ulong a, ulong b, ulong c)
32 b2cfc4e2 2003-09-30 devnull {
33 b2cfc4e2 2003-09-30 devnull double d;
34 b2cfc4e2 2003-09-30 devnull
35 b2cfc4e2 2003-09-30 devnull d = ((double)a * (double)b) / (double)c;
36 b2cfc4e2 2003-09-30 devnull if(d >= 4294967295.)
37 b2cfc4e2 2003-09-30 devnull d = 4294967295.;
38 b2cfc4e2 2003-09-30 devnull return (ulong)d;
39 b2cfc4e2 2003-09-30 devnull }
40 b2cfc4e2 2003-09-30 devnull
41 b2cfc4e2 2003-09-30 devnull /*
42 b2cfc4e2 2003-09-30 devnull * This routine will convert to arbitrary precision
43 b2cfc4e2 2003-09-30 devnull * floating point entirely in multi-precision fixed.
44 b2cfc4e2 2003-09-30 devnull * The answer is the closest floating point number to
45 b2cfc4e2 2003-09-30 devnull * the given decimal number. Exactly half way are
46 b2cfc4e2 2003-09-30 devnull * rounded ala ieee rules.
47 b2cfc4e2 2003-09-30 devnull * Method is to scale input decimal between .500 and .999...
48 b2cfc4e2 2003-09-30 devnull * with external power of 2, then binary search for the
49 b2cfc4e2 2003-09-30 devnull * closest mantissa to this decimal number.
50 b2cfc4e2 2003-09-30 devnull * Nmant is is the required precision. (53 for ieee dp)
51 b2cfc4e2 2003-09-30 devnull * Nbits is the max number of bits/word. (must be <= 28)
52 b2cfc4e2 2003-09-30 devnull * Prec is calculated - the number of words of fixed mantissa.
53 b2cfc4e2 2003-09-30 devnull */
54 b2cfc4e2 2003-09-30 devnull enum
55 b2cfc4e2 2003-09-30 devnull {
56 b2cfc4e2 2003-09-30 devnull Nbits = 28, /* bits safely represented in a ulong */
57 b2cfc4e2 2003-09-30 devnull Nmant = 53, /* bits of precision required */
58 b2cfc4e2 2003-09-30 devnull Prec = (Nmant+Nbits+1)/Nbits, /* words of Nbits each to represent mantissa */
59 b2cfc4e2 2003-09-30 devnull Sigbit = 1<<(Prec*Nbits-Nmant), /* first significant bit of Prec-th word */
60 b2cfc4e2 2003-09-30 devnull Ndig = 1500,
61 b2cfc4e2 2003-09-30 devnull One = (ulong)(1<<Nbits),
62 b2cfc4e2 2003-09-30 devnull Half = (ulong)(One>>1),
63 b2cfc4e2 2003-09-30 devnull Maxe = 310,
64 b2cfc4e2 2003-09-30 devnull
65 b2cfc4e2 2003-09-30 devnull Fsign = 1<<0, /* found - */
66 b2cfc4e2 2003-09-30 devnull Fesign = 1<<1, /* found e- */
67 b2cfc4e2 2003-09-30 devnull Fdpoint = 1<<2, /* found . */
68 b2cfc4e2 2003-09-30 devnull
69 b2cfc4e2 2003-09-30 devnull S0 = 0, /* _ _S0 +S1 #S2 .S3 */
70 b2cfc4e2 2003-09-30 devnull S1, /* _+ #S2 .S3 */
71 b2cfc4e2 2003-09-30 devnull S2, /* _+# #S2 .S4 eS5 */
72 b2cfc4e2 2003-09-30 devnull S3, /* _+. #S4 */
73 b2cfc4e2 2003-09-30 devnull S4, /* _+#.# #S4 eS5 */
74 b2cfc4e2 2003-09-30 devnull S5, /* _+#.#e +S6 #S7 */
75 b2cfc4e2 2003-09-30 devnull S6, /* _+#.#e+ #S7 */
76 b2cfc4e2 2003-09-30 devnull S7, /* _+#.#e+# #S7 */
77 b2cfc4e2 2003-09-30 devnull };
78 b2cfc4e2 2003-09-30 devnull
79 b2cfc4e2 2003-09-30 devnull static int xcmp(char*, char*);
80 b2cfc4e2 2003-09-30 devnull static int fpcmp(char*, ulong*);
81 b2cfc4e2 2003-09-30 devnull static void frnorm(ulong*);
82 b2cfc4e2 2003-09-30 devnull static void divascii(char*, int*, int*, int*);
83 b2cfc4e2 2003-09-30 devnull static void mulascii(char*, int*, int*, int*);
84 b2cfc4e2 2003-09-30 devnull
85 b2cfc4e2 2003-09-30 devnull typedef struct Tab Tab;
86 b2cfc4e2 2003-09-30 devnull struct Tab
87 b2cfc4e2 2003-09-30 devnull {
88 b2cfc4e2 2003-09-30 devnull int bp;
89 b2cfc4e2 2003-09-30 devnull int siz;
90 b2cfc4e2 2003-09-30 devnull char* cmp;
91 b2cfc4e2 2003-09-30 devnull };
92 b2cfc4e2 2003-09-30 devnull
93 b2cfc4e2 2003-09-30 devnull double
94 b2cfc4e2 2003-09-30 devnull fmtstrtod(const char *as, char **aas)
95 b2cfc4e2 2003-09-30 devnull {
96 b2cfc4e2 2003-09-30 devnull int na, ex, dp, bp, c, i, flag, state;
97 b2cfc4e2 2003-09-30 devnull ulong low[Prec], hig[Prec], mid[Prec];
98 b2cfc4e2 2003-09-30 devnull double d;
99 b2cfc4e2 2003-09-30 devnull char *s, a[Ndig];
100 b2cfc4e2 2003-09-30 devnull
101 b2cfc4e2 2003-09-30 devnull flag = 0; /* Fsign, Fesign, Fdpoint */
102 b2cfc4e2 2003-09-30 devnull na = 0; /* number of digits of a[] */
103 b2cfc4e2 2003-09-30 devnull dp = 0; /* na of decimal point */
104 b2cfc4e2 2003-09-30 devnull ex = 0; /* exonent */
105 b2cfc4e2 2003-09-30 devnull
106 b2cfc4e2 2003-09-30 devnull state = S0;
107 b2cfc4e2 2003-09-30 devnull for(s=(char*)as;; s++) {
108 b2cfc4e2 2003-09-30 devnull c = *s;
109 b2cfc4e2 2003-09-30 devnull if(c >= '0' && c <= '9') {
110 b2cfc4e2 2003-09-30 devnull switch(state) {
111 b2cfc4e2 2003-09-30 devnull case S0:
112 b2cfc4e2 2003-09-30 devnull case S1:
113 b2cfc4e2 2003-09-30 devnull case S2:
114 b2cfc4e2 2003-09-30 devnull state = S2;
115 b2cfc4e2 2003-09-30 devnull break;
116 b2cfc4e2 2003-09-30 devnull case S3:
117 b2cfc4e2 2003-09-30 devnull case S4:
118 b2cfc4e2 2003-09-30 devnull state = S4;
119 b2cfc4e2 2003-09-30 devnull break;
120 b2cfc4e2 2003-09-30 devnull
121 b2cfc4e2 2003-09-30 devnull case S5:
122 b2cfc4e2 2003-09-30 devnull case S6:
123 b2cfc4e2 2003-09-30 devnull case S7:
124 b2cfc4e2 2003-09-30 devnull state = S7;
125 b2cfc4e2 2003-09-30 devnull ex = ex*10 + (c-'0');
126 b2cfc4e2 2003-09-30 devnull continue;
127 b2cfc4e2 2003-09-30 devnull }
128 b2cfc4e2 2003-09-30 devnull if(na == 0 && c == '0') {
129 b2cfc4e2 2003-09-30 devnull dp--;
130 b2cfc4e2 2003-09-30 devnull continue;
131 b2cfc4e2 2003-09-30 devnull }
132 b2cfc4e2 2003-09-30 devnull if(na < Ndig-50)
133 b2cfc4e2 2003-09-30 devnull a[na++] = c;
134 b2cfc4e2 2003-09-30 devnull continue;
135 b2cfc4e2 2003-09-30 devnull }
136 b2cfc4e2 2003-09-30 devnull switch(c) {
137 b2cfc4e2 2003-09-30 devnull case '\t':
138 b2cfc4e2 2003-09-30 devnull case '\n':
139 b2cfc4e2 2003-09-30 devnull case '\v':
140 b2cfc4e2 2003-09-30 devnull case '\f':
141 b2cfc4e2 2003-09-30 devnull case '\r':
142 b2cfc4e2 2003-09-30 devnull case ' ':
143 b2cfc4e2 2003-09-30 devnull if(state == S0)
144 b2cfc4e2 2003-09-30 devnull continue;
145 b2cfc4e2 2003-09-30 devnull break;
146 b2cfc4e2 2003-09-30 devnull case '-':
147 b2cfc4e2 2003-09-30 devnull if(state == S0)
148 b2cfc4e2 2003-09-30 devnull flag |= Fsign;
149 b2cfc4e2 2003-09-30 devnull else
150 b2cfc4e2 2003-09-30 devnull flag |= Fesign;
151 b2cfc4e2 2003-09-30 devnull case '+':
152 b2cfc4e2 2003-09-30 devnull if(state == S0)
153 b2cfc4e2 2003-09-30 devnull state = S1;
154 b2cfc4e2 2003-09-30 devnull else
155 b2cfc4e2 2003-09-30 devnull if(state == S5)
156 b2cfc4e2 2003-09-30 devnull state = S6;
157 b2cfc4e2 2003-09-30 devnull else
158 b2cfc4e2 2003-09-30 devnull break; /* syntax */
159 b2cfc4e2 2003-09-30 devnull continue;
160 b2cfc4e2 2003-09-30 devnull case '.':
161 b2cfc4e2 2003-09-30 devnull flag |= Fdpoint;
162 b2cfc4e2 2003-09-30 devnull dp = na;
163 b2cfc4e2 2003-09-30 devnull if(state == S0 || state == S1) {
164 b2cfc4e2 2003-09-30 devnull state = S3;
165 b2cfc4e2 2003-09-30 devnull continue;
166 b2cfc4e2 2003-09-30 devnull }
167 b2cfc4e2 2003-09-30 devnull if(state == S2) {
168 b2cfc4e2 2003-09-30 devnull state = S4;
169 b2cfc4e2 2003-09-30 devnull continue;
170 b2cfc4e2 2003-09-30 devnull }
171 b2cfc4e2 2003-09-30 devnull break;
172 b2cfc4e2 2003-09-30 devnull case 'e':
173 b2cfc4e2 2003-09-30 devnull case 'E':
174 b2cfc4e2 2003-09-30 devnull if(state == S2 || state == S4) {
175 b2cfc4e2 2003-09-30 devnull state = S5;
176 b2cfc4e2 2003-09-30 devnull continue;
177 b2cfc4e2 2003-09-30 devnull }
178 b2cfc4e2 2003-09-30 devnull break;
179 b2cfc4e2 2003-09-30 devnull }
180 b2cfc4e2 2003-09-30 devnull break;
181 b2cfc4e2 2003-09-30 devnull }
182 b2cfc4e2 2003-09-30 devnull
183 b2cfc4e2 2003-09-30 devnull /*
184 b2cfc4e2 2003-09-30 devnull * clean up return char-pointer
185 b2cfc4e2 2003-09-30 devnull */
186 b2cfc4e2 2003-09-30 devnull switch(state) {
187 b2cfc4e2 2003-09-30 devnull case S0:
188 b2cfc4e2 2003-09-30 devnull if(xcmp(s, "nan") == 0) {
189 b2cfc4e2 2003-09-30 devnull if(aas != nil)
190 b2cfc4e2 2003-09-30 devnull *aas = s+3;
191 b2cfc4e2 2003-09-30 devnull goto retnan;
192 b2cfc4e2 2003-09-30 devnull }
193 b2cfc4e2 2003-09-30 devnull case S1:
194 b2cfc4e2 2003-09-30 devnull if(xcmp(s, "infinity") == 0) {
195 b2cfc4e2 2003-09-30 devnull if(aas != nil)
196 b2cfc4e2 2003-09-30 devnull *aas = s+8;
197 b2cfc4e2 2003-09-30 devnull goto retinf;
198 b2cfc4e2 2003-09-30 devnull }
199 b2cfc4e2 2003-09-30 devnull if(xcmp(s, "inf") == 0) {
200 b2cfc4e2 2003-09-30 devnull if(aas != nil)
201 b2cfc4e2 2003-09-30 devnull *aas = s+3;
202 b2cfc4e2 2003-09-30 devnull goto retinf;
203 b2cfc4e2 2003-09-30 devnull }
204 b2cfc4e2 2003-09-30 devnull case S3:
205 b2cfc4e2 2003-09-30 devnull if(aas != nil)
206 b2cfc4e2 2003-09-30 devnull *aas = (char*)as;
207 b2cfc4e2 2003-09-30 devnull goto ret0; /* no digits found */
208 b2cfc4e2 2003-09-30 devnull case S6:
209 b2cfc4e2 2003-09-30 devnull s--; /* back over +- */
210 b2cfc4e2 2003-09-30 devnull case S5:
211 b2cfc4e2 2003-09-30 devnull s--; /* back over e */
212 b2cfc4e2 2003-09-30 devnull break;
213 b2cfc4e2 2003-09-30 devnull }
214 b2cfc4e2 2003-09-30 devnull if(aas != nil)
215 b2cfc4e2 2003-09-30 devnull *aas = s;
216 b2cfc4e2 2003-09-30 devnull
217 b2cfc4e2 2003-09-30 devnull if(flag & Fdpoint)
218 b2cfc4e2 2003-09-30 devnull while(na > 0 && a[na-1] == '0')
219 b2cfc4e2 2003-09-30 devnull na--;
220 b2cfc4e2 2003-09-30 devnull if(na == 0)
221 b2cfc4e2 2003-09-30 devnull goto ret0; /* zero */
222 b2cfc4e2 2003-09-30 devnull a[na] = 0;
223 b2cfc4e2 2003-09-30 devnull if(!(flag & Fdpoint))
224 b2cfc4e2 2003-09-30 devnull dp = na;
225 b2cfc4e2 2003-09-30 devnull if(flag & Fesign)
226 b2cfc4e2 2003-09-30 devnull ex = -ex;
227 b2cfc4e2 2003-09-30 devnull dp += ex;
228 b2cfc4e2 2003-09-30 devnull if(dp < -Maxe){
229 b2cfc4e2 2003-09-30 devnull errno = ERANGE;
230 b2cfc4e2 2003-09-30 devnull goto ret0; /* underflow by exp */
231 b2cfc4e2 2003-09-30 devnull } else
232 b2cfc4e2 2003-09-30 devnull if(dp > +Maxe)
233 b2cfc4e2 2003-09-30 devnull goto retinf; /* overflow by exp */
234 b2cfc4e2 2003-09-30 devnull
235 b2cfc4e2 2003-09-30 devnull /*
236 b2cfc4e2 2003-09-30 devnull * normalize the decimal ascii number
237 b2cfc4e2 2003-09-30 devnull * to range .[5-9][0-9]* e0
238 b2cfc4e2 2003-09-30 devnull */
239 b2cfc4e2 2003-09-30 devnull bp = 0; /* binary exponent */
240 b2cfc4e2 2003-09-30 devnull while(dp > 0)
241 b2cfc4e2 2003-09-30 devnull divascii(a, &na, &dp, &bp);
242 b2cfc4e2 2003-09-30 devnull while(dp < 0 || a[0] < '5')
243 b2cfc4e2 2003-09-30 devnull mulascii(a, &na, &dp, &bp);
244 b2cfc4e2 2003-09-30 devnull
245 b2cfc4e2 2003-09-30 devnull /* close approx by naive conversion */
246 b2cfc4e2 2003-09-30 devnull mid[0] = 0;
247 b2cfc4e2 2003-09-30 devnull mid[1] = 1;
248 b2cfc4e2 2003-09-30 devnull for(i=0; c=a[i]; i++) {
249 b2cfc4e2 2003-09-30 devnull mid[0] = mid[0]*10 + (c-'0');
250 b2cfc4e2 2003-09-30 devnull mid[1] = mid[1]*10;
251 b2cfc4e2 2003-09-30 devnull if(i >= 8)
252 b2cfc4e2 2003-09-30 devnull break;
253 b2cfc4e2 2003-09-30 devnull }
254 b2cfc4e2 2003-09-30 devnull low[0] = umuldiv(mid[0], One, mid[1]);
255 b2cfc4e2 2003-09-30 devnull hig[0] = umuldiv(mid[0]+1, One, mid[1]);
256 b2cfc4e2 2003-09-30 devnull for(i=1; i<Prec; i++) {
257 b2cfc4e2 2003-09-30 devnull low[i] = 0;
258 b2cfc4e2 2003-09-30 devnull hig[i] = One-1;
259 b2cfc4e2 2003-09-30 devnull }
260 b2cfc4e2 2003-09-30 devnull
261 b2cfc4e2 2003-09-30 devnull /* binary search for closest mantissa */
262 b2cfc4e2 2003-09-30 devnull for(;;) {
263 b2cfc4e2 2003-09-30 devnull /* mid = (hig + low) / 2 */
264 b2cfc4e2 2003-09-30 devnull c = 0;
265 b2cfc4e2 2003-09-30 devnull for(i=0; i<Prec; i++) {
266 b2cfc4e2 2003-09-30 devnull mid[i] = hig[i] + low[i];
267 b2cfc4e2 2003-09-30 devnull if(c)
268 b2cfc4e2 2003-09-30 devnull mid[i] += One;
269 b2cfc4e2 2003-09-30 devnull c = mid[i] & 1;
270 b2cfc4e2 2003-09-30 devnull mid[i] >>= 1;
271 b2cfc4e2 2003-09-30 devnull }
272 b2cfc4e2 2003-09-30 devnull frnorm(mid);
273 b2cfc4e2 2003-09-30 devnull
274 b2cfc4e2 2003-09-30 devnull /* compare */
275 b2cfc4e2 2003-09-30 devnull c = fpcmp(a, mid);
276 b2cfc4e2 2003-09-30 devnull if(c > 0) {
277 b2cfc4e2 2003-09-30 devnull c = 1;
278 b2cfc4e2 2003-09-30 devnull for(i=0; i<Prec; i++)
279 b2cfc4e2 2003-09-30 devnull if(low[i] != mid[i]) {
280 b2cfc4e2 2003-09-30 devnull c = 0;
281 b2cfc4e2 2003-09-30 devnull low[i] = mid[i];
282 b2cfc4e2 2003-09-30 devnull }
283 b2cfc4e2 2003-09-30 devnull if(c)
284 b2cfc4e2 2003-09-30 devnull break; /* between mid and hig */
285 b2cfc4e2 2003-09-30 devnull continue;
286 b2cfc4e2 2003-09-30 devnull }
287 b2cfc4e2 2003-09-30 devnull if(c < 0) {
288 b2cfc4e2 2003-09-30 devnull for(i=0; i<Prec; i++)
289 b2cfc4e2 2003-09-30 devnull hig[i] = mid[i];
290 b2cfc4e2 2003-09-30 devnull continue;
291 b2cfc4e2 2003-09-30 devnull }
292 b2cfc4e2 2003-09-30 devnull
293 b2cfc4e2 2003-09-30 devnull /* only hard part is if even/odd roundings wants to go up */
294 b2cfc4e2 2003-09-30 devnull c = mid[Prec-1] & (Sigbit-1);
295 b2cfc4e2 2003-09-30 devnull if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0)
296 b2cfc4e2 2003-09-30 devnull mid[Prec-1] -= c;
297 b2cfc4e2 2003-09-30 devnull break; /* exactly mid */
298 b2cfc4e2 2003-09-30 devnull }
299 b2cfc4e2 2003-09-30 devnull
300 b2cfc4e2 2003-09-30 devnull /* normal rounding applies */
301 b2cfc4e2 2003-09-30 devnull c = mid[Prec-1] & (Sigbit-1);
302 b2cfc4e2 2003-09-30 devnull mid[Prec-1] -= c;
303 b2cfc4e2 2003-09-30 devnull if(c >= Sigbit/2) {
304 b2cfc4e2 2003-09-30 devnull mid[Prec-1] += Sigbit;
305 b2cfc4e2 2003-09-30 devnull frnorm(mid);
306 b2cfc4e2 2003-09-30 devnull }
307 b2cfc4e2 2003-09-30 devnull goto out;
308 b2cfc4e2 2003-09-30 devnull
309 b2cfc4e2 2003-09-30 devnull ret0:
310 b2cfc4e2 2003-09-30 devnull return 0;
311 b2cfc4e2 2003-09-30 devnull
312 b2cfc4e2 2003-09-30 devnull retnan:
313 b2cfc4e2 2003-09-30 devnull return __NaN();
314 b2cfc4e2 2003-09-30 devnull
315 b2cfc4e2 2003-09-30 devnull retinf:
316 b2cfc4e2 2003-09-30 devnull /*
317 b2cfc4e2 2003-09-30 devnull * Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */
318 b2cfc4e2 2003-09-30 devnull errno = ERANGE;
319 b2cfc4e2 2003-09-30 devnull if(flag & Fsign)
320 b2cfc4e2 2003-09-30 devnull return -HUGE_VAL;
321 b2cfc4e2 2003-09-30 devnull return HUGE_VAL;
322 b2cfc4e2 2003-09-30 devnull
323 b2cfc4e2 2003-09-30 devnull out:
324 b2cfc4e2 2003-09-30 devnull d = 0;
325 b2cfc4e2 2003-09-30 devnull for(i=0; i<Prec; i++)
326 b2cfc4e2 2003-09-30 devnull d = d*One + mid[i];
327 b2cfc4e2 2003-09-30 devnull if(flag & Fsign)
328 b2cfc4e2 2003-09-30 devnull d = -d;
329 b2cfc4e2 2003-09-30 devnull d = ldexp(d, bp - Prec*Nbits);
330 b2cfc4e2 2003-09-30 devnull if(d == 0){ /* underflow */
331 b2cfc4e2 2003-09-30 devnull errno = ERANGE;
332 b2cfc4e2 2003-09-30 devnull }
333 b2cfc4e2 2003-09-30 devnull return d;
334 b2cfc4e2 2003-09-30 devnull }
335 b2cfc4e2 2003-09-30 devnull
336 b2cfc4e2 2003-09-30 devnull static void
337 b2cfc4e2 2003-09-30 devnull frnorm(ulong *f)
338 b2cfc4e2 2003-09-30 devnull {
339 b2cfc4e2 2003-09-30 devnull int i, c;
340 b2cfc4e2 2003-09-30 devnull
341 b2cfc4e2 2003-09-30 devnull c = 0;
342 b2cfc4e2 2003-09-30 devnull for(i=Prec-1; i>0; i--) {
343 b2cfc4e2 2003-09-30 devnull f[i] += c;
344 b2cfc4e2 2003-09-30 devnull c = f[i] >> Nbits;
345 b2cfc4e2 2003-09-30 devnull f[i] &= One-1;
346 b2cfc4e2 2003-09-30 devnull }
347 b2cfc4e2 2003-09-30 devnull f[0] += c;
348 b2cfc4e2 2003-09-30 devnull }
349 b2cfc4e2 2003-09-30 devnull
350 b2cfc4e2 2003-09-30 devnull static int
351 b2cfc4e2 2003-09-30 devnull fpcmp(char *a, ulong* f)
352 b2cfc4e2 2003-09-30 devnull {
353 b2cfc4e2 2003-09-30 devnull ulong tf[Prec];
354 b2cfc4e2 2003-09-30 devnull int i, d, c;
355 b2cfc4e2 2003-09-30 devnull
356 b2cfc4e2 2003-09-30 devnull for(i=0; i<Prec; i++)
357 b2cfc4e2 2003-09-30 devnull tf[i] = f[i];
358 b2cfc4e2 2003-09-30 devnull
359 b2cfc4e2 2003-09-30 devnull for(;;) {
360 b2cfc4e2 2003-09-30 devnull /* tf *= 10 */
361 b2cfc4e2 2003-09-30 devnull for(i=0; i<Prec; i++)
362 b2cfc4e2 2003-09-30 devnull tf[i] = tf[i]*10;
363 b2cfc4e2 2003-09-30 devnull frnorm(tf);
364 b2cfc4e2 2003-09-30 devnull d = (tf[0] >> Nbits) + '0';
365 b2cfc4e2 2003-09-30 devnull tf[0] &= One-1;
366 b2cfc4e2 2003-09-30 devnull
367 b2cfc4e2 2003-09-30 devnull /* compare next digit */
368 b2cfc4e2 2003-09-30 devnull c = *a;
369 b2cfc4e2 2003-09-30 devnull if(c == 0) {
370 b2cfc4e2 2003-09-30 devnull if('0' < d)
371 b2cfc4e2 2003-09-30 devnull return -1;
372 b2cfc4e2 2003-09-30 devnull if(tf[0] != 0)
373 b2cfc4e2 2003-09-30 devnull goto cont;
374 b2cfc4e2 2003-09-30 devnull for(i=1; i<Prec; i++)
375 b2cfc4e2 2003-09-30 devnull if(tf[i] != 0)
376 b2cfc4e2 2003-09-30 devnull goto cont;
377 b2cfc4e2 2003-09-30 devnull return 0;
378 b2cfc4e2 2003-09-30 devnull }
379 b2cfc4e2 2003-09-30 devnull if(c > d)
380 b2cfc4e2 2003-09-30 devnull return +1;
381 b2cfc4e2 2003-09-30 devnull if(c < d)
382 b2cfc4e2 2003-09-30 devnull return -1;
383 b2cfc4e2 2003-09-30 devnull a++;
384 b2cfc4e2 2003-09-30 devnull cont:;
385 b2cfc4e2 2003-09-30 devnull }
386 b2cfc4e2 2003-09-30 devnull return 0;
387 b2cfc4e2 2003-09-30 devnull }
388 b2cfc4e2 2003-09-30 devnull
389 b2cfc4e2 2003-09-30 devnull static void
390 b2cfc4e2 2003-09-30 devnull divby(char *a, int *na, int b)
391 b2cfc4e2 2003-09-30 devnull {
392 b2cfc4e2 2003-09-30 devnull int n, c;
393 b2cfc4e2 2003-09-30 devnull char *p;
394 b2cfc4e2 2003-09-30 devnull
395 b2cfc4e2 2003-09-30 devnull p = a;
396 b2cfc4e2 2003-09-30 devnull n = 0;
397 b2cfc4e2 2003-09-30 devnull while(n>>b == 0) {
398 b2cfc4e2 2003-09-30 devnull c = *a++;
399 b2cfc4e2 2003-09-30 devnull if(c == 0) {
400 b2cfc4e2 2003-09-30 devnull while(n) {
401 b2cfc4e2 2003-09-30 devnull c = n*10;
402 b2cfc4e2 2003-09-30 devnull if(c>>b)
403 b2cfc4e2 2003-09-30 devnull break;
404 b2cfc4e2 2003-09-30 devnull n = c;
405 b2cfc4e2 2003-09-30 devnull }
406 b2cfc4e2 2003-09-30 devnull goto xx;
407 b2cfc4e2 2003-09-30 devnull }
408 b2cfc4e2 2003-09-30 devnull n = n*10 + c-'0';
409 b2cfc4e2 2003-09-30 devnull (*na)--;
410 b2cfc4e2 2003-09-30 devnull }
411 b2cfc4e2 2003-09-30 devnull for(;;) {
412 b2cfc4e2 2003-09-30 devnull c = n>>b;
413 b2cfc4e2 2003-09-30 devnull n -= c<<b;
414 b2cfc4e2 2003-09-30 devnull *p++ = c + '0';
415 b2cfc4e2 2003-09-30 devnull c = *a++;
416 b2cfc4e2 2003-09-30 devnull if(c == 0)
417 b2cfc4e2 2003-09-30 devnull break;
418 b2cfc4e2 2003-09-30 devnull n = n*10 + c-'0';
419 b2cfc4e2 2003-09-30 devnull }
420 b2cfc4e2 2003-09-30 devnull (*na)++;
421 b2cfc4e2 2003-09-30 devnull xx:
422 b2cfc4e2 2003-09-30 devnull while(n) {
423 b2cfc4e2 2003-09-30 devnull n = n*10;
424 b2cfc4e2 2003-09-30 devnull c = n>>b;
425 b2cfc4e2 2003-09-30 devnull n -= c<<b;
426 b2cfc4e2 2003-09-30 devnull *p++ = c + '0';
427 b2cfc4e2 2003-09-30 devnull (*na)++;
428 b2cfc4e2 2003-09-30 devnull }
429 b2cfc4e2 2003-09-30 devnull *p = 0;
430 b2cfc4e2 2003-09-30 devnull }
431 b2cfc4e2 2003-09-30 devnull
432 b2cfc4e2 2003-09-30 devnull static Tab tab1[] =
433 b2cfc4e2 2003-09-30 devnull {
434 b2cfc4e2 2003-09-30 devnull 1, 0, "",
435 b2cfc4e2 2003-09-30 devnull 3, 1, "7",
436 b2cfc4e2 2003-09-30 devnull 6, 2, "63",
437 b2cfc4e2 2003-09-30 devnull 9, 3, "511",
438 b2cfc4e2 2003-09-30 devnull 13, 4, "8191",
439 b2cfc4e2 2003-09-30 devnull 16, 5, "65535",
440 b2cfc4e2 2003-09-30 devnull 19, 6, "524287",
441 b2cfc4e2 2003-09-30 devnull 23, 7, "8388607",
442 b2cfc4e2 2003-09-30 devnull 26, 8, "67108863",
443 b2cfc4e2 2003-09-30 devnull 27, 9, "134217727",
444 b2cfc4e2 2003-09-30 devnull };
445 b2cfc4e2 2003-09-30 devnull
446 b2cfc4e2 2003-09-30 devnull static void
447 b2cfc4e2 2003-09-30 devnull divascii(char *a, int *na, int *dp, int *bp)
448 b2cfc4e2 2003-09-30 devnull {
449 b2cfc4e2 2003-09-30 devnull int b, d;
450 b2cfc4e2 2003-09-30 devnull Tab *t;
451 b2cfc4e2 2003-09-30 devnull
452 b2cfc4e2 2003-09-30 devnull d = *dp;
453 b2cfc4e2 2003-09-30 devnull if(d >= (int)(nelem(tab1)))
454 b2cfc4e2 2003-09-30 devnull d = (int)(nelem(tab1))-1;
455 b2cfc4e2 2003-09-30 devnull t = tab1 + d;
456 b2cfc4e2 2003-09-30 devnull b = t->bp;
457 b2cfc4e2 2003-09-30 devnull if(memcmp(a, t->cmp, t->siz) > 0)
458 b2cfc4e2 2003-09-30 devnull d--;
459 b2cfc4e2 2003-09-30 devnull *dp -= d;
460 b2cfc4e2 2003-09-30 devnull *bp += b;
461 b2cfc4e2 2003-09-30 devnull divby(a, na, b);
462 b2cfc4e2 2003-09-30 devnull }
463 b2cfc4e2 2003-09-30 devnull
464 b2cfc4e2 2003-09-30 devnull static void
465 b2cfc4e2 2003-09-30 devnull mulby(char *a, char *p, char *q, int b)
466 b2cfc4e2 2003-09-30 devnull {
467 b2cfc4e2 2003-09-30 devnull int n, c;
468 b2cfc4e2 2003-09-30 devnull
469 b2cfc4e2 2003-09-30 devnull n = 0;
470 b2cfc4e2 2003-09-30 devnull *p = 0;
471 b2cfc4e2 2003-09-30 devnull for(;;) {
472 b2cfc4e2 2003-09-30 devnull q--;
473 b2cfc4e2 2003-09-30 devnull if(q < a)
474 b2cfc4e2 2003-09-30 devnull break;
475 b2cfc4e2 2003-09-30 devnull c = *q - '0';
476 b2cfc4e2 2003-09-30 devnull c = (c<<b) + n;
477 b2cfc4e2 2003-09-30 devnull n = c/10;
478 b2cfc4e2 2003-09-30 devnull c -= n*10;
479 b2cfc4e2 2003-09-30 devnull p--;
480 b2cfc4e2 2003-09-30 devnull *p = c + '0';
481 b2cfc4e2 2003-09-30 devnull }
482 b2cfc4e2 2003-09-30 devnull while(n) {
483 b2cfc4e2 2003-09-30 devnull c = n;
484 b2cfc4e2 2003-09-30 devnull n = c/10;
485 b2cfc4e2 2003-09-30 devnull c -= n*10;
486 b2cfc4e2 2003-09-30 devnull p--;
487 b2cfc4e2 2003-09-30 devnull *p = c + '0';
488 b2cfc4e2 2003-09-30 devnull }
489 b2cfc4e2 2003-09-30 devnull }
490 b2cfc4e2 2003-09-30 devnull
491 b2cfc4e2 2003-09-30 devnull static Tab tab2[] =
492 b2cfc4e2 2003-09-30 devnull {
493 b2cfc4e2 2003-09-30 devnull 1, 1, "", /* dp = 0-0 */
494 b2cfc4e2 2003-09-30 devnull 3, 3, "125",
495 b2cfc4e2 2003-09-30 devnull 6, 5, "15625",
496 b2cfc4e2 2003-09-30 devnull 9, 7, "1953125",
497 b2cfc4e2 2003-09-30 devnull 13, 10, "1220703125",
498 b2cfc4e2 2003-09-30 devnull 16, 12, "152587890625",
499 b2cfc4e2 2003-09-30 devnull 19, 14, "19073486328125",
500 b2cfc4e2 2003-09-30 devnull 23, 17, "11920928955078125",
501 b2cfc4e2 2003-09-30 devnull 26, 19, "1490116119384765625",
502 b2cfc4e2 2003-09-30 devnull 27, 19, "7450580596923828125", /* dp 8-9 */
503 b2cfc4e2 2003-09-30 devnull };
504 b2cfc4e2 2003-09-30 devnull
505 b2cfc4e2 2003-09-30 devnull static void
506 b2cfc4e2 2003-09-30 devnull mulascii(char *a, int *na, int *dp, int *bp)
507 b2cfc4e2 2003-09-30 devnull {
508 b2cfc4e2 2003-09-30 devnull char *p;
509 b2cfc4e2 2003-09-30 devnull int d, b;
510 b2cfc4e2 2003-09-30 devnull Tab *t;
511 b2cfc4e2 2003-09-30 devnull
512 b2cfc4e2 2003-09-30 devnull d = -*dp;
513 b2cfc4e2 2003-09-30 devnull if(d >= (int)(nelem(tab2)))
514 b2cfc4e2 2003-09-30 devnull d = (int)(nelem(tab2))-1;
515 b2cfc4e2 2003-09-30 devnull t = tab2 + d;
516 b2cfc4e2 2003-09-30 devnull b = t->bp;
517 b2cfc4e2 2003-09-30 devnull if(memcmp(a, t->cmp, t->siz) < 0)
518 b2cfc4e2 2003-09-30 devnull d--;
519 b2cfc4e2 2003-09-30 devnull p = a + *na;
520 b2cfc4e2 2003-09-30 devnull *bp -= b;
521 b2cfc4e2 2003-09-30 devnull *dp += d;
522 b2cfc4e2 2003-09-30 devnull *na += d;
523 b2cfc4e2 2003-09-30 devnull mulby(a, p+d, p, b);
524 b2cfc4e2 2003-09-30 devnull }
525 b2cfc4e2 2003-09-30 devnull
526 b2cfc4e2 2003-09-30 devnull static int
527 b2cfc4e2 2003-09-30 devnull xcmp(char *a, char *b)
528 b2cfc4e2 2003-09-30 devnull {
529 b2cfc4e2 2003-09-30 devnull int c1, c2;
530 b2cfc4e2 2003-09-30 devnull
531 b2cfc4e2 2003-09-30 devnull while(c1 = *b++) {
532 b2cfc4e2 2003-09-30 devnull c2 = *a++;
533 b2cfc4e2 2003-09-30 devnull if(isupper(c2))
534 b2cfc4e2 2003-09-30 devnull c2 = tolower(c2);
535 b2cfc4e2 2003-09-30 devnull if(c1 != c2)
536 b2cfc4e2 2003-09-30 devnull return 1;
537 b2cfc4e2 2003-09-30 devnull }
538 b2cfc4e2 2003-09-30 devnull return 0;
539 b2cfc4e2 2003-09-30 devnull }