OpenCores
URL https://opencores.org/ocsvn/openrisc_me/openrisc_me/trunk

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [newlib-1.18.0/] [newlib/] [libc/] [stdlib/] [strtod.c] - Blame information for rev 309

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 207 jeremybenn
/*
2
FUNCTION
3
        <<strtod>>, <<strtof>>---string to double or float
4
 
5
INDEX
6
        strtod
7
INDEX
8
        _strtod_r
9
INDEX
10
        strtof
11
 
12
ANSI_SYNOPSIS
13
        #include <stdlib.h>
14
        double strtod(const char *<[str]>, char **<[tail]>);
15
        float strtof(const char *<[str]>, char **<[tail]>);
16
 
17
        double _strtod_r(void *<[reent]>,
18
                         const char *<[str]>, char **<[tail]>);
19
 
20
TRAD_SYNOPSIS
21
        #include <stdlib.h>
22
        double strtod(<[str]>,<[tail]>)
23
        char *<[str]>;
24
        char **<[tail]>;
25
 
26
        float strtof(<[str]>,<[tail]>)
27
        char *<[str]>;
28
        char **<[tail]>;
29
 
30
        double _strtod_r(<[reent]>,<[str]>,<[tail]>)
31
        char *<[reent]>;
32
        char *<[str]>;
33
        char **<[tail]>;
34
 
35
DESCRIPTION
36
        The function <<strtod>> parses the character string <[str]>,
37
        producing a substring which can be converted to a double
38
        value.  The substring converted is the longest initial
39
        subsequence of <[str]>, beginning with the first
40
        non-whitespace character, that has one of these formats:
41
        .[+|-]<[digits]>[.[<[digits]>]][(e|E)[+|-]<[digits]>]
42
        .[+|-].<[digits]>[(e|E)[+|-]<[digits]>]
43
        .[+|-](i|I)(n|N)(f|F)[(i|I)(n|N)(i|I)(t|T)(y|Y)]
44
        .[+|-](n|N)(a|A)(n|N)[<(>[<[hexdigits]>]<)>]
45
        .[+|-]0(x|X)<[hexdigits]>[.[<[hexdigits]>]][(p|P)[+|-]<[digits]>]
46
        .[+|-]0(x|X).<[hexdigits]>[(p|P)[+|-]<[digits]>]
47
        The substring contains no characters if <[str]> is empty, consists
48
        entirely of whitespace, or if the first non-whitespace
49
        character is something other than <<+>>, <<->>, <<.>>, or a
50
        digit, and cannot be parsed as infinity or NaN. If the platform
51
        does not support NaN, then NaN is treated as an empty substring.
52
        If the substring is empty, no conversion is done, and
53
        the value of <[str]> is stored in <<*<[tail]>>>.  Otherwise,
54
        the substring is converted, and a pointer to the final string
55
        (which will contain at least the terminating null character of
56
        <[str]>) is stored in <<*<[tail]>>>.  If you want no
57
        assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>.
58
        <<strtof>> is identical to <<strtod>> except for its return type.
59
 
60
        This implementation returns the nearest machine number to the
61
        input decimal string.  Ties are broken by using the IEEE
62
        round-even rule.  However, <<strtof>> is currently subject to
63
        double rounding errors.
64
 
65
        The alternate function <<_strtod_r>> is a reentrant version.
66
        The extra argument <[reent]> is a pointer to a reentrancy structure.
67
 
68
RETURNS
69
        <<strtod>> returns the converted substring value, if any.  If
70
        no conversion could be performed, 0 is returned.  If the
71
        correct value is out of the range of representable values,
72
        plus or minus <<HUGE_VAL>> is returned, and <<ERANGE>> is
73
        stored in errno. If the correct value would cause underflow, 0
74
        is returned and <<ERANGE>> is stored in errno.
75
 
76
Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
77
<<lseek>>, <<read>>, <<sbrk>>, <<write>>.
78
*/
79
 
80
/****************************************************************
81
 
82
The author of this software is David M. Gay.
83
 
84
Copyright (C) 1998-2001 by Lucent Technologies
85
All Rights Reserved
86
 
87
Permission to use, copy, modify, and distribute this software and
88
its documentation for any purpose and without fee is hereby
89
granted, provided that the above copyright notice appear in all
90
copies and that both that the copyright notice and this
91
permission notice and warranty disclaimer appear in supporting
92
documentation, and that the name of Lucent or any of its entities
93
not be used in advertising or publicity pertaining to
94
distribution of the software without specific, written prior
95
permission.
96
 
97
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
98
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
99
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
100
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
101
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
102
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
103
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
104
THIS SOFTWARE.
105
 
106
****************************************************************/
107
 
108
/* Please send bug reports to David M. Gay (dmg at acm dot org,
109
 * with " at " changed at "@" and " dot " changed to ".").      */
110
 
111
/* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib.  */
112
 
113
#include <_ansi.h>
114
#include <errno.h>
115
#include <stdlib.h>
116
#include <string.h>
117
#include "mprec.h"
118
#include "gdtoa.h"
119
#include "gd_qnan.h"
120
 
121
/* #ifndef NO_FENV_H */
122
/* #include <fenv.h> */
123
/* #endif */
124
 
125
#include "locale.h"
126
 
