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

Subversion Repositories openrisc_me

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

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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