1 /*
2  * Copyright (c) 1999
3  * Silicon Graphics Computer Systems, Inc.
4  *
5  * Copyright (c) 1999
6  * Boris Fomitchev
7  *
8  * This material is provided "as is", with absolutely no warranty expressed
9  * or implied. Any use is at your own risk.
10  *
11  * Permission to use or copy this software for any purpose is hereby granted
12  * without fee, provided the above notices are retained on all copies.
13  * Permission to modify the code and to distribute modified code is granted,
14  * provided the above notices are retained, and a notice that the code was
15  * modified is included with the above copyright notice.
16  *
17  */
18 
19 #include "stlport_prefix.h"
20 
21 #include <limits>
22 #include <locale>
23 #include <istream>
24 
25 #if (defined (__GNUC__) && !defined (__sun) && !defined (__hpux)) || \
26     defined (__DMC__)
27 #  include <stdint.h>
28 #endif
29 
30 #if defined (__linux__) || defined (__MINGW32__) || defined (__CYGWIN__) || \
31     defined (__BORLANDC__) || defined (__DMC__) || defined (__HP_aCC)
32 
33 #  if defined (__BORLANDC__)
34 typedef unsigned int uint32_t;
35 typedef unsigned __int64 uint64_t;
36 #  endif
37 
38 union _ll {
39   uint64_t i64;
40   struct {
41 #  if defined (_STLP_BIG_ENDIAN)
42     uint32_t hi;
43     uint32_t lo;
44 #  elif defined (_STLP_LITTLE_ENDIAN)
45     uint32_t lo;
46     uint32_t hi;
47 #  else
48 #    error Unknown endianess
49 #  endif
50   } i32;
51 };
52 
53 #  if defined (__linux__)
54 #    include <ieee754.h>
55 #  else
56 union ieee854_long_double {
57   long double d;
58 
59   /* This is the IEEE 854 double-extended-precision format.  */
60   struct {
61     unsigned int mantissa1:32;
62     unsigned int mantissa0:32;
63     unsigned int exponent:15;
64     unsigned int negative:1;
65     unsigned int empty:16;
66   } ieee;
67 };
68 
69 #    define IEEE854_LONG_DOUBLE_BIAS 0x3fff
70 #  endif
71 #endif
72 
73 _STLP_BEGIN_NAMESPACE
74 _STLP_MOVE_TO_PRIV_NAMESPACE
75 
76 //----------------------------------------------------------------------
77 // num_get
78 
79 // Helper functions for _M_do_get_float.
80 
81 #if !defined (_STLP_NO_WCHAR_T)
82 void  _STLP_CALL
_Initialize_get_float(const ctype<wchar_t> & ct,wchar_t & Plus,wchar_t & Minus,wchar_t & pow_e,wchar_t & pow_E,wchar_t * digits)83 _Initialize_get_float( const ctype<wchar_t>& ct,
84                        wchar_t& Plus, wchar_t& Minus,
85                        wchar_t& pow_e, wchar_t& pow_E,
86                        wchar_t* digits) {
87   char ndigits[11] = "0123456789";
88   Plus  = ct.widen('+');
89   Minus = ct.widen('-');
90   pow_e = ct.widen('e');
91   pow_E = ct.widen('E');
92   ct.widen(ndigits + 0, ndigits + 10, digits);
93 }
94 #endif /* WCHAR_T */
95 
96 /*
97  * __string_to_double is just lifted from atof, the difference being
98  * that we just use '.' for the decimal point, rather than let it
99  * be taken from the current C locale, which of course is not accessible
100  * to us.
101  */
102 #if defined (_STLP_MSVC) || defined (__BORLANDC__) || defined (__ICL)
103 typedef unsigned long uint32;
104 typedef unsigned __int64 uint64;
105 #  define ULL(x) x##Ui64
106 #elif defined (__unix) || defined (__MINGW32__) || \
107       (defined (__DMC__) && (__LONGLONG)) || defined (__WATCOMC__)
108 typedef uint32_t uint32;
109 typedef uint64_t uint64;
110 #  define ULL(x) x##ULL
111 #else
112 #  error There should be some unsigned 64-bit integer on the system!
113 #endif
114 
115 // Multiplication of two 64-bit integers, giving a 128-bit result.
116 // Taken from Algorithm M in Knuth section 4.3.1, with the loop
117 // hand-unrolled.
_Stl_mult64(const uint64 u,const uint64 v,uint64 & high,uint64 & low)118 static void _Stl_mult64(const uint64 u, const uint64 v,
119                         uint64& high, uint64& low) {
120   const uint64 low_mask = ULL(0xffffffff);
121   const uint64 u0 = u & low_mask;
122   const uint64 u1 = u >> 32;
123   const uint64 v0 = v & low_mask;
124   const uint64 v1 = v >> 32;
125 
126   uint64 t = u0 * v0;
127   low = t & low_mask;
128 
129   t = u1 * v0 + (t >> 32);
130   uint64 w1 = t & low_mask;
131   uint64 w2 = t >> 32;
132 
133   uint64 x = u0 * v1 + w1;
134   low += (x & low_mask) << 32;
135   high = u1 * v1 + w2 + (x >> 32);
136 }
137 
138 #ifndef __linux__
139 
140 #  define bit11 ULL(0x7ff)
141 #  define exponent_mask (bit11 << 52)
142 
143 #  if !defined (__GNUC__) || (__GNUC__ != 3) || (__GNUC_MINOR__ != 4) || \
144       (!defined (__CYGWIN__) && !defined (__MINGW32__))
145 //Generate bad code when compiled with -O2 option.
146 inline
147 #  endif
_Stl_set_exponent(uint64 & val,uint64 exp)148 void _Stl_set_exponent(uint64 &val, uint64 exp)
149 { val = (val & ~exponent_mask) | ((exp & bit11) << 52); }
150 
151 #endif // __linux__
152 
153 /* Power of ten fractions for tenscale*/
154 /* The constants are factored so that at most two constants
155  * and two multiplies are needed. Furthermore, one of the constants
156  * is represented exactly - 10**n where 1<= n <= 27.
157  */
158 
159 static const uint64 _Stl_tenpow[80] = {
160 ULL(0xa000000000000000), /* _Stl_tenpow[0]=(10**1)/(2**4) */
161 ULL(0xc800000000000000), /* _Stl_tenpow[1]=(10**2)/(2**7) */
162 ULL(0xfa00000000000000), /* _Stl_tenpow[2]=(10**3)/(2**10) */
163 ULL(0x9c40000000000000), /* _Stl_tenpow[3]=(10**4)/(2**14) */
164 ULL(0xc350000000000000), /* _Stl_tenpow[4]=(10**5)/(2**17) */
165 ULL(0xf424000000000000), /* _Stl_tenpow[5]=(10**6)/(2**20) */
166 ULL(0x9896800000000000), /* _Stl_tenpow[6]=(10**7)/(2**24) */
167 ULL(0xbebc200000000000), /* _Stl_tenpow[7]=(10**8)/(2**27) */
168 ULL(0xee6b280000000000), /* _Stl_tenpow[8]=(10**9)/(2**30) */
169 ULL(0x9502f90000000000), /* _Stl_tenpow[9]=(10**10)/(2**34) */
170 ULL(0xba43b74000000000), /* _Stl_tenpow[10]=(10**11)/(2**37) */
171 ULL(0xe8d4a51000000000), /* _Stl_tenpow[11]=(10**12)/(2**40) */
172 ULL(0x9184e72a00000000), /* _Stl_tenpow[12]=(10**13)/(2**44) */
173 ULL(0xb5e620f480000000), /* _Stl_tenpow[13]=(10**14)/(2**47) */
174 ULL(0xe35fa931a0000000), /* _Stl_tenpow[14]=(10**15)/(2**50) */
175 ULL(0x8e1bc9bf04000000), /* _Stl_tenpow[15]=(10**16)/(2**54) */
176 ULL(0xb1a2bc2ec5000000), /* _Stl_tenpow[16]=(10**17)/(2**57) */
177 ULL(0xde0b6b3a76400000), /* _Stl_tenpow[17]=(10**18)/(2**60) */
178 ULL(0x8ac7230489e80000), /* _Stl_tenpow[18]=(10**19)/(2**64) */
179 ULL(0xad78ebc5ac620000), /* _Stl_tenpow[19]=(10**20)/(2**67) */
180 ULL(0xd8d726b7177a8000), /* _Stl_tenpow[20]=(10**21)/(2**70) */
181 ULL(0x878678326eac9000), /* _Stl_tenpow[21]=(10**22)/(2**74) */
182 ULL(0xa968163f0a57b400), /* _Stl_tenpow[22]=(10**23)/(2**77) */
183 ULL(0xd3c21bcecceda100), /* _Stl_tenpow[23]=(10**24)/(2**80) */
184 ULL(0x84595161401484a0), /* _Stl_tenpow[24]=(10**25)/(2**84) */
185 ULL(0xa56fa5b99019a5c8), /* _Stl_tenpow[25]=(10**26)/(2**87) */
186 ULL(0xcecb8f27f4200f3a), /* _Stl_tenpow[26]=(10**27)/(2**90) */
187 
188 ULL(0xd0cf4b50cfe20766), /* _Stl_tenpow[27]=(10**55)/(2**183) */
189 ULL(0xd2d80db02aabd62c), /* _Stl_tenpow[28]=(10**83)/(2**276) */
190 ULL(0xd4e5e2cdc1d1ea96), /* _Stl_tenpow[29]=(10**111)/(2**369) */
191 ULL(0xd6f8d7509292d603), /* _Stl_tenpow[30]=(10**139)/(2**462) */
192 ULL(0xd910f7ff28069da4), /* _Stl_tenpow[31]=(10**167)/(2**555) */
193 ULL(0xdb2e51bfe9d0696a), /* _Stl_tenpow[32]=(10**195)/(2**648) */
194 ULL(0xdd50f1996b947519), /* _Stl_tenpow[33]=(10**223)/(2**741) */
195 ULL(0xdf78e4b2bd342cf7), /* _Stl_tenpow[34]=(10**251)/(2**834) */
196 ULL(0xe1a63853bbd26451), /* _Stl_tenpow[35]=(10**279)/(2**927) */
197 ULL(0xe3d8f9e563a198e5), /* _Stl_tenpow[36]=(10**307)/(2**1020) */
198 
199 // /* _Stl_tenpow[36]=(10**335)/(2**) */
200 // /* _Stl_tenpow[36]=(10**335)/(2**) */
201 
202 ULL(0xfd87b5f28300ca0e), /* _Stl_tenpow[37]=(10**-28)/(2**-93) */
203 ULL(0xfb158592be068d2f), /* _Stl_tenpow[38]=(10**-56)/(2**-186) */
204 ULL(0xf8a95fcf88747d94), /* _Stl_tenpow[39]=(10**-84)/(2**-279) */
205 ULL(0xf64335bcf065d37d), /* _Stl_tenpow[40]=(10**-112)/(2**-372) */
206 ULL(0xf3e2f893dec3f126), /* _Stl_tenpow[41]=(10**-140)/(2**-465) */
207 ULL(0xf18899b1bc3f8ca2), /* _Stl_tenpow[42]=(10**-168)/(2**-558) */
208 ULL(0xef340a98172aace5), /* _Stl_tenpow[43]=(10**-196)/(2**-651) */
209 ULL(0xece53cec4a314ebe), /* _Stl_tenpow[44]=(10**-224)/(2**-744) */
210 ULL(0xea9c227723ee8bcb), /* _Stl_tenpow[45]=(10**-252)/(2**-837)     */
211 ULL(0xe858ad248f5c22ca), /* _Stl_tenpow[46]=(10**-280)/(2**-930) */
212 ULL(0xe61acf033d1a45df), /* _Stl_tenpow[47]=(10**-308)/(2**-1023)    */
213 ULL(0xe3e27a444d8d98b8), /* _Stl_tenpow[48]=(10**-336)/(2**-1116) */
214 ULL(0xe1afa13afbd14d6e)  /* _Stl_tenpow[49]=(10**-364)/(2**-1209) */
215 };
216 
217 static const short _Stl_twoexp[80] = {
218 4,7,10,14,17,20,24,27,30,34,37,40,44,47,50,54,57,60,64,67,70,74,77,80,84,87,90,
219 183,276,369,462,555,648,741,834,927,1020,
220 -93,-186,-279,-372,-465,-558,-651,-744,-837,-930,-1023,-1116,-1209
221 };
222 
223 #define  TEN_1  0           /* offset to 10 **   1 */
224 #define  TEN_27   26        /* offset to 10 **  27 */
225 #define  TEN_M28  37        /* offset to 10 ** -28 */
226 #define  NUM_HI_P 11
227 #define  NUM_HI_N 13
228 
229 #define _Stl_HIBITULL (ULL(1) << 63)
230 
_Stl_norm_and_round(uint64 & p,int & norm,uint64 prodhi,uint64 prodlo)231 static void _Stl_norm_and_round(uint64& p, int& norm, uint64 prodhi, uint64 prodlo) {
232   norm = 0;
233   if ((prodhi & _Stl_HIBITULL) == 0) {
234                                 /* leading bit is a zero
235                                  * may have to normalize
236                                  */
237     if ((prodhi == ~_Stl_HIBITULL) &&
238         ((prodlo >> 62) == 0x3)) {  /* normalization followed by round
239                                      * would cause carry to create
240                                      * extra bit, so don't normalize
241                                      */
242       p = _Stl_HIBITULL;
243       return;
244     }
245     p = (prodhi << 1) | (prodlo >> 63); /* normalize */
246     norm = 1;
247     prodlo <<= 1;
248   }
249   else {
250     p = prodhi;
251   }
252 
253   if ((prodlo & _Stl_HIBITULL) != 0) {     /* first guard bit a one */
254     if (((p & 0x1) != 0) ||
255         prodlo != _Stl_HIBITULL ) {    /* not borderline for round to even */
256       /* round */
257       ++p;
258       if (p == 0)
259         ++p;
260     }
261   }
262 }
263 
264 // Convert a 64-bitb fraction * 10^exp to a 64-bit fraction * 2^bexp.
265 // p:    64-bit fraction
266 // exp:  base-10 exponent
267 // bexp: base-2 exponent (output parameter)
_Stl_tenscale(uint64 & p,int exp,int & bexp)268 static void _Stl_tenscale(uint64& p, int exp, int& bexp) {
269   bexp = 0;
270 
271   if ( exp == 0 ) {              /* no scaling needed */
272     return;
273   }
274 
275   int exp_hi = 0, exp_lo = exp; /* exp = exp_hi*32 + exp_lo */
276   int tlo = TEN_1, thi;         /* offsets in power of ten table */
277   int num_hi;                   /* number of high exponent powers */
278 
279   if (exp > 0) {                /* split exponent */
280     if (exp_lo > 27) {
281       exp_lo++;
282       while (exp_lo > 27) {
283         exp_hi++;
284         exp_lo -= 28;
285       }
286     }
287     thi = TEN_27;
288     num_hi = NUM_HI_P;
289   } else { // exp < 0
290     while (exp_lo < 0) {
291       exp_hi++;
292       exp_lo += 28;
293     }
294     thi = TEN_M28;
295     num_hi = NUM_HI_N;
296   }
297 
298   uint64 prodhi, prodlo;        /* 128b product */
299   int norm;                     /* number of bits of normalization */
300 
301   int hi, lo;                   /* offsets in power of ten table */
302   while (exp_hi) {              /* scale */
303     hi = (min) (exp_hi, num_hi);    /* only a few large powers of 10 */
304     exp_hi -= hi;               /* could iterate in extreme case */
305     hi += thi-1;
306     _Stl_mult64(p, _Stl_tenpow[hi], prodhi, prodlo);
307     _Stl_norm_and_round(p, norm, prodhi, prodlo);
308     bexp += _Stl_twoexp[hi] - norm;
309   }
310 
311   if (exp_lo) {
312     lo = tlo + exp_lo -1;
313     _Stl_mult64(p, _Stl_tenpow[lo], prodhi, prodlo);
314     _Stl_norm_and_round(p, norm, prodhi, prodlo);
315     bexp += _Stl_twoexp[lo] - norm;
316   }
317 
318   return;
319 }
320 
321 // First argument is a buffer of values from 0 to 9, NOT ascii.
322 // Second argument is number of digits in buffer, 1 <= digits <= 17.
323 // Third argument is base-10 exponent.
324 
325 /* IEEE representation */
326 #if !defined (__linux__)
327 
328 union _Double_rep {
329   uint64 ival;
330   double val;
331 };
332 
_Stl_atod(char * buffer,ptrdiff_t ndigit,int dexp)333 static double _Stl_atod(char *buffer, ptrdiff_t ndigit, int dexp) {
334   typedef numeric_limits<double> limits;
335   _Double_rep drep;
336   uint64 &value = drep.ival;  /* Value develops as follows:
337                                  * 1) decimal digits as an integer
338                                  * 2) left adjusted fraction
339                                  * 3) right adjusted fraction
340                                  * 4) exponent and fraction
341                                  */
342 
343   uint32 guard;         /* First guard bit */
344   uint64 rest;          /* Remaining guard bits */
345 
346   int bexp;             /* binary exponent */
347   int nzero;            /* number of non-zero bits */
348   int sexp;             /* scaling exponent */
349 
350   char *bufferend;              /* pointer to char after last digit */
351 
352   /* Convert the decimal digits to a binary integer. */
353   bufferend = buffer + ndigit;
354   value = 0;
355 
356   while (buffer < bufferend) {
357     value *= 10;
358     value += *buffer++;
359   }
360 
361   /* Check for zero and treat it as a special case */
362   if (value == 0) {
363     return 0.0;
364   }
365 
366   /* Normalize value */
367   bexp = 64;                    /* convert from 64b int to fraction */
368 
369   /* Count number of non-zeroes in value */
370   nzero = 0;
371   if ((value >> 32) != 0) { nzero  = 32; }    //*TY 03/25/2000 - added explicit comparison to zero to avoid uint64 to bool conversion operator
372   if ((value >> (16 + nzero)) != 0) { nzero += 16; }
373   if ((value >> ( 8 + nzero)) != 0) { nzero +=  8; }
374   if ((value >> ( 4 + nzero)) != 0) { nzero +=  4; }
375   if ((value >> ( 2 + nzero)) != 0) { nzero +=  2; }
376   if ((value >> ( 1 + nzero)) != 0) { nzero +=  1; }
377   if ((value >> (     nzero)) != 0) { nzero +=  1; }
378 
379   /* Normalize */
380   value <<= /*(uint64)*/ (64 - nzero);    //*TY 03/25/2000 - removed extraneous cast to uint64
381   bexp -= 64 - nzero;
382 
383   /* At this point we have a 64b fraction and a binary exponent
384    * but have yet to incorporate the decimal exponent.
385    */
386 
387   /* multiply by 10^dexp */
388   _Stl_tenscale(value, dexp, sexp);
389   bexp += sexp;
390 
391   if (bexp <= -1022) {          /* HI denorm or underflow */
392     bexp += 1022;
393     if (bexp < -53) {          /* guaranteed underflow */
394       value = 0;
395     }
396     else {                      /* denorm or possible underflow */
397       int lead0 = 12 - bexp;          /* 12 sign and exponent bits */
398 
399       /* we must special case right shifts of more than 63 */
400       if (lead0 > 64) {
401         rest = value;
402         guard = 0;
403         value = 0;
404       }
405       else if (lead0 == 64) {
406         rest = value & ((ULL(1)<< 63)-1);
407         guard = (uint32) ((value>> 63) & 1 );
408         value = 0;
409       }
410       else {
411         rest = value & (((ULL(1) << lead0)-1)-1);
412         guard = (uint32) (((value>> lead0)-1) & 1);
413         value >>= /*(uint64)*/ lead0; /* exponent is zero */
414       }
415 
416       /* Round */
417       if (guard && ((value & 1) || rest) ) {
418         ++value;
419         if (value == (ULL(1) << (limits::digits - 1))) { /* carry created normal number */
420           value = 0;
421           _Stl_set_exponent(value, 1);
422         }
423       }
424     }
425   }
426   else {                        /* not zero or denorm */
427     /* Round to 53 bits */
428     rest = value & ((1 << 10) - 1);
429     value >>= 10;
430     guard = (uint32) value & 1;
431     value >>= 1;
432 
433     /*  value&1 guard   rest    Action
434      *
435      *  dc      0       dc      none
436      *  1       1       dc      round
437      *  0       1       0       none
438      *  0       1       !=0     round
439      */
440     if (guard) {
441       if (((value&1)!=0) || (rest!=0)) {
442         ++value;                        /* round */
443         if ((value >> 53) != 0) {       /* carry all the way across */
444           value >>= 1;          /* renormalize */
445           ++bexp;
446         }
447       }
448     }
449     /*
450      * Check for overflow
451      * IEEE Double Precision Format
452      * (From Table 7-8 of Kane and Heinrich)
453      *
454      * Fraction bits               52
455      * Emax                     +1023
456      * Emin                     -1022
457      * Exponent bias            +1023
458      * Exponent bits               11
459      * Integer bit             hidden
460      * Total width in bits         64
461      */
462 
463     if (bexp > limits::max_exponent) {          /* overflow */
464       return limits::infinity();
465     }
466     else {                      /* value is normal */
467       value &= ~(ULL(1) << (limits::digits - 1));   /* hide hidden bit */
468       _Stl_set_exponent(value, bexp + 1022); /* add bias */
469     }
470   }
471 
472   _STLP_STATIC_ASSERT(sizeof(uint64) >= sizeof(double))
473   return drep.val;
474 }
475 
476 #endif
477 
478 #if defined (__linux__) || defined (__MINGW32__) || defined (__CYGWIN__) || \
479     defined (__BORLANDC__) || defined (__DMC__) || defined (__HP_aCC)
480 
481 template <class D, class IEEE, int M, int BIAS>
482 D _Stl_atodT(char *buffer, ptrdiff_t ndigit, int dexp)
483 {
484   typedef numeric_limits<D> limits;
485 
486   /* Convert the decimal digits to a binary integer. */
487   char *bufferend = buffer + ndigit; /* pointer to char after last digit */
488   _ll vv;
489   vv.i64 = 0L;
490 
491   while ( buffer < bufferend ) {
492     vv.i64 *= 10;
493     vv.i64 += *buffer++;
494   }
495 
496   if ( vv.i64 == ULL(0) ) { /* Check for zero and treat it as a special case */
497     return D(0.0);
498   }
499 
500   /* Normalize value */
501 
502   int bexp = 64; /* convert from 64b int to fraction */
503 
504   /* Count number of non-zeroes in value */
505   int nzero = 0;
506   if ((vv.i64 >> 32) != 0) { nzero = 32; }
507   if ((vv.i64 >> (16 + nzero)) != 0) { nzero += 16; }
508   if ((vv.i64 >> ( 8 + nzero)) != 0) { nzero +=  8; }
509   if ((vv.i64 >> ( 4 + nzero)) != 0) { nzero +=  4; }
510   if ((vv.i64 >> ( 2 + nzero)) != 0) { nzero +=  2; }
511   if ((vv.i64 >> ( 1 + nzero)) != 0) { nzero +=  1; }
512   if ((vv.i64 >> (     nzero)) != 0) { nzero +=  1; }
513 
514   /* Normalize */
515   nzero = 64 - nzero;
516   vv.i64 <<= nzero;    // * TY 03/25/2000 - removed extraneous cast to uint64
517   bexp -= nzero;
518 
519   /* At this point we have a 64b fraction and a binary exponent
520    * but have yet to incorporate the decimal exponent.
521    */
522 
523   /* multiply by 10^dexp */
524   int sexp;
525   _Stl_tenscale(vv.i64, dexp, sexp);
526   bexp += sexp;
527 
528   if ( bexp >= limits::min_exponent ) { /* not zero or denorm */
529     if ( limits::digits < 64 ) {
530       /* Round to (64 - M + 1) bits */
531       uint64_t rest = vv.i64 & ((~ULL(0) / ULL(2)) >> (limits::digits - 1));
532       vv.i64 >>= M - 2;
533       uint32_t guard = (uint32) vv.i64 & 1;
534       vv.i64 >>= 1;
535 
536       /*  value&1 guard   rest    Action
537        *
538        *  dc      0       dc      none
539        *  1       1       dc      round
540        *  0       1       0       none
541        *  0       1       !=0     round
542        */
543 
544       if (guard) {
545         if ( ((vv.i64 & 1) != 0) || (rest != 0) ) {
546           vv.i64++;       /* round */
547           if ( (vv.i64 >> (limits::digits < 64 ? limits::digits : 0)) != 0 ) { /* carry all the way across */
548             vv.i64 >>= 1; /* renormalize */
549             ++bexp;
550           }
551         }
552       }
553 
554       vv.i64 &= ~(ULL(1) << (limits::digits - 1)); /* hide hidden bit */
555     }
556     /*
557      * Check for overflow
558      * IEEE Double Precision Format
559      * (From Table 7-8 of Kane and Heinrich)
560      *
561      * Fraction bits               52
562      * Emax                     +1023
563      * Emin                     -1022
564      * Exponent bias            +1023
565      * Exponent bits               11
566      * Integer bit             hidden
567      * Total width in bits         64
568      */
569 
570     if (bexp > limits::max_exponent) { /* overflow */
571       return limits::infinity();
572     }
573 
574     /* value is normal */
575 
576     IEEE v;
577 
578     v.ieee.mantissa0 = vv.i32.hi;
579     v.ieee.mantissa1 = vv.i32.lo;
580     v.ieee.negative = 0;
581     v.ieee.exponent = bexp + BIAS - 1;
582 
583     return v.d;
584   }
585 
586   /* HI denorm or underflow */
587   bexp += BIAS - 1;
588   if (bexp < -limits::digits) { /* guaranteed underflow */
589     vv.i64 = 0;
590   } else {  /* denorm or possible underflow */
591 
592     /*
593      * Problem point for long double: looks like this code reflect shareing of mantissa
594      * and exponent in 64b int; not so for long double
595      */
596 
597     int lead0 = M - bexp; /* M = 12 sign and exponent bits */
598     uint64_t rest;
599     uint32_t guard;
600 
601     /* we must special case right shifts of more than 63 */
602 
603     if (lead0 > 64) {
604       rest = vv.i64;
605       guard = 0;
606       vv.i64 = 0;
607     } else if (lead0 == 64) {
608       rest = vv.i64 & ((ULL(1) << 63)-1);
609       guard = (uint32) ((vv.i64 >> 63) & 1 );
610       vv.i64 = 0;
611     } else {
612       rest = vv.i64 & (((ULL(1) << lead0)-1)-1);
613       guard = (uint32) (((vv.i64 >> lead0)-1) & 1);
614       vv.i64 >>= /*(uint64)*/ lead0; /* exponent is zero */
615     }
616 
617     /* Round */
618     if (guard && ( (vv.i64 & 1) || rest)) {
619       vv.i64++;
620       if (vv.i64 == (ULL(1) << (limits::digits - 1))) { /* carry created normal number */
621         IEEE v;
622 
623         v.ieee.mantissa0 = 0;
624         v.ieee.mantissa1 = 0;
625         v.ieee.negative = 0;
626         v.ieee.exponent = 1;
627         return v.d;
628       }
629     }
630   }
631 
632   IEEE v;
633 
634   v.ieee.mantissa0 = vv.i32.hi;
635   v.ieee.mantissa1 = vv.i32.lo;
636   v.ieee.negative = 0;
637   v.ieee.exponent = 0;
638 
639   return v.d;
640 }
641 #endif // __linux__
642 
643 #ifndef __linux__
_Stl_string_to_double(const char * s)644 static double _Stl_string_to_double(const char *s) {
645   typedef numeric_limits<double> limits;
646   const int max_digits = limits::digits10 + 2;
647   unsigned c;
648   unsigned Negate, decimal_point;
649   char *d;
650   int exp;
651   int dpchar;
652   char digits[max_digits];
653 
654   c = *s++;
655 
656   /* process sign */
657   Negate = 0;
658   if (c == '+') {
659     c = *s++;
660   } else if (c == '-') {
661     Negate = 1;
662     c = *s++;
663   }
664 
665   d = digits;
666   dpchar = '.' - '0';
667   decimal_point = 0;
668   exp = 0;
669 
670   for (;;) {
671     c -= '0';
672     if (c < 10) {
673       if (d == digits + max_digits) {
674         /* ignore more than max_digits digits, but adjust exponent */
675         exp += (decimal_point ^ 1);
676       } else {
677         if (c == 0 && d == digits) {
678           /* ignore leading zeros */
679         } else {
680           *d++ = (char) c;
681         }
682         exp -= decimal_point;
683       }
684     } else if (c == (unsigned int) dpchar && !decimal_point) { /* INTERNATIONAL */
685       decimal_point = 1;
686     } else {
687       break;
688     }
689     c = *s++;
690   }
691 
692   /* strtod cant return until it finds the end of the exponent */
693   if (d == digits) {
694     return 0.0;
695   }
696 
697   if (c == 'e' - '0' || c == 'E' - '0') {
698     register unsigned negate_exp = 0;
699     register int e = 0;
700     c = *s++;
701     if (c == '+' || c == ' ') {
702       c = *s++;
703     } else if (c == '-') {
704       negate_exp = 1;
705       c = *s++;
706     }
707     if (c -= '0', c < 10) {
708       do {
709         e = e * 10 + (int)c;
710         c = *s++;
711       } while (c -= '0', c < 10);
712 
713       if (negate_exp) {
714         e = -e;
715       }
716       exp += e;
717     }
718   }
719 
720   double x;
721   ptrdiff_t n = d - digits;
722   if ((exp + n - 1) < limits::min_exponent10) {
723     x = 0;
724   }
725   else if ((exp + n - 1) > limits::max_exponent10) {
726     x = limits::infinity();
727   }
728   else {
729     /* Let _Stl_atod diagnose under- and over-flows.
730      * If the input was == 0.0, we have already returned,
731      * so retval of +-Inf signals OVERFLOW, 0.0 UNDERFLOW */
732     x = _Stl_atod(digits, n, exp);
733   }
734 
735   if (Negate) {
736     x = -x;
737   }
738 
739   return x;
740 }
741 
742 #endif
743 
744 #if defined (__linux__) || defined (__MINGW32__) || defined (__CYGWIN__) || \
745     defined (__BORLANDC__) || defined (__DMC__) || defined (__HP_aCC)
746 
747 template <class D, class IEEE, int M, int BIAS>
748 D _Stl_string_to_doubleT(const char *s)
749 {
750   typedef numeric_limits<D> limits;
751   const int max_digits = limits::digits10; /* + 2 17 */;
752   unsigned c;
753   unsigned decimal_point;
754   char *d;
755   int exp;
756   D x;
757   int dpchar;
758   char digits[max_digits];
759 
760   c = *s++;
761 
762   /* process sign */
763   bool Negate = false;
764   if (c == '+') {
765     c = *s++;
766   } else if (c == '-') {
767     Negate = true;
768     c = *s++;
769   }
770 
771   d = digits;
772   dpchar = '.' - '0';
773   decimal_point = 0;
774   exp = 0;
775 
776   for (;;) {
777     c -= '0';
778     if (c < 10) {
779       if (d == digits + max_digits) {
780         /* ignore more than max_digits digits, but adjust exponent */
781         exp += (decimal_point ^ 1);
782       } else {
783         if (c == 0 && d == digits) {
784           /* ignore leading zeros */
785         } else {
786           *d++ = (char) c;
787         }
788         exp -= decimal_point;
789       }
790     } else if (c == (unsigned int) dpchar && !decimal_point) {    /* INTERNATIONAL */
791       decimal_point = 1;
792     } else {
793       break;
794     }
795     c = *s++;
796   }
797   /* strtod cant return until it finds the end of the exponent */
798   if (d == digits) {
799     return D(0.0);
800   }
801 
802   if (c == 'e'-'0' || c == 'E'-'0') {
803     bool negate_exp = false;
804     register int e = 0;
805     c = *s++;
806     if (c == '+' || c == ' ') {
807       c = *s++;
808     } else if (c == '-') {
809       negate_exp = true;
810       c = *s++;
811     }
812     if (c -= '0', c < 10) {
813       do {
814         e = e * 10 + (int)c;
815         c = *s++;
816       } while (c -= '0', c < 10);
817 
818       if (negate_exp) {
819         e = -e;
820       }
821       exp += e;
822     }
823   }
824 
825   ptrdiff_t n = d - digits;
826   if ((exp + n - 1) < limits::min_exponent10) {
827     return D(0.0); // +0.0 is the same as -0.0
828   } else if ((exp + n - 1) > limits::max_exponent10 ) {
829     // not good, because of x = -x below; this may lead to portability problems
830     x = limits::infinity();
831   } else {
832     /* let _Stl_atod diagnose under- and over-flows */
833     /* if the input was == 0.0, we have already returned,
834        so retval of +-Inf signals OVERFLOW, 0.0 UNDERFLOW
835     */
836     x = _Stl_atodT<D,IEEE,M,BIAS>(digits, n, exp);
837   }
838 
839   return Negate ? -x : x;
840 }
841 
842 #endif // __linux__
843 
844 void _STLP_CALL
__string_to_float(const __iostring & v,float & val)845 __string_to_float(const __iostring& v, float& val)
846 {
847 #if !defined (__linux__)
848   val = (float)_Stl_string_to_double(v.c_str());
849 #else
850   val = (float)_Stl_string_to_doubleT<double,ieee754_double,12,IEEE754_DOUBLE_BIAS>(v.c_str());
851 #endif
852 }
853 
854 void _STLP_CALL
__string_to_float(const __iostring & v,double & val)855 __string_to_float(const __iostring& v, double& val)
856 {
857 #if !defined (__linux__)
858   val = _Stl_string_to_double(v.c_str());
859 #else
860   val = _Stl_string_to_doubleT<double,ieee754_double,12,IEEE754_DOUBLE_BIAS>(v.c_str());
861 #endif
862 }
863 
864 #if !defined (_STLP_NO_LONG_DOUBLE)
865 void _STLP_CALL
__string_to_float(const __iostring & v,long double & val)866 __string_to_float(const __iostring& v, long double& val) {
867 #if !defined (__linux__) && !defined (__MINGW32__) && !defined (__CYGWIN__) && \
868     !defined (__BORLANDC__) && !defined (__DMC__) && !defined (__HP_aCC)
869   //The following function is valid only if long double is an alias for double.
870   _STLP_STATIC_ASSERT( sizeof(long double) <= sizeof(double) )
871   val = _Stl_string_to_double(v.c_str());
872 #else
873   val = _Stl_string_to_doubleT<long double,ieee854_long_double,16,IEEE854_LONG_DOUBLE_BIAS>(v.c_str());
874 #endif
875 }
876 #endif
877 
878 _STLP_MOVE_TO_STD_NAMESPACE
879 _STLP_END_NAMESPACE
880 
881 // Local Variables:
882 // mode:C++
883 // End:
884