127
#ifdef IEEE_Arith
128
#ifndef NO_IEEE_Scale
129
#define Avoid_Underflow
130
#undef tinytens
131
/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
132
/* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
133
static _CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
134
                9007199254740992.e-256
135
                };
136
#endif
137
#endif
138
 
139
#ifdef Honor_FLT_ROUNDS
140
#define Rounding rounding
141
#undef Check_FLT_ROUNDS
142
#define Check_FLT_ROUNDS
143
#else
144
#define Rounding Flt_Rounds
145
#endif
146
 
147
#ifndef NO_HEX_FP
148
 
149
static void
150
_DEFUN (ULtod, (L, bits, exp, k),
151
        __ULong *L _AND
152
        __ULong *bits _AND
153
        Long exp _AND
154
        int k)
155
{
156
        switch(k & STRTOG_Retmask) {
157
          case STRTOG_NoNumber:
158
          case STRTOG_Zero:
159
                L[0] = L[1] = 0;
160
                break;
161
 
162
          case STRTOG_Denormal:
163
                L[_1] = bits[0];
164
                L[_0] = bits[1];
165
                break;
166
 
167
          case STRTOG_Normal:
168
          case STRTOG_NaNbits:
169
                L[_1] = bits[0];
170
                L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
171
                break;
172
 
173
          case STRTOG_Infinite:
174
                L[_0] = 0x7ff00000;
175
                L[_1] = 0;
176
                break;
177
 
178
          case STRTOG_NaN:
179
                L[_0] = 0x7fffffff;
180
                L[_1] = (__ULong)-1;
181
          }
182
        if (k & STRTOG_Neg)
183
                L[_0] |= 0x80000000L;
184
}
185
#endif /* !NO_HEX_FP */
186
 
187
#ifdef INFNAN_CHECK
188
static int
189
_DEFUN (match, (sp, t),
190
        _CONST char **sp _AND
191
        char *t)
192
{
193
        int c, d;
194
        _CONST char *s = *sp;
195
 
196
        while( (d = *t++) !=0) {
197
                if ((c = *++s) >= 'A' && c <= 'Z')
198
                        c += 'a' - 'A';
199
                if (c != d)
200
                        return 0;
201
                }
202
        *sp = s + 1;
203
        return 1;
204
}
205
#endif /* INFNAN_CHECK */
206
 
207
 
208
double
209
_DEFUN (_strtod_r, (ptr, s00, se),
210
        struct _reent *ptr _AND
211
        _CONST char *s00 _AND
212
        char **se)
