WebSVN - shark - Blame - Rev 1618 - /shark/branches/xen/libc/stdlib/strtod.c

Rev	Author	Line No.	Line
2	pj	1
		2	/* MODIFIED to remove warnings */
		3
		4	/*-
		5	* Copyright (c) 1993
		6	* The Regents of the University of California. All rights reserved.
		7	*
		8	* Redistribution and use in source and binary forms, with or without
		9	* modification, are permitted provided that the following conditions
		10	* are met:
		11	* 1. Redistributions of source code must retain the above copyright
		12	* notice, this list of conditions and the following disclaimer.
		13	* 2. Redistributions in binary form must reproduce the above copyright
		14	* notice, this list of conditions and the following disclaimer in the
		15	* documentation and/or other materials provided with the distribution.
		16	* 3. All advertising materials mentioning features or use of this software
		17	* must display the following acknowledgement:
		18	* This product includes software developed by the University of
		19	* California, Berkeley and its contributors.
		20	* 4. Neither the name of the University nor the names of its contributors
		21	* may be used to endorse or promote products derived from this software
		22	* without specific prior written permission.
		23	*
		24	* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
		25	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		26	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		27	* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
		28	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		29	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		30	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		31	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		32	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		33	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		34	* SUCH DAMAGE.
		35	*/
		36
		37	#if defined(LIBC_SCCS) && !defined(lint)
		38	static char sccsid[] = "@(#)strtod.c 8.1 (Berkeley) 6/4/93";
		39	#endif /* LIBC_SCCS and not lint */
		40
		41	/****************************************************************
		42	*
		43	* The author of this software is David M. Gay.
		44	*
		45	* Copyright (c) 1991 by AT&T.
		46	*
		47	* Permission to use, copy, modify, and distribute this software for any
		48	* purpose without fee is hereby granted, provided that this entire notice
		49	* is included in all copies of any software which is or includes a copy
		50	* or modification of this software and in all copies of the supporting
		51	* documentation for such software.
		52	*
		53	* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
		54	* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
		55	* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
		56	* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
		57	*
		58	***************************************************************/
		59
		60	/* Please send bug reports to
		61	David M. Gay
		62	AT&T Bell Laboratories, Room 2C-463
		63	600 Mountain Avenue
		64	Murray Hill, NJ 07974-2070
		65	U.S.A.
		66	dmg@research.att.com or research!dmg
		67	*/
		68
		69	/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
		70	*
		71	* This strtod returns a nearest machine number to the input decimal
		72	* string (or sets errno to ERANGE). With IEEE arithmetic, ties are
		73	* broken by the IEEE round-even rule. Otherwise ties are broken by
		74	* biased rounding (add half and chop).
		75	*
		76	* Inspired loosely by William D. Clinger's paper "How to Read Floating
		77	* Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
		78	*
		79	* Modifications:
		80	*
		81	* 1. We only require IEEE, IBM, or VAX double-precision
		82	* arithmetic (not IEEE double-extended).
		83	* 2. We get by with floating-point arithmetic in a case that
		84	* Clinger missed -- when we're computing d * 10^n
		85	* for a small integer d and the integer n is not too
		86	* much larger than 22 (the maximum integer k for which
		87	* we can represent 10^k exactly), we may be able to
		88	* compute (d10^k) 10^(e-k) with just one roundoff.
		89	* 3. Rather than a bit-at-a-time adjustment of the binary
		90	* result in the hard case, we use floating-point
		91	* arithmetic to determine the adjustment to within
		92	* one bit; only in really hard cases do we need to
		93	* compute a second residual.
		94	* 4. Because of 3., we don't need a large table of powers of 10
		95	* for ten-to-e (just some small tables, e.g. of 10^k
		96	* for 0 <= k <= 22).
		97	*/
		98
		99	/*
		100	* #define IEEE_8087 for IEEE-arithmetic machines where the least
		101	* significant byte has the lowest address.
		102	* #define IEEE_MC68k for IEEE-arithmetic machines where the most
		103	* significant byte has the lowest address.
		104	* #define Sudden_Underflow for IEEE-format machines without gradual
		105	* underflow (i.e., that flush to zero on underflow).
		106	* #define IBM for IBM mainframe-style floating-point arithmetic.
		107	* #define VAX for VAX-style floating-point arithmetic.
		108	* #define Unsigned_Shifts if >> does treats its left operand as unsigned.
		109	* #define No_leftright to omit left-right logic in fast floating-point
		110	* computation of dtoa.
		111	* #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
		112	* #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
		113	* that use extended-precision instructions to compute rounded
		114	* products and quotients) with IBM.
		115	* #define ROUND_BIASED for IEEE-format with biased rounding.
		116	* #define Inaccurate_Divide for IEEE-format with correctly rounded
		117	* products but inaccurate quotients, e.g., for Intel i860.
		118	* #define Just_16 to store 16 bits per 32-bit long when doing high-precision
		119	* integer arithmetic. Whether this speeds things up or slows things
		120	* down depends on the machine and the number being converted.
		121	* #define KR_headers for old-style C function headers.
		122	* #define Bad_float_h if your system lacks a float.h or if it does not
		123	* define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
		124	* FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
		125	*/
		126
		127	#if defined(i386) \|\| defined(mips) && defined(MIPSEL)
		128	#define IEEE_8087
		129	#else
		130	#define IEEE_MC68k
		131	#endif
		132
		133	#ifdef DEBUG
		134	#include "stdio.h"
		135	#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
		136	#endif
		137
		138	#ifdef __cplusplus
		139	#include "malloc.h"
		140	#include "memory.h"
		141	#else
		142	#ifndef KR_headers
		143	#include "stdlib.h"
		144	#include "string.h"
		145	#else
		146	#include "malloc.h"
		147	#include "memory.h"
		148	#endif
		149	#endif
		150
		151	#include "errno.h"
		152	#include <ctype.h>
		153	#ifdef Bad_float_h
		154	#undef __STDC__
		155	#ifdef IEEE_MC68k
		156	#define IEEE_ARITHMETIC
		157	#endif
		158	#ifdef IEEE_8087
		159	#define IEEE_ARITHMETIC
		160	#endif
		161	#ifdef IEEE_ARITHMETIC
		162	#define DBL_DIG 15
		163	#define DBL_MAX_10_EXP 308
		164	#define DBL_MAX_EXP 1024
		165	#define FLT_RADIX 2
		166	#define FLT_ROUNDS 1
		167	#define DBL_MAX 1.7976931348623157e+308
		168	#endif
		169
		170	#ifdef IBM
		171	#define DBL_DIG 16
		172	#define DBL_MAX_10_EXP 75
		173	#define DBL_MAX_EXP 63
		174	#define FLT_RADIX 16
		175	#define FLT_ROUNDS 0
		176	#define DBL_MAX 7.2370055773322621e+75
		177	#endif
		178
		179	#ifdef VAX
		180	#define DBL_DIG 16
		181	#define DBL_MAX_10_EXP 38
		182	#define DBL_MAX_EXP 127
		183	#define FLT_RADIX 2
		184	#define FLT_ROUNDS 1
		185	#define DBL_MAX 1.7014118346046923e+38
		186	#endif
		187
		188	#ifndef LONG_MAX
		189	#define LONG_MAX 2147483647
		190	#endif
		191	#else
		192	#include "float.h"
		193	#endif
		194	#ifndef __MATH_H__
		195	#include "math.h"
		196	#endif
		197
		198	#ifdef __cplusplus
		199	extern "C" {
		200	#endif
		201
		202	#ifndef CONST
		203	#ifdef KR_headers
		204	#define CONST /* blank */
		205	#else
		206	#define CONST const
		207	#endif
		208	#endif
		209
		210	#ifdef Unsigned_Shifts
		211	#define Sign_Extend(a,b) if (b < 0) a \|= 0xffff0000;
		212	#else
		213	#define Sign_Extend(a,b) /no-op/
		214	#endif
		215
		216	#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
		217	Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
		218	#endif
		219
		220	#ifdef IEEE_8087
		221	#define word0(x) ((unsigned long *)&x)[1]
		222	#define word1(x) ((unsigned long *)&x)[0]
		223	#else
		224	#define word0(x) ((unsigned long *)&x)[0]
		225	#define word1(x) ((unsigned long *)&x)[1]
		226	#endif
		227
		228	/* The following definition of Storeinc is appropriate for MIPS processors.
		229	* An alternative that might be better on some machines is
		230	* #define Storeinc(a,b,c) (*a++ = b << 16 \| c & 0xffff)
		231	*/
		232	#if defined(IEEE_8087) + defined(VAX)
		233	#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
		234	((unsigned short *)a)[0] = (unsigned short)c, a++)
		235	#else
		236	#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
		237	((unsigned short *)a)[1] = (unsigned short)c, a++)
		238	#endif
		239
		240	/* #define P DBL_MANT_DIG */
		241	/* Ten_pmax = floor(Plog(2)/log(5)) /
		242	/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
		243	/* Quick_max = floor((P-1)log(FLT_RADIX)/log(10) - 1) /
		244	/* Int_max = floor(Plog(FLT_RADIX)/log(10) - 1) /
		245
		246	#if defined(IEEE_8087) + defined(IEEE_MC68k)
		247	#define Exp_shift 20
		248	#define Exp_shift1 20
		249	#define Exp_msk1 0x100000
		250	#define Exp_msk11 0x100000
		251	#define Exp_mask 0x7ff00000
		252	#define P 53
		253	#define Bias 1023
		254	#define IEEE_Arith
		255	#define Emin (-1022)
		256	#define Exp_1 0x3ff00000
		257	#define Exp_11 0x3ff00000
		258	#define Ebits 11
		259	#define Frac_mask 0xfffff
		260	#define Frac_mask1 0xfffff
		261	#define Ten_pmax 22
		262	#define Bletch 0x10
		263	#define Bndry_mask 0xfffff
		264	#define Bndry_mask1 0xfffff
		265	#define LSB 1
		266	#define Sign_bit 0x80000000
		267	#define Log2P 1
		268	#define Tiny0 0
		269	#define Tiny1 1
		270	#define Quick_max 14
		271	#define Int_max 14
		272	#define Infinite(x) (word0(x) == 0x7ff00000) /* sufficient test for here */
		273	#else
		274	#undef Sudden_Underflow
		275	#define Sudden_Underflow
		276	#ifdef IBM
		277	#define Exp_shift 24
		278	#define Exp_shift1 24
		279	#define Exp_msk1 0x1000000
		280	#define Exp_msk11 0x1000000
		281	#define Exp_mask 0x7f000000
		282	#define P 14
		283	#define Bias 65
		284	#define Exp_1 0x41000000
		285	#define Exp_11 0x41000000
		286	#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
		287	#define Frac_mask 0xffffff
		288	#define Frac_mask1 0xffffff
		289	#define Bletch 4
		290	#define Ten_pmax 22
		291	#define Bndry_mask 0xefffff
		292	#define Bndry_mask1 0xffffff
		293	#define LSB 1
		294	#define Sign_bit 0x80000000
		295	#define Log2P 4
		296	#define Tiny0 0x100000
		297	#define Tiny1 0
		298	#define Quick_max 14
		299	#define Int_max 15
		300	#else /* VAX */
		301	#define Exp_shift 23
		302	#define Exp_shift1 7
		303	#define Exp_msk1 0x80
		304	#define Exp_msk11 0x800000
		305	#define Exp_mask 0x7f80
		306	#define P 56
		307	#define Bias 129
		308	#define Exp_1 0x40800000
		309	#define Exp_11 0x4080
		310	#define Ebits 8
		311	#define Frac_mask 0x7fffff
		312	#define Frac_mask1 0xffff007f
		313	#define Ten_pmax 24
		314	#define Bletch 2
		315	#define Bndry_mask 0xffff007f
		316	#define Bndry_mask1 0xffff007f
		317	#define LSB 0x10000
		318	#define Sign_bit 0x8000
		319	#define Log2P 1
		320	#define Tiny0 0x80
		321	#define Tiny1 0
		322	#define Quick_max 15
		323	#define Int_max 15
		324	#endif
		325	#endif
		326
		327	#ifndef IEEE_Arith
		328	#define ROUND_BIASED
		329	#endif
		330
		331	#ifdef RND_PRODQUOT
		332	#define rounded_product(a,b) a = rnd_prod(a, b)
		333	#define rounded_quotient(a,b) a = rnd_quot(a, b)
		334	#ifdef KR_headers
		335	extern double rnd_prod(), rnd_quot();
		336	#else
		337	extern double rnd_prod(double, double), rnd_quot(double, double);
		338	#endif
		339	#else
		340	#define rounded_product(a,b) a *= b
		341	#define rounded_quotient(a,b) a /= b
		342	#endif
		343
		344	#define Big0 (Frac_mask1 \| Exp_msk1*(DBL_MAX_EXP+Bias-1))
		345	#define Big1 0xffffffff
		346
		347	#ifndef Just_16
		348	/* When Pack_32 is not defined, we store 16 bits per 32-bit long.
		349	* This makes some inner loops simpler and sometimes saves work
		350	* during multiplications, but it often seems to make things slightly
		351	* slower. Hence the default is now to store 32 bits per long.
		352	*/
		353	#ifndef Pack_32
		354	#define Pack_32
		355	#endif
		356	#endif
		357
		358	#define Kmax 15
		359
		360	#ifdef __cplusplus
		361	extern "C" double strtod(const char s00, char *se);
		362	extern "C" char *dtoa(double d, int mode, int ndigits,
		363	int decpt, int sign, char **rve);
		364	#endif
		365
		366	struct
		367	Bigint {
		368	struct Bigint *next;
		369	int k, maxwds, sign, wds;
		370	unsigned long x[1];
		371	};
		372
		373	typedef struct Bigint Bigint;
		374
		375	static Bigint *freelist[Kmax+1];
		376
		377	static Bigint *
		378	Balloc
		379	#ifdef KR_headers
		380	(k) int k;
		381	#else
		382	(int k)
		383	#endif
		384	{
		385	int x;
		386	Bigint *rv;
		387
		388	if ( (rv = freelist[k]) ) {
		389	freelist[k] = rv->next;
		390	} else {
		391	x = 1 << k;
		392	rv = (Bigint )malloc(sizeof(Bigint) + (x-1)sizeof(long));
		393	rv->k = k;
		394	rv->maxwds = x;
		395	}
		396	rv->sign = rv->wds = 0;
		397	return rv;
		398	}
		399
		400	static void
		401	Bfree
		402	#ifdef KR_headers
		403	(v) Bigint *v;
		404	#else
		405	(Bigint *v)
		406	#endif
		407	{
		408	if (v) {
		409	v->next = freelist[v->k];
		410	freelist[v->k] = v;
		411	}
		412	}
		413
		414	#define Bcopy(x,y) memcpy((char )&x->sign, (char )&y->sign, \
		415	y->wdssizeof(long) + 2sizeof(int))
		416
		417	static Bigint *
		418	multadd
		419	#ifdef KR_headers
		420	(b, m, a) Bigint *b; int m, a;
		421	#else
		422	(Bigint b, int m, int a) / multiply by m and add a */
		423	#endif
		424	{
		425	int i, wds;
		426	unsigned long *x, y;
		427	#ifdef Pack_32
		428	unsigned long xi, z;
		429	#endif
		430	Bigint *b1;
		431
		432	wds = b->wds;
		433	x = b->x;
		434	i = 0;
		435	do {
		436	#ifdef Pack_32
		437	xi = *x;
		438	y = (xi & 0xffff) * m + a;
		439	z = (xi >> 16) * m + (y >> 16);
		440	a = (int)(z >> 16);
		441	*x++ = (z << 16) + (y & 0xffff);
		442	#else
		443	y = x m + a;
		444	a = (int)(y >> 16);
		445	*x++ = y & 0xffff;
		446	#endif
		447	} while (++i < wds);
		448	if (a) {
		449	if (wds >= b->maxwds) {
		450	b1 = Balloc(b->k+1);
		451	Bcopy(b1, b);
		452	Bfree(b);
		453	b = b1;
		454	}
		455	b->x[wds++] = a;
		456	b->wds = wds;
		457	}
		458	return b;
		459	}
		460
		461	static Bigint *
		462	s2b
		463	#ifdef KR_headers
		464	(s, nd0, nd, y9) CONST char *s; int nd0, nd; unsigned long y9;
		465	#else
		466	(CONST char *s, int nd0, int nd, unsigned long y9)
		467	#endif
		468	{
		469	Bigint *b;
		470	int i, k;
		471	long x, y;
		472
		473	x = (nd + 8) / 9;
		474	for (k = 0, y = 1; x > y; y <<= 1, k++) ;
		475	#ifdef Pack_32
		476	b = Balloc(k);
		477	b->x[0] = y9;
		478	b->wds = 1;
		479	#else
		480	b = Balloc(k+1);
		481	b->x[0] = y9 & 0xffff;
		482	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
		483	#endif
		484
		485	i = 9;
		486	if (9 < nd0) {
		487	s += 9;
		488	do
		489	b = multadd(b, 10, *s++ - '0');
		490	while (++i < nd0);
		491	s++;
		492	} else
		493	s += 10;
		494	for (; i < nd; i++)
		495	b = multadd(b, 10, *s++ - '0');
		496	return b;
		497	}
		498
		499	static int
		500	hi0bits
		501	#ifdef KR_headers
		502	(x) register unsigned long x;
		503	#else
		504	(register unsigned long x)
		505	#endif
		506	{
		507	register int k = 0;
		508
		509	if (!(x & 0xffff0000)) {
		510	k = 16;
		511	x <<= 16;
		512	}
		513	if (!(x & 0xff000000)) {
		514	k += 8;
		515	x <<= 8;
		516	}
		517	if (!(x & 0xf0000000)) {
		518	k += 4;
		519	x <<= 4;
		520	}
		521	if (!(x & 0xc0000000)) {
		522	k += 2;
		523	x <<= 2;
		524	}
		525	if (!(x & 0x80000000)) {
		526	k++;
		527	if (!(x & 0x40000000))
		528	return 32;
		529	}
		530	return k;
		531	}
		532
		533	static int
		534	lo0bits
		535	#ifdef KR_headers
		536	(y) unsigned long *y;
		537	#else
		538	(unsigned long *y)
		539	#endif
		540	{
		541	register int k;
		542	register unsigned long x = *y;
		543
		544	if (x & 7) {
		545	if (x & 1)
		546	return 0;
		547	if (x & 2) {
		548	*y = x >> 1;
		549	return 1;
		550	}
		551	*y = x >> 2;
		552	return 2;
		553	}
		554	k = 0;
		555	if (!(x & 0xffff)) {
		556	k = 16;
		557	x >>= 16;
		558	}
		559	if (!(x & 0xff)) {
		560	k += 8;
		561	x >>= 8;
		562	}
		563	if (!(x & 0xf)) {
		564	k += 4;
		565	x >>= 4;
		566	}
		567	if (!(x & 0x3)) {
		568	k += 2;
		569	x >>= 2;
		570	}
		571	if (!(x & 1)) {
		572	k++;
		573	x >>= 1;
		574	if (!x & 1)
		575	return 32;
		576	}
		577	*y = x;
		578	return k;
		579	}
		580
		581	static Bigint *
		582	i2b
		583	#ifdef KR_headers
		584	(i) int i;
		585	#else
		586	(int i)
		587	#endif
		588	{
		589	Bigint *b;
		590
		591	b = Balloc(1);
		592	b->x[0] = i;
		593	b->wds = 1;
		594	return b;
		595	}
		596
		597	static Bigint *
		598	mult
		599	#ifdef KR_headers
		600	(a, b) Bigint a, b;
		601	#else
		602	(Bigint a, Bigint b)
		603	#endif
		604	{
		605	Bigint *c;
		606	int k, wa, wb, wc;
		607	unsigned long carry, y, z;
		608	unsigned long x, xa, xae, xb, xbe, xc, *xc0;
		609	#ifdef Pack_32
		610	unsigned long z2;
		611	#endif
		612
		613	if (a->wds < b->wds) {
		614	c = a;
		615	a = b;
		616	b = c;
		617	}
		618	k = a->k;
		619	wa = a->wds;
		620	wb = b->wds;
		621	wc = wa + wb;
		622	if (wc > a->maxwds)
		623	k++;
		624	c = Balloc(k);
		625	for (x = c->x, xa = x + wc; x < xa; x++)
		626	*x = 0;
		627	xa = a->x;
		628	xae = xa + wa;
		629	xb = b->x;
		630	xbe = xb + wb;
		631	xc0 = c->x;
		632	#ifdef Pack_32
		633	for (; xb < xbe; xb++, xc0++) {
		634	if ( (y = *xb & 0xffff) ) {
		635	x = xa;
		636	xc = xc0;
		637	carry = 0;
		638	do {
		639	z = (x & 0xffff) y + (*xc & 0xffff) + carry;
		640	carry = z >> 16;
		641	z2 = (x++ >> 16) y + (*xc >> 16) + carry;
		642	carry = z2 >> 16;
		643	Storeinc(xc, z2, z);
		644	} while (x < xae);
		645	*xc = carry;
		646	}
		647	if ( (y = *xb >> 16) ) {
		648	x = xa;
		649	xc = xc0;
		650	carry = 0;
		651	z2 = *xc;
		652	do {
		653	z = (x & 0xffff) y + (*xc >> 16) + carry;
		654	carry = z >> 16;
		655	Storeinc(xc, z, z2);
		656	z2 = (x++ >> 16) y + (*xc & 0xffff) + carry;
		657	carry = z2 >> 16;
		658	} while (x < xae);
		659	*xc = z2;
		660	}
		661	}
		662	#else
		663	for (; xb < xbe; xc0++) {
		664	if (y = *xb++) {
		665	x = xa;
		666	xc = xc0;
		667	carry = 0;
		668	do {
		669	z = x++ y + *xc + carry;
		670	carry = z >> 16;
		671	*xc++ = z & 0xffff;
		672	} while (x < xae);
		673	*xc = carry;
		674	}
		675	}
		676	#endif
		677	for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
		678	c->wds = wc;
		679	return c;
		680	}
		681
		682	static Bigint *p5s;
		683
		684	static Bigint *
		685	pow5mult
		686	#ifdef KR_headers
		687	(b, k) Bigint *b; int k;
		688	#else
		689	(Bigint *b, int k)
		690	#endif
		691	{
		692	Bigint b1, p5, *p51;
		693	int i;
		694	static int p05[3] = { 5, 25, 125 };
		695
		696	if ( (i = k & 3) )
		697	b = multadd(b, p05[i-1], 0);
		698
		699	if (!(k >>= 2))
		700	return b;
		701	if (!(p5 = p5s)) {
		702	/* first time */
		703	p5 = p5s = i2b(625);
		704	p5->next = 0;
		705	}
		706	for (;;) {
		707	if (k & 1) {
		708	b1 = mult(b, p5);
		709	Bfree(b);
		710	b = b1;
		711	}
		712	if (!(k >>= 1))
		713	break;
		714	if (!(p51 = p5->next)) {
		715	p51 = p5->next = mult(p5,p5);
		716	p51->next = 0;
		717	}
		718	p5 = p51;
		719	}
		720	return b;
		721	}
		722
		723	static Bigint *
		724	lshift
		725	#ifdef KR_headers
		726	(b, k) Bigint *b; int k;
		727	#else
		728	(Bigint *b, int k)
		729	#endif
		730	{
		731	int i, k1, n, n1;
		732	Bigint *b1;
		733	unsigned long x, x1, *xe, z;
		734
		735	#ifdef Pack_32
		736	n = k >> 5;
		737	#else
		738	n = k >> 4;
		739	#endif
		740	k1 = b->k;
		741	n1 = n + b->wds + 1;
		742	for (i = b->maxwds; n1 > i; i <<= 1)
		743	k1++;
		744	b1 = Balloc(k1);
		745	x1 = b1->x;
		746	for (i = 0; i < n; i++)
		747	*x1++ = 0;
		748	x = b->x;
		749	xe = x + b->wds;
		750	#ifdef Pack_32
		751	if (k &= 0x1f) {
		752	k1 = 32 - k;
		753	z = 0;
		754	do {
		755	x1++ = x << k \| z;
		756	z = *x++ >> k1;
		757	} while (x < xe);
		758	if ( (*x1 = z) )
		759	++n1;
		760	}
		761	#else
		762	if (k &= 0xf) {
		763	k1 = 16 - k;
		764	z = 0;
		765	do {
		766	x1++ = x << k & 0xffff \| z;
		767	z = *x++ >> k1;
		768	} while (x < xe);
		769	if (*x1 = z)
		770	++n1;
		771	}
		772	#endif
		773	else
		774	do
		775	x1++ = x++;
		776	while (x < xe);
		777	b1->wds = n1 - 1;
		778	Bfree(b);
		779	return b1;
		780	}
		781
		782	static int
		783	cmp
		784	#ifdef KR_headers
		785	(a, b) Bigint a, b;
		786	#else
		787	(Bigint a, Bigint b)
		788	#endif
		789	{
		790	unsigned long xa, xa0, xb, xb0;
		791	int i, j;
		792
		793	i = a->wds;
		794	j = b->wds;
		795	#ifdef DEBUG
		796	if (i > 1 && !a->x[i-1])
		797	Bug("cmp called with a->x[a->wds-1] == 0");
		798	if (j > 1 && !b->x[j-1])
		799	Bug("cmp called with b->x[b->wds-1] == 0");
		800	#endif
		801	if (i -= j)
		802	return i;
		803	xa0 = a->x;
		804	xa = xa0 + j;
		805	xb0 = b->x;
		806	xb = xb0 + j;
		807	for (;;) {
		808	if (--xa != --xb)
		809	return xa < xb ? -1 : 1;
		810	if (xa <= xa0)
		811	break;
		812	}
		813	return 0;
		814	}
		815
		816	static Bigint *
		817	diff
		818	#ifdef KR_headers
		819	(a, b) Bigint a, b;
		820	#else
		821	(Bigint a, Bigint b)
		822	#endif
		823	{
		824	Bigint *c;
		825	int i, wa, wb;
		826	long borrow, y; /* We need signed shifts here. */
		827	unsigned long xa, xae, xb, xbe, *xc;
		828	#ifdef Pack_32
		829	long z;
		830	#endif
		831
		832	i = cmp(a,b);
		833	if (!i) {
		834	c = Balloc(0);
		835	c->wds = 1;
		836	c->x[0] = 0;
		837	return c;
		838	}
		839	if (i < 0) {
		840	c = a;
		841	a = b;
		842	b = c;
		843	i = 1;
		844	} else
		845	i = 0;
		846	c = Balloc(a->k);
		847	c->sign = i;
		848	wa = a->wds;
		849	xa = a->x;
		850	xae = xa + wa;
		851	wb = b->wds;
		852	xb = b->x;
		853	xbe = xb + wb;
		854	xc = c->x;
		855	borrow = 0;
		856	#ifdef Pack_32
		857	do {
		858	y = (xa & 0xffff) - (xb & 0xffff) + borrow;
		859	borrow = y >> 16;
		860	Sign_Extend(borrow, y);
		861	z = (xa++ >> 16) - (xb++ >> 16) + borrow;
		862	borrow = z >> 16;
		863	Sign_Extend(borrow, z);
		864	Storeinc(xc, z, y);
		865	} while (xb < xbe);
		866	while (xa < xae) {
		867	y = (*xa & 0xffff) + borrow;
		868	borrow = y >> 16;
		869	Sign_Extend(borrow, y);
		870	z = (*xa++ >> 16) + borrow;
		871	borrow = z >> 16;
		872	Sign_Extend(borrow, z);
		873	Storeinc(xc, z, y);
		874	}
		875	#else
		876	do {
		877	y = xa++ - xb++ + borrow;
		878	borrow = y >> 16;
		879	Sign_Extend(borrow, y);
		880	*xc++ = y & 0xffff;
		881	} while (xb < xbe);
		882	while (xa < xae) {
		883	y = *xa++ + borrow;
		884	borrow = y >> 16;
		885	Sign_Extend(borrow, y);
		886	*xc++ = y & 0xffff;
		887	}
		888	#endif
		889	while (!*--xc)
		890	wa--;
		891	c->wds = wa;
		892	return c;
		893	}
		894
		895	static double
		896	ulp
		897	#ifdef KR_headers
		898	(x) double x;
		899	#else
		900	(double x)
		901	#endif
		902	{
		903	register long L;
		904	double a;
		905
		906	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
		907	#ifndef Sudden_Underflow
		908	if (L > 0) {
		909	#endif
		910	#ifdef IBM
		911	L \|= Exp_msk1 >> 4;
		912	#endif
		913	word0(a) = L;
		914	word1(a) = 0;
		915	#ifndef Sudden_Underflow
		916	} else {
		917	L = -L >> Exp_shift;
		918	if (L < Exp_shift) {
		919	word0(a) = 0x80000 >> L;
		920	word1(a) = 0;
		921	} else {
		922	word0(a) = 0;
		923	L -= Exp_shift;
		924	word1(a) = L >= 31 ? 1 : 1 << (31 - L);
		925	}
		926	}
		927	#endif
		928	return a;
		929	}
		930
		931	static double
		932	b2d
		933	#ifdef KR_headers
		934	(a, e) Bigint a; int e;
		935	#else
		936	(Bigint a, int e)
		937	#endif
		938	{
		939	unsigned long xa, xa0, w, y, z;
		940	int k;
		941	double d;
		942	#ifdef VAX
		943	unsigned long d0, d1;
		944	#else
		945	#define d0 word0(d)
		946	#define d1 word1(d)
		947	#endif
		948
		949	xa0 = a->x;
		950	xa = xa0 + a->wds;
		951	y = *--xa;
		952	#ifdef DEBUG
		953	if (!y) Bug("zero y in b2d");
		954	#endif
		955	k = hi0bits(y);
		956	*e = 32 - k;
		957	#ifdef Pack_32
		958	if (k < Ebits) {
		959	d0 = Exp_1 \| (y >> (Ebits - k));
		960	w = xa > xa0 ? *--xa : 0;
		961	d1 = (y << ((32-Ebits) + k)) \| (w >> (Ebits - k));
		962	goto ret_d;
		963	}
		964	z = xa > xa0 ? *--xa : 0;
		965	if (k -= Ebits) {
		966	d0 = Exp_1 \| (y << k) \| (z >> (32 - k));
		967	y = xa > xa0 ? *--xa : 0;
		968	d1 = (z << k) \| (y >> (32 - k));
		969	} else {
		970	d0 = Exp_1 \| y;
		971	d1 = z;
		972	}
		973	#else
		974	if (k < Ebits + 16) {
		975	z = xa > xa0 ? *--xa : 0;
		976	d0 = Exp_1 \| y << k - Ebits \| z >> Ebits + 16 - k;
		977	w = xa > xa0 ? *--xa : 0;
		978	y = xa > xa0 ? *--xa : 0;
		979	d1 = z << k + 16 - Ebits \| w << k - Ebits \| y >> 16 + Ebits - k;
		980	goto ret_d;
		981	}
		982	z = xa > xa0 ? *--xa : 0;
		983	w = xa > xa0 ? *--xa : 0;
		984	k -= Ebits + 16;
		985	d0 = Exp_1 \| y << k + 16 \| z << k \| w >> 16 - k;
		986	y = xa > xa0 ? *--xa : 0;
		987	d1 = w << k + 16 \| y << k;
		988	#endif
		989	ret_d:
		990	#ifdef VAX
		991	word0(d) = d0 >> 16 \| d0 << 16;
		992	word1(d) = d1 >> 16 \| d1 << 16;
		993	#else
		994	#undef d0
		995	#undef d1
		996	#endif
		997	return d;
		998	}
		999
		1000	static Bigint *
		1001	d2b
		1002	#ifdef KR_headers
		1003	(d, e, bits) double d; int e, bits;
		1004	#else
		1005	(double d, int e, int bits)
		1006	#endif
		1007	{
		1008	Bigint *b;
		1009	int de, i, k;
		1010	unsigned long *x, y, z;
		1011	#ifdef VAX
		1012	unsigned long d0, d1;
		1013	d0 = word0(d) >> 16 \| word0(d) << 16;
		1014	d1 = word1(d) >> 16 \| word1(d) << 16;
		1015	#else
		1016	#define d0 word0(d)
		1017	#define d1 word1(d)
		1018	#endif
		1019
		1020	#ifdef Pack_32
		1021	b = Balloc(1);
		1022	#else
		1023	b = Balloc(2);
		1024	#endif
		1025	x = b->x;
		1026
		1027	z = d0 & Frac_mask;
		1028	d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
		1029	#ifdef Sudden_Underflow
		1030	de = (int)(d0 >> Exp_shift);
		1031	#ifndef IBM
		1032	z \|= Exp_msk11;
		1033	#endif
		1034	#else
		1035	if ( (de = (int)(d0 >> Exp_shift)) )
		1036	z \|= Exp_msk1;
		1037	#endif
		1038	#ifdef Pack_32
		1039	if ( (y = d1) ) {
		1040	if ( (k = lo0bits(&y)) ) {
		1041	x[0] = y \| (z << (32 - k));
		1042	z >>= k;
		1043	}
		1044	else
		1045	x[0] = y;
		1046	i = b->wds = (x[1] = z) ? 2 : 1;
		1047	} else {
		1048	#ifdef DEBUG
		1049	if (!z)
		1050	Bug("Zero passed to d2b");
		1051	#endif
		1052	k = lo0bits(&z);
		1053	x[0] = z;
		1054	i = b->wds = 1;
		1055	k += 32;
		1056	}
		1057	#else
		1058	if (y = d1) {
		1059	if (k = lo0bits(&y))
		1060	if (k >= 16) {
		1061	x[0] = y \| z << 32 - k & 0xffff;
		1062	x[1] = z >> k - 16 & 0xffff;
		1063	x[2] = z >> k;
		1064	i = 2;
		1065	} else {
		1066	x[0] = y & 0xffff;
		1067	x[1] = y >> 16 \| z << 16 - k & 0xffff;
		1068	x[2] = z >> k & 0xffff;
		1069	x[3] = z >> k+16;
		1070	i = 3;
		1071	}
		1072	else {
		1073	x[0] = y & 0xffff;
		1074	x[1] = y >> 16;
		1075	x[2] = z & 0xffff;
		1076	x[3] = z >> 16;
		1077	i = 3;
		1078	}
		1079	} else {
		1080	#ifdef DEBUG
		1081	if (!z)
		1082	Bug("Zero passed to d2b");
		1083	#endif
		1084	k = lo0bits(&z);
		1085	if (k >= 16) {
		1086	x[0] = z;
		1087	i = 0;
		1088	} else {
		1089	x[0] = z & 0xffff;
		1090	x[1] = z >> 16;
		1091	i = 1;
		1092	}
		1093	k += 32;
		1094	}
		1095	while (!x[i])
		1096	--i;
		1097	b->wds = i + 1;
		1098	#endif
		1099	#ifndef Sudden_Underflow
		1100	if (de) {
		1101	#endif
		1102	#ifdef IBM
		1103	*e = (de - Bias - (P-1) << 2) + k;
		1104	bits = 4P + 8 - k - hi0bits(word0(d) & Frac_mask);
		1105	#else
		1106	*e = de - Bias - (P-1) + k;
		1107	*bits = P - k;
		1108	#endif
		1109	#ifndef Sudden_Underflow
		1110	} else {
		1111	*e = de - Bias - (P-1) + 1 + k;
		1112	#ifdef Pack_32
		1113	bits = 32i - hi0bits(x[i-1]);
		1114	#else
		1115	bits = (i+2)16 - hi0bits(x[i]);
		1116	#endif
		1117	}
		1118	#endif
		1119	return b;
		1120	}
		1121	#undef d0
		1122	#undef d1
		1123
		1124	static double
		1125	ratio
		1126	#ifdef KR_headers
		1127	(a, b) Bigint a, b;
		1128	#else
		1129	(Bigint a, Bigint b)
		1130	#endif
		1131	{
		1132	double da, db;
		1133	int k, ka, kb;
		1134
		1135	da = b2d(a, &ka);
		1136	db = b2d(b, &kb);
		1137	#ifdef Pack_32
		1138	k = ka - kb + 32*(a->wds - b->wds);
		1139	#else
		1140	k = ka - kb + 16*(a->wds - b->wds);
		1141	#endif
		1142	#ifdef IBM
		1143	if (k > 0) {
		1144	word0(da) += (k >> 2)*Exp_msk1;
		1145	if (k &= 3)
		1146	da *= 1 << k;
		1147	} else {
		1148	k = -k;
		1149	word0(db) += (k >> 2)*Exp_msk1;
		1150	if (k &= 3)
		1151	db *= 1 << k;
		1152	}
		1153	#else
		1154	if (k > 0)
		1155	word0(da) += k*Exp_msk1;
		1156	else {
		1157	k = -k;
		1158	word0(db) += k*Exp_msk1;
		1159	}
		1160	#endif
		1161	return da / db;
		1162	}
		1163
		1164	static double
		1165	tens[] = {
		1166	1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
		1167	1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
		1168	1e20, 1e21, 1e22
		1169	#ifdef VAX
		1170	, 1e23, 1e24
		1171	#endif
		1172	};
		1173
		1174	static double
		1175	#ifdef IEEE_Arith
		1176	bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
		1177	static double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };
		1178	#define n_bigtens 5
		1179	#else
		1180	#ifdef IBM
		1181	bigtens[] = { 1e16, 1e32, 1e64 };
		1182	static double tinytens[] = { 1e-16, 1e-32, 1e-64 };
		1183	#define n_bigtens 3
		1184	#else
		1185	bigtens[] = { 1e16, 1e32 };
		1186	static double tinytens[] = { 1e-16, 1e-32 };
		1187	#define n_bigtens 2
		1188	#endif
		1189	#endif
		1190
		1191	/* MG */
		1192	double
		1193	torelink__strtod
		1194	#ifdef KR_headers
		1195	(s00, se) CONST char s00; char *se;
		1196	#else
		1197	(CONST char s00, char *se)
		1198	#endif
		1199	{
		1200	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign,
		1201	e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
		1202	CONST char s, s0, *s1;
		1203	double aadj, aadj1, adj, rv, rv0;
		1204	long L;
		1205	unsigned long y, z;
		1206	Bigint bb, bb1, bd, bd0, bs, delta;
		1207	sign = nz0 = nz = 0;
		1208	rv = 0.;
		1209	for (s = s00;;s++) switch(*s) {
		1210	case '-':
		1211	sign = 1;
		1212	/* no break */
		1213	case '+':
		1214	if (*++s)
		1215	goto break2;
		1216	/* no break */
		1217	case 0:
		1218	s = s00;
		1219	goto ret;
		1220	default:
		1221	if (isspace((unsigned char)*s))
		1222	continue;
		1223	goto break2;
		1224	}
		1225	break2:
		1226	if (*s == '0') {
		1227	nz0 = 1;
		1228	while (*++s == '0') ;
		1229	if (!*s)
		1230	goto ret;
		1231	}
		1232	s0 = s;
		1233	y = z = 0;
		1234	for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
		1235	if (nd < 9)
		1236	y = 10*y + c - '0';
		1237	else if (nd < 16)
		1238	z = 10*z + c - '0';
		1239	nd0 = nd;
		1240	if (c == '.') {
		1241	c = *++s;
		1242	if (!nd) {
		1243	for (; c == '0'; c = *++s)
		1244	nz++;
		1245	if (c > '0' && c <= '9') {
		1246	s0 = s;
		1247	nf += nz;
		1248	nz = 0;
		1249	goto have_dig;
		1250	}
		1251	goto dig_done;
		1252	}
		1253	for (; c >= '0' && c <= '9'; c = *++s) {
		1254	have_dig:
		1255	nz++;
		1256	if (c -= '0') {
		1257	nf += nz;
		1258	for (i = 1; i < nz; i++)
		1259	if (nd++ < 9)
		1260	y *= 10;
		1261	else if (nd <= DBL_DIG + 1)
		1262	z *= 10;
		1263	if (nd++ < 9)
		1264	y = 10*y + c;
		1265	else if (nd <= DBL_DIG + 1)
		1266	z = 10*z + c;
		1267	nz = 0;
		1268	}
		1269	}
		1270	}
		1271	dig_done:
		1272	e = 0;
		1273	if (c == 'e' \|\| c == 'E') {
		1274	if (!nd && !nz && !nz0) {
		1275	s = s00;
		1276	goto ret;
		1277	}
		1278	s00 = s;
		1279	esign = 0;
		1280	switch(c = *++s) {
		1281	case '-':
		1282	esign = 1;
		1283	case '+':
		1284	c = *++s;
		1285	}
		1286	if (c >= '0' && c <= '9') {
		1287	while (c == '0')
		1288	c = *++s;
		1289	if (c > '0' && c <= '9') {
		1290	L = c - '0';
		1291	s1 = s;
		1292	while ((c = *++s) >= '0' && c <= '9')
		1293	L = 10*L + c - '0';
		1294	if (s - s1 > 8 \|\| L > 19999)
		1295	/* Avoid confusion from exponents
		1296	* so large that e might overflow.
		1297	*/
		1298	e = 19999; /* safe for 16 bit ints */
		1299	else
		1300	e = (int)L;
		1301	if (esign)
		1302	e = -e;
		1303	} else
		1304	e = 0;
		1305	} else
		1306	s = s00;
		1307	}
		1308	if (!nd) {
		1309	if (!nz && !nz0)
		1310	s = s00;
		1311	goto ret;
		1312	}
		1313	e1 = e -= nf;
		1314
		1315	/* Now we have nd0 digits, starting at s0, followed by a
		1316	* decimal point, followed by nd-nd0 digits. The number we're
		1317	* after is the integer represented by those digits times
		1318	* 10*e /
		1319
		1320	if (!nd0)
		1321	nd0 = nd;
		1322	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
		1323	rv = y;
		1324	if (k > 9)
		1325	rv = tens[k - 9] * rv + z;
		1326	if (nd <= DBL_DIG
		1327	#ifndef RND_PRODQUOT
		1328	&& FLT_ROUNDS == 1
		1329	#endif
		1330	) {
		1331	if (!e)
		1332	goto ret;
		1333	if (e > 0) {
		1334	if (e <= Ten_pmax) {
		1335	#ifdef VAX
		1336	goto vax_ovfl_check;
		1337	#else
		1338	/* rv = */ rounded_product(rv, tens[e]);
		1339	goto ret;
		1340	#endif
		1341	}
		1342	i = DBL_DIG - nd;
		1343	if (e <= Ten_pmax + i) {
		1344	/* A fancier test would sometimes let us do
		1345	* this for larger i values.
		1346	*/
		1347	e -= i;
		1348	rv *= tens[i];
		1349	#ifdef VAX
		1350	/* VAX exponent range is so narrow we must
		1351	* worry about overflow here...
		1352	*/
		1353	vax_ovfl_check:
		1354	word0(rv) -= P*Exp_msk1;
		1355	/* rv = */ rounded_product(rv, tens[e]);
		1356	if ((word0(rv) & Exp_mask)
		1357	> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
		1358	goto ovfl;
		1359	word0(rv) += P*Exp_msk1;
		1360	#else
		1361	/* rv = */ rounded_product(rv, tens[e]);
		1362	#endif
		1363	goto ret;
		1364	}
		1365	}
		1366	#ifndef Inaccurate_Divide
		1367	else if (e >= -Ten_pmax) {
		1368	/* rv = */ rounded_quotient(rv, tens[-e]);
		1369	goto ret;
		1370	}
		1371	#endif
		1372	}
		1373	e1 += nd - k;
		1374
		1375	/* Get starting approximation = rv * 10*e1 /
		1376
		1377	if (e1 > 0) {
		1378	if ( (i = e1 & 15) )
		1379	rv *= tens[i];
		1380	if ( (e1 &= ~15) ) {
		1381	if (e1 > DBL_MAX_10_EXP) {
		1382	ovfl:
		1383	errno = ERANGE;
		1384	#ifdef __STDC__
		1385	rv = HUGE_VAL;
		1386	#else
		1387	/* Can't trust HUGE_VAL */
		1388	#ifdef IEEE_Arith
		1389	word0(rv) = Exp_mask;
		1390	word1(rv) = 0;
		1391	#else
		1392	word0(rv) = Big0;
		1393	word1(rv) = Big1;
		1394	#endif
		1395	#endif
		1396	goto ret;
		1397	}
		1398	if (e1 >>= 4) {
		1399	for (j = 0; e1 > 1; j++, e1 >>= 1)
		1400	if (e1 & 1)
		1401	rv *= bigtens[j];
		1402	/* The last multiplication could overflow. */
		1403	word0(rv) -= P*Exp_msk1;
		1404	rv *= bigtens[j];
		1405	if ((z = word0(rv) & Exp_mask)
		1406	> Exp_msk1*(DBL_MAX_EXP+Bias-P))
		1407	goto ovfl;
		1408	if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
		1409	/* set to largest number */
		1410	/* (Can't trust DBL_MAX) */
		1411	word0(rv) = Big0;
		1412	word1(rv) = Big1;
		1413	}
		1414	else
		1415	word0(rv) += P*Exp_msk1;
		1416	}
		1417	}
		1418	} else if (e1 < 0) {
		1419	e1 = -e1;
		1420	if ( (i = e1 & 15) )
		1421	rv /= tens[i];
		1422	if ( (e1 &= ~15) ) {
		1423	e1 >>= 4;
		1424	for (j = 0; e1 > 1; j++, e1 >>= 1)
		1425	if (e1 & 1)
		1426	rv *= tinytens[j];
		1427	/* The last multiplication could underflow. */
		1428	rv0 = rv;
		1429	rv *= tinytens[j];
		1430	if (!rv) {
		1431	rv = 2.*rv0;
		1432	rv *= tinytens[j];
		1433	if (!rv) {
		1434	undfl:
		1435	rv = 0.;
		1436	errno = ERANGE;
		1437	goto ret;
		1438	}
		1439	word0(rv) = Tiny0;
		1440	word1(rv) = Tiny1;
		1441	/* The refinement below will clean
		1442	* this approximation up.
		1443	*/
		1444	}
		1445	}
		1446	}
		1447
		1448	/* Now the hard part -- adjusting rv to the correct value.*/
		1449
		1450	/* Put digits into bd: true value = bd * 10^e */
		1451
		1452	bd0 = s2b(s0, nd0, nd, y);
		1453
		1454	for (;;) {
		1455	bd = Balloc(bd0->k);
		1456	Bcopy(bd, bd0);
		1457	bb = d2b(rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
		1458	bs = i2b(1);
		1459
		1460	if (e >= 0) {
		1461	bb2 = bb5 = 0;
		1462	bd2 = bd5 = e;
		1463	} else {
		1464	bb2 = bb5 = -e;
		1465	bd2 = bd5 = 0;
		1466	}
		1467	if (bbe >= 0)
		1468	bb2 += bbe;
		1469	else
		1470	bd2 -= bbe;
		1471	bs2 = bb2;
		1472	#ifdef Sudden_Underflow
		1473	#ifdef IBM
		1474	j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
		1475	#else
		1476	j = P + 1 - bbbits;
		1477	#endif
		1478	#else
		1479	i = bbe + bbbits - 1; /* logb(rv) */
		1480	if (i < Emin) /* denormal */
		1481	j = bbe + (P-Emin);
		1482	else
		1483	j = P + 1 - bbbits;
		1484	#endif
		1485	bb2 += j;
		1486	bd2 += j;
		1487	i = bb2 < bd2 ? bb2 : bd2;
		1488	if (i > bs2)
		1489	i = bs2;
		1490	if (i > 0) {
		1491	bb2 -= i;
		1492	bd2 -= i;
		1493	bs2 -= i;
		1494	}
		1495	if (bb5 > 0) {
		1496	bs = pow5mult(bs, bb5);
		1497	bb1 = mult(bs, bb);
		1498	Bfree(bb);
		1499	bb = bb1;
		1500	}
		1501	if (bb2 > 0)
		1502	bb = lshift(bb, bb2);
		1503	if (bd5 > 0)
		1504	bd = pow5mult(bd, bd5);
		1505	if (bd2 > 0)
		1506	bd = lshift(bd, bd2);
		1507	if (bs2 > 0)
		1508	bs = lshift(bs, bs2);
		1509	delta = diff(bb, bd);
		1510	dsign = delta->sign;
		1511	delta->sign = 0;
		1512	i = cmp(delta, bs);
		1513	if (i < 0) {
		1514	/* Error is less than half an ulp -- check for
		1515	* special case of mantissa a power of two.
		1516	*/
		1517	if (dsign \|\| word1(rv) \|\| word0(rv) & Bndry_mask)
		1518	break;
		1519	delta = lshift(delta,Log2P);
		1520	if (cmp(delta, bs) > 0)
		1521	goto drop_down;
		1522	break;
		1523	}
		1524	if (i == 0) {
		1525	/* exactly half-way between */
		1526	if (dsign) {
		1527	if ((word0(rv) & Bndry_mask1) == Bndry_mask1
		1528	&& word1(rv) == 0xffffffff) {
		1529	/boundary case -- increment exponent/
		1530	word0(rv) = (word0(rv) & Exp_mask)
		1531	+ Exp_msk1
		1532	#ifdef IBM
		1533	\| Exp_msk1 >> 4
		1534	#endif
		1535	;
		1536	word1(rv) = 0;
		1537	break;
		1538	}
		1539	} else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
		1540	drop_down:
		1541	/* boundary case -- decrement exponent */
		1542	#ifdef Sudden_Underflow
		1543	L = word0(rv) & Exp_mask;
		1544	#ifdef IBM
		1545	if (L < Exp_msk1)
		1546	#else
		1547	if (L <= Exp_msk1)
		1548	#endif
		1549	goto undfl;
		1550	L -= Exp_msk1;
		1551	#else
		1552	L = (word0(rv) & Exp_mask) - Exp_msk1;
		1553	#endif
		1554	word0(rv) = L \| Bndry_mask1;
		1555	word1(rv) = 0xffffffff;
		1556	#ifdef IBM
		1557	goto cont;
		1558	#else
		1559	break;
		1560	#endif
		1561	}
		1562	#ifndef ROUND_BIASED
		1563	if (!(word1(rv) & LSB))
		1564	break;
		1565	#endif
		1566	if (dsign)
		1567	rv += ulp(rv);
		1568	#ifndef ROUND_BIASED
		1569	else {
		1570	rv -= ulp(rv);
		1571	#ifndef Sudden_Underflow
		1572	if (!rv)
		1573	goto undfl;
		1574	#endif
		1575	}
		1576	#endif
		1577	break;
		1578	}
		1579	if ((aadj = ratio(delta, bs)) <= 2.) {
		1580	if (dsign)
		1581	aadj = aadj1 = 1.;
		1582	else if (word1(rv) \|\| word0(rv) & Bndry_mask) {
		1583	#ifndef Sudden_Underflow
		1584	if (word1(rv) == Tiny1 && !word0(rv))
		1585	goto undfl;
		1586	#endif
		1587	aadj = 1.;
		1588	aadj1 = -1.;
		1589	} else {
		1590	/* special case -- power of FLT_RADIX to be */
		1591	/* rounded down... */
		1592
		1593	if (aadj < 2./FLT_RADIX)
		1594	aadj = 1./FLT_RADIX;
		1595	else
		1596	aadj *= 0.5;
		1597	aadj1 = -aadj;
		1598	}
		1599	} else {
		1600	aadj *= 0.5;
		1601	aadj1 = dsign ? aadj : -aadj;
		1602	#ifdef Check_FLT_ROUNDS
		1603	switch(FLT_ROUNDS) {
		1604	case 2: /* towards +infinity */
		1605	aadj1 -= 0.5;
		1606	break;
		1607	case 0: /* towards 0 */
		1608	case 3: /* towards -infinity */
		1609	aadj1 += 0.5;
		1610	}
		1611	#else
		1612	if (FLT_ROUNDS == 0)
		1613	aadj1 += 0.5;
		1614	#endif
		1615	}
		1616	y = word0(rv) & Exp_mask;
		1617
		1618	/* Check for overflow */
		1619
		1620	if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
		1621	rv0 = rv;
		1622	word0(rv) -= P*Exp_msk1;
		1623	adj = aadj1 * ulp(rv);
		1624	rv += adj;
		1625	if ((word0(rv) & Exp_mask) >=
		1626	Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
		1627	if (word0(rv0) == Big0 && word1(rv0) == Big1)
		1628	goto ovfl;
		1629	word0(rv) = Big0;
		1630	word1(rv) = Big1;
		1631	goto cont;
		1632	} else
		1633	word0(rv) += P*Exp_msk1;
		1634	} else {
		1635	#ifdef Sudden_Underflow
		1636	if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
		1637	rv0 = rv;
		1638	word0(rv) += P*Exp_msk1;
		1639	adj = aadj1 * ulp(rv);
		1640	rv += adj;
		1641	#ifdef IBM
		1642	if ((word0(rv) & Exp_mask) < P*Exp_msk1)
		1643	#else
		1644	if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
		1645	#endif
		1646	{
		1647	if (word0(rv0) == Tiny0
		1648	&& word1(rv0) == Tiny1)
		1649	goto undfl;
		1650	word0(rv) = Tiny0;
		1651	word1(rv) = Tiny1;
		1652	goto cont;
		1653	} else
		1654	word0(rv) -= P*Exp_msk1;
		1655	} else {
		1656	adj = aadj1 * ulp(rv);
		1657	rv += adj;
		1658	}
		1659	#else
		1660	/* Compute adj so that the IEEE rounding rules will
		1661	* correctly round rv + adj in some half-way cases.
		1662	* If rv * ulp(rv) is denormalized (i.e.,
		1663	* y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
		1664	* trouble from bits lost to denormalization;
		1665	* example: 1.2e-307 .
		1666	*/
		1667	if (y <= (P-1)*Exp_msk1 && aadj >= 1.) {
		1668	aadj1 = (double)(int)(aadj + 0.5);
		1669	if (!dsign)
		1670	aadj1 = -aadj1;
		1671	}
		1672	adj = aadj1 * ulp(rv);
		1673	rv += adj;
		1674	#endif
		1675	}
		1676	z = word0(rv) & Exp_mask;
		1677	if (y == z) {
		1678	/* Can we stop now? */
		1679	L = aadj;
		1680	aadj -= L;
		1681	/* The tolerances below are conservative. */
		1682	if (dsign \|\| word1(rv) \|\| word0(rv) & Bndry_mask) {
		1683	if (aadj < .4999999 \|\| aadj > .5000001)
		1684	break;
		1685	} else if (aadj < .4999999/FLT_RADIX)
		1686	break;
		1687	}
		1688	cont:
		1689	Bfree(bb);
		1690	Bfree(bd);
		1691	Bfree(bs);
		1692	Bfree(delta);
		1693	}
		1694	Bfree(bb);
		1695	Bfree(bd);
		1696	Bfree(bs);
		1697	Bfree(bd0);
		1698	Bfree(delta);
		1699	ret:
		1700	if (se)
		1701	se = (char )s;
		1702	return sign ? -rv : rv;
		1703	}
		1704
		1705	static int
		1706	quorem
		1707	#ifdef KR_headers
		1708	(b, S) Bigint b, S;
		1709	#else
		1710	(Bigint b, Bigint S)
		1711	#endif
		1712	{
		1713	int n;
		1714	long borrow, y;
		1715	unsigned long carry, q, ys;
		1716	unsigned long bx, bxe, sx, sxe;
		1717	#ifdef Pack_32
		1718	long z;
		1719	unsigned long si, zs;
		1720	#endif
		1721
		1722	n = S->wds;
		1723	#ifdef DEBUG
		1724	/debug/ if (b->wds > n)
		1725	/debug/ Bug("oversize b in quorem");
		1726	#endif
		1727	if (b->wds < n)
		1728	return 0;
		1729	sx = S->x;
		1730	sxe = sx + --n;
		1731	bx = b->x;
		1732	bxe = bx + n;
		1733	q = bxe / (sxe + 1); /* ensure q <= true quotient */
		1734	#ifdef DEBUG
		1735	/debug/ if (q > 9)
		1736	/debug/ Bug("oversized quotient in quorem");
		1737	#endif
		1738	if (q) {
		1739	borrow = 0;
		1740	carry = 0;
		1741	do {
		1742	#ifdef Pack_32
		1743	si = *sx++;
		1744	ys = (si & 0xffff) * q + carry;
		1745	zs = (si >> 16) * q + (ys >> 16);
		1746	carry = zs >> 16;
		1747	y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
		1748	borrow = y >> 16;
		1749	Sign_Extend(borrow, y);
		1750	z = (*bx >> 16) - (zs & 0xffff) + borrow;
		1751	borrow = z >> 16;
		1752	Sign_Extend(borrow, z);
		1753	Storeinc(bx, z, y);
		1754	#else
		1755	ys = sx++ q + carry;
		1756	carry = ys >> 16;
		1757	y = *bx - (ys & 0xffff) + borrow;
		1758	borrow = y >> 16;
		1759	Sign_Extend(borrow, y);
		1760	*bx++ = y & 0xffff;
		1761	#endif
		1762	} while (sx <= sxe);
		1763	if (!*bxe) {
		1764	bx = b->x;
		1765	while (--bxe > bx && !*bxe)
		1766	--n;
		1767	b->wds = n;
		1768	}
		1769	}
		1770	if (cmp(b, S) >= 0) {
		1771	q++;
		1772	borrow = 0;
		1773	carry = 0;
		1774	bx = b->x;
		1775	sx = S->x;
		1776	do {
		1777	#ifdef Pack_32
		1778	si = *sx++;
		1779	ys = (si & 0xffff) + carry;
		1780	zs = (si >> 16) + (ys >> 16);
		1781	carry = zs >> 16;
		1782	y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
		1783	borrow = y >> 16;
		1784	Sign_Extend(borrow, y);
		1785	z = (*bx >> 16) - (zs & 0xffff) + borrow;
		1786	borrow = z >> 16;
		1787	Sign_Extend(borrow, z);
		1788	Storeinc(bx, z, y);
		1789	#else
		1790	ys = *sx++ + carry;
		1791	carry = ys >> 16;
		1792	y = *bx - (ys & 0xffff) + borrow;
		1793	borrow = y >> 16;
		1794	Sign_Extend(borrow, y);
		1795	*bx++ = y & 0xffff;
		1796	#endif
		1797	} while (sx <= sxe);
		1798	bx = b->x;
		1799	bxe = bx + n;
		1800	if (!*bxe) {
		1801	while (--bxe > bx && !*bxe)
		1802	--n;
		1803	b->wds = n;
		1804	}
		1805	}
		1806	return q;
		1807	}
		1808
		1809	/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
		1810	*
		1811	* Inspired by "How to Print Floating-Point Numbers Accurately" by
		1812	* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
		1813	*
		1814	* Modifications:
		1815	* 1. Rather than iterating, we use a simple numeric overestimate
		1816	* to determine k = floor(log10(d)). We scale relevant
		1817	* quantities using O(log2(k)) rather than O(k) multiplications.
		1818	* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
		1819	* try to generate digits strictly left to right. Instead, we
		1820	* compute with fewer bits and propagate the carry if necessary
		1821	* when rounding the final digit up. This is often faster.
		1822	* 3. Under the assumption that input will be rounded nearest,
		1823	* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
		1824	* That is, we allow equality in stopping tests when the
		1825	* round-nearest rule will give the same floating-point value
		1826	* as would satisfaction of the stopping test with strict
		1827	* inequality.
		1828	* 4. We remove common factors of powers of 2 from relevant
		1829	* quantities.
		1830	* 5. When converting floating-point integers less than 1e16,
		1831	* we use floating-point arithmetic rather than resorting
		1832	* to multiple-precision integers.
		1833	* 6. When asked to produce fewer than 15 digits, we first try
		1834	* to get by with floating-point arithmetic; we resort to
		1835	* multiple-precision integer arithmetic only if we cannot
		1836	* guarantee that the floating-point calculation has given
		1837	* the correctly rounded result. For k requested digits and
		1838	* "uniformly" distributed input, the probability is
		1839	* something like 10^(k-15) that we must resort to the long
		1840	* calculation.
		1841	*/
		1842
		1843	char *
		1844	__dtoa
		1845	#ifdef KR_headers
		1846	(d, mode, ndigits, decpt, sign, rve)
		1847	double d; int mode, ndigits, decpt, sign; char **rve;
		1848	#else
		1849	(double d, int mode, int ndigits, int decpt, int sign, char **rve)
		1850	#endif
		1851	{
		1852	/* Arguments ndigits, decpt, sign are similar to those
		1853	of ecvt and fcvt; trailing zeros are suppressed from
		1854	the returned string. If not null, *rve is set to point
		1855	to the end of the return value. If d is +-Infinity or NaN,
		1856	then *decpt is set to 9999.
		1857
		1858	mode:
		1859
		1860	and rounded to nearest.
		1861	1 ==> like 0, but with Steele & White stopping rule;
		1862	e.g. with IEEE P754 arithmetic , mode 0 gives
		1863	1e23 whereas mode 1 gives 9.999999999999999e22.
		1864	2 ==> max(1,ndigits) significant digits. This gives a
		1865	return value similar to that of ecvt, except
		1866	that trailing zeros are suppressed.
		1867	3 ==> through ndigits past the decimal point. This
		1868	gives a return value similar to that from fcvt,
		1869	except that trailing zeros are suppressed, and
		1870	ndigits can be negative.
		1871	4-9 should give the same return values as 2-3, i.e.,
		1872	4 <= mode <= 9 ==> same return as mode
		1873	2 + (mode & 1). These modes are mainly for
		1874	debugging; often they run slower but sometimes
		1875	faster than modes 2-3.
		1876	4,5,8,9 ==> left-to-right digit generation.
		1877	6-9 ==> don't try fast floating-point estimate
		1878	(if applicable).
		1879
		1880	Values of mode other than 0-9 are treated as mode 0.
		1881
		1882	Sufficient space is allocated to the return value
		1883	to hold the suppressed trailing zeros.
		1884	*/
		1885
		1886	int bbits, b2, b5, be, dig, i, ieps, ilim=0, ilim0=0, ilim1=0,
		1887	j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
		1888	spec_case=0, try_quick;
		1889	long L;
		1890	#ifndef Sudden_Underflow
		1891	int denorm;
		1892	unsigned long x;
		1893	#endif
		1894	Bigint b, b1, delta, mlo=NULL, mhi, S;
		1895	double d2, ds, eps;
		1896	char s, s0;
		1897	static Bigint *result;
		1898	static int result_k;
		1899
		1900	if (result) {
		1901	result->k = result_k;
		1902	result->maxwds = 1 << result_k;
		1903	Bfree(result);
		1904	result = 0;
		1905	}
		1906
		1907	if (word0(d) & Sign_bit) {
		1908	/* set sign for everything, including 0's and NaNs */
		1909	*sign = 1;
		1910	word0(d) &= ~Sign_bit; /* clear sign bit */
		1911	}
		1912	else
		1913	*sign = 0;
		1914
		1915	#if defined(IEEE_Arith) + defined(VAX)
		1916	#ifdef IEEE_Arith
		1917	if ((word0(d) & Exp_mask) == Exp_mask)
		1918	#else
		1919	if (word0(d) == 0x8000)
		1920	#endif
		1921	{
		1922	/* Infinity or NaN */
		1923	*decpt = 9999;
		1924	s =
		1925	#ifdef IEEE_Arith
		1926	!word1(d) && !(word0(d) & 0xfffff) ? "Infinity" :
		1927	#endif
		1928	"NaN";
		1929	if (rve)
		1930	*rve =
		1931	#ifdef IEEE_Arith
		1932	s[3] ? s + 8 :
		1933	#endif
		1934	s + 3;
		1935	return s;
		1936	}
		1937	#endif
		1938	#ifdef IBM
		1939	d += 0; /* normalize */
		1940	#endif
		1941	if (!d) {
		1942	*decpt = 1;
		1943	s = "0";
		1944	if (rve)
		1945	*rve = s + 1;
		1946	return s;
		1947	}
		1948
		1949	b = d2b(d, &be, &bbits);
		1950	#ifdef Sudden_Underflow
		1951	i = (int)(word0(d) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
		1952	#else
		1953	if ( (i = (int)((word0(d) >> Exp_shift1) & (Exp_mask>>Exp_shift1))) ) {
		1954	#endif
		1955	d2 = d;
		1956	word0(d2) &= Frac_mask1;
		1957	word0(d2) \|= Exp_11;
		1958	#ifdef IBM
		1959	if ( (j = 11 - hi0bits(word0(d2) & Frac_mask)) )
		1960	d2 /= 1 << j;
		1961	#endif
		1962
		1963	/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
		1964	* log10(x) = log(x) / log(10)
		1965	* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
		1966	* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
		1967	*
		1968	* This suggests computing an approximation k to log10(d) by
		1969	*
		1970	* k = (i - Bias)*0.301029995663981
		1971	* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
		1972	*
		1973	* We want k to be too large rather than too small.
		1974	* The error in the first-order Taylor series approximation
		1975	* is in our favor, so we just round up the constant enough
		1976	* to compensate for any error in the multiplication of
		1977	* (i - Bias) by 0.301029995663981; since \|i - Bias\| <= 1077,
		1978	* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
		1979	* adding 1e-13 to the constant term more than suffices.
		1980	* Hence we adjust the constant term to 0.1760912590558.
		1981	* (We could get a more accurate k by invoking log10,
		1982	* but this is probably not worthwhile.)
		1983	*/
		1984
		1985	i -= Bias;
		1986	#ifdef IBM
		1987	i <<= 2;
		1988	i += j;
		1989	#endif
		1990	#ifndef Sudden_Underflow
		1991	denorm = 0;
		1992	} else {
		1993	/* d is denormalized */
		1994
		1995	i = bbits + be + (Bias + (P-1) - 1);
		1996	x = i > 32 ? ((word0(d) << (64 - i)) \| (word1(d) >> (i - 32)))
		1997	: (word1(d) << (32 - i));
		1998	d2 = x;
		1999	word0(d2) -= 31Exp_msk1; / adjust exponent */
		2000	i -= (Bias + (P-1) - 1) + 1;
		2001	denorm = 1;
		2002	}
		2003	#endif
		2004	ds = (d2-1.5)0.289529654602168 + 0.1760912590558 + i0.301029995663981;
		2005	k = (int)ds;
		2006	if (ds < 0. && ds != k)
		2007	k--; /* want k = floor(ds) */
		2008	k_check = 1;
		2009	if (k >= 0 && k <= Ten_pmax) {
		2010	if (d < tens[k])
		2011	k--;
		2012	k_check = 0;
		2013	}
		2014	j = bbits - i - 1;
		2015	if (j >= 0) {
		2016	b2 = 0;
		2017	s2 = j;
		2018	} else {
		2019	b2 = -j;
		2020	s2 = 0;
		2021	}
		2022	if (k >= 0) {
		2023	b5 = 0;
		2024	s5 = k;
		2025	s2 += k;
		2026	} else {
		2027	b2 -= k;
		2028	b5 = -k;
		2029	s5 = 0;
		2030	}
		2031	if (mode < 0 \|\| mode > 9)
		2032	mode = 0;
		2033	try_quick = 1;
		2034	if (mode > 5) {
		2035	mode -= 4;
		2036	try_quick = 0;
		2037	}
		2038	leftright = 1;
		2039	switch(mode) {
		2040	case 0:
		2041	case 1:
		2042	ilim = ilim1 = -1;
		2043	i = 18;
		2044	ndigits = 0;
		2045	break;
		2046	case 2:
		2047	leftright = 0;
		2048	/* no break */
		2049	case 4:
		2050	if (ndigits <= 0)
		2051	ndigits = 1;
		2052	ilim = ilim1 = i = ndigits;
		2053	break;
		2054	case 3:
		2055	leftright = 0;
		2056	/* no break */
		2057	case 5:
		2058	i = ndigits + k + 1;
		2059	ilim = i;
		2060	ilim1 = i - 1;
		2061	if (i <= 0)
		2062	i = 1;
		2063	}
		2064	j = sizeof(unsigned long);
		2065	for (result_k = 0; sizeof(Bigint) - sizeof(unsigned long) + j < i;
		2066	j <<= 1) result_k++;
		2067	result = Balloc(result_k);
		2068	s = s0 = (char *)result;
		2069
		2070	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
		2071
		2072	/* Try to get by with floating-point arithmetic. */
		2073
		2074	i = 0;
		2075	d2 = d;
		2076	k0 = k;
		2077	ilim0 = ilim;
		2078	ieps = 2; /* conservative */
		2079	if (k > 0) {
		2080	ds = tens[k&0xf];
		2081	j = k >> 4;
		2082	if (j & Bletch) {
		2083	/* prevent overflows */
		2084	j &= Bletch - 1;
		2085	d /= bigtens[n_bigtens-1];
		2086	ieps++;
		2087	}
		2088	for (; j; j >>= 1, i++)
		2089	if (j & 1) {
		2090	ieps++;
		2091	ds *= bigtens[i];
		2092	}
		2093	d /= ds;
		2094	} else if ( (j1 = -k) ) {
		2095	d *= tens[j1 & 0xf];
		2096	for (j = j1 >> 4; j; j >>= 1, i++)
		2097	if (j & 1) {
		2098	ieps++;
		2099	d *= bigtens[i];
		2100	}
		2101	}
		2102	if (k_check && d < 1. && ilim > 0) {
		2103	if (ilim1 <= 0)
		2104	goto fast_failed;
		2105	ilim = ilim1;
		2106	k--;
		2107	d *= 10.;
		2108	ieps++;
		2109	}
		2110	eps = ieps*d + 7.;
		2111	word0(eps) -= (P-1)*Exp_msk1;
		2112	if (ilim == 0) {
		2113	S = mhi = 0;
		2114	d -= 5.;
		2115	if (d > eps)
		2116	goto one_digit;
		2117	if (d < -eps)
		2118	goto no_digits;
		2119	goto fast_failed;
		2120	}
		2121	#ifndef No_leftright
		2122	if (leftright) {
		2123	/* Use Steele & White method of only
		2124	* generating digits needed.
		2125	*/
		2126	eps = 0.5/tens[ilim-1] - eps;
		2127	for (i = 0;;) {
		2128	L = d;
		2129	d -= L;
		2130	*s++ = '0' + (int)L;
		2131	if (d < eps)
		2132	goto ret1;
		2133	if (1. - d < eps)
		2134	goto bump_up;
		2135	if (++i >= ilim)
		2136	break;
		2137	eps *= 10.;
		2138	d *= 10.;
		2139	}
		2140	} else {
		2141	#endif
		2142	/* Generate ilim digits, then fix them up. */
		2143	eps *= tens[ilim-1];
		2144	for (i = 1;; i++, d *= 10.) {
		2145	L = d;
		2146	d -= L;
		2147	*s++ = '0' + (int)L;
		2148	if (i == ilim) {
		2149	if (d > 0.5 + eps)
		2150	goto bump_up;
		2151	else if (d < 0.5 - eps) {
		2152	while (*--s == '0');
		2153	s++;
		2154	goto ret1;
		2155	}
		2156	break;
		2157	}
		2158	}
		2159	#ifndef No_leftright
		2160	}
		2161	#endif
		2162	fast_failed:
		2163	s = s0;
		2164	d = d2;
		2165	k = k0;
		2166	ilim = ilim0;
		2167	}
		2168
		2169	/* Do we have a "small" integer? */
		2170
		2171	if (be >= 0 && k <= Int_max) {
		2172	/* Yes. */
		2173	ds = tens[k];
		2174	if (ndigits < 0 && ilim <= 0) {
		2175	S = mhi = 0;
		2176	if (ilim < 0 \|\| d <= 5*ds)
		2177	goto no_digits;
		2178	goto one_digit;
		2179	}
		2180	for (i = 1;; i++) {
		2181	L = d / ds;
		2182	d -= L*ds;
		2183	#ifdef Check_FLT_ROUNDS
		2184	/* If FLT_ROUNDS == 2, L will usually be high by 1 */
		2185	if (d < 0) {
		2186	L--;
		2187	d += ds;
		2188	}
		2189	#endif
		2190	*s++ = '0' + (int)L;
		2191	if (i == ilim) {
		2192	d += d;
		2193	if (d > ds \|\| (d == ds && L & 1)) {
		2194	bump_up:
		2195	while (*--s == '9')
		2196	if (s == s0) {
		2197	k++;
		2198	*s = '0';
		2199	break;
		2200	}
		2201	++*s++;
		2202	}
		2203	break;
		2204	}
		2205	if (!(d *= 10.))
		2206	break;
		2207	}
		2208	goto ret1;
		2209	}
		2210
		2211	m2 = b2;
		2212	m5 = b5;
		2213	mhi = mlo = 0;
		2214	if (leftright) {
		2215	if (mode < 2) {
		2216	i =
		2217	#ifndef Sudden_Underflow
		2218	denorm ? be + (Bias + (P-1) - 1 + 1) :
		2219	#endif
		2220	#ifdef IBM
		2221	1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
		2222	#else
		2223	1 + P - bbits;
		2224	#endif
		2225	} else {
		2226	j = ilim - 1;
		2227	if (m5 >= j)
		2228	m5 -= j;
		2229	else {
		2230	s5 += j -= m5;
		2231	b5 += j;
		2232	m5 = 0;
		2233	}
		2234	if ((i = ilim) < 0) {
		2235	m2 -= i;
		2236	i = 0;
		2237	}
		2238	}
		2239	b2 += i;
		2240	s2 += i;
		2241	mhi = i2b(1);
		2242	}
		2243	if (m2 > 0 && s2 > 0) {
		2244	i = m2 < s2 ? m2 : s2;
		2245	b2 -= i;
		2246	m2 -= i;
		2247	s2 -= i;
		2248	}
		2249	if (b5 > 0) {
		2250	if (leftright) {
		2251	if (m5 > 0) {
		2252	mhi = pow5mult(mhi, m5);
		2253	b1 = mult(mhi, b);
		2254	Bfree(b);
		2255	b = b1;
		2256	}
		2257	if ( (j = b5 - m5) )
		2258	b = pow5mult(b, j);
		2259	} else
		2260	b = pow5mult(b, b5);
		2261	}
		2262	S = i2b(1);
		2263	if (s5 > 0)
		2264	S = pow5mult(S, s5);
		2265
		2266	/* Check for special case that d is a normalized power of 2. */
		2267
		2268	if (mode < 2) {
		2269	if (!word1(d) && !(word0(d) & Bndry_mask)
		2270	#ifndef Sudden_Underflow
		2271	&& word0(d) & Exp_mask
		2272	#endif
		2273	) {
		2274	/* The special case */
		2275	b2 += Log2P;
		2276	s2 += Log2P;
		2277	spec_case = 1;
		2278	} else
		2279	spec_case = 0;
		2280	}
		2281
		2282	/* Arrange for convenient computation of quotients:
		2283	* shift left if necessary so divisor has 4 leading 0 bits.
		2284	*
		2285	* Perhaps we should just compute leading 28 bits of S once
		2286	* and for all and pass them and a shift to quorem, so it
		2287	* can do shifts and ors to compute the numerator for q.
		2288	*/
		2289	#ifdef Pack_32
		2290	if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f) )
		2291	i = 32 - i;
		2292	#else
		2293	if ( (i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0xf) )
		2294	i = 16 - i;
		2295	#endif
		2296	if (i > 4) {
		2297	i -= 4;
		2298	b2 += i;
		2299	m2 += i;
		2300	s2 += i;
		2301	} else if (i < 4) {
		2302	i += 28;
		2303	b2 += i;
		2304	m2 += i;
		2305	s2 += i;
		2306	}
		2307	if (b2 > 0)
		2308	b = lshift(b, b2);
		2309	if (s2 > 0)
		2310	S = lshift(S, s2);
		2311	if (k_check) {
		2312	if (cmp(b,S) < 0) {
		2313	k--;
		2314	b = multadd(b, 10, 0); /* we botched the k estimate */
		2315	if (leftright)
		2316	mhi = multadd(mhi, 10, 0);
		2317	ilim = ilim1;
		2318	}
		2319	}
		2320	if (ilim <= 0 && mode > 2) {
		2321	if (ilim < 0 \|\| cmp(b,S = multadd(S,5,0)) <= 0) {
		2322	/* no digits, fcvt style */
		2323	no_digits:
		2324	k = -1 - ndigits;
		2325	goto ret;
		2326	}
		2327	one_digit:
		2328	*s++ = '1';
		2329	k++;
		2330	goto ret;
		2331	}
		2332	if (leftright) {
		2333	if (m2 > 0)
		2334	mhi = lshift(mhi, m2);
		2335
		2336	/* Compute mlo -- check for special case
		2337	* that d is a normalized power of 2.
		2338	*/
		2339
		2340	mlo = mhi;
		2341	if (spec_case) {
		2342	mhi = Balloc(mhi->k);
		2343	Bcopy(mhi, mlo);
		2344	mhi = lshift(mhi, Log2P);
		2345	}
		2346
		2347	for (i = 1;;i++) {
		2348	dig = quorem(b,S) + '0';
		2349	/* Do we yet have the shortest decimal string
		2350	* that will round to d?
		2351	*/
		2352	j = cmp(b, mlo);
		2353	delta = diff(S, mhi);
		2354	j1 = delta->sign ? 1 : cmp(b, delta);
		2355	Bfree(delta);
		2356	#ifndef ROUND_BIASED
		2357	if (j1 == 0 && !mode && !(word1(d) & 1)) {
		2358	if (dig == '9')
		2359	goto round_9_up;
		2360	if (j > 0)
		2361	dig++;
		2362	*s++ = dig;
		2363	goto ret;
		2364	}
		2365	#endif
		2366	if (j < 0 \|\| (j == 0 && !mode
		2367	#ifndef ROUND_BIASED
		2368	&& !(word1(d) & 1)
		2369	#endif
		2370	)) {
		2371	if (j1 > 0) {
		2372	b = lshift(b, 1);
		2373	j1 = cmp(b, S);
		2374	if ((j1 > 0 \|\| (j1 == 0 && dig & 1))
		2375	&& dig++ == '9')
		2376	goto round_9_up;
		2377	}
		2378	*s++ = dig;
		2379	goto ret;
		2380	}
		2381	if (j1 > 0) {
		2382	if (dig == '9') { /* possible if i == 1 */
		2383	round_9_up:
		2384	*s++ = '9';
		2385	goto roundoff;
		2386	}
		2387	*s++ = dig + 1;
		2388	goto ret;
		2389	}
		2390	*s++ = dig;
		2391	if (i == ilim)
		2392	break;
		2393	b = multadd(b, 10, 0);
		2394	if (mlo == mhi)
		2395	mlo = mhi = multadd(mhi, 10, 0);
		2396	else {
		2397	mlo = multadd(mlo, 10, 0);
		2398	mhi = multadd(mhi, 10, 0);
		2399	}
		2400	}
		2401	} else
		2402	for (i = 1;; i++) {
		2403	*s++ = dig = quorem(b,S) + '0';
		2404	if (i >= ilim)
		2405	break;
		2406	b = multadd(b, 10, 0);
		2407	}
		2408
		2409	/* Round off last digit */
		2410
		2411	b = lshift(b, 1);
		2412	j = cmp(b, S);
		2413	if (j > 0 \|\| (j == 0 && dig & 1)) {
		2414	roundoff:
		2415	while (*--s == '9')
		2416	if (s == s0) {
		2417	k++;
		2418	*s++ = '1';
		2419	goto ret;
		2420	}
		2421	++*s++;
		2422	} else {
		2423	while (*--s == '0');
		2424	s++;
		2425	}
		2426	ret:
		2427	Bfree(S);
		2428	if (mhi) {
		2429	if (mlo && mlo != mhi)
		2430	Bfree(mlo);
		2431	Bfree(mhi);
		2432	}
		2433	ret1:
		2434	Bfree(b);
		2435	if (s == s0) { /* don't return empty string */
		2436	*s++ = '0';
		2437	k = 0;
		2438	}
		2439	*s = 0;
		2440	*decpt = k + 1;
		2441	if (rve)
		2442	*rve = s;
		2443	return s0;
		2444	}
		2445	#ifdef __cplusplus
		2446	}
		2447	#endif

Subversion Repositories shark

(root)/shark/branches/xen/libc/stdlib/strtod.c @ 1680 - Rev 1618