Subversion Repositories Kolibri OS

Rev

Rev 1693 | 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;
302
	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
1906 serge 303
			if (nd < 9)
304
				y = 10*y + c - '0';
1693 serge 305
		else if (nd < 16)
1906 serge 306
				z = 10*z + c - '0';
1693 serge 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)
1906 serge 331
							y *= 10;
1693 serge 332
					else if (nd <= DBL_DIG + 1)
1906 serge 333
							z *= 10;
1693 serge 334
				if (nd++ < 9)
1906 serge 335
						y = 10*y + c;
1693 serge 336
				else if (nd <= DBL_DIG + 1)
1906 serge 337
						z = 10*z + c;
1693 serge 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