213
{
214
#ifdef Avoid_Underflow
215
        int scale;
216
#endif
217
        int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
218
                 e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
219
        _CONST char *s, *s0, *s1;
220
        double aadj, adj;
221
        U aadj1, rv, rv0;
222
        Long L;
223
        __ULong y, z;
224
        _Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
225
#ifdef SET_INEXACT
226
        int inexact, oldinexact;
227
#endif
228
#ifdef Honor_FLT_ROUNDS
229
        int rounding;
230
#endif
231
 
232
        delta = bs = bd = NULL;
233
        sign = nz0 = nz = decpt = 0;
234
        dval(rv) = 0.;
235
        for(s = s00;;s++) switch(*s) {
236
                case '-':
237
                        sign = 1;
238
                        /* no break */
239
                case '+':
240
                        if (*++s)
241
                                goto break2;
242
                        /* no break */
243
                case 0:
244
                        goto ret0;
245
                case '\t':
246
                case '\n':
247
                case '\v':
248
                case '\f':
249
                case '\r':
250
                case ' ':
251
                        continue;
252
                default:
253
                        goto break2;
254
                }
255
 break2:
256
        if (*s == '0') {
257
#ifndef NO_HEX_FP
258
                {
259
                static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
260
                Long exp;
261
                __ULong bits[2];
262
                switch(s[1]) {
263
                  case 'x':
264
                  case 'X':
265
                        /* If the number is not hex, then the parse of
266
 
267
                        s00 = s + 1;
268
                        {
269
#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
270
                        FPI fpi1 = fpi;
271
                        switch(fegetround()) {
272
                          case FE_TOWARDZERO:   fpi1.rounding = 0; break;
273
                          case FE_UPWARD:       fpi1.rounding = 2; break;
274
                          case FE_DOWNWARD:     fpi1.rounding = 3;
275
                          }
276
#else
277
#define fpi1 fpi
278
#endif
279
                        switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) {
280
                          case STRTOG_NoNumber:
281
                                s = s00;
282
                          case STRTOG_Zero:
283
                                break;
284
                          default:
285
                                if (bb) {
286
                                        copybits(bits, fpi.nbits, bb);
287
                                        Bfree(ptr,bb);
288
                                        }
289
                                ULtod(rv.i, bits, exp, i);
290
                          }}
291
                        goto ret;
292
                  }
293
                }
294
#endif
295
                nz0 = 1;
296
                while(*++s == '0') ;
297
                if (!*s)
298
                        goto ret;
299
                }
300
        s0 = s;
301
        y = z = 0;
302
        for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
303
                if (nd < 9)
304
                        y = 10*y + c - '0';
305
                else if (nd < 16)
306
                        z = 10*z + c - '0';
307
        nd0 = nd;
308
        if (strncmp (s, _localeconv_r (ptr)->decimal_point,
309
                     strlen (_localeconv_r (ptr)->decimal_point)) == 0)
310
                {
311
                decpt = 1;
312
                c = *(s += strlen (_localeconv_r (ptr)->decimal_point));
313
                if (!nd) {
314
                        for(; c == '0'; c = *++s)
315
                                nz++;
316
                        if (c > '0' && c <= '9') {
317
                                s0 = s;
318
                                nf += nz;
319
                                nz = 0;
320
                                goto have_dig;
321
                                }
322
                        goto dig_done;
323
                        }
324
                for(; c >= '0' && c <= '9'; c = *++s) {
325
 have_dig:
326
                        nz++;
327
                        if (c -= '0') {
328
                                nf += nz;
329
                                for(i = 1; i < nz; i++)
330
                                        if (nd++ < 9)
331
                                                y *= 10;
332
                                        else if (nd <= DBL_DIG + 1)
333
                                                z *= 10;
334
                                if (nd++ < 9)
335
                                        y = 10*y + c;
336
                                else if (nd <= DBL_DIG + 1)
337
                                        z = 10*z + c;
338
                                nz = 0;
339
                                }
340
                        }
341
                }
342
 dig_done:
343
        e = 0;
344
        if (c == 'e' || c == 'E') {
345
                if (!nd && !nz && !nz0) {
346
                        goto ret0;
347
                        }
348
                s00 = s;
349
                esign = 0;
350
                switch(c = *++s) {
351
                        case '-':
352
                                esign = 1;
353
                        case '+':
354
                                c = *++s;
355
                        }
356
                if (c >= '0' && c <= '9') {
357
                        while(c == '0')
358
                                c = *++s;
359
                        if (c > '0' && c <= '9') {
360
                                L = c - '0';
361
                                s1 = s;
362
                                while((c = *++s) >= '0' && c <= '9')
363
                                        L = 10*L + c - '0';
364
                                if (s - s1 > 8 || L > 19999)
365
                                        /* Avoid confusion from exponents
366
                                         * so large that e might overflow.
367
                                         */
368
                                        e = 19999; /* safe for 16 bit ints */
369
                                else
370
                                        e = (int)L;
371
                                if (esign)
372
                                        e = -e;
373
                                }
374
                        else
375
                                e = 0;
376
                        }
377
                else
378
                        s = s00;
379
                }
380
        if (!nd) {
381
                if (!nz && !nz0) {
382
#ifdef INFNAN_CHECK
383
                        /* Check for Nan and Infinity */
384
                        __ULong bits[2];
385
                        static FPI fpinan =     /* only 52 explicit bits */
386
                                { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
387
                        if (!decpt)
388
                         switch(c) {
389
                          case 'i':
390
                          case 'I':
391
                                if (match(&s,"nf")) {
392
                                        --s;
393
                                        if (!match(&s,"inity"))
394
                                                ++s;
395
                                        dword0(rv) = 0x7ff00000;
396
#ifndef _DOUBLE_IS_32BITS
397
                                        dword1(rv) = 0;
398
#endif /*!_DOUBLE_IS_32BITS*/
399
                                        goto ret;
400
                                        }
401
                                break;
402
                          case 'n':
403
                          case 'N':
404
                                if (match(&s, "an")) {
405
#ifndef No_Hex_NaN
406
                                        if (*s == '(' /*)*/
407
                                         && hexnan(&s, &fpinan, bits)
408
                                                        == STRTOG_NaNbits) {
409
                                                dword0(rv) = 0x7ff00000 | bits[1];
410
#ifndef _DOUBLE_IS_32BITS
411
                                                dword1(rv) = bits[0];
412
#endif /*!_DOUBLE_IS_32BITS*/
413
                                                }
414
                                        else {
415
#endif
416
                                                dword0(rv) = NAN_WORD0;
417
#ifndef _DOUBLE_IS_32BITS
418
                                                dword1(rv) = NAN_WORD1;
419
#endif /*!_DOUBLE_IS_32BITS*/
420
#ifndef No_Hex_NaN
421
                                                }
422
#endif
423
                                        goto ret;
424
                                        }
425
                          }
426
#endif /* INFNAN_CHECK */
427
 ret0:
428
                        s = s00;
429
                        sign = 0;
430
                        }
431
                goto ret;
432
                }
433
        e1 = e -= nf;
434
 
435
        /* Now we have nd0 digits, starting at s0, followed by a
436
         * decimal point, followed by nd-nd0 digits.  The number we're
437
         * after is the integer represented by those digits times
438
         * 10**e */
439
 
440
        if (!nd0)
441
                nd0 = nd;
442
        k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
443
        dval(rv) = y;
444
        if (k > 9) {
445
#ifdef SET_INEXACT
446
                if (k > DBL_DIG)
447
                        oldinexact = get_inexact();
448
#endif
449
                dval(rv) = tens[k - 9] * dval(rv) + z;
450
                }
451
        bd0 = 0;
452
        if (nd <= DBL_DIG
453
#ifndef RND_PRODQUOT
454
#ifndef Honor_FLT_ROUNDS
455
                && Flt_Rounds == 1
456
#endif
457
#endif
458
                        ) {
459
                if (!e)
460
                        goto ret;
461
                if (e > 0) {
462
                        if (e <= Ten_pmax) {
463
#ifdef VAX
464
                                goto vax_ovfl_check;
465
#else
466
#ifdef Honor_FLT_ROUNDS
467
                                /* round correctly FLT_ROUNDS = 2 or 3 */
468
                                if (sign) {
469
                                        dval(rv) = -dval(rv);
470
                                        sign = 0;
471
                                        }
472
#endif
473
                                /* rv = */ rounded_product(dval(rv), tens[e]);
474
                                goto ret;
475
#endif
476
                                }
477
                        i = DBL_DIG - nd;
478
                        if (e <= Ten_pmax + i) {
479
                                /* A fancier test would sometimes let us do
480
                                 * this for larger i values.
481
                                 */
482
#ifdef Honor_FLT_ROUNDS
483
                                /* round correctly FLT_ROUNDS = 2 or 3 */
484
                                if (sign) {
485
                                        dval(rv) = -dval(rv);
486
                                        sign = 0;
487
                                        }
488
#endif
489
                                e -= i;
490
                                dval(rv) *= tens[i];
491
#ifdef VAX
492
                                /* VAX exponent range is so narrow we must
493
                                 * worry about overflow here...
494
                                 */
495
 vax_ovfl_check:
496
                                dword0(rv) -= P*Exp_msk1;
497
                                /* rv = */ rounded_product(dval(rv), tens[e]);
498
                                if ((dword0(rv) & Exp_mask)
499
                                 > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
500
                                        goto ovfl;
501
                                dword0(rv) += P*Exp_msk1;
502
#else
503
                                /* rv = */ rounded_product(dval(rv), tens[e]);
504
#endif
505
                                goto ret;
506
                                }
507
                        }
508
#ifndef Inaccurate_Divide
509
                else if (e >= -Ten_pmax) {
510
#ifdef Honor_FLT_ROUNDS
511
                        /* round correctly FLT_ROUNDS = 2 or 3 */
512
                        if (sign) {
513
                                dval(rv) = -dval(rv);
514
                                sign = 0;
515
                                }
516
#endif
517
                        /* rv = */ rounded_quotient(dval(rv), tens[-e]);
518
                        goto ret;
519
                        }
520
#endif
521
                }
522
        e1 += nd - k;
523
 
524
#ifdef IEEE_Arith
525
#ifdef SET_INEXACT
526
        inexact = 1;
527
        if (k <= DBL_DIG)
528
                oldinexact = get_inexact();
529
#endif
530
#ifdef Avoid_Underflow
531
        scale = 0;
532
#endif
533
#ifdef Honor_FLT_ROUNDS
534
        if ((rounding = Flt_Rounds) >= 2) {
535
                if (sign)
536
                        rounding = rounding == 2 ? 0 : 2;
537
                else
538
                        if (rounding != 2)
539
                                rounding = 0;
540
                }
541
#endif
542
#endif /*IEEE_Arith*/
543
 
544
        /* Get starting approximation = rv * 10**e1 */
545
 
546
        if (e1 > 0) {
547
                if ( (i = e1 & 15) !=0)
548
                        dval(rv) *= tens[i];
549
                if (e1 &= ~15) {
550
                        if (e1 > DBL_MAX_10_EXP) {
551
 ovfl:
552
#ifndef NO_ERRNO
553
                                ptr->_errno = ERANGE;
554
#endif
555
                                /* Can't trust HUGE_VAL */
556
#ifdef IEEE_Arith
557
#ifdef Honor_FLT_ROUNDS
558
                                switch(rounding) {
559
                                  case 0: /* toward 0 */
560
                                  case 3: /* toward -infinity */
561
                                        dword0(rv) = Big0;
562
#ifndef _DOUBLE_IS_32BITS
563
                                        dword1(rv) = Big1;
564
#endif /*!_DOUBLE_IS_32BITS*/
565
                                        break;
566
                                  default:
567
                                        dword0(rv) = Exp_mask;
568
#ifndef _DOUBLE_IS_32BITS
569
                                        dword1(rv) = 0;
570
#endif /*!_DOUBLE_IS_32BITS*/
571
                                  }
572
#else /*Honor_FLT_ROUNDS*/
573
                                dword0(rv) = Exp_mask;
574
#ifndef _DOUBLE_IS_32BITS
575
                                dword1(rv) = 0;
576
#endif /*!_DOUBLE_IS_32BITS*/
577
#endif /*Honor_FLT_ROUNDS*/
578
#ifdef SET_INEXACT
579
                                /* set overflow bit */
580
                                dval(rv0) = 1e300;
581
                                dval(rv0) *= dval(rv0);
582
#endif
583
#else /*IEEE_Arith*/
584
                                dword0(rv) = Big0;
585
#ifndef _DOUBLE_IS_32BITS
586
                                dword1(rv) = Big1;
587
#endif /*!_DOUBLE_IS_32BITS*/
588
#endif /*IEEE_Arith*/
589
                                if (bd0)
590
                                        goto retfree;
591
                                goto ret;
592
                                }
593
                        e1 >>= 4;
594
                        for(j = 0; e1 > 1; j++, e1 >>= 1)
595
                                if (e1 & 1)
596
                                        dval(rv) *= bigtens[j];
597
                /* The last multiplication could overflow. */
598
                        dword0(rv) -= P*Exp_msk1;
599
                        dval(rv) *= bigtens[j];
600
                        if ((z = dword0(rv) & Exp_mask)
601
                         > Exp_msk1*(DBL_MAX_EXP+Bias-P))
602
                                goto ovfl;
603
                        if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
604
                                /* set to largest number */
605
                                /* (Can't trust DBL_MAX) */
606
                                dword0(rv) = Big0;
607
#ifndef _DOUBLE_IS_32BITS
608
                                dword1(rv) = Big1;
609
#endif /*!_DOUBLE_IS_32BITS*/
610
                                }
611
                        else
612
                                dword0(rv) += P*Exp_msk1;
613
                        }
614
                }
615
        else if (e1 < 0) {
616
                e1 = -e1;
617
                if ( (i = e1 & 15) !=0)
618
                        dval(rv) /= tens[i];
619
                if (e1 >>= 4) {
620
                        if (e1 >= 1 << n_bigtens)
621
                                goto undfl;
622
#ifdef Avoid_Underflow
623
                        if (e1 & Scale_Bit)
624
                                scale = 2*P;
625
                        for(j = 0; e1 > 0; j++, e1 >>= 1)
626
                                if (e1 & 1)
627
                                        dval(rv) *= tinytens[j];
628
                        if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask)
629
                                                >> Exp_shift)) > 0) {
630
                                /* scaled rv is denormal; zap j low bits */
631
                                if (j >= 32) {
632
#ifndef _DOUBLE_IS_32BITS
633
                                        dword1(rv) = 0;
634
#endif /*!_DOUBLE_IS_32BITS*/
635
                                        if (j >= 53)
636
                                         dword0(rv) = (P+2)*Exp_msk1;
637
                                        else
638
                                         dword0(rv) &= 0xffffffff << (j-32);
639
                                        }
640
#ifndef _DOUBLE_IS_32BITS
641
                                else
642
                                        dword1(rv) &= 0xffffffff << j;
643
#endif /*!_DOUBLE_IS_32BITS*/
644
                                }
