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
7 #include "math/cephes_emath.h"
8
9 #if NE == 10
10 /* 1.0E0 */
11 static const unsigned short __eone[NE] = {
12 0x0000, 0x0000, 0x0000, 0x0000,
13 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,
14 };
15 #else
16 static const unsigned short __eone[NE] = {
17 0, 0000000,0000000,0000000,0100000,0x3fff,
18 };
19 #endif
20
21 #if NE == 10
22 static const unsigned short __etens[NTEN + 1][NE] =
23 {
24 {0x6576, 0x4a92, 0x804a, 0x153f,
25 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
26 {0x6a32, 0xce52, 0x329a, 0x28ce,
27 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
28 {0x526c, 0x50ce, 0xf18b, 0x3d28,
29 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
30 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
31 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
32 {0x851e, 0xeab7, 0x98fe, 0x901b,
33 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
34 {0x0235, 0x0137, 0x36b1, 0x336c,
35 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
36 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
37 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
38 {0x0000, 0x0000, 0x0000, 0x0000,
39 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
40 {0x0000, 0x0000, 0x0000, 0x0000,
41 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
42 {0x0000, 0x0000, 0x0000, 0x0000,
43 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
44 {0x0000, 0x0000, 0x0000, 0x0000,
45 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
46 {0x0000, 0x0000, 0x0000, 0x0000,
47 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
48 {0x0000, 0x0000, 0x0000, 0x0000,
49 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
50 };
51 #else
52 static const unsigned short __etens[NTEN+1][NE] = {
53 {0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */
54 {0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */
55 {0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,},
56 {0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,},
57 {0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,},
58 {0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,},
59 {0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,},
60 {0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,},
61 {0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,},
62 {0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,},
63 {0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,},
64 {0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,},
65 {0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */
66 };
67 #endif
68
__asctoe64(const char * __restrict__ ss,short unsigned int * __restrict__ y)69 int __asctoe64(const char * __restrict__ ss, short unsigned int * __restrict__ y)
70 {
71 unsigned short yy[NI], xt[NI], tt[NI];
72 int esign, decflg, nexp, expo, lost;
73 int k, c;
74 int valid_lead_string = 0;
75 int have_non_zero_mant = 0;
76 int prec = 0;
77 /* int trail = 0; */
78 int lexp;
79 unsigned short nsign = 0;
80 const unsigned short *p;
81 char *sp, *lstr;
82 char *s;
83
84 const char dec_sym = *(localeconv ()->decimal_point);
85
86 int lenldstr = 0;
87
88 /* Copy the input string. */
89 c = strlen (ss) + 2;
90 lstr = (char *) alloca (c);
91 s = (char *) ss;
92 while( isspace ((int)(unsigned char)*s)) /* skip leading spaces */
93 {
94 ++s;
95 ++lenldstr;
96 }
97 sp = lstr;
98 for (k = 0; k < c; k++)
99 {
100 if ((*sp++ = *s++) == '\0')
101 break;
102 }
103 *sp = '\0';
104 s = lstr;
105
106 if (*s == '-')
107 {
108 nsign = 0xffff;
109 ++s;
110 }
111 else if (*s == '+')
112 {
113 ++s;
114 }
115
116 if (_strnicmp("INF", s , 3) == 0)
117 {
118 valid_lead_string = 1;
119 s += 3;
120 if ( _strnicmp ("INITY", s, 5) == 0)
121 s += 5;
122 __ecleaz(yy);
123 yy[E] = 0x7fff; /* infinity */
124 goto aexit;
125 }
126 else if(_strnicmp ("NAN", s, 3) == 0)
127 {
128 valid_lead_string = 1;
129 s += 3;
130 __enan_NI16( yy );
131 goto aexit;
132 }
133
134 /* FIXME: Handle case of strtold ("NAN(n_char_seq)",endptr) */
135
136 /* Now get some digits. */
137 lost = 0;
138 decflg = 0;
139 nexp = 0;
140 expo = 0;
141 __ecleaz( yy );
142
143 /* Ignore leading zeros */
144 while (*s == '0')
145 {
146 valid_lead_string = 1;
147 s++;
148 }
149
150 nxtcom:
151
152 k = *s - '0';
153 if ((k >= 0) && (k <= 9))
154 {
155 #if 0
156 /* The use of a special char as a flag for trailing zeroes causes problems when input
157 actually contains the char */
158 /* Identify and strip trailing zeros after the decimal point. */
159 if ((trail == 0) && (decflg != 0))
160 {
161 sp = s;
162 while ((*sp >= '0') && (*sp <= '9'))
163 ++sp;
164 --sp;
165 while (*sp == '0')
166 {
167 *sp-- = (char)-1;
168 trail++;
169 }
170 if( *s == (char)-1 )
171 goto donchr;
172 }
173 #endif
174
175 /* If enough digits were given to more than fill up the yy register,
176 * continuing until overflow into the high guard word yy[2]
177 * guarantees that there will be a roundoff bit at the top
178 * of the low guard word after normalization.
179 */
180 if (yy[2] == 0)
181 {
182 if( decflg )
183 nexp += 1; /* count digits after decimal point */
184 __eshup1( yy ); /* multiply current number by 10 */
185 __emovz( yy, xt );
186 __eshup1( xt );
187 __eshup1( xt );
188 __eaddm( xt, yy );
189 __ecleaz( xt );
190 xt[NI-2] = (unsigned short )k;
191 __eaddm( xt, yy );
192 }
193 else
194 {
195 /* Mark any lost non-zero digit. */
196 lost |= k;
197 /* Count lost digits before the decimal point. */
198 if (decflg == 0)
199 nexp -= 1;
200 }
201 have_non_zero_mant |= k;
202 prec ++;
203 /* goto donchr; */
204 }
205 else if (*s == dec_sym)
206 {
207 if( decflg )
208 goto daldone;
209 ++decflg;
210 }
211 else if ((*s == 'E') || (*s == 'e') )
212 {
213 if (prec || valid_lead_string)
214 goto expnt;
215 else
216 goto daldone;
217 }
218 #if 0
219 else if (*s == (char)-1)
220 goto donchr;
221 #endif
222 else /* an invalid char */
223 goto daldone;
224
225 /* donchr: */
226 ++s;
227 goto nxtcom;
228
229 /* Exponent interpretation */
230 expnt:
231
232 esign = 1;
233 expo = 0;
234 /* Save position in case we need to fall back. */
235 sp = s;
236 ++s;
237 /* check for + or - */
238 if (*s == '-')
239 {
240 esign = -1;
241 ++s;
242 }
243 if (*s == '+')
244 ++s;
245
246 /* Check for valid exponent. */
247 if (!(*s >= '0' && *s <= '9'))
248 {
249 s = sp;
250 goto daldone;
251 }
252
253 while ((*s >= '0') && (*s <= '9'))
254 {
255 /* Stop modifying exp if we are going to overflow anyway,
256 but keep parsing the string. */
257 if (expo < 4978)
258 {
259 expo *= 10;
260 expo += *s - '0';
261 }
262 s++;
263 }
264
265 if (esign < 0)
266 expo = -expo;
267
268 if (expo > 4977) /* maybe overflow */
269 {
270 __ecleaz(yy);
271 if (have_non_zero_mant)
272 yy[E] = 0x7fff;
273 goto aexit;
274 }
275 else if (expo < -4977) /* underflow */
276 {
277 __ecleaz(yy);
278 goto aexit;
279 }
280
281 daldone:
282
283 nexp = expo - nexp;
284
285 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
286 while ((nexp > 0) && (yy[2] == 0))
287 {
288 __emovz( yy, xt );
289 __eshup1( xt );
290 __eshup1( xt );
291 __eaddm( yy, xt );
292 __eshup1( xt );
293 if (xt[2] != 0)
294 break;
295 nexp -= 1;
296 __emovz( xt, yy );
297 }
298 if ((k = __enormlz(yy)) > NBITS)
299 {
300 __ecleaz(yy);
301 goto aexit;
302 }
303 lexp = (EXONE - 1 + NBITS) - k;
304 __emdnorm( yy, lost, 0, lexp, 64, NBITS );
305 /* convert to external format */
306
307 /* Multiply by 10**nexp. If precision is 64 bits,
308 * the maximum relative error incurred in forming 10**n
309 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
310 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
311 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
312 */
313 lexp = yy[E];
314 if (nexp == 0)
315 {
316 k = 0;
317 goto expdon;
318 }
319 esign = 1;
320 if (nexp < 0)
321 {
322 nexp = -nexp;
323 esign = -1;
324 if (nexp > 4096)
325 { /* Punt. Can't handle this without 2 divides. */
326 __emovi( __etens[0], tt );
327 lexp -= tt[E];
328 k = __edivm( tt, yy );
329 lexp += EXONE;
330 nexp -= 4096;
331 }
332 }
333 p = &__etens[NTEN][0];
334 __emov( __eone, xt );
335 expo = 1;
336 do
337 {
338 if (expo & nexp)
339 __emul( p, xt, xt );
340 p -= NE;
341 expo = expo + expo;
342 }
343 while (expo <= MAXP);
344
345 __emovi( xt, tt );
346 if (esign < 0)
347 {
348 lexp -= tt[E];
349 k = __edivm( tt, yy );
350 lexp += EXONE;
351 }
352 else
353 {
354 lexp += tt[E];
355 k = __emulm( tt, yy );
356 lexp -= EXONE - 1;
357 }
358
359 expdon:
360
361 /* Round and convert directly to the destination type */
362
363 __emdnorm( yy, k, 0, lexp, 64, 64 );
364
365 aexit:
366
367 yy[0] = nsign;
368
369 __toe64( yy, y );
370
371 /* Check for overflow, undeflow */
372 if (have_non_zero_mant &&
373 (*((long double*) y) == 0.0L || isinf (*((long double*) y))))
374 errno = ERANGE;
375
376 if (prec || valid_lead_string)
377 return (lenldstr + (s - lstr));
378
379 return 0;
380 }
381
382
strtold(const char * __restrict__ s,char ** __restrict__ se)383 long double strtold (const char * __restrict__ s, char ** __restrict__ se)
384 {
385 int lenldstr;
386 union
387 {
388 unsigned short int us[6];
389 long double ld;
390 } xx = {{0}};
391
392 lenldstr = __asctoe64( s, xx.us);
393 if (se)
394 *se = (char*)s + lenldstr;
395
396 return xx.ld;
397 }
398
399