Subversion Repositories Kolibri OS

Rev

Rev 1906 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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