645
#else
646
                        for(j = 0; e1 > 1; j++, e1 >>= 1)
647
                                if (e1 & 1)
648
                                        dval(rv) *= tinytens[j];
649
                        /* The last multiplication could underflow. */
650
                        dval(rv0) = dval(rv);
651
                        dval(rv) *= tinytens[j];
652
                        if (!dval(rv)) {
653
                                dval(rv) = 2.*dval(rv0);
654
                                dval(rv) *= tinytens[j];
655
#endif
656
                                if (!dval(rv)) {
657
 undfl:
658
                                        dval(rv) = 0.;
659
#ifndef NO_ERRNO
660
                                        ptr->_errno = ERANGE;
661
#endif
662
                                        if (bd0)
663
                                                goto retfree;
664
                                        goto ret;
665
                                        }
666
#ifndef Avoid_Underflow
667
#ifndef _DOUBLE_IS_32BITS
668
                                dword0(rv) = Tiny0;
669
                                dword1(rv) = Tiny1;
670
#else
671
                                dword0(rv) = Tiny1;
672
#endif /*_DOUBLE_IS_32BITS*/
673
                                /* The refinement below will clean
674
                                 * this approximation up.
675
                                 */
676
                                }
677
#endif
678
                        }
679
                }
680
 
681
        /* Now the hard part -- adjusting rv to the correct value.*/
682
 
683
        /* Put digits into bd: true value = bd * 10^e */
684
 
685
        bd0 = s2b(ptr, s0, nd0, nd, y);
686
 
687
        for(;;) {
688
                bd = Balloc(ptr,bd0->_k);
689
                Bcopy(bd, bd0);
690
                bb = d2b(ptr,dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
691
                bs = i2b(ptr,1);
692
 
693
                if (e >= 0) {
694
                        bb2 = bb5 = 0;
695
                        bd2 = bd5 = e;
696
                        }
697
                else {
698
                        bb2 = bb5 = -e;
699
                        bd2 = bd5 = 0;
700
                        }
701
                if (bbe >= 0)
702
                        bb2 += bbe;
703
                else
704
                        bd2 -= bbe;
705
                bs2 = bb2;
706
#ifdef Honor_FLT_ROUNDS
707
                if (rounding != 1)
708
                        bs2++;
709
#endif
710
#ifdef Avoid_Underflow
711
                j = bbe - scale;
712
                i = j + bbbits - 1;     /* logb(rv) */
713
                if (i < Emin)   /* denormal */
714
                        j += P - Emin;
715
                else
716
                        j = P + 1 - bbbits;
717
#else /*Avoid_Underflow*/
718
#ifdef Sudden_Underflow
719
#ifdef IBM
720
                j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
721
#else
722
                j = P + 1 - bbbits;
723
#endif
724
#else /*Sudden_Underflow*/
725
                j = bbe;
726
                i = j + bbbits - 1;     /* logb(rv) */
727
                if (i < Emin)   /* denormal */
728
                        j += P - Emin;
729
                else
730
                        j = P + 1 - bbbits;
731
#endif /*Sudden_Underflow*/
732
#endif /*Avoid_Underflow*/
733
                bb2 += j;
734
                bd2 += j;
735
#ifdef Avoid_Underflow
736
                bd2 += scale;
737
#endif
738
                i = bb2 < bd2 ? bb2 : bd2;
739
                if (i > bs2)
740
                        i = bs2;
741
                if (i > 0) {
742
                        bb2 -= i;
743
                        bd2 -= i;
744
                        bs2 -= i;
745
                        }
746
                if (bb5 > 0) {
747
                        bs = pow5mult(ptr, bs, bb5);
748
                        bb1 = mult(ptr, bs, bb);
749
                        Bfree(ptr, bb);
750
                        bb = bb1;
751
                        }
752
                if (bb2 > 0)
753
                        bb = lshift(ptr, bb, bb2);
754
                if (bd5 > 0)
755
                        bd = pow5mult(ptr, bd, bd5);
756
                if (bd2 > 0)
757
                        bd = lshift(ptr, bd, bd2);
758
                if (bs2 > 0)
759
                        bs = lshift(ptr, bs, bs2);
760
                delta = diff(ptr, bb, bd);
761
                dsign = delta->_sign;
762
                delta->_sign = 0;
763
                i = cmp(delta, bs);
764
#ifdef Honor_FLT_ROUNDS
765
                if (rounding != 1) {
766
                        if (i < 0) {
767
                                /* Error is less than an ulp */
768
                                if (!delta->_x[0] && delta->_wds <= 1) {
769
                                        /* exact */
770
#ifdef SET_INEXACT
771
                                        inexact = 0;
772
#endif
773
                                        break;
774
                                        }
775
                                if (rounding) {
776
                                        if (dsign) {
777
                                                adj = 1.;
778
                                                goto apply_adj;
779
                                                }
780
                                        }
781
                                else if (!dsign) {
782
                                        adj = -1.;
783
                                        if (!dword1(rv)
784
                                         && !(dword0(rv) & Frac_mask)) {
785
                                                y = dword0(rv) & Exp_mask;
786
#ifdef Avoid_Underflow
787
                                                if (!scale || y > 2*P*Exp_msk1)
788
#else
789
                                                if (y)
790
#endif
791
                                                  {
792
                                                  delta = lshift(ptr, delta,Log2P);
793
                                                  if (cmp(delta, bs) <= 0)
794
                                                        adj = -0.5;
795
                                                  }
796
                                                }
797
 apply_adj:
798
#ifdef Avoid_Underflow
799
                                        if (scale && (y = dword0(rv) & Exp_mask)
800
                                                <= 2*P*Exp_msk1)
801
                                          dword0(adj) += (2*P+1)*Exp_msk1 - y;
802
#else
803
#ifdef Sudden_Underflow
804
                                        if ((dword0(rv) & Exp_mask) <=
805
                                                        P*Exp_msk1) {
806
                                                dword0(rv) += P*Exp_msk1;
807
                                                dval(rv) += adj*ulp(dval(rv));
808
                                                dword0(rv) -= P*Exp_msk1;
809
                                                }
810
                                        else
811
#endif /*Sudden_Underflow*/
812
#endif /*Avoid_Underflow*/
813
                                        dval(rv) += adj*ulp(dval(rv));
814
                                        }
815
                                break;
816
                                }
817
                        adj = ratio(delta, bs);
818
                        if (adj < 1.)
819
                                adj = 1.;
820
                        if (adj <= 0x7ffffffe) {
821
                                /* adj = rounding ? ceil(adj) : floor(adj); */
822
                                y = adj;
823
                                if (y != adj) {
824
                                        if (!((rounding>>1) ^ dsign))
825
                                                y++;
826
                                        adj = y;
827
                                        }
828
                                }
829
#ifdef Avoid_Underflow
830
                        if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
831
                                dword0(adj) += (2*P+1)*Exp_msk1 - y;
832
#else
833
#ifdef Sudden_Underflow
834
                        if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
835
                                dword0(rv) += P*Exp_msk1;
836
                                adj *= ulp(dval(rv));
837
                                if (dsign)
838
                                        dval(rv) += adj;
839
                                else
840
                                        dval(rv) -= adj;
841
                                dword0(rv) -= P*Exp_msk1;
842
                                goto cont;
843
                                }
844
#endif /*Sudden_Underflow*/
845
#endif /*Avoid_Underflow*/
846
                        adj *= ulp(dval(rv));
847
                        if (dsign)
848
                                dval(rv) += adj;
849
                        else
850
                                dval(rv) -= adj;
851
                        goto cont;
852
                        }
853
#endif /*Honor_FLT_ROUNDS*/
854
 
855
                if (i < 0) {
856
                        /* Error is less than half an ulp -- check for
857
                         * special case of mantissa a power of two.
858
                         */
859
                        if (dsign || dword1(rv) || dword0(rv) & Bndry_mask
860
#ifdef IEEE_Arith
861
#ifdef Avoid_Underflow
862
                         || (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
863
#else
864
                         || (dword0(rv) & Exp_mask) <= Exp_msk1
865
#endif
866
#endif
867
                                ) {
868
#ifdef SET_INEXACT
869
                                if (!delta->x[0] && delta->wds <= 1)
870
                                        inexact = 0;
871
#endif
872
                                break;
873
                                }
874
                        if (!delta->_x[0] && delta->_wds <= 1) {
875
                                /* exact result */
876
#ifdef SET_INEXACT
877
                                inexact = 0;
878
#endif
879
                                break;
880
                                }
881
                        delta = lshift(ptr,delta,Log2P);
882
                        if (cmp(delta, bs) > 0)
883
                                goto drop_down;
884
                        break;
885
                        }
886
                if (i == 0) {
887
                        /* exactly half-way between */
888
                        if (dsign) {
889
                                if ((dword0(rv) & Bndry_mask1) == Bndry_mask1
890
                                 &&  dword1(rv) == (
891
#ifdef Avoid_Underflow
892
                        (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1)
893
                ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
894
#endif
895
                                                   0xffffffff)) {
896
                                        /*boundary case -- increment exponent*/
897
                                        dword0(rv) = (dword0(rv) & Exp_mask)
898
                                                + Exp_msk1
899
#ifdef IBM
900
                                                | Exp_msk1 >> 4
901
#endif
902
                                                ;
903
#ifndef _DOUBLE_IS_32BITS
904
                                        dword1(rv) = 0;
905
#endif /*!_DOUBLE_IS_32BITS*/
906
#ifdef Avoid_Underflow
907
                                        dsign = 0;
908
#endif
909
                                        break;
910
                                        }
911
                                }
912
                        else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) {
913
 drop_down:
914
                                /* boundary case -- decrement exponent */
915
#ifdef Sudden_Underflow /*{{*/
916
                                L = dword0(rv) & Exp_mask;
917
#ifdef IBM
918
                                if (L <  Exp_msk1)
919
#else
920
#ifdef Avoid_Underflow
921
                                if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
922
#else
923
                                if (L <= Exp_msk1)
924
#endif /*Avoid_Underflow*/
925
#endif /*IBM*/
926
                                        goto undfl;
927
                                L -= Exp_msk1;
928
#else /*Sudden_Underflow}{*/
929
#ifdef Avoid_Underflow
930
                                if (scale) {
931
                                        L = dword0(rv) & Exp_mask;
932
                                        if (L <= (2*P+1)*Exp_msk1) {
933
                                                if (L > (P+2)*Exp_msk1)
934
                                                        /* round even ==> */
935
                                                        /* accept rv */
936
                                                        break;
937
                                                /* rv = smallest denormal */
938
                                                goto undfl;
939
                                                }
940
                                        }
941
#endif /*Avoid_Underflow*/
942
                                L = (dword0(rv) & Exp_mask) - Exp_msk1;
943
#endif /*Sudden_Underflow}*/
944
                                dword0(rv) = L | Bndry_mask1;
945
#ifndef _DOUBLE_IS_32BITS
946
                                dword1(rv) = 0xffffffff;
947
#endif /*!_DOUBLE_IS_32BITS*/
948
#ifdef IBM
949
                                goto cont;
950
#else
951
                                break;
952
#endif
953
                                }
954
#ifndef ROUND_BIASED
955
                        if (!(dword1(rv) & LSB))
956
                                break;
957
#endif
958
                        if (dsign)
959
                                dval(rv) += ulp(dval(rv));
960
#ifndef ROUND_BIASED
961
                        else {
962
                                dval(rv) -= ulp(dval(rv));
963
#ifndef Sudden_Underflow
964
                                if (!dval(rv))
965
                                        goto undfl;
966
#endif
967
                                }
968
#ifdef Avoid_Underflow
969
                        dsign = 1 - dsign;
970
#endif
971
#endif
972
                        break;
973
                        }
974
                if ((aadj = ratio(delta, bs)) <= 2.) {
975
                        if (dsign)
976
                                aadj = dval(aadj1) = 1.;
977
                        else if (dword1(rv) || dword0(rv) & Bndry_mask) {
978
#ifndef Sudden_Underflow
979
                                if (dword1(rv) == Tiny1 && !dword0(rv))
980
                                        goto undfl;
981
#endif
982
                                aadj = 1.;
983
                                dval(aadj1) = -1.;
984
                                }
985
                        else {
986
                                /* special case -- power of FLT_RADIX to be */
987
                                /* rounded down... */
988
 
989
                                if (aadj < 2./FLT_RADIX)
990
                                        aadj = 1./FLT_RADIX;
991
                                else
992
                                        aadj *= 0.5;
993
                                dval(aadj1) = -aadj;
994
                                }
995
                        }
996
                else {
997
                        aadj *= 0.5;
998
                        dval(aadj1) = dsign ? aadj : -aadj;
999
#ifdef Check_FLT_ROUNDS
1000
                        switch(Rounding) {
1001
                                case 2: /* towards +infinity */
1002
                                        dval(aadj1) -= 0.5;
1003
                                        break;
1004
                                case 0: /* towards 0 */
1005
                                case 3: /* towards -infinity */
1006
                                        dval(aadj1) += 0.5;
1007
                                }
1008
#else
1009
                        if (Flt_Rounds == 0)
1010
                                dval(aadj1) += 0.5;
1011
#endif /*Check_FLT_ROUNDS*/
1012
                        }
1013
                y = dword0(rv) & Exp_mask;
1014
 
1015
                /* Check for overflow */
1016
 
1017
                if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
1018
                        dval(rv0) = dval(rv);
1019
                        dword0(rv) -= P*Exp_msk1;
1020
                        adj = dval(aadj1) * ulp(dval(rv));
1021
                        dval(rv) += adj;
1022
                        if ((dword0(rv) & Exp_mask) >=
1023
                                        Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
1024
                                if (dword0(rv0) == Big0 && dword1(rv0) == Big1)
1025
                                        goto ovfl;
1026
                                dword0(rv) = Big0;
1027
#ifndef _DOUBLE_IS_32BITS
1028
                                dword1(rv) = Big1;
1029
#endif /*!_DOUBLE_IS_32BITS*/
1030
                                goto cont;
1031
                                }
1032
                        else
1033
                                dword0(rv) += P*Exp_msk1;
1034
                        }
1035
                else {
1036
#ifdef Avoid_Underflow
1037
                        if (scale && y <= 2*P*Exp_msk1) {
1038
                                if (aadj <= 0x7fffffff) {
1039
                                        if ((z = aadj) <= 0)
1040
                                                z = 1;
1041
                                        aadj = z;
1042
                                        dval(aadj1) = dsign ? aadj : -aadj;
1043
                                        }
1044
                                dword0(aadj1) += (2*P+1)*Exp_msk1 - y;
1045
                                }
1046
                        adj = dval(aadj1) * ulp(dval(rv));
1047
                        dval(rv) += adj;
1048
#else
1049
#ifdef Sudden_Underflow
1050
                        if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) {
1051
                                dval(rv0) = dval(rv);
1052
                                dword0(rv) += P*Exp_msk1;
1053
                                adj = dval(aadj1) * ulp(dval(rv));
1054
                                dval(rv) += adj;
1055
#ifdef IBM
1056
                                if ((dword0(rv) & Exp_mask) <  P*Exp_msk1)
1057
#else
1058
                                if ((dword0(rv) & Exp_mask) <= P*Exp_msk1)
1059
#endif
1060
                                        {
1061
                                        if (dword0(rv0) == Tiny0
1062
                                         && dword1(rv0) == Tiny1)
1063
                                                goto undfl;
1064
#ifndef _DOUBLE_IS_32BITS
1065
                                        dword0(rv) = Tiny0;
1066
                                        dword1(rv) = Tiny1;
1067
#else
1068
                                        dword0(rv) = Tiny1;
1069
#endif /*_DOUBLE_IS_32BITS*/
1070
                                        goto cont;
1071
                                        }
1072
                                else
1073
                                        dword0(rv) -= P*Exp_msk1;
1074
                                }
1075
                        else {
1076
                                adj = dval(aadj1) * ulp(dval(rv));
1077
                                dval(rv) += adj;
1078
                                }
1079
#else /*Sudden_Underflow*/
1080
                        /* Compute adj so that the IEEE rounding rules will
1081
                         * correctly round rv + adj in some half-way cases.
1082
                         * If rv * ulp(rv) is denormalized (i.e.,
1083
                         * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
1084
                         * trouble from bits lost to denormalization;
1085
                         * example: 1.2e-307 .
1086
                         */
1087
                        if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
1088
                                dval(aadj1) = (double)(int)(aadj + 0.5);
1089
                                if (!dsign)
1090
                                        dval(aadj1) = -dval(aadj1);
1091
                                }
1092
                        adj = dval(aadj1) * ulp(dval(rv));
1093
                        dval(rv) += adj;
1094
#endif /*Sudden_Underflow*/
1095
#endif /*Avoid_Underflow*/
1096
                        }
1097
                z = dword0(rv) & Exp_mask;
1098
#ifndef SET_INEXACT
1099
#ifdef Avoid_Underflow
1100
                if (!scale)
1101
#endif
1102
                if (y == z) {
1103
                        /* Can we stop now? */
1104
                        L = (Long)aadj;
1105
                        aadj -= L;
1106
                        /* The tolerances below are conservative. */
1107
                        if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) {
1108
                                if (aadj < .4999999 || aadj > .5000001)
1109
                                        break;
1110
                                }
1111
                        else if (aadj < .4999999/FLT_RADIX)
1112
                                break;
1113
                        }
1114
#endif
1115
 cont:
1116
                Bfree(ptr,bb);
1117
                Bfree(ptr,bd);
1118
                Bfree(ptr,bs);
1119
                Bfree(ptr,delta);
1120
                }
1121
#ifdef SET_INEXACT
1122
        if (inexact) {
1123
                if (!oldinexact) {
1124
                        dword0(rv0) = Exp_1 + (70 << Exp_shift);
1125
#ifndef _DOUBLE_IS_32BITS
1126
                        dword1(rv0) = 0;
1127
#endif /*!_DOUBLE_IS_32BITS*/
1128
                        dval(rv0) += 1.;
1129
                        }
1130
                }
1131
        else if (!oldinexact)
1132
                clear_inexact();
1133
#endif
1134
#ifdef Avoid_Underflow
1135
        if (scale) {
1136
                dword0(rv0) = Exp_1 - 2*P*Exp_msk1;
1137
#ifndef _DOUBLE_IS_32BITS
1138
                dword1(rv0) = 0;
1139
#endif /*!_DOUBLE_IS_32BITS*/
1140
                dval(rv) *= dval(rv0);
1141
#ifndef NO_ERRNO
1142
                /* try to avoid the bug of testing an 8087 register value */
1143
                if (dword0(rv) == 0 && dword1(rv) == 0)
1144
                        ptr->_errno = ERANGE;
1145
#endif
1146
                }
1147
#endif /* Avoid_Underflow */
1148
#ifdef SET_INEXACT
1149
        if (inexact && !(dword0(rv) & Exp_mask)) {
1150
                /* set underflow bit */
1151
                dval(rv0) = 1e-300;
1152
                dval(rv0) *= dval(rv0);
1153
                }
1154
#endif
1155
 retfree:
1156
        Bfree(ptr,bb);
1157
        Bfree(ptr,bd);
1158
        Bfree(ptr,bs);
1159
        Bfree(ptr,bd0);
1160
        Bfree(ptr,delta);
1161
 ret:
1162
        if (se)
1163
                *se = (char *)s;
1164
        return sign ? -dval(rv) : dval(rv);
1165
}
1166
 
1167
#ifndef _REENT_ONLY
1168
 
1169
double
1170
_DEFUN (strtod, (s00, se),
1171
        _CONST char *s00 _AND char **se)
1172
{
1173
  return _strtod_r (_REENT, s00, se);
1174
}
1175
 
1176
float
1177
_DEFUN (strtof, (s00, se),
1178
        _CONST char *s00 _AND
1179
        char **se)
1180
{
1181
  double retval = _strtod_r (_REENT, s00, se);
1182
  if (isnan (retval))
1183
    return nanf (NULL);
1184
  return (float)retval;
1185
}
1186
 
1187
#endif

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.