• file
• revisions
• annotate
• diff
• comparison
• raw

nsprpub/pr/src/misc/dtoa.c@6474c204b198 (annotated)

nsprpub/pr/src/misc/dtoa.c

Wed, 31 Dec 2014 06:09:35 +0100

author

Michael Schloh von Bennewitz <michael@schloh.com>

date

Wed, 31 Dec 2014 06:09:35 +0100

changeset 0

6474c204b198

permissions

-rw-r--r--

Cloned upstream origin tor-browser at tor-browser-31.3.0esr-4.5-1-build1
revision ID fc1c9ff7c1b2defdbc039f12214767608f46423f for hacking purpose.

michael@0	1	/****************************************************************
michael@0	2	*
michael@0	3	* The author of this software is David M. Gay.
michael@0	4	*
michael@0	5	* Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
michael@0	6	*
michael@0	7	* Permission to use, copy, modify, and distribute this software for any
michael@0	8	* purpose without fee is hereby granted, provided that this entire notice
michael@0	9	* is included in all copies of any software which is or includes a copy
michael@0	10	* or modification of this software and in all copies of the supporting
michael@0	11	* documentation for such software.
michael@0	12	*
michael@0	13	* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
michael@0	14	* WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
michael@0	15	* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
michael@0	16	* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
michael@0	17	*
michael@0	18	***************************************************************/
michael@0	19
michael@0	20	/* Please send bug reports to David M. Gay (dmg at acm dot org,
michael@0	21	* with " at " changed at "@" and " dot " changed to "."). */
michael@0	22
michael@0	23	/* On a machine with IEEE extended-precision registers, it is
michael@0	24	* necessary to specify double-precision (53-bit) rounding precision
michael@0	25	* before invoking strtod or dtoa. If the machine uses (the equivalent
michael@0	26	* of) Intel 80x87 arithmetic, the call
michael@0	27	* _control87(PC_53, MCW_PC);
michael@0	28	* does this with many compilers. Whether this or another call is
michael@0	29	* appropriate depends on the compiler; for this to work, it may be
michael@0	30	* necessary to #include "float.h" or another system-dependent header
michael@0	31	* file.
michael@0	32	*/
michael@0	33
michael@0	34	/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
michael@0	35	*
michael@0	36	* This strtod returns a nearest machine number to the input decimal
michael@0	37	* string (or sets errno to ERANGE). With IEEE arithmetic, ties are
michael@0	38	* broken by the IEEE round-even rule. Otherwise ties are broken by
michael@0	39	* biased rounding (add half and chop).
michael@0	40	*
michael@0	41	* Inspired loosely by William D. Clinger's paper "How to Read Floating
michael@0	42	* Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
michael@0	43	*
michael@0	44	* Modifications:
michael@0	45	*
michael@0	46	* 1. We only require IEEE, IBM, or VAX double-precision
michael@0	47	* arithmetic (not IEEE double-extended).
michael@0	48	* 2. We get by with floating-point arithmetic in a case that
michael@0	49	* Clinger missed -- when we're computing d * 10^n
michael@0	50	* for a small integer d and the integer n is not too
michael@0	51	* much larger than 22 (the maximum integer k for which
michael@0	52	* we can represent 10^k exactly), we may be able to
michael@0	53	* compute (d*10^k) * 10^(e-k) with just one roundoff.
michael@0	54	* 3. Rather than a bit-at-a-time adjustment of the binary
michael@0	55	* result in the hard case, we use floating-point
michael@0	56	* arithmetic to determine the adjustment to within
michael@0	57	* one bit; only in really hard cases do we need to
michael@0	58	* compute a second residual.
michael@0	59	* 4. Because of 3., we don't need a large table of powers of 10
michael@0	60	* for ten-to-e (just some small tables, e.g. of 10^k
michael@0	61	* for 0 <= k <= 22).
michael@0	62	*/
michael@0	63
michael@0	64	/*
michael@0	65	* #define IEEE_8087 for IEEE-arithmetic machines where the least
michael@0	66	* significant byte has the lowest address.
michael@0	67	* #define IEEE_MC68k for IEEE-arithmetic machines where the most
michael@0	68	* significant byte has the lowest address.
michael@0	69	* #define Long int on machines with 32-bit ints and 64-bit longs.
michael@0	70	* #define IBM for IBM mainframe-style floating-point arithmetic.
michael@0	71	* #define VAX for VAX-style floating-point arithmetic (D_floating).
michael@0	72	* #define No_leftright to omit left-right logic in fast floating-point
michael@0	73	* computation of dtoa. This will cause dtoa modes 4 and 5 to be
michael@0	74	* treated the same as modes 2 and 3 for some inputs.
michael@0	75	* #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
michael@0	76	* and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS
michael@0	77	* is also #defined, fegetround() will be queried for the rounding mode.
michael@0	78	* Note that both FLT_ROUNDS and fegetround() are specified by the C99
michael@0	79	* standard (and are specified to be consistent, with fesetround()
michael@0	80	* affecting the value of FLT_ROUNDS), but that some (Linux) systems
michael@0	81	* do not work correctly in this regard, so using fegetround() is more
michael@0	82	* portable than using FLT_ROUNDS directly.
michael@0	83	* #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
michael@0	84	* and Honor_FLT_ROUNDS is not #defined.
michael@0	85	* #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
michael@0	86	* that use extended-precision instructions to compute rounded
michael@0	87	* products and quotients) with IBM.
michael@0	88	* #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
michael@0	89	* that rounds toward +Infinity.
michael@0	90	* #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
michael@0	91	* rounding when the underlying floating-point arithmetic uses
michael@0	92	* unbiased rounding. This prevent using ordinary floating-point
michael@0	93	* arithmetic when the result could be computed with one rounding error.
michael@0	94	* #define Inaccurate_Divide for IEEE-format with correctly rounded
michael@0	95	* products but inaccurate quotients, e.g., for Intel i860.
michael@0	96	* #define NO_LONG_LONG on machines that do not have a "long long"
michael@0	97	* integer type (of >= 64 bits). On such machines, you can
michael@0	98	* #define Just_16 to store 16 bits per 32-bit Long when doing
michael@0	99	* high-precision integer arithmetic. Whether this speeds things
michael@0	100	* up or slows things down depends on the machine and the number
michael@0	101	* being converted. If long long is available and the name is
michael@0	102	* something other than "long long", #define Llong to be the name,
michael@0	103	* and if "unsigned Llong" does not work as an unsigned version of
michael@0	104	* Llong, #define #ULLong to be the corresponding unsigned type.
michael@0	105	* #define KR_headers for old-style C function headers.
michael@0	106	* #define Bad_float_h if your system lacks a float.h or if it does not
michael@0	107	* define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
michael@0	108	* FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
michael@0	109	* #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
michael@0	110	* if memory is available and otherwise does something you deem
michael@0	111	* appropriate. If MALLOC is undefined, malloc will be invoked
michael@0	112	* directly -- and assumed always to succeed. Similarly, if you
michael@0	113	* want something other than the system's free() to be called to
michael@0	114	* recycle memory acquired from MALLOC, #define FREE to be the
michael@0	115	* name of the alternate routine. (FREE or free is only called in
michael@0	116	* pathological cases, e.g., in a dtoa call after a dtoa return in
michael@0	117	* mode 3 with thousands of digits requested.)
michael@0	118	* #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
michael@0	119	* memory allocations from a private pool of memory when possible.
michael@0	120	* When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
michael@0	121	* unless #defined to be a different length. This default length
michael@0	122	* suffices to get rid of MALLOC calls except for unusual cases,
michael@0	123	* such as decimal-to-binary conversion of a very long string of
michael@0	124	* digits. The longest string dtoa can return is about 751 bytes
michael@0	125	* long. For conversions by strtod of strings of 800 digits and
michael@0	126	* all dtoa conversions in single-threaded executions with 8-byte
michael@0	127	* pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
michael@0	128	* pointers, PRIVATE_MEM >= 7112 appears adequate.
michael@0	129	* #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
michael@0	130	* #defined automatically on IEEE systems. On such systems,
michael@0	131	* when INFNAN_CHECK is #defined, strtod checks
michael@0	132	* for Infinity and NaN (case insensitively). On some systems
michael@0	133	* (e.g., some HP systems), it may be necessary to #define NAN_WORD0
michael@0	134	* appropriately -- to the most significant word of a quiet NaN.
michael@0	135	* (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
michael@0	136	* When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
michael@0	137	* strtod also accepts (case insensitively) strings of the form
michael@0	138	* NaN(x), where x is a string of hexadecimal digits and spaces;
michael@0	139	* if there is only one string of hexadecimal digits, it is taken
michael@0	140	* for the 52 fraction bits of the resulting NaN; if there are two
michael@0	141	* or more strings of hex digits, the first is for the high 20 bits,
michael@0	142	* the second and subsequent for the low 32 bits, with intervening
michael@0	143	* white space ignored; but if this results in none of the 52
michael@0	144	* fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
michael@0	145	* and NAN_WORD1 are used instead.
michael@0	146	* #define MULTIPLE_THREADS if the system offers preemptively scheduled
michael@0	147	* multiple threads. In this case, you must provide (or suitably
michael@0	148	* #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
michael@0	149	* by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
michael@0	150	* in pow5mult, ensures lazy evaluation of only one copy of high
michael@0	151	* powers of 5; omitting this lock would introduce a small
michael@0	152	* probability of wasting memory, but would otherwise be harmless.)
michael@0	153	* You must also invoke freedtoa(s) to free the value s returned by
michael@0	154	* dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
michael@0	155	* #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
michael@0	156	* avoids underflows on inputs whose result does not underflow.
michael@0	157	* If you #define NO_IEEE_Scale on a machine that uses IEEE-format
michael@0	158	* floating-point numbers and flushes underflows to zero rather
michael@0	159	* than implementing gradual underflow, then you must also #define
michael@0	160	* Sudden_Underflow.
michael@0	161	* #define USE_LOCALE to use the current locale's decimal_point value.
michael@0	162	* #define SET_INEXACT if IEEE arithmetic is being used and extra
michael@0	163	* computation should be done to set the inexact flag when the
michael@0	164	* result is inexact and avoid setting inexact when the result
michael@0	165	* is exact. In this case, dtoa.c must be compiled in
michael@0	166	* an environment, perhaps provided by #include "dtoa.c" in a
michael@0	167	* suitable wrapper, that defines two functions,
michael@0	168	* int get_inexact(void);
michael@0	169	* void clear_inexact(void);
michael@0	170	* such that get_inexact() returns a nonzero value if the
michael@0	171	* inexact bit is already set, and clear_inexact() sets the
michael@0	172	* inexact bit to 0. When SET_INEXACT is #defined, strtod
michael@0	173	* also does extra computations to set the underflow and overflow
michael@0	174	* flags when appropriate (i.e., when the result is tiny and
michael@0	175	* inexact or when it is a numeric value rounded to +-infinity).
michael@0	176	* #define NO_ERRNO if strtod should not assign errno = ERANGE when
michael@0	177	* the result overflows to +-Infinity or underflows to 0.
michael@0	178	* #define NO_HEX_FP to omit recognition of hexadecimal floating-point
michael@0	179	* values by strtod.
michael@0	180	* #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
michael@0	181	* to disable logic for "fast" testing of very long input strings
michael@0	182	* to strtod. This testing proceeds by initially truncating the
michael@0	183	* input string, then if necessary comparing the whole string with
michael@0	184	* a decimal expansion to decide close cases. This logic is only
michael@0	185	* used for input more than STRTOD_DIGLIM digits long (default 40).
michael@0	186	*/
michael@0	187
michael@0	188	#ifndef Long
michael@0	189	#define Long long
michael@0	190	#endif
michael@0	191	#ifndef ULong
michael@0	192	typedef unsigned Long ULong;
michael@0	193	#endif
michael@0	194
michael@0	195	#ifdef DEBUG
michael@0	196	#include "stdio.h"
michael@0	197	#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
michael@0	198	#endif
michael@0	199
michael@0	200	#include "stdlib.h"
michael@0	201	#include "string.h"
michael@0	202
michael@0	203	#ifdef USE_LOCALE
michael@0	204	#include "locale.h"
michael@0	205	#endif
michael@0	206
michael@0	207	#ifdef Honor_FLT_ROUNDS
michael@0	208	#ifndef Trust_FLT_ROUNDS
michael@0	209	#include <fenv.h>
michael@0	210	#endif
michael@0	211	#endif
michael@0	212
michael@0	213	#ifdef MALLOC
michael@0	214	#ifdef KR_headers
michael@0	215	extern char *MALLOC();
michael@0	216	#else
michael@0	217	extern void *MALLOC(size_t);
michael@0	218	#endif
michael@0	219	#else
michael@0	220	#define MALLOC malloc
michael@0	221	#endif
michael@0	222
michael@0	223	#ifndef Omit_Private_Memory
michael@0	224	#ifndef PRIVATE_MEM
michael@0	225	#define PRIVATE_MEM 2304
michael@0	226	#endif
michael@0	227	#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
michael@0	228	static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
michael@0	229	#endif
michael@0	230
michael@0	231	#undef IEEE_Arith
michael@0	232	#undef Avoid_Underflow
michael@0	233	#ifdef IEEE_MC68k
michael@0	234	#define IEEE_Arith
michael@0	235	#endif
michael@0	236	#ifdef IEEE_8087
michael@0	237	#define IEEE_Arith
michael@0	238	#endif
michael@0	239
michael@0	240	#ifdef IEEE_Arith
michael@0	241	#ifndef NO_INFNAN_CHECK
michael@0	242	#undef INFNAN_CHECK
michael@0	243	#define INFNAN_CHECK
michael@0	244	#endif
michael@0	245	#else
michael@0	246	#undef INFNAN_CHECK
michael@0	247	#define NO_STRTOD_BIGCOMP
michael@0	248	#endif
michael@0	249
michael@0	250	#include "errno.h"
michael@0	251
michael@0	252	#ifdef Bad_float_h
michael@0	253
michael@0	254	#ifdef IEEE_Arith
michael@0	255	#define DBL_DIG 15
michael@0	256	#define DBL_MAX_10_EXP 308
michael@0	257	#define DBL_MAX_EXP 1024
michael@0	258	#define FLT_RADIX 2
michael@0	259	#endif /*IEEE_Arith*/
michael@0	260
michael@0	261	#ifdef IBM
michael@0	262	#define DBL_DIG 16
michael@0	263	#define DBL_MAX_10_EXP 75
michael@0	264	#define DBL_MAX_EXP 63
michael@0	265	#define FLT_RADIX 16
michael@0	266	#define DBL_MAX 7.2370055773322621e+75
michael@0	267	#endif
michael@0	268
michael@0	269	#ifdef VAX
michael@0	270	#define DBL_DIG 16
michael@0	271	#define DBL_MAX_10_EXP 38
michael@0	272	#define DBL_MAX_EXP 127
michael@0	273	#define FLT_RADIX 2
michael@0	274	#define DBL_MAX 1.7014118346046923e+38
michael@0	275	#endif
michael@0	276
michael@0	277	#ifndef LONG_MAX
michael@0	278	#define LONG_MAX 2147483647
michael@0	279	#endif
michael@0	280
michael@0	281	#else /* ifndef Bad_float_h */
michael@0	282	#include "float.h"
michael@0	283	#endif /* Bad_float_h */
michael@0	284
michael@0	285	#ifndef __MATH_H__
michael@0	286	#include "math.h"
michael@0	287	#endif
michael@0	288
michael@0	289	#ifdef __cplusplus
michael@0	290	extern "C" {
michael@0	291	#endif
michael@0	292
michael@0	293	#ifndef CONST
michael@0	294	#ifdef KR_headers
michael@0	295	#define CONST /* blank */
michael@0	296	#else
michael@0	297	#define CONST const
michael@0	298	#endif
michael@0	299	#endif
michael@0	300
michael@0	301	#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
michael@0	302	Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
michael@0	303	#endif
michael@0	304
michael@0	305	typedef union { double d; ULong L[2]; } U;
michael@0	306
michael@0	307	#ifdef IEEE_8087
michael@0	308	#define word0(x) (x)->L[1]
michael@0	309	#define word1(x) (x)->L[0]
michael@0	310	#else
michael@0	311	#define word0(x) (x)->L[0]
michael@0	312	#define word1(x) (x)->L[1]
michael@0	313	#endif
michael@0	314	#define dval(x) (x)->d
michael@0	315
michael@0	316	#ifndef STRTOD_DIGLIM
michael@0	317	#define STRTOD_DIGLIM 40
michael@0	318	#endif
michael@0	319
michael@0	320	#ifdef DIGLIM_DEBUG
michael@0	321	extern int strtod_diglim;
michael@0	322	#else
michael@0	323	#define strtod_diglim STRTOD_DIGLIM
michael@0	324	#endif
michael@0	325
michael@0	326	/* The following definition of Storeinc is appropriate for MIPS processors.
michael@0	327	* An alternative that might be better on some machines is
michael@0	328	* #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
michael@0	329	*/
michael@0	330	#if defined(IEEE_8087) + defined(VAX)
michael@0	331	#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
michael@0	332	((unsigned short *)a)[0] = (unsigned short)c, a++)
michael@0	333	#else
michael@0	334	#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
michael@0	335	((unsigned short *)a)[1] = (unsigned short)c, a++)
michael@0	336	#endif
michael@0	337
michael@0	338	/* #define P DBL_MANT_DIG */
michael@0	339	/* Ten_pmax = floor(P*log(2)/log(5)) */
michael@0	340	/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
michael@0	341	/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
michael@0	342	/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
michael@0	343
michael@0	344	#ifdef IEEE_Arith
michael@0	345	#define Exp_shift 20
michael@0	346	#define Exp_shift1 20
michael@0	347	#define Exp_msk1 0x100000
michael@0	348	#define Exp_msk11 0x100000
michael@0	349	#define Exp_mask 0x7ff00000
michael@0	350	#define P 53
michael@0	351	#define Nbits 53
michael@0	352	#define Bias 1023
michael@0	353	#define Emax 1023
michael@0	354	#define Emin (-1022)
michael@0	355	#define Exp_1 0x3ff00000
michael@0	356	#define Exp_11 0x3ff00000
michael@0	357	#define Ebits 11
michael@0	358	#define Frac_mask 0xfffff
michael@0	359	#define Frac_mask1 0xfffff
michael@0	360	#define Ten_pmax 22
michael@0	361	#define Bletch 0x10
michael@0	362	#define Bndry_mask 0xfffff
michael@0	363	#define Bndry_mask1 0xfffff
michael@0	364	#define LSB 1
michael@0	365	#define Sign_bit 0x80000000
michael@0	366	#define Log2P 1
michael@0	367	#define Tiny0 0
michael@0	368	#define Tiny1 1
michael@0	369	#define Quick_max 14
michael@0	370	#define Int_max 14
michael@0	371	#ifndef NO_IEEE_Scale
michael@0	372	#define Avoid_Underflow
michael@0	373	#ifdef Flush_Denorm /* debugging option */
michael@0	374	#undef Sudden_Underflow
michael@0	375	#endif
michael@0	376	#endif
michael@0	377
michael@0	378	#ifndef Flt_Rounds
michael@0	379	#ifdef FLT_ROUNDS
michael@0	380	#define Flt_Rounds FLT_ROUNDS
michael@0	381	#else
michael@0	382	#define Flt_Rounds 1
michael@0	383	#endif
michael@0	384	#endif /*Flt_Rounds*/
michael@0	385
michael@0	386	#ifdef Honor_FLT_ROUNDS
michael@0	387	#undef Check_FLT_ROUNDS
michael@0	388	#define Check_FLT_ROUNDS
michael@0	389	#else
michael@0	390	#define Rounding Flt_Rounds
michael@0	391	#endif
michael@0	392
michael@0	393	#else /* ifndef IEEE_Arith */
michael@0	394	#undef Check_FLT_ROUNDS
michael@0	395	#undef Honor_FLT_ROUNDS
michael@0	396	#undef SET_INEXACT
michael@0	397	#undef Sudden_Underflow
michael@0	398	#define Sudden_Underflow
michael@0	399	#ifdef IBM
michael@0	400	#undef Flt_Rounds
michael@0	401	#define Flt_Rounds 0
michael@0	402	#define Exp_shift 24
michael@0	403	#define Exp_shift1 24
michael@0	404	#define Exp_msk1 0x1000000
michael@0	405	#define Exp_msk11 0x1000000
michael@0	406	#define Exp_mask 0x7f000000
michael@0	407	#define P 14
michael@0	408	#define Nbits 56
michael@0	409	#define Bias 65
michael@0	410	#define Emax 248
michael@0	411	#define Emin (-260)
michael@0	412	#define Exp_1 0x41000000
michael@0	413	#define Exp_11 0x41000000
michael@0	414	#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
michael@0	415	#define Frac_mask 0xffffff
michael@0	416	#define Frac_mask1 0xffffff
michael@0	417	#define Bletch 4
michael@0	418	#define Ten_pmax 22
michael@0	419	#define Bndry_mask 0xefffff
michael@0	420	#define Bndry_mask1 0xffffff
michael@0	421	#define LSB 1
michael@0	422	#define Sign_bit 0x80000000
michael@0	423	#define Log2P 4
michael@0	424	#define Tiny0 0x100000
michael@0	425	#define Tiny1 0
michael@0	426	#define Quick_max 14
michael@0	427	#define Int_max 15
michael@0	428	#else /* VAX */
michael@0	429	#undef Flt_Rounds
michael@0	430	#define Flt_Rounds 1
michael@0	431	#define Exp_shift 23
michael@0	432	#define Exp_shift1 7
michael@0	433	#define Exp_msk1 0x80
michael@0	434	#define Exp_msk11 0x800000
michael@0	435	#define Exp_mask 0x7f80
michael@0	436	#define P 56
michael@0	437	#define Nbits 56
michael@0	438	#define Bias 129
michael@0	439	#define Emax 126
michael@0	440	#define Emin (-129)
michael@0	441	#define Exp_1 0x40800000
michael@0	442	#define Exp_11 0x4080
michael@0	443	#define Ebits 8
michael@0	444	#define Frac_mask 0x7fffff
michael@0	445	#define Frac_mask1 0xffff007f
michael@0	446	#define Ten_pmax 24
michael@0	447	#define Bletch 2
michael@0	448	#define Bndry_mask 0xffff007f
michael@0	449	#define Bndry_mask1 0xffff007f
michael@0	450	#define LSB 0x10000
michael@0	451	#define Sign_bit 0x8000
michael@0	452	#define Log2P 1
michael@0	453	#define Tiny0 0x80
michael@0	454	#define Tiny1 0
michael@0	455	#define Quick_max 15
michael@0	456	#define Int_max 15
michael@0	457	#endif /* IBM, VAX */
michael@0	458	#endif /* IEEE_Arith */
michael@0	459
michael@0	460	#ifndef IEEE_Arith
michael@0	461	#define ROUND_BIASED
michael@0	462	#else
michael@0	463	#ifdef ROUND_BIASED_without_Round_Up
michael@0	464	#undef ROUND_BIASED
michael@0	465	#define ROUND_BIASED
michael@0	466	#endif
michael@0	467	#endif
michael@0	468
michael@0	469	#ifdef RND_PRODQUOT
michael@0	470	#define rounded_product(a,b) a = rnd_prod(a, b)
michael@0	471	#define rounded_quotient(a,b) a = rnd_quot(a, b)
michael@0	472	#ifdef KR_headers
michael@0	473	extern double rnd_prod(), rnd_quot();
michael@0	474	#else
michael@0	475	extern double rnd_prod(double, double), rnd_quot(double, double);
michael@0	476	#endif
michael@0	477	#else
michael@0	478	#define rounded_product(a,b) a *= b
michael@0	479	#define rounded_quotient(a,b) a /= b
michael@0	480	#endif
michael@0	481
michael@0	482	#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
michael@0	483	#define Big1 0xffffffff
michael@0	484
michael@0	485	#ifndef Pack_32
michael@0	486	#define Pack_32
michael@0	487	#endif
michael@0	488
michael@0	489	typedef struct BCinfo BCinfo;
michael@0	490	struct
michael@0	491	BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
michael@0	492
michael@0	493	#ifdef KR_headers
michael@0	494	#define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
michael@0	495	#else
michael@0	496	#define FFFFFFFF 0xffffffffUL
michael@0	497	#endif
michael@0	498
michael@0	499	#ifdef NO_LONG_LONG
michael@0	500	#undef ULLong
michael@0	501	#ifdef Just_16
michael@0	502	#undef Pack_32
michael@0	503	/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
michael@0	504	* This makes some inner loops simpler and sometimes saves work
michael@0	505	* during multiplications, but it often seems to make things slightly
michael@0	506	* slower. Hence the default is now to store 32 bits per Long.
michael@0	507	*/
michael@0	508	#endif
michael@0	509	#else /* long long available */
michael@0	510	#ifndef Llong
michael@0	511	#define Llong long long
michael@0	512	#endif
michael@0	513	#ifndef ULLong
michael@0	514	#define ULLong unsigned Llong
michael@0	515	#endif
michael@0	516	#endif /* NO_LONG_LONG */
michael@0	517
michael@0	518	#ifndef MULTIPLE_THREADS
michael@0	519	#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
michael@0	520	#define FREE_DTOA_LOCK(n) /*nothing*/
michael@0	521	#endif
michael@0	522
michael@0	523	#define Kmax 7
michael@0	524
michael@0	525	#ifdef __cplusplus
michael@0	526	extern "C" double strtod(const char *s00, char **se);
michael@0	527	extern "C" char *dtoa(double d, int mode, int ndigits,
michael@0	528	int *decpt, int *sign, char **rve);
michael@0	529	#endif
michael@0	530
michael@0	531	struct
michael@0	532	Bigint {
michael@0	533	struct Bigint *next;
michael@0	534	int k, maxwds, sign, wds;
michael@0	535	ULong x[1];
michael@0	536	};
michael@0	537
michael@0	538	typedef struct Bigint Bigint;
michael@0	539
michael@0	540	static Bigint *freelist[Kmax+1];
michael@0	541
michael@0	542	static Bigint *
michael@0	543	Balloc
michael@0	544	#ifdef KR_headers
michael@0	545	(k) int k;
michael@0	546	#else
michael@0	547	(int k)
michael@0	548	#endif
michael@0	549	{
michael@0	550	int x;
michael@0	551	Bigint *rv;
michael@0	552	#ifndef Omit_Private_Memory
michael@0	553	unsigned int len;
michael@0	554	#endif
michael@0	555
michael@0	556	ACQUIRE_DTOA_LOCK(0);
michael@0	557	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
michael@0	558	/* but this case seems very unlikely. */
michael@0	559	if (k <= Kmax && (rv = freelist[k]))
michael@0	560	freelist[k] = rv->next;
michael@0	561	else {
michael@0	562	x = 1 << k;
michael@0	563	#ifdef Omit_Private_Memory
michael@0	564	rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
michael@0	565	#else
michael@0	566	len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
michael@0	567	/sizeof(double);
michael@0	568	if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
michael@0	569	rv = (Bigint*)pmem_next;
michael@0	570	pmem_next += len;
michael@0	571	}
michael@0	572	else
michael@0	573	rv = (Bigint*)MALLOC(len*sizeof(double));
michael@0	574	#endif
michael@0	575	rv->k = k;
michael@0	576	rv->maxwds = x;
michael@0	577	}
michael@0	578	FREE_DTOA_LOCK(0);
michael@0	579	rv->sign = rv->wds = 0;
michael@0	580	return rv;
michael@0	581	}
michael@0	582
michael@0	583	static void
michael@0	584	Bfree
michael@0	585	#ifdef KR_headers
michael@0	586	(v) Bigint *v;
michael@0	587	#else
michael@0	588	(Bigint *v)
michael@0	589	#endif
michael@0	590	{
michael@0	591	if (v) {
michael@0	592	if (v->k > Kmax)
michael@0	593	#ifdef FREE
michael@0	594	FREE((void*)v);
michael@0	595	#else
michael@0	596	free((void*)v);
michael@0	597	#endif
michael@0	598	else {
michael@0	599	ACQUIRE_DTOA_LOCK(0);
michael@0	600	v->next = freelist[v->k];
michael@0	601	freelist[v->k] = v;
michael@0	602	FREE_DTOA_LOCK(0);
michael@0	603	}
michael@0	604	}
michael@0	605	}
michael@0	606
michael@0	607	#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
michael@0	608	y->wds*sizeof(Long) + 2*sizeof(int))
michael@0	609
michael@0	610	static Bigint *
michael@0	611	multadd
michael@0	612	#ifdef KR_headers
michael@0	613	(b, m, a) Bigint *b; int m, a;
michael@0	614	#else
michael@0	615	(Bigint *b, int m, int a) /* multiply by m and add a */
michael@0	616	#endif
michael@0	617	{
michael@0	618	int i, wds;
michael@0	619	#ifdef ULLong
michael@0	620	ULong *x;
michael@0	621	ULLong carry, y;
michael@0	622	#else
michael@0	623	ULong carry, *x, y;
michael@0	624	#ifdef Pack_32
michael@0	625	ULong xi, z;
michael@0	626	#endif
michael@0	627	#endif
michael@0	628	Bigint *b1;
michael@0	629
michael@0	630	wds = b->wds;
michael@0	631	x = b->x;
michael@0	632	i = 0;
michael@0	633	carry = a;
michael@0	634	do {
michael@0	635	#ifdef ULLong
michael@0	636	y = *x * (ULLong)m + carry;
michael@0	637	carry = y >> 32;
michael@0	638	*x++ = y & FFFFFFFF;
michael@0	639	#else
michael@0	640	#ifdef Pack_32
michael@0	641	xi = *x;
michael@0	642	y = (xi & 0xffff) * m + carry;
michael@0	643	z = (xi >> 16) * m + (y >> 16);
michael@0	644	carry = z >> 16;
michael@0	645	*x++ = (z << 16) + (y & 0xffff);
michael@0	646	#else
michael@0	647	y = *x * m + carry;
michael@0	648	carry = y >> 16;
michael@0	649	*x++ = y & 0xffff;
michael@0	650	#endif
michael@0	651	#endif
michael@0	652	}
michael@0	653	while(++i < wds);
michael@0	654	if (carry) {
michael@0	655	if (wds >= b->maxwds) {
michael@0	656	b1 = Balloc(b->k+1);
michael@0	657	Bcopy(b1, b);
michael@0	658	Bfree(b);
michael@0	659	b = b1;
michael@0	660	}
michael@0	661	b->x[wds++] = carry;
michael@0	662	b->wds = wds;
michael@0	663	}
michael@0	664	return b;
michael@0	665	}
michael@0	666
michael@0	667	static Bigint *
michael@0	668	s2b
michael@0	669	#ifdef KR_headers
michael@0	670	(s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
michael@0	671	#else
michael@0	672	(const char *s, int nd0, int nd, ULong y9, int dplen)
michael@0	673	#endif
michael@0	674	{
michael@0	675	Bigint *b;
michael@0	676	int i, k;
michael@0	677	Long x, y;
michael@0	678
michael@0	679	x = (nd + 8) / 9;
michael@0	680	for(k = 0, y = 1; x > y; y <<= 1, k++) ;
michael@0	681	#ifdef Pack_32
michael@0	682	b = Balloc(k);
michael@0	683	b->x[0] = y9;
michael@0	684	b->wds = 1;
michael@0	685	#else
michael@0	686	b = Balloc(k+1);
michael@0	687	b->x[0] = y9 & 0xffff;
michael@0	688	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
michael@0	689	#endif
michael@0	690
michael@0	691	i = 9;
michael@0	692	if (9 < nd0) {
michael@0	693	s += 9;
michael@0	694	do b = multadd(b, 10, *s++ - '0');
michael@0	695	while(++i < nd0);
michael@0	696	s += dplen;
michael@0	697	}
michael@0	698	else
michael@0	699	s += dplen + 9;
michael@0	700	for(; i < nd; i++)
michael@0	701	b = multadd(b, 10, *s++ - '0');
michael@0	702	return b;
michael@0	703	}
michael@0	704
michael@0	705	static int
michael@0	706	hi0bits
michael@0	707	#ifdef KR_headers
michael@0	708	(x) ULong x;
michael@0	709	#else
michael@0	710	(ULong x)
michael@0	711	#endif
michael@0	712	{
michael@0	713	int k = 0;
michael@0	714
michael@0	715	if (!(x & 0xffff0000)) {
michael@0	716	k = 16;
michael@0	717	x <<= 16;
michael@0	718	}
michael@0	719	if (!(x & 0xff000000)) {
michael@0	720	k += 8;
michael@0	721	x <<= 8;
michael@0	722	}
michael@0	723	if (!(x & 0xf0000000)) {
michael@0	724	k += 4;
michael@0	725	x <<= 4;
michael@0	726	}
michael@0	727	if (!(x & 0xc0000000)) {
michael@0	728	k += 2;
michael@0	729	x <<= 2;
michael@0	730	}
michael@0	731	if (!(x & 0x80000000)) {
michael@0	732	k++;
michael@0	733	if (!(x & 0x40000000))
michael@0	734	return 32;
michael@0	735	}
michael@0	736	return k;
michael@0	737	}
michael@0	738
michael@0	739	static int
michael@0	740	lo0bits
michael@0	741	#ifdef KR_headers
michael@0	742	(y) ULong *y;
michael@0	743	#else
michael@0	744	(ULong *y)
michael@0	745	#endif
michael@0	746	{
michael@0	747	int k;
michael@0	748	ULong x = *y;
michael@0	749
michael@0	750	if (x & 7) {
michael@0	751	if (x & 1)
michael@0	752	return 0;
michael@0	753	if (x & 2) {
michael@0	754	*y = x >> 1;
michael@0	755	return 1;
michael@0	756	}
michael@0	757	*y = x >> 2;
michael@0	758	return 2;
michael@0	759	}
michael@0	760	k = 0;
michael@0	761	if (!(x & 0xffff)) {
michael@0	762	k = 16;
michael@0	763	x >>= 16;
michael@0	764	}
michael@0	765	if (!(x & 0xff)) {
michael@0	766	k += 8;
michael@0	767	x >>= 8;
michael@0	768	}
michael@0	769	if (!(x & 0xf)) {
michael@0	770	k += 4;
michael@0	771	x >>= 4;
michael@0	772	}
michael@0	773	if (!(x & 0x3)) {
michael@0	774	k += 2;
michael@0	775	x >>= 2;
michael@0	776	}
michael@0	777	if (!(x & 1)) {
michael@0	778	k++;
michael@0	779	x >>= 1;
michael@0	780	if (!x)
michael@0	781	return 32;
michael@0	782	}
michael@0	783	*y = x;
michael@0	784	return k;
michael@0	785	}
michael@0	786
michael@0	787	static Bigint *
michael@0	788	i2b
michael@0	789	#ifdef KR_headers
michael@0	790	(i) int i;
michael@0	791	#else
michael@0	792	(int i)
michael@0	793	#endif
michael@0	794	{
michael@0	795	Bigint *b;
michael@0	796
michael@0	797	b = Balloc(1);
michael@0	798	b->x[0] = i;
michael@0	799	b->wds = 1;
michael@0	800	return b;
michael@0	801	}
michael@0	802
michael@0	803	static Bigint *
michael@0	804	mult
michael@0	805	#ifdef KR_headers
michael@0	806	(a, b) Bigint *a, *b;
michael@0	807	#else
michael@0	808	(Bigint *a, Bigint *b)
michael@0	809	#endif
michael@0	810	{
michael@0	811	Bigint *c;
michael@0	812	int k, wa, wb, wc;
michael@0	813	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
michael@0	814	ULong y;
michael@0	815	#ifdef ULLong
michael@0	816	ULLong carry, z;
michael@0	817	#else
michael@0	818	ULong carry, z;
michael@0	819	#ifdef Pack_32
michael@0	820	ULong z2;
michael@0	821	#endif
michael@0	822	#endif
michael@0	823
michael@0	824	if (a->wds < b->wds) {
michael@0	825	c = a;
michael@0	826	a = b;
michael@0	827	b = c;
michael@0	828	}
michael@0	829	k = a->k;
michael@0	830	wa = a->wds;
michael@0	831	wb = b->wds;
michael@0	832	wc = wa + wb;
michael@0	833	if (wc > a->maxwds)
michael@0	834	k++;
michael@0	835	c = Balloc(k);
michael@0	836	for(x = c->x, xa = x + wc; x < xa; x++)
michael@0	837	*x = 0;
michael@0	838	xa = a->x;
michael@0	839	xae = xa + wa;
michael@0	840	xb = b->x;
michael@0	841	xbe = xb + wb;
michael@0	842	xc0 = c->x;
michael@0	843	#ifdef ULLong
michael@0	844	for(; xb < xbe; xc0++) {
michael@0	845	if ((y = *xb++)) {
michael@0	846	x = xa;
michael@0	847	xc = xc0;
michael@0	848	carry = 0;
michael@0	849	do {
michael@0	850	z = *x++ * (ULLong)y + *xc + carry;
michael@0	851	carry = z >> 32;
michael@0	852	*xc++ = z & FFFFFFFF;
michael@0	853	}
michael@0	854	while(x < xae);
michael@0	855	*xc = carry;
michael@0	856	}
michael@0	857	}
michael@0	858	#else
michael@0	859	#ifdef Pack_32
michael@0	860	for(; xb < xbe; xb++, xc0++) {
michael@0	861	if (y = *xb & 0xffff) {
michael@0	862	x = xa;
michael@0	863	xc = xc0;
michael@0	864	carry = 0;
michael@0	865	do {
michael@0	866	z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
michael@0	867	carry = z >> 16;
michael@0	868	z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
michael@0	869	carry = z2 >> 16;
michael@0	870	Storeinc(xc, z2, z);
michael@0	871	}
michael@0	872	while(x < xae);
michael@0	873	*xc = carry;
michael@0	874	}
michael@0	875	if (y = *xb >> 16) {
michael@0	876	x = xa;
michael@0	877	xc = xc0;
michael@0	878	carry = 0;
michael@0	879	z2 = *xc;
michael@0	880	do {
michael@0	881	z = (*x & 0xffff) * y + (*xc >> 16) + carry;
michael@0	882	carry = z >> 16;
michael@0	883	Storeinc(xc, z, z2);
michael@0	884	z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
michael@0	885	carry = z2 >> 16;
michael@0	886	}
michael@0	887	while(x < xae);
michael@0	888	*xc = z2;
michael@0	889	}
michael@0	890	}
michael@0	891	#else
michael@0	892	for(; xb < xbe; xc0++) {
michael@0	893	if (y = *xb++) {
michael@0	894	x = xa;
michael@0	895	xc = xc0;
michael@0	896	carry = 0;
michael@0	897	do {
michael@0	898	z = *x++ * y + *xc + carry;
michael@0	899	carry = z >> 16;
michael@0	900	*xc++ = z & 0xffff;
michael@0	901	}
michael@0	902	while(x < xae);
michael@0	903	*xc = carry;
michael@0	904	}
michael@0	905	}
michael@0	906	#endif
michael@0	907	#endif
michael@0	908	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
michael@0	909	c->wds = wc;
michael@0	910	return c;
michael@0	911	}
michael@0	912
michael@0	913	static Bigint *p5s;
michael@0	914
michael@0	915	static Bigint *
michael@0	916	pow5mult
michael@0	917	#ifdef KR_headers
michael@0	918	(b, k) Bigint *b; int k;
michael@0	919	#else
michael@0	920	(Bigint *b, int k)
michael@0	921	#endif
michael@0	922	{
michael@0	923	Bigint *b1, *p5, *p51;
michael@0	924	int i;
michael@0	925	static int p05[3] = { 5, 25, 125 };
michael@0	926
michael@0	927	if ((i = k & 3))
michael@0	928	b = multadd(b, p05[i-1], 0);
michael@0	929
michael@0	930	if (!(k >>= 2))
michael@0	931	return b;
michael@0	932	if (!(p5 = p5s)) {
michael@0	933	/* first time */
michael@0	934	#ifdef MULTIPLE_THREADS
michael@0	935	ACQUIRE_DTOA_LOCK(1);
michael@0	936	if (!(p5 = p5s)) {
michael@0	937	p5 = p5s = i2b(625);
michael@0	938	p5->next = 0;
michael@0	939	}
michael@0	940	FREE_DTOA_LOCK(1);
michael@0	941	#else
michael@0	942	p5 = p5s = i2b(625);
michael@0	943	p5->next = 0;
michael@0	944	#endif
michael@0	945	}
michael@0	946	for(;;) {
michael@0	947	if (k & 1) {
michael@0	948	b1 = mult(b, p5);
michael@0	949	Bfree(b);
michael@0	950	b = b1;
michael@0	951	}
michael@0	952	if (!(k >>= 1))
michael@0	953	break;
michael@0	954	if (!(p51 = p5->next)) {
michael@0	955	#ifdef MULTIPLE_THREADS
michael@0	956	ACQUIRE_DTOA_LOCK(1);
michael@0	957	if (!(p51 = p5->next)) {
michael@0	958	p51 = p5->next = mult(p5,p5);
michael@0	959	p51->next = 0;
michael@0	960	}
michael@0	961	FREE_DTOA_LOCK(1);
michael@0	962	#else
michael@0	963	p51 = p5->next = mult(p5,p5);
michael@0	964	p51->next = 0;
michael@0	965	#endif
michael@0	966	}
michael@0	967	p5 = p51;
michael@0	968	}
michael@0	969	return b;
michael@0	970	}
michael@0	971
michael@0	972	static Bigint *
michael@0	973	lshift
michael@0	974	#ifdef KR_headers
michael@0	975	(b, k) Bigint *b; int k;
michael@0	976	#else
michael@0	977	(Bigint *b, int k)
michael@0	978	#endif
michael@0	979	{
michael@0	980	int i, k1, n, n1;
michael@0	981	Bigint *b1;
michael@0	982	ULong *x, *x1, *xe, z;
michael@0	983
michael@0	984	#ifdef Pack_32
michael@0	985	n = k >> 5;
michael@0	986	#else
michael@0	987	n = k >> 4;
michael@0	988	#endif
michael@0	989	k1 = b->k;
michael@0	990	n1 = n + b->wds + 1;
michael@0	991	for(i = b->maxwds; n1 > i; i <<= 1)
michael@0	992	k1++;
michael@0	993	b1 = Balloc(k1);
michael@0	994	x1 = b1->x;
michael@0	995	for(i = 0; i < n; i++)
michael@0	996	*x1++ = 0;
michael@0	997	x = b->x;
michael@0	998	xe = x + b->wds;
michael@0	999	#ifdef Pack_32
michael@0	1000	if (k &= 0x1f) {
michael@0	1001	k1 = 32 - k;
michael@0	1002	z = 0;
michael@0	1003	do {
michael@0	1004	*x1++ = *x << k | z;
michael@0	1005	z = *x++ >> k1;
michael@0	1006	}
michael@0	1007	while(x < xe);
michael@0	1008	if ((*x1 = z))
michael@0	1009	++n1;
michael@0	1010	}
michael@0	1011	#else
michael@0	1012	if (k &= 0xf) {
michael@0	1013	k1 = 16 - k;
michael@0	1014	z = 0;
michael@0	1015	do {
michael@0	1016	*x1++ = *x << k & 0xffff | z;
michael@0	1017	z = *x++ >> k1;
michael@0	1018	}
michael@0	1019	while(x < xe);
michael@0	1020	if (*x1 = z)
michael@0	1021	++n1;
michael@0	1022	}
michael@0	1023	#endif
michael@0	1024	else do
michael@0	1025	*x1++ = *x++;
michael@0	1026	while(x < xe);
michael@0	1027	b1->wds = n1 - 1;
michael@0	1028	Bfree(b);
michael@0	1029	return b1;
michael@0	1030	}
michael@0	1031
michael@0	1032	static int
michael@0	1033	cmp
michael@0	1034	#ifdef KR_headers
michael@0	1035	(a, b) Bigint *a, *b;
michael@0	1036	#else
michael@0	1037	(Bigint *a, Bigint *b)
michael@0	1038	#endif
michael@0	1039	{
michael@0	1040	ULong *xa, *xa0, *xb, *xb0;
michael@0	1041	int i, j;
michael@0	1042
michael@0	1043	i = a->wds;
michael@0	1044	j = b->wds;
michael@0	1045	#ifdef DEBUG
michael@0	1046	if (i > 1 && !a->x[i-1])
michael@0	1047	Bug("cmp called with a->x[a->wds-1] == 0");
michael@0	1048	if (j > 1 && !b->x[j-1])
michael@0	1049	Bug("cmp called with b->x[b->wds-1] == 0");
michael@0	1050	#endif
michael@0	1051	if (i -= j)
michael@0	1052	return i;
michael@0	1053	xa0 = a->x;
michael@0	1054	xa = xa0 + j;
michael@0	1055	xb0 = b->x;
michael@0	1056	xb = xb0 + j;
michael@0	1057	for(;;) {
michael@0	1058	if (*--xa != *--xb)
michael@0	1059	return *xa < *xb ? -1 : 1;
michael@0	1060	if (xa <= xa0)
michael@0	1061	break;
michael@0	1062	}
michael@0	1063	return 0;
michael@0	1064	}
michael@0	1065
michael@0	1066	static Bigint *
michael@0	1067	diff
michael@0	1068	#ifdef KR_headers
michael@0	1069	(a, b) Bigint *a, *b;
michael@0	1070	#else
michael@0	1071	(Bigint *a, Bigint *b)
michael@0	1072	#endif
michael@0	1073	{
michael@0	1074	Bigint *c;
michael@0	1075	int i, wa, wb;
michael@0	1076	ULong *xa, *xae, *xb, *xbe, *xc;
michael@0	1077	#ifdef ULLong
michael@0	1078	ULLong borrow, y;
michael@0	1079	#else
michael@0	1080	ULong borrow, y;
michael@0	1081	#ifdef Pack_32
michael@0	1082	ULong z;
michael@0	1083	#endif
michael@0	1084	#endif
michael@0	1085
michael@0	1086	i = cmp(a,b);
michael@0	1087	if (!i) {
michael@0	1088	c = Balloc(0);
michael@0	1089	c->wds = 1;
michael@0	1090	c->x[0] = 0;
michael@0	1091	return c;
michael@0	1092	}
michael@0	1093	if (i < 0) {
michael@0	1094	c = a;
michael@0	1095	a = b;
michael@0	1096	b = c;
michael@0	1097	i = 1;
michael@0	1098	}
michael@0	1099	else
michael@0	1100	i = 0;
michael@0	1101	c = Balloc(a->k);
michael@0	1102	c->sign = i;
michael@0	1103	wa = a->wds;
michael@0	1104	xa = a->x;
michael@0	1105	xae = xa + wa;
michael@0	1106	wb = b->wds;
michael@0	1107	xb = b->x;
michael@0	1108	xbe = xb + wb;
michael@0	1109	xc = c->x;
michael@0	1110	borrow = 0;
michael@0	1111	#ifdef ULLong
michael@0	1112	do {
michael@0	1113	y = (ULLong)*xa++ - *xb++ - borrow;
michael@0	1114	borrow = y >> 32 & (ULong)1;
michael@0	1115	*xc++ = y & FFFFFFFF;
michael@0	1116	}
michael@0	1117	while(xb < xbe);
michael@0	1118	while(xa < xae) {
michael@0	1119	y = *xa++ - borrow;
michael@0	1120	borrow = y >> 32 & (ULong)1;
michael@0	1121	*xc++ = y & FFFFFFFF;
michael@0	1122	}
michael@0	1123	#else
michael@0	1124	#ifdef Pack_32
michael@0	1125	do {
michael@0	1126	y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
michael@0	1127	borrow = (y & 0x10000) >> 16;
michael@0	1128	z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
michael@0	1129	borrow = (z & 0x10000) >> 16;
michael@0	1130	Storeinc(xc, z, y);
michael@0	1131	}
michael@0	1132	while(xb < xbe);
michael@0	1133	while(xa < xae) {
michael@0	1134	y = (*xa & 0xffff) - borrow;
michael@0	1135	borrow = (y & 0x10000) >> 16;
michael@0	1136	z = (*xa++ >> 16) - borrow;
michael@0	1137	borrow = (z & 0x10000) >> 16;
michael@0	1138	Storeinc(xc, z, y);
michael@0	1139	}
michael@0	1140	#else
michael@0	1141	do {
michael@0	1142	y = *xa++ - *xb++ - borrow;
michael@0	1143	borrow = (y & 0x10000) >> 16;
michael@0	1144	*xc++ = y & 0xffff;
michael@0	1145	}
michael@0	1146	while(xb < xbe);
michael@0	1147	while(xa < xae) {
michael@0	1148	y = *xa++ - borrow;
michael@0	1149	borrow = (y & 0x10000) >> 16;
michael@0	1150	*xc++ = y & 0xffff;
michael@0	1151	}
michael@0	1152	#endif
michael@0	1153	#endif
michael@0	1154	while(!*--xc)
michael@0	1155	wa--;
michael@0	1156	c->wds = wa;
michael@0	1157	return c;
michael@0	1158	}
michael@0	1159
michael@0	1160	static double
michael@0	1161	ulp
michael@0	1162	#ifdef KR_headers
michael@0	1163	(x) U *x;
michael@0	1164	#else
michael@0	1165	(U *x)
michael@0	1166	#endif
michael@0	1167	{
michael@0	1168	Long L;
michael@0	1169	U u;
michael@0	1170
michael@0	1171	L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
michael@0	1172	#ifndef Avoid_Underflow
michael@0	1173	#ifndef Sudden_Underflow
michael@0	1174	if (L > 0) {
michael@0	1175	#endif
michael@0	1176	#endif
michael@0	1177	#ifdef IBM
michael@0	1178	L |= Exp_msk1 >> 4;
michael@0	1179	#endif
michael@0	1180	word0(&u) = L;
michael@0	1181	word1(&u) = 0;
michael@0	1182	#ifndef Avoid_Underflow
michael@0	1183	#ifndef Sudden_Underflow
michael@0	1184	}
michael@0	1185	else {
michael@0	1186	L = -L >> Exp_shift;
michael@0	1187	if (L < Exp_shift) {
michael@0	1188	word0(&u) = 0x80000 >> L;
michael@0	1189	word1(&u) = 0;
michael@0	1190	}
michael@0	1191	else {
michael@0	1192	word0(&u) = 0;
michael@0	1193	L -= Exp_shift;
michael@0	1194	word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
michael@0	1195	}
michael@0	1196	}
michael@0	1197	#endif
michael@0	1198	#endif
michael@0	1199	return dval(&u);
michael@0	1200	}
michael@0	1201
michael@0	1202	static double
michael@0	1203	b2d
michael@0	1204	#ifdef KR_headers
michael@0	1205	(a, e) Bigint *a; int *e;
michael@0	1206	#else
michael@0	1207	(Bigint *a, int *e)
michael@0	1208	#endif
michael@0	1209	{
michael@0	1210	ULong *xa, *xa0, w, y, z;
michael@0	1211	int k;
michael@0	1212	U d;
michael@0	1213	#ifdef VAX
michael@0	1214	ULong d0, d1;
michael@0	1215	#else
michael@0	1216	#define d0 word0(&d)
michael@0	1217	#define d1 word1(&d)
michael@0	1218	#endif
michael@0	1219
michael@0	1220	xa0 = a->x;
michael@0	1221	xa = xa0 + a->wds;
michael@0	1222	y = *--xa;
michael@0	1223	#ifdef DEBUG
michael@0	1224	if (!y) Bug("zero y in b2d");
michael@0	1225	#endif
michael@0	1226	k = hi0bits(y);
michael@0	1227	*e = 32 - k;
michael@0	1228	#ifdef Pack_32
michael@0	1229	if (k < Ebits) {
michael@0	1230	d0 = Exp_1 | y >> (Ebits - k);
michael@0	1231	w = xa > xa0 ? *--xa : 0;
michael@0	1232	d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
michael@0	1233	goto ret_d;
michael@0	1234	}
michael@0	1235	z = xa > xa0 ? *--xa : 0;
michael@0	1236	if (k -= Ebits) {
michael@0	1237	d0 = Exp_1 | y << k | z >> (32 - k);
michael@0	1238	y = xa > xa0 ? *--xa : 0;
michael@0	1239	d1 = z << k | y >> (32 - k);
michael@0	1240	}
michael@0	1241	else {
michael@0	1242	d0 = Exp_1 | y;
michael@0	1243	d1 = z;
michael@0	1244	}
michael@0	1245	#else
michael@0	1246	if (k < Ebits + 16) {
michael@0	1247	z = xa > xa0 ? *--xa : 0;
michael@0	1248	d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
michael@0	1249	w = xa > xa0 ? *--xa : 0;
michael@0	1250	y = xa > xa0 ? *--xa : 0;
michael@0	1251	d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
michael@0	1252	goto ret_d;
michael@0	1253	}
michael@0	1254	z = xa > xa0 ? *--xa : 0;
michael@0	1255	w = xa > xa0 ? *--xa : 0;
michael@0	1256	k -= Ebits + 16;
michael@0	1257	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
michael@0	1258	y = xa > xa0 ? *--xa : 0;
michael@0	1259	d1 = w << k + 16 | y << k;
michael@0	1260	#endif
michael@0	1261	ret_d:
michael@0	1262	#ifdef VAX
michael@0	1263	word0(&d) = d0 >> 16 | d0 << 16;
michael@0	1264	word1(&d) = d1 >> 16 | d1 << 16;
michael@0	1265	#else
michael@0	1266	#undef d0
michael@0	1267	#undef d1
michael@0	1268	#endif
michael@0	1269	return dval(&d);
michael@0	1270	}
michael@0	1271
michael@0	1272	static Bigint *
michael@0	1273	d2b
michael@0	1274	#ifdef KR_headers
michael@0	1275	(d, e, bits) U *d; int *e, *bits;
michael@0	1276	#else
michael@0	1277	(U *d, int *e, int *bits)
michael@0	1278	#endif
michael@0	1279	{
michael@0	1280	Bigint *b;
michael@0	1281	int de, k;
michael@0	1282	ULong *x, y, z;
michael@0	1283	#ifndef Sudden_Underflow
michael@0	1284	int i;
michael@0	1285	#endif
michael@0	1286	#ifdef VAX
michael@0	1287	ULong d0, d1;
michael@0	1288	d0 = word0(d) >> 16 | word0(d) << 16;
michael@0	1289	d1 = word1(d) >> 16 | word1(d) << 16;
michael@0	1290	#else
michael@0	1291	#define d0 word0(d)
michael@0	1292	#define d1 word1(d)
michael@0	1293	#endif
michael@0	1294
michael@0	1295	#ifdef Pack_32
michael@0	1296	b = Balloc(1);
michael@0	1297	#else
michael@0	1298	b = Balloc(2);
michael@0	1299	#endif
michael@0	1300	x = b->x;
michael@0	1301
michael@0	1302	z = d0 & Frac_mask;
michael@0	1303	d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
michael@0	1304	#ifdef Sudden_Underflow
michael@0	1305	de = (int)(d0 >> Exp_shift);
michael@0	1306	#ifndef IBM
michael@0	1307	z |= Exp_msk11;
michael@0	1308	#endif
michael@0	1309	#else
michael@0	1310	if ((de = (int)(d0 >> Exp_shift)))
michael@0	1311	z |= Exp_msk1;
michael@0	1312	#endif
michael@0	1313	#ifdef Pack_32
michael@0	1314	if ((y = d1)) {
michael@0	1315	if ((k = lo0bits(&y))) {
michael@0	1316	x[0] = y | z << (32 - k);
michael@0	1317	z >>= k;
michael@0	1318	}
michael@0	1319	else
michael@0	1320	x[0] = y;
michael@0	1321	#ifndef Sudden_Underflow
michael@0	1322	i =
michael@0	1323	#endif
michael@0	1324	b->wds = (x[1] = z) ? 2 : 1;
michael@0	1325	}
michael@0	1326	else {
michael@0	1327	k = lo0bits(&z);
michael@0	1328	x[0] = z;
michael@0	1329	#ifndef Sudden_Underflow
michael@0	1330	i =
michael@0	1331	#endif
michael@0	1332	b->wds = 1;
michael@0	1333	k += 32;
michael@0	1334	}
michael@0	1335	#else
michael@0	1336	if (y = d1) {
michael@0	1337	if (k = lo0bits(&y))
michael@0	1338	if (k >= 16) {
michael@0	1339	x[0] = y | z << 32 - k & 0xffff;
michael@0	1340	x[1] = z >> k - 16 & 0xffff;
michael@0	1341	x[2] = z >> k;
michael@0	1342	i = 2;
michael@0	1343	}
michael@0	1344	else {
michael@0	1345	x[0] = y & 0xffff;
michael@0	1346	x[1] = y >> 16 | z << 16 - k & 0xffff;
michael@0	1347	x[2] = z >> k & 0xffff;
michael@0	1348	x[3] = z >> k+16;
michael@0	1349	i = 3;
michael@0	1350	}
michael@0	1351	else {
michael@0	1352	x[0] = y & 0xffff;
michael@0	1353	x[1] = y >> 16;
michael@0	1354	x[2] = z & 0xffff;
michael@0	1355	x[3] = z >> 16;
michael@0	1356	i = 3;
michael@0	1357	}
michael@0	1358	}
michael@0	1359	else {
michael@0	1360	#ifdef DEBUG
michael@0	1361	if (!z)
michael@0	1362	Bug("Zero passed to d2b");
michael@0	1363	#endif
michael@0	1364	k = lo0bits(&z);
michael@0	1365	if (k >= 16) {
michael@0	1366	x[0] = z;
michael@0	1367	i = 0;
michael@0	1368	}
michael@0	1369	else {
michael@0	1370	x[0] = z & 0xffff;
michael@0	1371	x[1] = z >> 16;
michael@0	1372	i = 1;
michael@0	1373	}
michael@0	1374	k += 32;
michael@0	1375	}
michael@0	1376	while(!x[i])
michael@0	1377	--i;
michael@0	1378	b->wds = i + 1;
michael@0	1379	#endif
michael@0	1380	#ifndef Sudden_Underflow
michael@0	1381	if (de) {
michael@0	1382	#endif
michael@0	1383	#ifdef IBM
michael@0	1384	*e = (de - Bias - (P-1) << 2) + k;
michael@0	1385	*bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
michael@0	1386	#else
michael@0	1387	*e = de - Bias - (P-1) + k;
michael@0	1388	*bits = P - k;
michael@0	1389	#endif
michael@0	1390	#ifndef Sudden_Underflow
michael@0	1391	}
michael@0	1392	else {
michael@0	1393	*e = de - Bias - (P-1) + 1 + k;
michael@0	1394	#ifdef Pack_32
michael@0	1395	*bits = 32*i - hi0bits(x[i-1]);
michael@0	1396	#else
michael@0	1397	*bits = (i+2)*16 - hi0bits(x[i]);
michael@0	1398	#endif
michael@0	1399	}
michael@0	1400	#endif
michael@0	1401	return b;
michael@0	1402	}
michael@0	1403	#undef d0
michael@0	1404	#undef d1
michael@0	1405
michael@0	1406	static double
michael@0	1407	ratio
michael@0	1408	#ifdef KR_headers
michael@0	1409	(a, b) Bigint *a, *b;
michael@0	1410	#else
michael@0	1411	(Bigint *a, Bigint *b)
michael@0	1412	#endif
michael@0	1413	{
michael@0	1414	U da, db;
michael@0	1415	int k, ka, kb;
michael@0	1416
michael@0	1417	dval(&da) = b2d(a, &ka);
michael@0	1418	dval(&db) = b2d(b, &kb);
michael@0	1419	#ifdef Pack_32
michael@0	1420	k = ka - kb + 32*(a->wds - b->wds);
michael@0	1421	#else
michael@0	1422	k = ka - kb + 16*(a->wds - b->wds);
michael@0	1423	#endif
michael@0	1424	#ifdef IBM
michael@0	1425	if (k > 0) {
michael@0	1426	word0(&da) += (k >> 2)*Exp_msk1;
michael@0	1427	if (k &= 3)
michael@0	1428	dval(&da) *= 1 << k;
michael@0	1429	}
michael@0	1430	else {
michael@0	1431	k = -k;
michael@0	1432	word0(&db) += (k >> 2)*Exp_msk1;
michael@0	1433	if (k &= 3)
michael@0	1434	dval(&db) *= 1 << k;
michael@0	1435	}
michael@0	1436	#else
michael@0	1437	if (k > 0)
michael@0	1438	word0(&da) += k*Exp_msk1;
michael@0	1439	else {
michael@0	1440	k = -k;
michael@0	1441	word0(&db) += k*Exp_msk1;
michael@0	1442	}
michael@0	1443	#endif
michael@0	1444	return dval(&da) / dval(&db);
michael@0	1445	}
michael@0	1446
michael@0	1447	static CONST double
michael@0	1448	tens[] = {
michael@0	1449	1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
michael@0	1450	1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
michael@0	1451	1e20, 1e21, 1e22
michael@0	1452	#ifdef VAX
michael@0	1453	, 1e23, 1e24
michael@0	1454	#endif
michael@0	1455	};
michael@0	1456
michael@0	1457	static CONST double
michael@0	1458	#ifdef IEEE_Arith
michael@0	1459	bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
michael@0	1460	static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
michael@0	1461	#ifdef Avoid_Underflow
michael@0	1462	9007199254740992.*9007199254740992.e-256
michael@0	1463	/* = 2^106 * 1e-256 */
michael@0	1464	#else
michael@0	1465	1e-256
michael@0	1466	#endif
michael@0	1467	};
michael@0	1468	/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
michael@0	1469	/* flag unnecessarily. It leads to a song and dance at the end of strtod. */
michael@0	1470	#define Scale_Bit 0x10
michael@0	1471	#define n_bigtens 5
michael@0	1472	#else
michael@0	1473	#ifdef IBM
michael@0	1474	bigtens[] = { 1e16, 1e32, 1e64 };
michael@0	1475	static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
michael@0	1476	#define n_bigtens 3
michael@0	1477	#else
michael@0	1478	bigtens[] = { 1e16, 1e32 };
michael@0	1479	static CONST double tinytens[] = { 1e-16, 1e-32 };
michael@0	1480	#define n_bigtens 2
michael@0	1481	#endif
michael@0	1482	#endif
michael@0	1483
michael@0	1484	#undef Need_Hexdig
michael@0	1485	#ifdef INFNAN_CHECK
michael@0	1486	#ifndef No_Hex_NaN
michael@0	1487	#define Need_Hexdig
michael@0	1488	#endif
michael@0	1489	#endif
michael@0	1490
michael@0	1491	#ifndef Need_Hexdig
michael@0	1492	#ifndef NO_HEX_FP
michael@0	1493	#define Need_Hexdig
michael@0	1494	#endif
michael@0	1495	#endif
michael@0	1496
michael@0	1497	#ifdef Need_Hexdig /*{*/
michael@0	1498	static unsigned char hexdig[256];
michael@0	1499
michael@0	1500	static void
michael@0	1501	#ifdef KR_headers
michael@0	1502	htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
michael@0	1503	#else
michael@0	1504	htinit(unsigned char *h, unsigned char *s, int inc)
michael@0	1505	#endif
michael@0	1506	{
michael@0	1507	int i, j;
michael@0	1508	for(i = 0; (j = s[i]) !=0; i++)
michael@0	1509	h[j] = i + inc;
michael@0	1510	}
michael@0	1511
michael@0	1512	static void
michael@0	1513	#ifdef KR_headers
michael@0	1514	hexdig_init()
michael@0	1515	#else
michael@0	1516	hexdig_init(void)
michael@0	1517	#endif
michael@0	1518	{
michael@0	1519	#define USC (unsigned char *)
michael@0	1520	htinit(hexdig, USC "0123456789", 0x10);
michael@0	1521	htinit(hexdig, USC "abcdef", 0x10 + 10);
michael@0	1522	htinit(hexdig, USC "ABCDEF", 0x10 + 10);
michael@0	1523	}
michael@0	1524	#endif /* } Need_Hexdig */
michael@0	1525
michael@0	1526	#ifdef INFNAN_CHECK
michael@0	1527
michael@0	1528	#ifndef NAN_WORD0
michael@0	1529	#define NAN_WORD0 0x7ff80000
michael@0	1530	#endif
michael@0	1531
michael@0	1532	#ifndef NAN_WORD1
michael@0	1533	#define NAN_WORD1 0
michael@0	1534	#endif
michael@0	1535
michael@0	1536	static int
michael@0	1537	match
michael@0	1538	#ifdef KR_headers
michael@0	1539	(sp, t) char **sp, *t;
michael@0	1540	#else
michael@0	1541	(const char **sp, const char *t)
michael@0	1542	#endif
michael@0	1543	{
michael@0	1544	int c, d;
michael@0	1545	CONST char *s = *sp;
michael@0	1546
michael@0	1547	while((d = *t++)) {
michael@0	1548	if ((c = *++s) >= 'A' && c <= 'Z')
michael@0	1549	c += 'a' - 'A';
michael@0	1550	if (c != d)
michael@0	1551	return 0;
michael@0	1552	}
michael@0	1553	*sp = s + 1;
michael@0	1554	return 1;
michael@0	1555	}
michael@0	1556
michael@0	1557	#ifndef No_Hex_NaN
michael@0	1558	static void
michael@0	1559	hexnan
michael@0	1560	#ifdef KR_headers
michael@0	1561	(rvp, sp) U *rvp; CONST char **sp;
michael@0	1562	#else
michael@0	1563	(U *rvp, const char **sp)
michael@0	1564	#endif
michael@0	1565	{
michael@0	1566	ULong c, x[2];
michael@0	1567	CONST char *s;
michael@0	1568	int c1, havedig, udx0, xshift;
michael@0	1569
michael@0	1570	if (!hexdig['0'])
michael@0	1571	hexdig_init();
michael@0	1572	x[0] = x[1] = 0;
michael@0	1573	havedig = xshift = 0;
michael@0	1574	udx0 = 1;
michael@0	1575	s = *sp;
michael@0	1576	/* allow optional initial 0x or 0X */
michael@0	1577	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
michael@0	1578	++s;
michael@0	1579	if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
michael@0	1580	s += 2;
michael@0	1581	while((c = *(CONST unsigned char*)++s)) {
michael@0	1582	if ((c1 = hexdig[c]))
michael@0	1583	c = c1 & 0xf;
michael@0	1584	else if (c <= ' ') {
michael@0	1585	if (udx0 && havedig) {
michael@0	1586	udx0 = 0;
michael@0	1587	xshift = 1;
michael@0	1588	}
michael@0	1589	continue;
michael@0	1590	}
michael@0	1591	#ifdef GDTOA_NON_PEDANTIC_NANCHECK
michael@0	1592	else if (/*(*/ c == ')' && havedig) {
michael@0	1593	*sp = s + 1;
michael@0	1594	break;
michael@0	1595	}
michael@0	1596	else
michael@0	1597	return; /* invalid form: don't change *sp */
michael@0	1598	#else
michael@0	1599	else {
michael@0	1600	do {
michael@0	1601	if (/*(*/ c == ')') {
michael@0	1602	*sp = s + 1;
michael@0	1603	break;
michael@0	1604	}
michael@0	1605	} while((c = *++s));
michael@0	1606	break;
michael@0	1607	}
michael@0	1608	#endif
michael@0	1609	havedig = 1;
michael@0	1610	if (xshift) {
michael@0	1611	xshift = 0;
michael@0	1612	x[0] = x[1];
michael@0	1613	x[1] = 0;
michael@0	1614	}
michael@0	1615	if (udx0)
michael@0	1616	x[0] = (x[0] << 4) | (x[1] >> 28);
michael@0	1617	x[1] = (x[1] << 4) | c;
michael@0	1618	}
michael@0	1619	if ((x[0] &= 0xfffff) || x[1]) {
michael@0	1620	word0(rvp) = Exp_mask | x[0];
michael@0	1621	word1(rvp) = x[1];
michael@0	1622	}
michael@0	1623	}
michael@0	1624	#endif /*No_Hex_NaN*/
michael@0	1625	#endif /* INFNAN_CHECK */
michael@0	1626
michael@0	1627	#ifdef Pack_32
michael@0	1628	#define ULbits 32
michael@0	1629	#define kshift 5
michael@0	1630	#define kmask 31
michael@0	1631	#else
michael@0	1632	#define ULbits 16
michael@0	1633	#define kshift 4
michael@0	1634	#define kmask 15
michael@0	1635	#endif
michael@0	1636
michael@0	1637	#if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
michael@0	1638	static Bigint *
michael@0	1639	#ifdef KR_headers
michael@0	1640	increment(b) Bigint *b;
michael@0	1641	#else
michael@0	1642	increment(Bigint *b)
michael@0	1643	#endif
michael@0	1644	{
michael@0	1645	ULong *x, *xe;
michael@0	1646	Bigint *b1;
michael@0	1647
michael@0	1648	x = b->x;
michael@0	1649	xe = x + b->wds;
michael@0	1650	do {
michael@0	1651	if (*x < (ULong)0xffffffffL) {
michael@0	1652	++*x;
michael@0	1653	return b;
michael@0	1654	}
michael@0	1655	*x++ = 0;
michael@0	1656	} while(x < xe);
michael@0	1657	{
michael@0	1658	if (b->wds >= b->maxwds) {
michael@0	1659	b1 = Balloc(b->k+1);
michael@0	1660	Bcopy(b1,b);
michael@0	1661	Bfree(b);
michael@0	1662	b = b1;
michael@0	1663	}
michael@0	1664	b->x[b->wds++] = 1;
michael@0	1665	}
michael@0	1666	return b;
michael@0	1667	}
michael@0	1668
michael@0	1669	#endif /*}*/
michael@0	1670
michael@0	1671	#ifndef NO_HEX_FP /*{*/
michael@0	1672
michael@0	1673	static void
michael@0	1674	#ifdef KR_headers
michael@0	1675	rshift(b, k) Bigint *b; int k;
michael@0	1676	#else
michael@0	1677	rshift(Bigint *b, int k)
michael@0	1678	#endif
michael@0	1679	{
michael@0	1680	ULong *x, *x1, *xe, y;
michael@0	1681	int n;
michael@0	1682
michael@0	1683	x = x1 = b->x;
michael@0	1684	n = k >> kshift;
michael@0	1685	if (n < b->wds) {
michael@0	1686	xe = x + b->wds;
michael@0	1687	x += n;
michael@0	1688	if (k &= kmask) {
michael@0	1689	n = 32 - k;
michael@0	1690	y = *x++ >> k;
michael@0	1691	while(x < xe) {
michael@0	1692	*x1++ = (y | (*x << n)) & 0xffffffff;
michael@0	1693	y = *x++ >> k;
michael@0	1694	}
michael@0	1695	if ((*x1 = y) !=0)
michael@0	1696	x1++;
michael@0	1697	}
michael@0	1698	else
michael@0	1699	while(x < xe)
michael@0	1700	*x1++ = *x++;
michael@0	1701	}
michael@0	1702	if ((b->wds = x1 - b->x) == 0)
michael@0	1703	b->x[0] = 0;
michael@0	1704	}
michael@0	1705
michael@0	1706	static ULong
michael@0	1707	#ifdef KR_headers
michael@0	1708	any_on(b, k) Bigint *b; int k;
michael@0	1709	#else
michael@0	1710	any_on(Bigint *b, int k)
michael@0	1711	#endif
michael@0	1712	{
michael@0	1713	int n, nwds;
michael@0	1714	ULong *x, *x0, x1, x2;
michael@0	1715
michael@0	1716	x = b->x;
michael@0	1717	nwds = b->wds;
michael@0	1718	n = k >> kshift;
michael@0	1719	if (n > nwds)
michael@0	1720	n = nwds;
michael@0	1721	else if (n < nwds && (k &= kmask)) {
michael@0	1722	x1 = x2 = x[n];
michael@0	1723	x1 >>= k;
michael@0	1724	x1 <<= k;
michael@0	1725	if (x1 != x2)
michael@0	1726	return 1;
michael@0	1727	}
michael@0	1728	x0 = x;
michael@0	1729	x += n;
michael@0	1730	while(x > x0)
michael@0	1731	if (*--x)
michael@0	1732	return 1;
michael@0	1733	return 0;
michael@0	1734	}
michael@0	1735
michael@0	1736	enum { /* rounding values: same as FLT_ROUNDS */
michael@0	1737	Round_zero = 0,
michael@0	1738	Round_near = 1,
michael@0	1739	Round_up = 2,
michael@0	1740	Round_down = 3
michael@0	1741	};
michael@0	1742
michael@0	1743	void
michael@0	1744	#ifdef KR_headers
michael@0	1745	gethex(sp, rvp, rounding, sign)
michael@0	1746	CONST char **sp; U *rvp; int rounding, sign;
michael@0	1747	#else
michael@0	1748	gethex( CONST char **sp, U *rvp, int rounding, int sign)
michael@0	1749	#endif
michael@0	1750	{
michael@0	1751	Bigint *b;
michael@0	1752	CONST unsigned char *decpt, *s0, *s, *s1;
michael@0	1753	Long e, e1;
michael@0	1754	ULong L, lostbits, *x;
michael@0	1755	int big, denorm, esign, havedig, k, n, nbits, up, zret;
michael@0	1756	#ifdef IBM
michael@0	1757	int j;
michael@0	1758	#endif
michael@0	1759	enum {
michael@0	1760	#ifdef IEEE_Arith /*{{*/
michael@0	1761	emax = 0x7fe - Bias - P + 1,
michael@0	1762	emin = Emin - P + 1
michael@0	1763	#else /*}{*/
michael@0	1764	emin = Emin - P,
michael@0	1765	#ifdef VAX
michael@0	1766	emax = 0x7ff - Bias - P + 1
michael@0	1767	#endif
michael@0	1768	#ifdef IBM
michael@0	1769	emax = 0x7f - Bias - P
michael@0	1770	#endif
michael@0	1771	#endif /*}}*/
michael@0	1772	};
michael@0	1773	#ifdef USE_LOCALE
michael@0	1774	int i;
michael@0	1775	#ifdef NO_LOCALE_CACHE
michael@0	1776	const unsigned char *decimalpoint = (unsigned char*)
michael@0	1777	localeconv()->decimal_point;
michael@0	1778	#else
michael@0	1779	const unsigned char *decimalpoint;
michael@0	1780	static unsigned char *decimalpoint_cache;
michael@0	1781	if (!(s0 = decimalpoint_cache)) {
michael@0	1782	s0 = (unsigned char*)localeconv()->decimal_point;
michael@0	1783	if ((decimalpoint_cache = (unsigned char*)
michael@0	1784	MALLOC(strlen((CONST char*)s0) + 1))) {
michael@0	1785	strcpy((char*)decimalpoint_cache, (CONST char*)s0);
michael@0	1786	s0 = decimalpoint_cache;
michael@0	1787	}
michael@0	1788	}
michael@0	1789	decimalpoint = s0;
michael@0	1790	#endif
michael@0	1791	#endif
michael@0	1792
michael@0	1793	if (!hexdig['0'])
michael@0	1794	hexdig_init();
michael@0	1795	havedig = 0;
michael@0	1796	s0 = *(CONST unsigned char **)sp + 2;
michael@0	1797	while(s0[havedig] == '0')
michael@0	1798	havedig++;
michael@0	1799	s0 += havedig;
michael@0	1800	s = s0;
michael@0	1801	decpt = 0;
michael@0	1802	zret = 0;
michael@0	1803	e = 0;
michael@0	1804	if (hexdig[*s])
michael@0	1805	havedig++;
michael@0	1806	else {
michael@0	1807	zret = 1;
michael@0	1808	#ifdef USE_LOCALE
michael@0	1809	for(i = 0; decimalpoint[i]; ++i) {
michael@0	1810	if (s[i] != decimalpoint[i])
michael@0	1811	goto pcheck;
michael@0	1812	}
michael@0	1813	decpt = s += i;
michael@0	1814	#else
michael@0	1815	if (*s != '.')
michael@0	1816	goto pcheck;
michael@0	1817	decpt = ++s;
michael@0	1818	#endif
michael@0	1819	if (!hexdig[*s])
michael@0	1820	goto pcheck;
michael@0	1821	while(*s == '0')
michael@0	1822	s++;
michael@0	1823	if (hexdig[*s])
michael@0	1824	zret = 0;
michael@0	1825	havedig = 1;
michael@0	1826	s0 = s;
michael@0	1827	}
michael@0	1828	while(hexdig[*s])
michael@0	1829	s++;
michael@0	1830	#ifdef USE_LOCALE
michael@0	1831	if (*s == *decimalpoint && !decpt) {
michael@0	1832	for(i = 1; decimalpoint[i]; ++i) {
michael@0	1833	if (s[i] != decimalpoint[i])
michael@0	1834	goto pcheck;
michael@0	1835	}
michael@0	1836	decpt = s += i;
michael@0	1837	#else
michael@0	1838	if (*s == '.' && !decpt) {
michael@0	1839	decpt = ++s;
michael@0	1840	#endif
michael@0	1841	while(hexdig[*s])
michael@0	1842	s++;
michael@0	1843	}/*}*/
michael@0	1844	if (decpt)
michael@0	1845	e = -(((Long)(s-decpt)) << 2);
michael@0	1846	pcheck:
michael@0	1847	s1 = s;
michael@0	1848	big = esign = 0;
michael@0	1849	switch(*s) {
michael@0	1850	case 'p':
michael@0	1851	case 'P':
michael@0	1852	switch(*++s) {
michael@0	1853	case '-':
michael@0	1854	esign = 1;
michael@0	1855	/* no break */
michael@0	1856	case '+':
michael@0	1857	s++;
michael@0	1858	}
michael@0	1859	if ((n = hexdig[*s]) == 0 || n > 0x19) {
michael@0	1860	s = s1;
michael@0	1861	break;
michael@0	1862	}
michael@0	1863	e1 = n - 0x10;
michael@0	1864	while((n = hexdig[*++s]) !=0 && n <= 0x19) {
michael@0	1865	if (e1 & 0xf8000000)
michael@0	1866	big = 1;
michael@0	1867	e1 = 10*e1 + n - 0x10;
michael@0	1868	}
michael@0	1869	if (esign)
michael@0	1870	e1 = -e1;
michael@0	1871	e += e1;
michael@0	1872	}
michael@0	1873	*sp = (char*)s;
michael@0	1874	if (!havedig)
michael@0	1875	*sp = (char*)s0 - 1;
michael@0	1876	if (zret)
michael@0	1877	goto retz1;
michael@0	1878	if (big) {
michael@0	1879	if (esign) {
michael@0	1880	#ifdef IEEE_Arith
michael@0	1881	switch(rounding) {
michael@0	1882	case Round_up:
michael@0	1883	if (sign)
michael@0	1884	break;
michael@0	1885	goto ret_tiny;
michael@0	1886	case Round_down:
michael@0	1887	if (!sign)
michael@0	1888	break;
michael@0	1889	goto ret_tiny;
michael@0	1890	}
michael@0	1891	#endif
michael@0	1892	goto retz;
michael@0	1893	#ifdef IEEE_Arith
michael@0	1894	ret_tiny:
michael@0	1895	#ifndef NO_ERRNO
michael@0	1896	errno = ERANGE;
michael@0	1897	#endif
michael@0	1898	word0(rvp) = 0;
michael@0	1899	word1(rvp) = 1;
michael@0	1900	return;
michael@0	1901	#endif /* IEEE_Arith */
michael@0	1902	}
michael@0	1903	switch(rounding) {
michael@0	1904	case Round_near:
michael@0	1905	goto ovfl1;
michael@0	1906	case Round_up:
michael@0	1907	if (!sign)
michael@0	1908	goto ovfl1;
michael@0	1909	goto ret_big;
michael@0	1910	case Round_down:
michael@0	1911	if (sign)
michael@0	1912	goto ovfl1;
michael@0	1913	goto ret_big;
michael@0	1914	}
michael@0	1915	ret_big:
michael@0	1916	word0(rvp) = Big0;
michael@0	1917	word1(rvp) = Big1;
michael@0	1918	return;
michael@0	1919	}
michael@0	1920	n = s1 - s0 - 1;
michael@0	1921	for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
michael@0	1922	k++;
michael@0	1923	b = Balloc(k);
michael@0	1924	x = b->x;
michael@0	1925	n = 0;
michael@0	1926	L = 0;
michael@0	1927	#ifdef USE_LOCALE
michael@0	1928	for(i = 0; decimalpoint[i+1]; ++i);
michael@0	1929	#endif
michael@0	1930	while(s1 > s0) {
michael@0	1931	#ifdef USE_LOCALE
michael@0	1932	if (*--s1 == decimalpoint[i]) {
michael@0	1933	s1 -= i;
michael@0	1934	continue;
michael@0	1935	}
michael@0	1936	#else
michael@0	1937	if (*--s1 == '.')
michael@0	1938	continue;
michael@0	1939	#endif
michael@0	1940	if (n == ULbits) {
michael@0	1941	*x++ = L;
michael@0	1942	L = 0;
michael@0	1943	n = 0;
michael@0	1944	}
michael@0	1945	L |= (hexdig[*s1] & 0x0f) << n;
michael@0	1946	n += 4;
michael@0	1947	}
michael@0	1948	*x++ = L;
michael@0	1949	b->wds = n = x - b->x;
michael@0	1950	n = ULbits*n - hi0bits(L);
michael@0	1951	nbits = Nbits;
michael@0	1952	lostbits = 0;
michael@0	1953	x = b->x;
michael@0	1954	if (n > nbits) {
michael@0	1955	n -= nbits;
michael@0	1956	if (any_on(b,n)) {
michael@0	1957	lostbits = 1;
michael@0	1958	k = n - 1;
michael@0	1959	if (x[k>>kshift] & 1 << (k & kmask)) {
michael@0	1960	lostbits = 2;
michael@0	1961	if (k > 0 && any_on(b,k))
michael@0	1962	lostbits = 3;
michael@0	1963	}
michael@0	1964	}
michael@0	1965	rshift(b, n);
michael@0	1966	e += n;
michael@0	1967	}
michael@0	1968	else if (n < nbits) {
michael@0	1969	n = nbits - n;
michael@0	1970	b = lshift(b, n);
michael@0	1971	e -= n;
michael@0	1972	x = b->x;
michael@0	1973	}
michael@0	1974	if (e > Emax) {
michael@0	1975	ovfl:
michael@0	1976	Bfree(b);
michael@0	1977	ovfl1:
michael@0	1978	#ifndef NO_ERRNO
michael@0	1979	errno = ERANGE;
michael@0	1980	#endif
michael@0	1981	word0(rvp) = Exp_mask;
michael@0	1982	word1(rvp) = 0;
michael@0	1983	return;
michael@0	1984	}
michael@0	1985	denorm = 0;
michael@0	1986	if (e < emin) {
michael@0	1987	denorm = 1;
michael@0	1988	n = emin - e;
michael@0	1989	if (n >= nbits) {
michael@0	1990	#ifdef IEEE_Arith /*{*/
michael@0	1991	switch (rounding) {
michael@0	1992	case Round_near:
michael@0	1993	if (n == nbits && (n < 2 || any_on(b,n-1)))
michael@0	1994	goto ret_tiny;
michael@0	1995	break;
michael@0	1996	case Round_up:
michael@0	1997	if (!sign)
michael@0	1998	goto ret_tiny;
michael@0	1999	break;
michael@0	2000	case Round_down:
michael@0	2001	if (sign)
michael@0	2002	goto ret_tiny;
michael@0	2003	}
michael@0	2004	#endif /* } IEEE_Arith */
michael@0	2005	Bfree(b);
michael@0	2006	retz:
michael@0	2007	#ifndef NO_ERRNO
michael@0	2008	errno = ERANGE;
michael@0	2009	#endif
michael@0	2010	retz1:
michael@0	2011	rvp->d = 0.;
michael@0	2012	return;
michael@0	2013	}
michael@0	2014	k = n - 1;
michael@0	2015	if (lostbits)
michael@0	2016	lostbits = 1;
michael@0	2017	else if (k > 0)
michael@0	2018	lostbits = any_on(b,k);
michael@0	2019	if (x[k>>kshift] & 1 << (k & kmask))
michael@0	2020	lostbits |= 2;
michael@0	2021	nbits -= n;
michael@0	2022	rshift(b,n);
michael@0	2023	e = emin;
michael@0	2024	}
michael@0	2025	if (lostbits) {
michael@0	2026	up = 0;
michael@0	2027	switch(rounding) {
michael@0	2028	case Round_zero:
michael@0	2029	break;
michael@0	2030	case Round_near:
michael@0	2031	if (lostbits & 2
michael@0	2032	&& (lostbits & 1) | (x[0] & 1))
michael@0	2033	up = 1;
michael@0	2034	break;
michael@0	2035	case Round_up:
michael@0	2036	up = 1 - sign;
michael@0	2037	break;
michael@0	2038	case Round_down:
michael@0	2039	up = sign;
michael@0	2040	}
michael@0	2041	if (up) {
michael@0	2042	k = b->wds;
michael@0	2043	b = increment(b);
michael@0	2044	x = b->x;
michael@0	2045	if (denorm) {
michael@0	2046	#if 0
michael@0	2047	if (nbits == Nbits - 1
michael@0	2048	&& x[nbits >> kshift] & 1 << (nbits & kmask))
michael@0	2049	denorm = 0; /* not currently used */
michael@0	2050	#endif
michael@0	2051	}
michael@0	2052	else if (b->wds > k
michael@0	2053	|| ((n = nbits & kmask) !=0
michael@0	2054	&& hi0bits(x[k-1]) < 32-n)) {
michael@0	2055	rshift(b,1);
michael@0	2056	if (++e > Emax)
michael@0	2057	goto ovfl;
michael@0	2058	}
michael@0	2059	}
michael@0	2060	}
michael@0	2061	#ifdef IEEE_Arith
michael@0	2062	if (denorm)
michael@0	2063	word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
michael@0	2064	else
michael@0	2065	word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
michael@0	2066	word1(rvp) = b->x[0];
michael@0	2067	#endif
michael@0	2068	#ifdef IBM
michael@0	2069	if ((j = e & 3)) {
michael@0	2070	k = b->x[0] & ((1 << j) - 1);
michael@0	2071	rshift(b,j);
michael@0	2072	if (k) {
michael@0	2073	switch(rounding) {
michael@0	2074	case Round_up:
michael@0	2075	if (!sign)
michael@0	2076	increment(b);
michael@0	2077	break;
michael@0	2078	case Round_down:
michael@0	2079	if (sign)
michael@0	2080	increment(b);
michael@0	2081	break;
michael@0	2082	case Round_near:
michael@0	2083	j = 1 << (j-1);
michael@0	2084	if (k & j && ((k & (j-1)) | lostbits))
michael@0	2085	increment(b);
michael@0	2086	}
michael@0	2087	}
michael@0	2088	}
michael@0	2089	e >>= 2;
michael@0	2090	word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
michael@0	2091	word1(rvp) = b->x[0];
michael@0	2092	#endif
michael@0	2093	#ifdef VAX
michael@0	2094	/* The next two lines ignore swap of low- and high-order 2 bytes. */
michael@0	2095	/* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
michael@0	2096	/* word1(rvp) = b->x[0]; */
michael@0	2097	word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
michael@0	2098	word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
michael@0	2099	#endif
michael@0	2100	Bfree(b);
michael@0	2101	}
michael@0	2102	#endif /*!NO_HEX_FP}*/
michael@0	2103
michael@0	2104	static int
michael@0	2105	#ifdef KR_headers
michael@0	2106	dshift(b, p2) Bigint *b; int p2;
michael@0	2107	#else
michael@0	2108	dshift(Bigint *b, int p2)
michael@0	2109	#endif
michael@0	2110	{
michael@0	2111	int rv = hi0bits(b->x[b->wds-1]) - 4;
michael@0	2112	if (p2 > 0)
michael@0	2113	rv -= p2;
michael@0	2114	return rv & kmask;
michael@0	2115	}
michael@0	2116
michael@0	2117	static int
michael@0	2118	quorem
michael@0	2119	#ifdef KR_headers
michael@0	2120	(b, S) Bigint *b, *S;
michael@0	2121	#else
michael@0	2122	(Bigint *b, Bigint *S)
michael@0	2123	#endif
michael@0	2124	{
michael@0	2125	int n;
michael@0	2126	ULong *bx, *bxe, q, *sx, *sxe;
michael@0	2127	#ifdef ULLong
michael@0	2128	ULLong borrow, carry, y, ys;
michael@0	2129	#else
michael@0	2130	ULong borrow, carry, y, ys;
michael@0	2131	#ifdef Pack_32
michael@0	2132	ULong si, z, zs;
michael@0	2133	#endif
michael@0	2134	#endif
michael@0	2135
michael@0	2136	n = S->wds;
michael@0	2137	#ifdef DEBUG
michael@0	2138	/*debug*/ if (b->wds > n)
michael@0	2139	/*debug*/ Bug("oversize b in quorem");
michael@0	2140	#endif
michael@0	2141	if (b->wds < n)
michael@0	2142	return 0;
michael@0	2143	sx = S->x;
michael@0	2144	sxe = sx + --n;
michael@0	2145	bx = b->x;
michael@0	2146	bxe = bx + n;
michael@0	2147	q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
michael@0	2148	#ifdef DEBUG
michael@0	2149	#ifdef NO_STRTOD_BIGCOMP
michael@0	2150	/*debug*/ if (q > 9)
michael@0	2151	#else
michael@0	2152	/* An oversized q is possible when quorem is called from bigcomp and */
michael@0	2153	/* the input is near, e.g., twice the smallest denormalized number. */
michael@0	2154	/*debug*/ if (q > 15)
michael@0	2155	#endif
michael@0	2156	/*debug*/ Bug("oversized quotient in quorem");
michael@0	2157	#endif
michael@0	2158	if (q) {
michael@0	2159	borrow = 0;
michael@0	2160	carry = 0;
michael@0	2161	do {
michael@0	2162	#ifdef ULLong
michael@0	2163	ys = *sx++ * (ULLong)q + carry;
michael@0	2164	carry = ys >> 32;
michael@0	2165	y = *bx - (ys & FFFFFFFF) - borrow;
michael@0	2166	borrow = y >> 32 & (ULong)1;
michael@0	2167	*bx++ = y & FFFFFFFF;
michael@0	2168	#else
michael@0	2169	#ifdef Pack_32
michael@0	2170	si = *sx++;
michael@0	2171	ys = (si & 0xffff) * q + carry;
michael@0	2172	zs = (si >> 16) * q + (ys >> 16);
michael@0	2173	carry = zs >> 16;
michael@0	2174	y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
michael@0	2175	borrow = (y & 0x10000) >> 16;
michael@0	2176	z = (*bx >> 16) - (zs & 0xffff) - borrow;
michael@0	2177	borrow = (z & 0x10000) >> 16;
michael@0	2178	Storeinc(bx, z, y);
michael@0	2179	#else
michael@0	2180	ys = *sx++ * q + carry;
michael@0	2181	carry = ys >> 16;
michael@0	2182	y = *bx - (ys & 0xffff) - borrow;
michael@0	2183	borrow = (y & 0x10000) >> 16;
michael@0	2184	*bx++ = y & 0xffff;
michael@0	2185	#endif
michael@0	2186	#endif
michael@0	2187	}
michael@0	2188	while(sx <= sxe);
michael@0	2189	if (!*bxe) {
michael@0	2190	bx = b->x;
michael@0	2191	while(--bxe > bx && !*bxe)
michael@0	2192	--n;
michael@0	2193	b->wds = n;
michael@0	2194	}
michael@0	2195	}
michael@0	2196	if (cmp(b, S) >= 0) {
michael@0	2197	q++;
michael@0	2198	borrow = 0;
michael@0	2199	carry = 0;
michael@0	2200	bx = b->x;
michael@0	2201	sx = S->x;
michael@0	2202	do {
michael@0	2203	#ifdef ULLong
michael@0	2204	ys = *sx++ + carry;
michael@0	2205	carry = ys >> 32;
michael@0	2206	y = *bx - (ys & FFFFFFFF) - borrow;
michael@0	2207	borrow = y >> 32 & (ULong)1;
michael@0	2208	*bx++ = y & FFFFFFFF;
michael@0	2209	#else
michael@0	2210	#ifdef Pack_32
michael@0	2211	si = *sx++;
michael@0	2212	ys = (si & 0xffff) + carry;
michael@0	2213	zs = (si >> 16) + (ys >> 16);
michael@0	2214	carry = zs >> 16;
michael@0	2215	y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
michael@0	2216	borrow = (y & 0x10000) >> 16;
michael@0	2217	z = (*bx >> 16) - (zs & 0xffff) - borrow;
michael@0	2218	borrow = (z & 0x10000) >> 16;
michael@0	2219	Storeinc(bx, z, y);
michael@0	2220	#else
michael@0	2221	ys = *sx++ + carry;
michael@0	2222	carry = ys >> 16;
michael@0	2223	y = *bx - (ys & 0xffff) - borrow;
michael@0	2224	borrow = (y & 0x10000) >> 16;
michael@0	2225	*bx++ = y & 0xffff;
michael@0	2226	#endif
michael@0	2227	#endif
michael@0	2228	}
michael@0	2229	while(sx <= sxe);
michael@0	2230	bx = b->x;
michael@0	2231	bxe = bx + n;
michael@0	2232	if (!*bxe) {
michael@0	2233	while(--bxe > bx && !*bxe)
michael@0	2234	--n;
michael@0	2235	b->wds = n;
michael@0	2236	}
michael@0	2237	}
michael@0	2238	return q;
michael@0	2239	}
michael@0	2240
michael@0	2241	#if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
michael@0	2242	static double
michael@0	2243	sulp
michael@0	2244	#ifdef KR_headers
michael@0	2245	(x, bc) U *x; BCinfo *bc;
michael@0	2246	#else
michael@0	2247	(U *x, BCinfo *bc)
michael@0	2248	#endif
michael@0	2249	{
michael@0	2250	U u;
michael@0	2251	double rv;
michael@0	2252	int i;
michael@0	2253
michael@0	2254	rv = ulp(x);
michael@0	2255	if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
michael@0	2256	return rv; /* Is there an example where i <= 0 ? */
michael@0	2257	word0(&u) = Exp_1 + (i << Exp_shift);
michael@0	2258	word1(&u) = 0;
michael@0	2259	return rv * u.d;
michael@0	2260	}
michael@0	2261	#endif /*}*/
michael@0	2262
michael@0	2263	#ifndef NO_STRTOD_BIGCOMP
michael@0	2264	static void
michael@0	2265	bigcomp
michael@0	2266	#ifdef KR_headers
michael@0	2267	(rv, s0, bc)
michael@0	2268	U *rv; CONST char *s0; BCinfo *bc;
michael@0	2269	#else
michael@0	2270	(U *rv, const char *s0, BCinfo *bc)
michael@0	2271	#endif
michael@0	2272	{
michael@0	2273	Bigint *b, *d;
michael@0	2274	int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
michael@0	2275
michael@0	2276	dsign = bc->dsign;
michael@0	2277	nd = bc->nd;
michael@0	2278	nd0 = bc->nd0;
michael@0	2279	p5 = nd + bc->e0 - 1;
michael@0	2280	speccase = 0;
michael@0	2281	#ifndef Sudden_Underflow
michael@0	2282	if (rv->d == 0.) { /* special case: value near underflow-to-zero */
michael@0	2283	/* threshold was rounded to zero */
michael@0	2284	b = i2b(1);
michael@0	2285	p2 = Emin - P + 1;
michael@0	2286	bbits = 1;
michael@0	2287	#ifdef Avoid_Underflow
michael@0	2288	word0(rv) = (P+2) << Exp_shift;
michael@0	2289	#else
michael@0	2290	word1(rv) = 1;
michael@0	2291	#endif
michael@0	2292	i = 0;
michael@0	2293	#ifdef Honor_FLT_ROUNDS
michael@0	2294	if (bc->rounding == 1)
michael@0	2295	#endif
michael@0	2296	{
michael@0	2297	speccase = 1;
michael@0	2298	--p2;
michael@0	2299	dsign = 0;
michael@0	2300	goto have_i;
michael@0	2301	}
michael@0	2302	}
michael@0	2303	else
michael@0	2304	#endif
michael@0	2305	b = d2b(rv, &p2, &bbits);
michael@0	2306	#ifdef Avoid_Underflow
michael@0	2307	p2 -= bc->scale;
michael@0	2308	#endif
michael@0	2309	/* floor(log2(rv)) == bbits - 1 + p2 */
michael@0	2310	/* Check for denormal case. */
michael@0	2311	i = P - bbits;
michael@0	2312	if (i > (j = P - Emin - 1 + p2)) {
michael@0	2313	#ifdef Sudden_Underflow
michael@0	2314	Bfree(b);
michael@0	2315	b = i2b(1);
michael@0	2316	p2 = Emin;
michael@0	2317	i = P - 1;
michael@0	2318	#ifdef Avoid_Underflow
michael@0	2319	word0(rv) = (1 + bc->scale) << Exp_shift;
michael@0	2320	#else
michael@0	2321	word0(rv) = Exp_msk1;
michael@0	2322	#endif
michael@0	2323	word1(rv) = 0;
michael@0	2324	#else
michael@0	2325	i = j;
michael@0	2326	#endif
michael@0	2327	}
michael@0	2328	#ifdef Honor_FLT_ROUNDS
michael@0	2329	if (bc->rounding != 1) {
michael@0	2330	if (i > 0)
michael@0	2331	b = lshift(b, i);
michael@0	2332	if (dsign)
michael@0	2333	b = increment(b);
michael@0	2334	}
michael@0	2335	else
michael@0	2336	#endif
michael@0	2337	{
michael@0	2338	b = lshift(b, ++i);
michael@0	2339	b->x[0] |= 1;
michael@0	2340	}
michael@0	2341	#ifndef Sudden_Underflow
michael@0	2342	have_i:
michael@0	2343	#endif
michael@0	2344	p2 -= p5 + i;
michael@0	2345	d = i2b(1);
michael@0	2346	/* Arrange for convenient computation of quotients:
michael@0	2347	* shift left if necessary so divisor has 4 leading 0 bits.
michael@0	2348	*/
michael@0	2349	if (p5 > 0)
michael@0	2350	d = pow5mult(d, p5);
michael@0	2351	else if (p5 < 0)
michael@0	2352	b = pow5mult(b, -p5);
michael@0	2353	if (p2 > 0) {
michael@0	2354	b2 = p2;
michael@0	2355	d2 = 0;
michael@0	2356	}
michael@0	2357	else {
michael@0	2358	b2 = 0;
michael@0	2359	d2 = -p2;
michael@0	2360	}
michael@0	2361	i = dshift(d, d2);
michael@0	2362	if ((b2 += i) > 0)
michael@0	2363	b = lshift(b, b2);
michael@0	2364	if ((d2 += i) > 0)
michael@0	2365	d = lshift(d, d2);
michael@0	2366
michael@0	2367	/* Now b/d = exactly half-way between the two floating-point values */
michael@0	2368	/* on either side of the input string. Compute first digit of b/d. */
michael@0	2369
michael@0	2370	if (!(dig = quorem(b,d))) {
michael@0	2371	b = multadd(b, 10, 0); /* very unlikely */
michael@0	2372	dig = quorem(b,d);
michael@0	2373	}
michael@0	2374
michael@0	2375	/* Compare b/d with s0 */
michael@0	2376
michael@0	2377	for(i = 0; i < nd0; ) {
michael@0	2378	if ((dd = s0[i++] - '0' - dig))
michael@0	2379	goto ret;
michael@0	2380	if (!b->x[0] && b->wds == 1) {
michael@0	2381	if (i < nd)
michael@0	2382	dd = 1;
michael@0	2383	goto ret;
michael@0	2384	}
michael@0	2385	b = multadd(b, 10, 0);
michael@0	2386	dig = quorem(b,d);
michael@0	2387	}
michael@0	2388	for(j = bc->dp1; i++ < nd;) {
michael@0	2389	if ((dd = s0[j++] - '0' - dig))
michael@0	2390	goto ret;
michael@0	2391	if (!b->x[0] && b->wds == 1) {
michael@0	2392	if (i < nd)
michael@0	2393	dd = 1;
michael@0	2394	goto ret;
michael@0	2395	}
michael@0	2396	b = multadd(b, 10, 0);
michael@0	2397	dig = quorem(b,d);
michael@0	2398	}
michael@0	2399	if (b->x[0] || b->wds > 1 || dig > 0)
michael@0	2400	dd = -1;
michael@0	2401	ret:
michael@0	2402	Bfree(b);
michael@0	2403	Bfree(d);
michael@0	2404	#ifdef Honor_FLT_ROUNDS
michael@0	2405	if (bc->rounding != 1) {
michael@0	2406	if (dd < 0) {
michael@0	2407	if (bc->rounding == 0) {
michael@0	2408	if (!dsign)
michael@0	2409	goto retlow1;
michael@0	2410	}
michael@0	2411	else if (dsign)
michael@0	2412	goto rethi1;
michael@0	2413	}
michael@0	2414	else if (dd > 0) {
michael@0	2415	if (bc->rounding == 0) {
michael@0	2416	if (dsign)
michael@0	2417	goto rethi1;
michael@0	2418	goto ret1;
michael@0	2419	}
michael@0	2420	if (!dsign)
michael@0	2421	goto rethi1;
michael@0	2422	dval(rv) += 2.*sulp(rv,bc);
michael@0	2423	}
michael@0	2424	else {
michael@0	2425	bc->inexact = 0;
michael@0	2426	if (dsign)
michael@0	2427	goto rethi1;
michael@0	2428	}
michael@0	2429	}
michael@0	2430	else
michael@0	2431	#endif
michael@0	2432	if (speccase) {
michael@0	2433	if (dd <= 0)
michael@0	2434	rv->d = 0.;
michael@0	2435	}
michael@0	2436	else if (dd < 0) {
michael@0	2437	if (!dsign) /* does not happen for round-near */
michael@0	2438	retlow1:
michael@0	2439	dval(rv) -= sulp(rv,bc);
michael@0	2440	}
michael@0	2441	else if (dd > 0) {
michael@0	2442	if (dsign) {
michael@0	2443	rethi1:
michael@0	2444	dval(rv) += sulp(rv,bc);
michael@0	2445	}
michael@0	2446	}
michael@0	2447	else {
michael@0	2448	/* Exact half-way case: apply round-even rule. */
michael@0	2449	if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) {
michael@0	2450	i = 1 - j;
michael@0	2451	if (i <= 31) {
michael@0	2452	if (word1(rv) & (0x1 << i))
michael@0	2453	goto odd;
michael@0	2454	}
michael@0	2455	else if (word0(rv) & (0x1 << (i-32)))
michael@0	2456	goto odd;
michael@0	2457	}
michael@0	2458	else if (word1(rv) & 1) {
michael@0	2459	odd:
michael@0	2460	if (dsign)
michael@0	2461	goto rethi1;
michael@0	2462	goto retlow1;
michael@0	2463	}
michael@0	2464	}
michael@0	2465
michael@0	2466	#ifdef Honor_FLT_ROUNDS
michael@0	2467	ret1:
michael@0	2468	#endif
michael@0	2469	return;
michael@0	2470	}
michael@0	2471	#endif /* NO_STRTOD_BIGCOMP */
michael@0	2472
michael@0	2473	double
michael@0	2474	strtod
michael@0	2475	#ifdef KR_headers
michael@0	2476	(s00, se) CONST char *s00; char **se;
michael@0	2477	#else
michael@0	2478	(const char *s00, char **se)
michael@0	2479	#endif
michael@0	2480	{
michael@0	2481	int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
michael@0	2482	int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign;
michael@0	2483	CONST char *s, *s0, *s1;
michael@0	2484	double aadj, aadj1;
michael@0	2485	Long L;
michael@0	2486	U aadj2, adj, rv, rv0;
michael@0	2487	ULong y, z;
michael@0	2488	BCinfo bc;
michael@0	2489	Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
michael@0	2490	#ifdef Avoid_Underflow
michael@0	2491	ULong Lsb, Lsb1;
michael@0	2492	#endif
michael@0	2493	#ifdef SET_INEXACT
michael@0	2494	int oldinexact;
michael@0	2495	#endif
michael@0	2496	#ifndef NO_STRTOD_BIGCOMP
michael@0	2497	int req_bigcomp = 0;
michael@0	2498	#endif
michael@0	2499	#ifdef Honor_FLT_ROUNDS /*{*/
michael@0	2500	#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
michael@0	2501	bc.rounding = Flt_Rounds;
michael@0	2502	#else /*}{*/
michael@0	2503	bc.rounding = 1;
michael@0	2504	switch(fegetround()) {
michael@0	2505	case FE_TOWARDZERO: bc.rounding = 0; break;
michael@0	2506	case FE_UPWARD: bc.rounding = 2; break;
michael@0	2507	case FE_DOWNWARD: bc.rounding = 3;
michael@0	2508	}
michael@0	2509	#endif /*}}*/
michael@0	2510	#endif /*}*/
michael@0	2511	#ifdef USE_LOCALE
michael@0	2512	CONST char *s2;
michael@0	2513	#endif
michael@0	2514
michael@0	2515	sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0;
michael@0	2516	dval(&rv) = 0.;
michael@0	2517	for(s = s00;;s++) switch(*s) {
michael@0	2518	case '-':
michael@0	2519	sign = 1;
michael@0	2520	/* no break */
michael@0	2521	case '+':
michael@0	2522	if (*++s)
michael@0	2523	goto break2;
michael@0	2524	/* no break */
michael@0	2525	case 0:
michael@0	2526	goto ret0;
michael@0	2527	case '\t':
michael@0	2528	case '\n':
michael@0	2529	case '\v':
michael@0	2530	case '\f':
michael@0	2531	case '\r':
michael@0	2532	case ' ':
michael@0	2533	continue;
michael@0	2534	default:
michael@0	2535	goto break2;
michael@0	2536	}
michael@0	2537	break2:
michael@0	2538	if (*s == '0') {
michael@0	2539	#ifndef NO_HEX_FP /*{*/
michael@0	2540	switch(s[1]) {
michael@0	2541	case 'x':
michael@0	2542	case 'X':
michael@0	2543	#ifdef Honor_FLT_ROUNDS
michael@0	2544	gethex(&s, &rv, bc.rounding, sign);
michael@0	2545	#else
michael@0	2546	gethex(&s, &rv, 1, sign);
michael@0	2547	#endif
michael@0	2548	goto ret;
michael@0	2549	}
michael@0	2550	#endif /*}*/
michael@0	2551	nz0 = 1;
michael@0	2552	while(*++s == '0') ;
michael@0	2553	if (!*s)
michael@0	2554	goto ret;
michael@0	2555	}
michael@0	2556	s0 = s;
michael@0	2557	y = z = 0;
michael@0	2558	for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
michael@0	2559	if (nd < 9)
michael@0	2560	y = 10*y + c - '0';
michael@0	2561	else if (nd < 16)
michael@0	2562	z = 10*z + c - '0';
michael@0	2563	nd0 = nd;
michael@0	2564	bc.dp0 = bc.dp1 = s - s0;
michael@0	2565	for(s1 = s; s1 > s0 && *--s1 == '0'; )
michael@0	2566	++nz1;
michael@0	2567	#ifdef USE_LOCALE
michael@0	2568	s1 = localeconv()->decimal_point;
michael@0	2569	if (c == *s1) {
michael@0	2570	c = '.';
michael@0	2571	if (*++s1) {
michael@0	2572	s2 = s;
michael@0	2573	for(;;) {
michael@0	2574	if (*++s2 != *s1) {
michael@0	2575	c = 0;
michael@0	2576	break;
michael@0	2577	}
michael@0	2578	if (!*++s1) {
michael@0	2579	s = s2;
michael@0	2580	break;
michael@0	2581	}
michael@0	2582	}
michael@0	2583	}
michael@0	2584	}
michael@0	2585	#endif
michael@0	2586	if (c == '.') {
michael@0	2587	c = *++s;
michael@0	2588	bc.dp1 = s - s0;
michael@0	2589	bc.dplen = bc.dp1 - bc.dp0;
michael@0	2590	if (!nd) {
michael@0	2591	for(; c == '0'; c = *++s)
michael@0	2592	nz++;
michael@0	2593	if (c > '0' && c <= '9') {
michael@0	2594	bc.dp0 = s0 - s;
michael@0	2595	bc.dp1 = bc.dp0 + bc.dplen;
michael@0	2596	s0 = s;
michael@0	2597	nf += nz;
michael@0	2598	nz = 0;
michael@0	2599	goto have_dig;
michael@0	2600	}
michael@0	2601	goto dig_done;
michael@0	2602	}
michael@0	2603	for(; c >= '0' && c <= '9'; c = *++s) {
michael@0	2604	have_dig:
michael@0	2605	nz++;
michael@0	2606	if (c -= '0') {
michael@0	2607	nf += nz;
michael@0	2608	for(i = 1; i < nz; i++)
michael@0	2609	if (nd++ < 9)
michael@0	2610	y *= 10;
michael@0	2611	else if (nd <= DBL_DIG + 1)
michael@0	2612	z *= 10;
michael@0	2613	if (nd++ < 9)
michael@0	2614	y = 10*y + c;
michael@0	2615	else if (nd <= DBL_DIG + 1)
michael@0	2616	z = 10*z + c;
michael@0	2617	nz = nz1 = 0;
michael@0	2618	}
michael@0	2619	}
michael@0	2620	}
michael@0	2621	dig_done:
michael@0	2622	e = 0;
michael@0	2623	if (c == 'e' || c == 'E') {
michael@0	2624	if (!nd && !nz && !nz0) {
michael@0	2625	goto ret0;
michael@0	2626	}
michael@0	2627	s00 = s;
michael@0	2628	esign = 0;
michael@0	2629	switch(c = *++s) {
michael@0	2630	case '-':
michael@0	2631	esign = 1;
michael@0	2632	case '+':
michael@0	2633	c = *++s;
michael@0	2634	}
michael@0	2635	if (c >= '0' && c <= '9') {
michael@0	2636	while(c == '0')
michael@0	2637	c = *++s;
michael@0	2638	if (c > '0' && c <= '9') {
michael@0	2639	L = c - '0';
michael@0	2640	s1 = s;
michael@0	2641	while((c = *++s) >= '0' && c <= '9')
michael@0	2642	L = 10*L + c - '0';
michael@0	2643	if (s - s1 > 8 || L > 19999)
michael@0	2644	/* Avoid confusion from exponents
michael@0	2645	* so large that e might overflow.
michael@0	2646	*/
michael@0	2647	e = 19999; /* safe for 16 bit ints */
michael@0	2648	else
michael@0	2649	e = (int)L;
michael@0	2650	if (esign)
michael@0	2651	e = -e;
michael@0	2652	}
michael@0	2653	else
michael@0	2654	e = 0;
michael@0	2655	}
michael@0	2656	else
michael@0	2657	s = s00;
michael@0	2658	}
michael@0	2659	if (!nd) {
michael@0	2660	if (!nz && !nz0) {
michael@0	2661	#ifdef INFNAN_CHECK
michael@0	2662	/* Check for Nan and Infinity */
michael@0	2663	if (!bc.dplen)
michael@0	2664	switch(c) {
michael@0	2665	case 'i':
michael@0	2666	case 'I':
michael@0	2667	if (match(&s,"nf")) {
michael@0	2668	--s;
michael@0	2669	if (!match(&s,"inity"))
michael@0	2670	++s;
michael@0	2671	word0(&rv) = 0x7ff00000;
michael@0	2672	word1(&rv) = 0;
michael@0	2673	goto ret;
michael@0	2674	}
michael@0	2675	break;
michael@0	2676	case 'n':
michael@0	2677	case 'N':
michael@0	2678	if (match(&s, "an")) {
michael@0	2679	word0(&rv) = NAN_WORD0;
michael@0	2680	word1(&rv) = NAN_WORD1;
michael@0	2681	#ifndef No_Hex_NaN
michael@0	2682	if (*s == '(') /*)*/
michael@0	2683	hexnan(&rv, &s);
michael@0	2684	#endif
michael@0	2685	goto ret;
michael@0	2686	}
michael@0	2687	}
michael@0	2688	#endif /* INFNAN_CHECK */
michael@0	2689	ret0:
michael@0	2690	s = s00;
michael@0	2691	sign = 0;
michael@0	2692	}
michael@0	2693	goto ret;
michael@0	2694	}
michael@0	2695	bc.e0 = e1 = e -= nf;
michael@0	2696
michael@0	2697	/* Now we have nd0 digits, starting at s0, followed by a
michael@0	2698	* decimal point, followed by nd-nd0 digits. The number we're
michael@0	2699	* after is the integer represented by those digits times
michael@0	2700	* 10**e */
michael@0	2701
michael@0	2702	if (!nd0)
michael@0	2703	nd0 = nd;
michael@0	2704	k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
michael@0	2705	dval(&rv) = y;
michael@0	2706	if (k > 9) {
michael@0	2707	#ifdef SET_INEXACT
michael@0	2708	if (k > DBL_DIG)
michael@0	2709	oldinexact = get_inexact();
michael@0	2710	#endif
michael@0	2711	dval(&rv) = tens[k - 9] * dval(&rv) + z;
michael@0	2712	}
michael@0	2713	bd0 = 0;
michael@0	2714	if (nd <= DBL_DIG
michael@0	2715	#ifndef RND_PRODQUOT
michael@0	2716	#ifndef Honor_FLT_ROUNDS
michael@0	2717	&& Flt_Rounds == 1
michael@0	2718	#endif
michael@0	2719	#endif
michael@0	2720	) {
michael@0	2721	if (!e)
michael@0	2722	goto ret;
michael@0	2723	#ifndef ROUND_BIASED_without_Round_Up
michael@0	2724	if (e > 0) {
michael@0	2725	if (e <= Ten_pmax) {
michael@0	2726	#ifdef VAX
michael@0	2727	goto vax_ovfl_check;
michael@0	2728	#else
michael@0	2729	#ifdef Honor_FLT_ROUNDS
michael@0	2730	/* round correctly FLT_ROUNDS = 2 or 3 */
michael@0	2731	if (sign) {
michael@0	2732	rv.d = -rv.d;
michael@0	2733	sign = 0;
michael@0	2734	}
michael@0	2735	#endif
michael@0	2736	/* rv = */ rounded_product(dval(&rv), tens[e]);
michael@0	2737	goto ret;
michael@0	2738	#endif
michael@0	2739	}
michael@0	2740	i = DBL_DIG - nd;
michael@0	2741	if (e <= Ten_pmax + i) {
michael@0	2742	/* A fancier test would sometimes let us do
michael@0	2743	* this for larger i values.
michael@0	2744	*/
michael@0	2745	#ifdef Honor_FLT_ROUNDS
michael@0	2746	/* round correctly FLT_ROUNDS = 2 or 3 */
michael@0	2747	if (sign) {
michael@0	2748	rv.d = -rv.d;
michael@0	2749	sign = 0;
michael@0	2750	}
michael@0	2751	#endif
michael@0	2752	e -= i;
michael@0	2753	dval(&rv) *= tens[i];
michael@0	2754	#ifdef VAX
michael@0	2755	/* VAX exponent range is so narrow we must
michael@0	2756	* worry about overflow here...
michael@0	2757	*/
michael@0	2758	vax_ovfl_check:
michael@0	2759	word0(&rv) -= P*Exp_msk1;
michael@0	2760	/* rv = */ rounded_product(dval(&rv), tens[e]);
michael@0	2761	if ((word0(&rv) & Exp_mask)
michael@0	2762	> Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
michael@0	2763	goto ovfl;
michael@0	2764	word0(&rv) += P*Exp_msk1;
michael@0	2765	#else
michael@0	2766	/* rv = */ rounded_product(dval(&rv), tens[e]);
michael@0	2767	#endif
michael@0	2768	goto ret;
michael@0	2769	}
michael@0	2770	}
michael@0	2771	#ifndef Inaccurate_Divide
michael@0	2772	else if (e >= -Ten_pmax) {
michael@0	2773	#ifdef Honor_FLT_ROUNDS
michael@0	2774	/* round correctly FLT_ROUNDS = 2 or 3 */
michael@0	2775	if (sign) {
michael@0	2776	rv.d = -rv.d;
michael@0	2777	sign = 0;
michael@0	2778	}
michael@0	2779	#endif
michael@0	2780	/* rv = */ rounded_quotient(dval(&rv), tens[-e]);
michael@0	2781	goto ret;
michael@0	2782	}
michael@0	2783	#endif
michael@0	2784	#endif /* ROUND_BIASED_without_Round_Up */
michael@0	2785	}
michael@0	2786	e1 += nd - k;
michael@0	2787
michael@0	2788	#ifdef IEEE_Arith
michael@0	2789	#ifdef SET_INEXACT
michael@0	2790	bc.inexact = 1;
michael@0	2791	if (k <= DBL_DIG)
michael@0	2792	oldinexact = get_inexact();
michael@0	2793	#endif
michael@0	2794	#ifdef Avoid_Underflow
michael@0	2795	bc.scale = 0;
michael@0	2796	#endif
michael@0	2797	#ifdef Honor_FLT_ROUNDS
michael@0	2798	if (bc.rounding >= 2) {
michael@0	2799	if (sign)
michael@0	2800	bc.rounding = bc.rounding == 2 ? 0 : 2;
michael@0	2801	else
michael@0	2802	if (bc.rounding != 2)
michael@0	2803	bc.rounding = 0;
michael@0	2804	}
michael@0	2805	#endif
michael@0	2806	#endif /*IEEE_Arith*/
michael@0	2807
michael@0	2808	/* Get starting approximation = rv * 10**e1 */
michael@0	2809
michael@0	2810	if (e1 > 0) {
michael@0	2811	if ((i = e1 & 15))
michael@0	2812	dval(&rv) *= tens[i];
michael@0	2813	if (e1 &= ~15) {
michael@0	2814	if (e1 > DBL_MAX_10_EXP) {
michael@0	2815	ovfl:
michael@0	2816	/* Can't trust HUGE_VAL */
michael@0	2817	#ifdef IEEE_Arith
michael@0	2818	#ifdef Honor_FLT_ROUNDS
michael@0	2819	switch(bc.rounding) {
michael@0	2820	case 0: /* toward 0 */
michael@0	2821	case 3: /* toward -infinity */
michael@0	2822	word0(&rv) = Big0;
michael@0	2823	word1(&rv) = Big1;
michael@0	2824	break;
michael@0	2825	default:
michael@0	2826	word0(&rv) = Exp_mask;
michael@0	2827	word1(&rv) = 0;
michael@0	2828	}
michael@0	2829	#else /*Honor_FLT_ROUNDS*/
michael@0	2830	word0(&rv) = Exp_mask;
michael@0	2831	word1(&rv) = 0;
michael@0	2832	#endif /*Honor_FLT_ROUNDS*/
michael@0	2833	#ifdef SET_INEXACT
michael@0	2834	/* set overflow bit */
michael@0	2835	dval(&rv0) = 1e300;
michael@0	2836	dval(&rv0) *= dval(&rv0);
michael@0	2837	#endif
michael@0	2838	#else /*IEEE_Arith*/
michael@0	2839	word0(&rv) = Big0;
michael@0	2840	word1(&rv) = Big1;
michael@0	2841	#endif /*IEEE_Arith*/
michael@0	2842	range_err:
michael@0	2843	if (bd0) {
michael@0	2844	Bfree(bb);
michael@0	2845	Bfree(bd);
michael@0	2846	Bfree(bs);
michael@0	2847	Bfree(bd0);
michael@0	2848	Bfree(delta);
michael@0	2849	}
michael@0	2850	#ifndef NO_ERRNO
michael@0	2851	errno = ERANGE;
michael@0	2852	#endif
michael@0	2853	goto ret;
michael@0	2854	}
michael@0	2855	e1 >>= 4;
michael@0	2856	for(j = 0; e1 > 1; j++, e1 >>= 1)
michael@0	2857	if (e1 & 1)
michael@0	2858	dval(&rv) *= bigtens[j];
michael@0	2859	/* The last multiplication could overflow. */
michael@0	2860	word0(&rv) -= P*Exp_msk1;
michael@0	2861	dval(&rv) *= bigtens[j];
michael@0	2862	if ((z = word0(&rv) & Exp_mask)
michael@0	2863	> Exp_msk1*(DBL_MAX_EXP+Bias-P))
michael@0	2864	goto ovfl;
michael@0	2865	if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
michael@0	2866	/* set to largest number */
michael@0	2867	/* (Can't trust DBL_MAX) */
michael@0	2868	word0(&rv) = Big0;
michael@0	2869	word1(&rv) = Big1;
michael@0	2870	}
michael@0	2871	else
michael@0	2872	word0(&rv) += P*Exp_msk1;
michael@0	2873	}
michael@0	2874	}
michael@0	2875	else if (e1 < 0) {
michael@0	2876	e1 = -e1;
michael@0	2877	if ((i = e1 & 15))
michael@0	2878	dval(&rv) /= tens[i];
michael@0	2879	if (e1 >>= 4) {
michael@0	2880	if (e1 >= 1 << n_bigtens)
michael@0	2881	goto undfl;
michael@0	2882	#ifdef Avoid_Underflow
michael@0	2883	if (e1 & Scale_Bit)
michael@0	2884	bc.scale = 2*P;
michael@0	2885	for(j = 0; e1 > 0; j++, e1 >>= 1)
michael@0	2886	if (e1 & 1)
michael@0	2887	dval(&rv) *= tinytens[j];
michael@0	2888	if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
michael@0	2889	>> Exp_shift)) > 0) {
michael@0	2890	/* scaled rv is denormal; clear j low bits */
michael@0	2891	if (j >= 32) {
michael@0	2892	if (j > 54)
michael@0	2893	goto undfl;
michael@0	2894	word1(&rv) = 0;
michael@0	2895	if (j >= 53)
michael@0	2896	word0(&rv) = (P+2)*Exp_msk1;
michael@0	2897	else
michael@0	2898	word0(&rv) &= 0xffffffff << (j-32);
michael@0	2899	}
michael@0	2900	else
michael@0	2901	word1(&rv) &= 0xffffffff << j;
michael@0	2902	}
michael@0	2903	#else
michael@0	2904	for(j = 0; e1 > 1; j++, e1 >>= 1)
michael@0	2905	if (e1 & 1)
michael@0	2906	dval(&rv) *= tinytens[j];
michael@0	2907	/* The last multiplication could underflow. */
michael@0	2908	dval(&rv0) = dval(&rv);
michael@0	2909	dval(&rv) *= tinytens[j];
michael@0	2910	if (!dval(&rv)) {
michael@0	2911	dval(&rv) = 2.*dval(&rv0);
michael@0	2912	dval(&rv) *= tinytens[j];
michael@0	2913	#endif
michael@0	2914	if (!dval(&rv)) {
michael@0	2915	undfl:
michael@0	2916	dval(&rv) = 0.;
michael@0	2917	goto range_err;
michael@0	2918	}
michael@0	2919	#ifndef Avoid_Underflow
michael@0	2920	word0(&rv) = Tiny0;
michael@0	2921	word1(&rv) = Tiny1;
michael@0	2922	/* The refinement below will clean
michael@0	2923	* this approximation up.
michael@0	2924	*/
michael@0	2925	}
michael@0	2926	#endif
michael@0	2927	}
michael@0	2928	}
michael@0	2929
michael@0	2930	/* Now the hard part -- adjusting rv to the correct value.*/
michael@0	2931
michael@0	2932	/* Put digits into bd: true value = bd * 10^e */
michael@0	2933
michael@0	2934	bc.nd = nd - nz1;
michael@0	2935	#ifndef NO_STRTOD_BIGCOMP
michael@0	2936	bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
michael@0	2937	/* to silence an erroneous warning about bc.nd0 */
michael@0	2938	/* possibly not being initialized. */
michael@0	2939	if (nd > strtod_diglim) {
michael@0	2940	/* ASSERT(strtod_diglim >= 18); 18 == one more than the */
michael@0	2941	/* minimum number of decimal digits to distinguish double values */
michael@0	2942	/* in IEEE arithmetic. */
michael@0	2943	i = j = 18;
michael@0	2944	if (i > nd0)
michael@0	2945	j += bc.dplen;
michael@0	2946	for(;;) {
michael@0	2947	if (--j < bc.dp1 && j >= bc.dp0)
michael@0	2948	j = bc.dp0 - 1;
michael@0	2949	if (s0[j] != '0')
michael@0	2950	break;
michael@0	2951	--i;
michael@0	2952	}
michael@0	2953	e += nd - i;
michael@0	2954	nd = i;
michael@0	2955	if (nd0 > nd)
michael@0	2956	nd0 = nd;
michael@0	2957	if (nd < 9) { /* must recompute y */
michael@0	2958	y = 0;
michael@0	2959	for(i = 0; i < nd0; ++i)
michael@0	2960	y = 10*y + s0[i] - '0';
michael@0	2961	for(j = bc.dp1; i < nd; ++i)
michael@0	2962	y = 10*y + s0[j++] - '0';
michael@0	2963	}
michael@0	2964	}
michael@0	2965	#endif
michael@0	2966	bd0 = s2b(s0, nd0, nd, y, bc.dplen);
michael@0	2967
michael@0	2968	for(;;) {
michael@0	2969	bd = Balloc(bd0->k);
michael@0	2970	Bcopy(bd, bd0);
michael@0	2971	bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
michael@0	2972	bs = i2b(1);
michael@0	2973
michael@0	2974	if (e >= 0) {
michael@0	2975	bb2 = bb5 = 0;
michael@0	2976	bd2 = bd5 = e;
michael@0	2977	}
michael@0	2978	else {
michael@0	2979	bb2 = bb5 = -e;
michael@0	2980	bd2 = bd5 = 0;
michael@0	2981	}
michael@0	2982	if (bbe >= 0)
michael@0	2983	bb2 += bbe;
michael@0	2984	else
michael@0	2985	bd2 -= bbe;
michael@0	2986	bs2 = bb2;
michael@0	2987	#ifdef Honor_FLT_ROUNDS
michael@0	2988	if (bc.rounding != 1)
michael@0	2989	bs2++;
michael@0	2990	#endif
michael@0	2991	#ifdef Avoid_Underflow
michael@0	2992	Lsb = LSB;
michael@0	2993	Lsb1 = 0;
michael@0	2994	j = bbe - bc.scale;
michael@0	2995	i = j + bbbits - 1; /* logb(rv) */
michael@0	2996	j = P + 1 - bbbits;
michael@0	2997	if (i < Emin) { /* denormal */
michael@0	2998	i = Emin - i;
michael@0	2999	j -= i;
michael@0	3000	if (i < 32)
michael@0	3001	Lsb <<= i;
michael@0	3002	else if (i < 52)
michael@0	3003	Lsb1 = Lsb << (i-32);
michael@0	3004	else
michael@0	3005	Lsb1 = Exp_mask;
michael@0	3006	}
michael@0	3007	#else /*Avoid_Underflow*/
michael@0	3008	#ifdef Sudden_Underflow
michael@0	3009	#ifdef IBM
michael@0	3010	j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
michael@0	3011	#else
michael@0	3012	j = P + 1 - bbbits;
michael@0	3013	#endif
michael@0	3014	#else /*Sudden_Underflow*/
michael@0	3015	j = bbe;
michael@0	3016	i = j + bbbits - 1; /* logb(rv) */
michael@0	3017	if (i < Emin) /* denormal */
michael@0	3018	j += P - Emin;
michael@0	3019	else
michael@0	3020	j = P + 1 - bbbits;
michael@0	3021	#endif /*Sudden_Underflow*/
michael@0	3022	#endif /*Avoid_Underflow*/
michael@0	3023	bb2 += j;
michael@0	3024	bd2 += j;
michael@0	3025	#ifdef Avoid_Underflow
michael@0	3026	bd2 += bc.scale;
michael@0	3027	#endif
michael@0	3028	i = bb2 < bd2 ? bb2 : bd2;
michael@0	3029	if (i > bs2)
michael@0	3030	i = bs2;
michael@0	3031	if (i > 0) {
michael@0	3032	bb2 -= i;
michael@0	3033	bd2 -= i;
michael@0	3034	bs2 -= i;
michael@0	3035	}
michael@0	3036	if (bb5 > 0) {
michael@0	3037	bs = pow5mult(bs, bb5);
michael@0	3038	bb1 = mult(bs, bb);
michael@0	3039	Bfree(bb);
michael@0	3040	bb = bb1;
michael@0	3041	}
michael@0	3042	if (bb2 > 0)
michael@0	3043	bb = lshift(bb, bb2);
michael@0	3044	if (bd5 > 0)
michael@0	3045	bd = pow5mult(bd, bd5);
michael@0	3046	if (bd2 > 0)
michael@0	3047	bd = lshift(bd, bd2);
michael@0	3048	if (bs2 > 0)
michael@0	3049	bs = lshift(bs, bs2);
michael@0	3050	delta = diff(bb, bd);
michael@0	3051	bc.dsign = delta->sign;
michael@0	3052	delta->sign = 0;
michael@0	3053	i = cmp(delta, bs);
michael@0	3054	#ifndef NO_STRTOD_BIGCOMP /*{*/
michael@0	3055	if (bc.nd > nd && i <= 0) {
michael@0	3056	if (bc.dsign) {
michael@0	3057	/* Must use bigcomp(). */
michael@0	3058	req_bigcomp = 1;
michael@0	3059	break;
michael@0	3060	}
michael@0	3061	#ifdef Honor_FLT_ROUNDS
michael@0	3062	if (bc.rounding != 1) {
michael@0	3063	if (i < 0) {
michael@0	3064	req_bigcomp = 1;
michael@0	3065	break;
michael@0	3066	}
michael@0	3067	}
michael@0	3068	else
michael@0	3069	#endif
michael@0	3070	i = -1; /* Discarded digits make delta smaller. */
michael@0	3071	}
michael@0	3072	#endif /*}*/
michael@0	3073	#ifdef Honor_FLT_ROUNDS /*{*/
michael@0	3074	if (bc.rounding != 1) {
michael@0	3075	if (i < 0) {
michael@0	3076	/* Error is less than an ulp */
michael@0	3077	if (!delta->x[0] && delta->wds <= 1) {
michael@0	3078	/* exact */
michael@0	3079	#ifdef SET_INEXACT
michael@0	3080	bc.inexact = 0;
michael@0	3081	#endif
michael@0	3082	break;
michael@0	3083	}
michael@0	3084	if (bc.rounding) {
michael@0	3085	if (bc.dsign) {
michael@0	3086	adj.d = 1.;
michael@0	3087	goto apply_adj;
michael@0	3088	}
michael@0	3089	}
michael@0	3090	else if (!bc.dsign) {
michael@0	3091	adj.d = -1.;
michael@0	3092	if (!word1(&rv)
michael@0	3093	&& !(word0(&rv) & Frac_mask)) {
michael@0	3094	y = word0(&rv) & Exp_mask;
michael@0	3095	#ifdef Avoid_Underflow
michael@0	3096	if (!bc.scale || y > 2*P*Exp_msk1)
michael@0	3097	#else
michael@0	3098	if (y)
michael@0	3099	#endif
michael@0	3100	{
michael@0	3101	delta = lshift(delta,Log2P);
michael@0	3102	if (cmp(delta, bs) <= 0)
michael@0	3103	adj.d = -0.5;
michael@0	3104	}
michael@0	3105	}
michael@0	3106	apply_adj:
michael@0	3107	#ifdef Avoid_Underflow /*{*/
michael@0	3108	if (bc.scale && (y = word0(&rv) & Exp_mask)
michael@0	3109	<= 2*P*Exp_msk1)
michael@0	3110	word0(&adj) += (2*P+1)*Exp_msk1 - y;
michael@0	3111	#else
michael@0	3112	#ifdef Sudden_Underflow
michael@0	3113	if ((word0(&rv) & Exp_mask) <=
michael@0	3114	P*Exp_msk1) {
michael@0	3115	word0(&rv) += P*Exp_msk1;
michael@0	3116	dval(&rv) += adj.d*ulp(dval(&rv));
michael@0	3117	word0(&rv) -= P*Exp_msk1;
michael@0	3118	}
michael@0	3119	else
michael@0	3120	#endif /*Sudden_Underflow*/
michael@0	3121	#endif /*Avoid_Underflow}*/
michael@0	3122	dval(&rv) += adj.d*ulp(&rv);
michael@0	3123	}
michael@0	3124	break;
michael@0	3125	}
michael@0	3126	adj.d = ratio(delta, bs);
michael@0	3127	if (adj.d < 1.)
michael@0	3128	adj.d = 1.;
michael@0	3129	if (adj.d <= 0x7ffffffe) {
michael@0	3130	/* adj = rounding ? ceil(adj) : floor(adj); */
michael@0	3131	y = adj.d;
michael@0	3132	if (y != adj.d) {
michael@0	3133	if (!((bc.rounding>>1) ^ bc.dsign))
michael@0	3134	y++;
michael@0	3135	adj.d = y;
michael@0	3136	}
michael@0	3137	}
michael@0	3138	#ifdef Avoid_Underflow /*{*/
michael@0	3139	if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
michael@0	3140	word0(&adj) += (2*P+1)*Exp_msk1 - y;
michael@0	3141	#else
michael@0	3142	#ifdef Sudden_Underflow
michael@0	3143	if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
michael@0	3144	word0(&rv) += P*Exp_msk1;
michael@0	3145	adj.d *= ulp(dval(&rv));
michael@0	3146	if (bc.dsign)
michael@0	3147	dval(&rv) += adj.d;
michael@0	3148	else
michael@0	3149	dval(&rv) -= adj.d;
michael@0	3150	word0(&rv) -= P*Exp_msk1;
michael@0	3151	goto cont;
michael@0	3152	}
michael@0	3153	#endif /*Sudden_Underflow*/
michael@0	3154	#endif /*Avoid_Underflow}*/
michael@0	3155	adj.d *= ulp(&rv);
michael@0	3156	if (bc.dsign) {
michael@0	3157	if (word0(&rv) == Big0 && word1(&rv) == Big1)
michael@0	3158	goto ovfl;
michael@0	3159	dval(&rv) += adj.d;
michael@0	3160	}
michael@0	3161	else
michael@0	3162	dval(&rv) -= adj.d;
michael@0	3163	goto cont;
michael@0	3164	}
michael@0	3165	#endif /*}Honor_FLT_ROUNDS*/
michael@0	3166
michael@0	3167	if (i < 0) {
michael@0	3168	/* Error is less than half an ulp -- check for
michael@0	3169	* special case of mantissa a power of two.
michael@0	3170	*/
michael@0	3171	if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
michael@0	3172	#ifdef IEEE_Arith /*{*/
michael@0	3173	#ifdef Avoid_Underflow
michael@0	3174	|| (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
michael@0	3175	#else
michael@0	3176	|| (word0(&rv) & Exp_mask) <= Exp_msk1
michael@0	3177	#endif
michael@0	3178	#endif /*}*/
michael@0	3179	) {
michael@0	3180	#ifdef SET_INEXACT
michael@0	3181	if (!delta->x[0] && delta->wds <= 1)
michael@0	3182	bc.inexact = 0;
michael@0	3183	#endif
michael@0	3184	break;
michael@0	3185	}
michael@0	3186	if (!delta->x[0] && delta->wds <= 1) {
michael@0	3187	/* exact result */
michael@0	3188	#ifdef SET_INEXACT
michael@0	3189	bc.inexact = 0;
michael@0	3190	#endif
michael@0	3191	break;
michael@0	3192	}
michael@0	3193	delta = lshift(delta,Log2P);
michael@0	3194	if (cmp(delta, bs) > 0)
michael@0	3195	goto drop_down;
michael@0	3196	break;
michael@0	3197	}
michael@0	3198	if (i == 0) {
michael@0	3199	/* exactly half-way between */
michael@0	3200	if (bc.dsign) {
michael@0	3201	if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
michael@0	3202	&& word1(&rv) == (
michael@0	3203	#ifdef Avoid_Underflow
michael@0	3204	(bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
michael@0	3205	? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
michael@0	3206	#endif
michael@0	3207	0xffffffff)) {
michael@0	3208	/*boundary case -- increment exponent*/
michael@0	3209	if (word0(&rv) == Big0 && word1(&rv) == Big1)
michael@0	3210	goto ovfl;
michael@0	3211	word0(&rv) = (word0(&rv) & Exp_mask)
michael@0	3212	+ Exp_msk1
michael@0	3213	#ifdef IBM
michael@0	3214	| Exp_msk1 >> 4
michael@0	3215	#endif
michael@0	3216	;
michael@0	3217	word1(&rv) = 0;
michael@0	3218	#ifdef Avoid_Underflow
michael@0	3219	bc.dsign = 0;
michael@0	3220	#endif
michael@0	3221	break;
michael@0	3222	}
michael@0	3223	}
michael@0	3224	else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
michael@0	3225	drop_down:
michael@0	3226	/* boundary case -- decrement exponent */
michael@0	3227	#ifdef Sudden_Underflow /*{{*/
michael@0	3228	L = word0(&rv) & Exp_mask;
michael@0	3229	#ifdef IBM
michael@0	3230	if (L < Exp_msk1)
michael@0	3231	#else
michael@0	3232	#ifdef Avoid_Underflow
michael@0	3233	if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
michael@0	3234	#else
michael@0	3235	if (L <= Exp_msk1)
michael@0	3236	#endif /*Avoid_Underflow*/
michael@0	3237	#endif /*IBM*/
michael@0	3238	{
michael@0	3239	if (bc.nd >nd) {
michael@0	3240	bc.uflchk = 1;
michael@0	3241	break;
michael@0	3242	}
michael@0	3243	goto undfl;
michael@0	3244	}
michael@0	3245	L -= Exp_msk1;
michael@0	3246	#else /*Sudden_Underflow}{*/
michael@0	3247	#ifdef Avoid_Underflow
michael@0	3248	if (bc.scale) {
michael@0	3249	L = word0(&rv) & Exp_mask;
michael@0	3250	if (L <= (2*P+1)*Exp_msk1) {
michael@0	3251	if (L > (P+2)*Exp_msk1)
michael@0	3252	/* round even ==> */
michael@0	3253	/* accept rv */
michael@0	3254	break;
michael@0	3255	/* rv = smallest denormal */
michael@0	3256	if (bc.nd >nd) {
michael@0	3257	bc.uflchk = 1;
michael@0	3258	break;
michael@0	3259	}
michael@0	3260	goto undfl;
michael@0	3261	}
michael@0	3262	}
michael@0	3263	#endif /*Avoid_Underflow*/
michael@0	3264	L = (word0(&rv) & Exp_mask) - Exp_msk1;
michael@0	3265	#endif /*Sudden_Underflow}}*/
michael@0	3266	word0(&rv) = L | Bndry_mask1;
michael@0	3267	word1(&rv) = 0xffffffff;
michael@0	3268	#ifdef IBM
michael@0	3269	goto cont;
michael@0	3270	#else
michael@0	3271	#ifndef NO_STRTOD_BIGCOMP
michael@0	3272	if (bc.nd > nd)
michael@0	3273	goto cont;
michael@0	3274	#endif
michael@0	3275	break;
michael@0	3276	#endif
michael@0	3277	}
michael@0	3278	#ifndef ROUND_BIASED
michael@0	3279	#ifdef Avoid_Underflow
michael@0	3280	if (Lsb1) {
michael@0	3281	if (!(word0(&rv) & Lsb1))
michael@0	3282	break;
michael@0	3283	}
michael@0	3284	else if (!(word1(&rv) & Lsb))
michael@0	3285	break;
michael@0	3286	#else
michael@0	3287	if (!(word1(&rv) & LSB))
michael@0	3288	break;
michael@0	3289	#endif
michael@0	3290	#endif
michael@0	3291	if (bc.dsign)
michael@0	3292	#ifdef Avoid_Underflow
michael@0	3293	dval(&rv) += sulp(&rv, &bc);
michael@0	3294	#else
michael@0	3295	dval(&rv) += ulp(&rv);
michael@0	3296	#endif
michael@0	3297	#ifndef ROUND_BIASED
michael@0	3298	else {
michael@0	3299	#ifdef Avoid_Underflow
michael@0	3300	dval(&rv) -= sulp(&rv, &bc);
michael@0	3301	#else
michael@0	3302	dval(&rv) -= ulp(&rv);
michael@0	3303	#endif
michael@0	3304	#ifndef Sudden_Underflow
michael@0	3305	if (!dval(&rv)) {
michael@0	3306	if (bc.nd >nd) {
michael@0	3307	bc.uflchk = 1;
michael@0	3308	break;
michael@0	3309	}
michael@0	3310	goto undfl;
michael@0	3311	}
michael@0	3312	#endif
michael@0	3313	}
michael@0	3314	#ifdef Avoid_Underflow
michael@0	3315	bc.dsign = 1 - bc.dsign;
michael@0	3316	#endif
michael@0	3317	#endif
michael@0	3318	break;
michael@0	3319	}
michael@0	3320	if ((aadj = ratio(delta, bs)) <= 2.) {
michael@0	3321	if (bc.dsign)
michael@0	3322	aadj = aadj1 = 1.;
michael@0	3323	else if (word1(&rv) || word0(&rv) & Bndry_mask) {
michael@0	3324	#ifndef Sudden_Underflow
michael@0	3325	if (word1(&rv) == Tiny1 && !word0(&rv)) {
michael@0	3326	if (bc.nd >nd) {
michael@0	3327	bc.uflchk = 1;
michael@0	3328	break;
michael@0	3329	}
michael@0	3330	goto undfl;
michael@0	3331	}
michael@0	3332	#endif
michael@0	3333	aadj = 1.;
michael@0	3334	aadj1 = -1.;
michael@0	3335	}
michael@0	3336	else {
michael@0	3337	/* special case -- power of FLT_RADIX to be */
michael@0	3338	/* rounded down... */
michael@0	3339
michael@0	3340	if (aadj < 2./FLT_RADIX)
michael@0	3341	aadj = 1./FLT_RADIX;
michael@0	3342	else
michael@0	3343	aadj *= 0.5;
michael@0	3344	aadj1 = -aadj;
michael@0	3345	}
michael@0	3346	}
michael@0	3347	else {
michael@0	3348	aadj *= 0.5;
michael@0	3349	aadj1 = bc.dsign ? aadj : -aadj;
michael@0	3350	#ifdef Check_FLT_ROUNDS
michael@0	3351	switch(bc.rounding) {
michael@0	3352	case 2: /* towards +infinity */
michael@0	3353	aadj1 -= 0.5;
michael@0	3354	break;
michael@0	3355	case 0: /* towards 0 */
michael@0	3356	case 3: /* towards -infinity */
michael@0	3357	aadj1 += 0.5;
michael@0	3358	}
michael@0	3359	#else
michael@0	3360	if (Flt_Rounds == 0)
michael@0	3361	aadj1 += 0.5;
michael@0	3362	#endif /*Check_FLT_ROUNDS*/
michael@0	3363	}
michael@0	3364	y = word0(&rv) & Exp_mask;
michael@0	3365
michael@0	3366	/* Check for overflow */
michael@0	3367
michael@0	3368	if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
michael@0	3369	dval(&rv0) = dval(&rv);
michael@0	3370	word0(&rv) -= P*Exp_msk1;
michael@0	3371	adj.d = aadj1 * ulp(&rv);
michael@0	3372	dval(&rv) += adj.d;
michael@0	3373	if ((word0(&rv) & Exp_mask) >=
michael@0	3374	Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
michael@0	3375	if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
michael@0	3376	goto ovfl;
michael@0	3377	word0(&rv) = Big0;
michael@0	3378	word1(&rv) = Big1;
michael@0	3379	goto cont;
michael@0	3380	}
michael@0	3381	else
michael@0	3382	word0(&rv) += P*Exp_msk1;
michael@0	3383	}
michael@0	3384	else {
michael@0	3385	#ifdef Avoid_Underflow
michael@0	3386	if (bc.scale && y <= 2*P*Exp_msk1) {
michael@0	3387	if (aadj <= 0x7fffffff) {
michael@0	3388	if ((z = aadj) <= 0)
michael@0	3389	z = 1;
michael@0	3390	aadj = z;
michael@0	3391	aadj1 = bc.dsign ? aadj : -aadj;
michael@0	3392	}
michael@0	3393	dval(&aadj2) = aadj1;
michael@0	3394	word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
michael@0	3395	aadj1 = dval(&aadj2);
michael@0	3396	adj.d = aadj1 * ulp(&rv);
michael@0	3397	dval(&rv) += adj.d;
michael@0	3398	if (rv.d == 0.)
michael@0	3399	#ifdef NO_STRTOD_BIGCOMP
michael@0	3400	goto undfl;
michael@0	3401	#else
michael@0	3402	{
michael@0	3403	if (bc.nd > nd)
michael@0	3404	bc.dsign = 1;
michael@0	3405	break;
michael@0	3406	}
michael@0	3407	#endif
michael@0	3408	}
michael@0	3409	else {
michael@0	3410	adj.d = aadj1 * ulp(&rv);
michael@0	3411	dval(&rv) += adj.d;
michael@0	3412	}
michael@0	3413	#else
michael@0	3414	#ifdef Sudden_Underflow
michael@0	3415	if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
michael@0	3416	dval(&rv0) = dval(&rv);
michael@0	3417	word0(&rv) += P*Exp_msk1;
michael@0	3418	adj.d = aadj1 * ulp(&rv);
michael@0	3419	dval(&rv) += adj.d;
michael@0	3420	#ifdef IBM
michael@0	3421	if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
michael@0	3422	#else
michael@0	3423	if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
michael@0	3424	#endif
michael@0	3425	{
michael@0	3426	if (word0(&rv0) == Tiny0
michael@0	3427	&& word1(&rv0) == Tiny1) {
michael@0	3428	if (bc.nd >nd) {
michael@0	3429	bc.uflchk = 1;
michael@0	3430	break;
michael@0	3431	}
michael@0	3432	goto undfl;
michael@0	3433	}
michael@0	3434	word0(&rv) = Tiny0;
michael@0	3435	word1(&rv) = Tiny1;
michael@0	3436	goto cont;
michael@0	3437	}
michael@0	3438	else
michael@0	3439	word0(&rv) -= P*Exp_msk1;
michael@0	3440	}
michael@0	3441	else {
michael@0	3442	adj.d = aadj1 * ulp(&rv);
michael@0	3443	dval(&rv) += adj.d;
michael@0	3444	}
michael@0	3445	#else /*Sudden_Underflow*/
michael@0	3446	/* Compute adj so that the IEEE rounding rules will
michael@0	3447	* correctly round rv + adj in some half-way cases.
michael@0	3448	* If rv * ulp(rv) is denormalized (i.e.,
michael@0	3449	* y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
michael@0	3450	* trouble from bits lost to denormalization;
michael@0	3451	* example: 1.2e-307 .
michael@0	3452	*/
michael@0	3453	if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
michael@0	3454	aadj1 = (double)(int)(aadj + 0.5);
michael@0	3455	if (!bc.dsign)
michael@0	3456	aadj1 = -aadj1;
michael@0	3457	}
michael@0	3458	adj.d = aadj1 * ulp(&rv);
michael@0	3459	dval(&rv) += adj.d;
michael@0	3460	#endif /*Sudden_Underflow*/
michael@0	3461	#endif /*Avoid_Underflow*/
michael@0	3462	}
michael@0	3463	z = word0(&rv) & Exp_mask;
michael@0	3464	#ifndef SET_INEXACT
michael@0	3465	if (bc.nd == nd) {
michael@0	3466	#ifdef Avoid_Underflow
michael@0	3467	if (!bc.scale)
michael@0	3468	#endif
michael@0	3469	if (y == z) {
michael@0	3470	/* Can we stop now? */
michael@0	3471	L = (Long)aadj;
michael@0	3472	aadj -= L;
michael@0	3473	/* The tolerances below are conservative. */
michael@0	3474	if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
michael@0	3475	if (aadj < .4999999 || aadj > .5000001)
michael@0	3476	break;
michael@0	3477	}
michael@0	3478	else if (aadj < .4999999/FLT_RADIX)
michael@0	3479	break;
michael@0	3480	}
michael@0	3481	}
michael@0	3482	#endif
michael@0	3483	cont:
michael@0	3484	Bfree(bb);
michael@0	3485	Bfree(bd);
michael@0	3486	Bfree(bs);
michael@0	3487	Bfree(delta);
michael@0	3488	}
michael@0	3489	Bfree(bb);
michael@0	3490	Bfree(bd);
michael@0	3491	Bfree(bs);
michael@0	3492	Bfree(bd0);
michael@0	3493	Bfree(delta);
michael@0	3494	#ifndef NO_STRTOD_BIGCOMP
michael@0	3495	if (req_bigcomp) {
michael@0	3496	bd0 = 0;
michael@0	3497	bc.e0 += nz1;
michael@0	3498	bigcomp(&rv, s0, &bc);
michael@0	3499	y = word0(&rv) & Exp_mask;
michael@0	3500	if (y == Exp_mask)
michael@0	3501	goto ovfl;
michael@0	3502	if (y == 0 && rv.d == 0.)
michael@0	3503	goto undfl;
michael@0	3504	}
michael@0	3505	#endif
michael@0	3506	#ifdef SET_INEXACT
michael@0	3507	if (bc.inexact) {
michael@0	3508	if (!oldinexact) {
michael@0	3509	word0(&rv0) = Exp_1 + (70 << Exp_shift);
michael@0	3510	word1(&rv0) = 0;
michael@0	3511	dval(&rv0) += 1.;
michael@0	3512	}
michael@0	3513	}
michael@0	3514	else if (!oldinexact)
michael@0	3515	clear_inexact();
michael@0	3516	#endif
michael@0	3517	#ifdef Avoid_Underflow
michael@0	3518	if (bc.scale) {
michael@0	3519	word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
michael@0	3520	word1(&rv0) = 0;
michael@0	3521	dval(&rv) *= dval(&rv0);
michael@0	3522	#ifndef NO_ERRNO
michael@0	3523	/* try to avoid the bug of testing an 8087 register value */
michael@0	3524	#ifdef IEEE_Arith
michael@0	3525	if (!(word0(&rv) & Exp_mask))
michael@0	3526	#else
michael@0	3527	if (word0(&rv) == 0 && word1(&rv) == 0)
michael@0	3528	#endif
michael@0	3529	errno = ERANGE;
michael@0	3530	#endif
michael@0	3531	}
michael@0	3532	#endif /* Avoid_Underflow */
michael@0	3533	#ifdef SET_INEXACT
michael@0	3534	if (bc.inexact && !(word0(&rv) & Exp_mask)) {
michael@0	3535	/* set underflow bit */
michael@0	3536	dval(&rv0) = 1e-300;
michael@0	3537	dval(&rv0) *= dval(&rv0);
michael@0	3538	}
michael@0	3539	#endif
michael@0	3540	ret:
michael@0	3541	if (se)
michael@0	3542	*se = (char *)s;
michael@0	3543	return sign ? -dval(&rv) : dval(&rv);
michael@0	3544	}
michael@0	3545
michael@0	3546	#ifndef MULTIPLE_THREADS
michael@0	3547	static char *dtoa_result;
michael@0	3548	#endif
michael@0	3549
michael@0	3550	static char *
michael@0	3551	#ifdef KR_headers
michael@0	3552	rv_alloc(i) int i;
michael@0	3553	#else
michael@0	3554	rv_alloc(int i)
michael@0	3555	#endif
michael@0	3556	{
michael@0	3557	int j, k, *r;
michael@0	3558
michael@0	3559	j = sizeof(ULong);
michael@0	3560	for(k = 0;
michael@0	3561	sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
michael@0	3562	j <<= 1)
michael@0	3563	k++;
michael@0	3564	r = (int*)Balloc(k);
michael@0	3565	*r = k;
michael@0	3566	return
michael@0	3567	#ifndef MULTIPLE_THREADS
michael@0	3568	dtoa_result =
michael@0	3569	#endif
michael@0	3570	(char *)(r+1);
michael@0	3571	}
michael@0	3572
michael@0	3573	static char *
michael@0	3574	#ifdef KR_headers
michael@0	3575	nrv_alloc(s, rve, n) char *s, **rve; int n;
michael@0	3576	#else
michael@0	3577	nrv_alloc(const char *s, char **rve, int n)
michael@0	3578	#endif
michael@0	3579	{
michael@0	3580	char *rv, *t;
michael@0	3581
michael@0	3582	t = rv = rv_alloc(n);
michael@0	3583	while((*t = *s++)) t++;
michael@0	3584	if (rve)
michael@0	3585	*rve = t;
michael@0	3586	return rv;
michael@0	3587	}
michael@0	3588
michael@0	3589	/* freedtoa(s) must be used to free values s returned by dtoa
michael@0	3590	* when MULTIPLE_THREADS is #defined. It should be used in all cases,
michael@0	3591	* but for consistency with earlier versions of dtoa, it is optional
michael@0	3592	* when MULTIPLE_THREADS is not defined.
michael@0	3593	*/
michael@0	3594
michael@0	3595	void
michael@0	3596	#ifdef KR_headers
michael@0	3597	freedtoa(s) char *s;
michael@0	3598	#else
michael@0	3599	freedtoa(char *s)
michael@0	3600	#endif
michael@0	3601	{
michael@0	3602	Bigint *b = (Bigint *)((int *)s - 1);
michael@0	3603	b->maxwds = 1 << (b->k = *(int*)b);
michael@0	3604	Bfree(b);
michael@0	3605	#ifndef MULTIPLE_THREADS
michael@0	3606	if (s == dtoa_result)
michael@0	3607	dtoa_result = 0;
michael@0	3608	#endif
michael@0	3609	}
michael@0	3610
michael@0	3611	/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
michael@0	3612	*
michael@0	3613	* Inspired by "How to Print Floating-Point Numbers Accurately" by
michael@0	3614	* Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
michael@0	3615	*
michael@0	3616	* Modifications:
michael@0	3617	* 1. Rather than iterating, we use a simple numeric overestimate
michael@0	3618	* to determine k = floor(log10(d)). We scale relevant
michael@0	3619	* quantities using O(log2(k)) rather than O(k) multiplications.
michael@0	3620	* 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
michael@0	3621	* try to generate digits strictly left to right. Instead, we
michael@0	3622	* compute with fewer bits and propagate the carry if necessary
michael@0	3623	* when rounding the final digit up. This is often faster.
michael@0	3624	* 3. Under the assumption that input will be rounded nearest,
michael@0	3625	* mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
michael@0	3626	* That is, we allow equality in stopping tests when the
michael@0	3627	* round-nearest rule will give the same floating-point value
michael@0	3628	* as would satisfaction of the stopping test with strict
michael@0	3629	* inequality.
michael@0	3630	* 4. We remove common factors of powers of 2 from relevant
michael@0	3631	* quantities.
michael@0	3632	* 5. When converting floating-point integers less than 1e16,
michael@0	3633	* we use floating-point arithmetic rather than resorting
michael@0	3634	* to multiple-precision integers.
michael@0	3635	* 6. When asked to produce fewer than 15 digits, we first try
michael@0	3636	* to get by with floating-point arithmetic; we resort to
michael@0	3637	* multiple-precision integer arithmetic only if we cannot
michael@0	3638	* guarantee that the floating-point calculation has given
michael@0	3639	* the correctly rounded result. For k requested digits and
michael@0	3640	* "uniformly" distributed input, the probability is
michael@0	3641	* something like 10^(k-15) that we must resort to the Long
michael@0	3642	* calculation.
michael@0	3643	*/
michael@0	3644
michael@0	3645	char *
michael@0	3646	dtoa
michael@0	3647	#ifdef KR_headers
michael@0	3648	(dd, mode, ndigits, decpt, sign, rve)
michael@0	3649	double dd; int mode, ndigits, *decpt, *sign; char **rve;
michael@0	3650	#else
michael@0	3651	(double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
michael@0	3652	#endif
michael@0	3653	{
michael@0	3654	/* Arguments ndigits, decpt, sign are similar to those
michael@0	3655	of ecvt and fcvt; trailing zeros are suppressed from
michael@0	3656	the returned string. If not null, *rve is set to point
michael@0	3657	to the end of the return value. If d is +-Infinity or NaN,
michael@0	3658	then *decpt is set to 9999.
michael@0	3659
michael@0	3660	mode:
michael@0	3661	0 ==> shortest string that yields d when read in
michael@0	3662	and rounded to nearest.
michael@0	3663	1 ==> like 0, but with Steele & White stopping rule;
michael@0	3664	e.g. with IEEE P754 arithmetic , mode 0 gives
michael@0	3665	1e23 whereas mode 1 gives 9.999999999999999e22.
michael@0	3666	2 ==> max(1,ndigits) significant digits. This gives a
michael@0	3667	return value similar to that of ecvt, except
michael@0	3668	that trailing zeros are suppressed.
michael@0	3669	3 ==> through ndigits past the decimal point. This
michael@0	3670	gives a return value similar to that from fcvt,
michael@0	3671	except that trailing zeros are suppressed, and
michael@0	3672	ndigits can be negative.
michael@0	3673	4,5 ==> similar to 2 and 3, respectively, but (in
michael@0	3674	round-nearest mode) with the tests of mode 0 to
michael@0	3675	possibly return a shorter string that rounds to d.
michael@0	3676	With IEEE arithmetic and compilation with
michael@0	3677	-DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
michael@0	3678	as modes 2 and 3 when FLT_ROUNDS != 1.
michael@0	3679	6-9 ==> Debugging modes similar to mode - 4: don't try
michael@0	3680	fast floating-point estimate (if applicable).
michael@0	3681
michael@0	3682	Values of mode other than 0-9 are treated as mode 0.
michael@0	3683
michael@0	3684	Sufficient space is allocated to the return value
michael@0	3685	to hold the suppressed trailing zeros.
michael@0	3686	*/
michael@0	3687
michael@0	3688	int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
michael@0	3689	j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
michael@0	3690	spec_case, try_quick;
michael@0	3691	Long L;
michael@0	3692	#ifndef Sudden_Underflow
michael@0	3693	int denorm;
michael@0	3694	ULong x;
michael@0	3695	#endif
michael@0	3696	Bigint *b, *b1, *delta, *mlo, *mhi, *S;
michael@0	3697	U d2, eps, u;
michael@0	3698	double ds;
michael@0	3699	char *s, *s0;
michael@0	3700	#ifndef No_leftright
michael@0	3701	#ifdef IEEE_Arith
michael@0	3702	U eps1;
michael@0	3703	#endif
michael@0	3704	#endif
michael@0	3705	#ifdef SET_INEXACT
michael@0	3706	int inexact, oldinexact;
michael@0	3707	#endif
michael@0	3708	#ifdef Honor_FLT_ROUNDS /*{*/
michael@0	3709	int Rounding;
michael@0	3710	#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
michael@0	3711	Rounding = Flt_Rounds;
michael@0	3712	#else /*}{*/
michael@0	3713	Rounding = 1;
michael@0	3714	switch(fegetround()) {
michael@0	3715	case FE_TOWARDZERO: Rounding = 0; break;
michael@0	3716	case FE_UPWARD: Rounding = 2; break;
michael@0	3717	case FE_DOWNWARD: Rounding = 3;
michael@0	3718	}
michael@0	3719	#endif /*}}*/
michael@0	3720	#endif /*}*/
michael@0	3721
michael@0	3722	#ifndef MULTIPLE_THREADS
michael@0	3723	if (dtoa_result) {
michael@0	3724	freedtoa(dtoa_result);
michael@0	3725	dtoa_result = 0;
michael@0	3726	}
michael@0	3727	#endif
michael@0	3728
michael@0	3729	u.d = dd;
michael@0	3730	if (word0(&u) & Sign_bit) {
michael@0	3731	/* set sign for everything, including 0's and NaNs */
michael@0	3732	*sign = 1;
michael@0	3733	word0(&u) &= ~Sign_bit; /* clear sign bit */
michael@0	3734	}
michael@0	3735	else
michael@0	3736	*sign = 0;
michael@0	3737
michael@0	3738	#if defined(IEEE_Arith) + defined(VAX)
michael@0	3739	#ifdef IEEE_Arith
michael@0	3740	if ((word0(&u) & Exp_mask) == Exp_mask)
michael@0	3741	#else
michael@0	3742	if (word0(&u) == 0x8000)
michael@0	3743	#endif
michael@0	3744	{
michael@0	3745	/* Infinity or NaN */
michael@0	3746	*decpt = 9999;
michael@0	3747	#ifdef IEEE_Arith
michael@0	3748	if (!word1(&u) && !(word0(&u) & 0xfffff))
michael@0	3749	return nrv_alloc("Infinity", rve, 8);
michael@0	3750	#endif
michael@0	3751	return nrv_alloc("NaN", rve, 3);
michael@0	3752	}
michael@0	3753	#endif
michael@0	3754	#ifdef IBM
michael@0	3755	dval(&u) += 0; /* normalize */
michael@0	3756	#endif
michael@0	3757	if (!dval(&u)) {
michael@0	3758	*decpt = 1;
michael@0	3759	return nrv_alloc("0", rve, 1);
michael@0	3760	}
michael@0	3761
michael@0	3762	#ifdef SET_INEXACT
michael@0	3763	try_quick = oldinexact = get_inexact();
michael@0	3764	inexact = 1;
michael@0	3765	#endif
michael@0	3766	#ifdef Honor_FLT_ROUNDS
michael@0	3767	if (Rounding >= 2) {
michael@0	3768	if (*sign)
michael@0	3769	Rounding = Rounding == 2 ? 0 : 2;
michael@0	3770	else
michael@0	3771	if (Rounding != 2)
michael@0	3772	Rounding = 0;
michael@0	3773	}
michael@0	3774	#endif
michael@0	3775
michael@0	3776	b = d2b(&u, &be, &bbits);
michael@0	3777	#ifdef Sudden_Underflow
michael@0	3778	i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
michael@0	3779	#else
michael@0	3780	if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
michael@0	3781	#endif
michael@0	3782	dval(&d2) = dval(&u);
michael@0	3783	word0(&d2) &= Frac_mask1;
michael@0	3784	word0(&d2) |= Exp_11;
michael@0	3785	#ifdef IBM
michael@0	3786	if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
michael@0	3787	dval(&d2) /= 1 << j;
michael@0	3788	#endif
michael@0	3789
michael@0	3790	/* log(x) ~=~ log(1.5) + (x-1.5)/1.5
michael@0	3791	* log10(x) = log(x) / log(10)
michael@0	3792	* ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
michael@0	3793	* log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
michael@0	3794	*
michael@0	3795	* This suggests computing an approximation k to log10(d) by
michael@0	3796	*
michael@0	3797	* k = (i - Bias)*0.301029995663981
michael@0	3798	* + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
michael@0	3799	*
michael@0	3800	* We want k to be too large rather than too small.
michael@0	3801	* The error in the first-order Taylor series approximation
michael@0	3802	* is in our favor, so we just round up the constant enough
michael@0	3803	* to compensate for any error in the multiplication of
michael@0	3804	* (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
michael@0	3805	* and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
michael@0	3806	* adding 1e-13 to the constant term more than suffices.
michael@0	3807	* Hence we adjust the constant term to 0.1760912590558.
michael@0	3808	* (We could get a more accurate k by invoking log10,
michael@0	3809	* but this is probably not worthwhile.)
michael@0	3810	*/
michael@0	3811
michael@0	3812	i -= Bias;
michael@0	3813	#ifdef IBM
michael@0	3814	i <<= 2;
michael@0	3815	i += j;
michael@0	3816	#endif
michael@0	3817	#ifndef Sudden_Underflow
michael@0	3818	denorm = 0;
michael@0	3819	}
michael@0	3820	else {
michael@0	3821	/* d is denormalized */
michael@0	3822
michael@0	3823	i = bbits + be + (Bias + (P-1) - 1);
michael@0	3824	x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
michael@0	3825	: word1(&u) << (32 - i);
michael@0	3826	dval(&d2) = x;
michael@0	3827	word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
michael@0	3828	i -= (Bias + (P-1) - 1) + 1;
michael@0	3829	denorm = 1;
michael@0	3830	}
michael@0	3831	#endif
michael@0	3832	ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
michael@0	3833	k = (int)ds;
michael@0	3834	if (ds < 0. && ds != k)
michael@0	3835	k--; /* want k = floor(ds) */
michael@0	3836	k_check = 1;
michael@0	3837	if (k >= 0 && k <= Ten_pmax) {
michael@0	3838	if (dval(&u) < tens[k])
michael@0	3839	k--;
michael@0	3840	k_check = 0;
michael@0	3841	}
michael@0	3842	j = bbits - i - 1;
michael@0	3843	if (j >= 0) {
michael@0	3844	b2 = 0;
michael@0	3845	s2 = j;
michael@0	3846	}
michael@0	3847	else {
michael@0	3848	b2 = -j;
michael@0	3849	s2 = 0;
michael@0	3850	}
michael@0	3851	if (k >= 0) {
michael@0	3852	b5 = 0;
michael@0	3853	s5 = k;
michael@0	3854	s2 += k;
michael@0	3855	}
michael@0	3856	else {
michael@0	3857	b2 -= k;
michael@0	3858	b5 = -k;
michael@0	3859	s5 = 0;
michael@0	3860	}
michael@0	3861	if (mode < 0 || mode > 9)
michael@0	3862	mode = 0;
michael@0	3863
michael@0	3864	#ifndef SET_INEXACT
michael@0	3865	#ifdef Check_FLT_ROUNDS
michael@0	3866	try_quick = Rounding == 1;
michael@0	3867	#else
michael@0	3868	try_quick = 1;
michael@0	3869	#endif
michael@0	3870	#endif /*SET_INEXACT*/
michael@0	3871
michael@0	3872	if (mode > 5) {
michael@0	3873	mode -= 4;
michael@0	3874	try_quick = 0;
michael@0	3875	}
michael@0	3876	leftright = 1;
michael@0	3877	ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
michael@0	3878	/* silence erroneous "gcc -Wall" warning. */
michael@0	3879	switch(mode) {
michael@0	3880	case 0:
michael@0	3881	case 1:
michael@0	3882	i = 18;
michael@0	3883	ndigits = 0;
michael@0	3884	break;
michael@0	3885	case 2:
michael@0	3886	leftright = 0;
michael@0	3887	/* no break */
michael@0	3888	case 4:
michael@0	3889	if (ndigits <= 0)
michael@0	3890	ndigits = 1;
michael@0	3891	ilim = ilim1 = i = ndigits;
michael@0	3892	break;
michael@0	3893	case 3:
michael@0	3894	leftright = 0;
michael@0	3895	/* no break */
michael@0	3896	case 5:
michael@0	3897	i = ndigits + k + 1;
michael@0	3898	ilim = i;
michael@0	3899	ilim1 = i - 1;
michael@0	3900	if (i <= 0)
michael@0	3901	i = 1;
michael@0	3902	}
michael@0	3903	s = s0 = rv_alloc(i);
michael@0	3904
michael@0	3905	#ifdef Honor_FLT_ROUNDS
michael@0	3906	if (mode > 1 && Rounding != 1)
michael@0	3907	leftright = 0;
michael@0	3908	#endif
michael@0	3909
michael@0	3910	if (ilim >= 0 && ilim <= Quick_max && try_quick) {
michael@0	3911
michael@0	3912	/* Try to get by with floating-point arithmetic. */
michael@0	3913
michael@0	3914	i = 0;
michael@0	3915	dval(&d2) = dval(&u);
michael@0	3916	k0 = k;
michael@0	3917	ilim0 = ilim;
michael@0	3918	ieps = 2; /* conservative */
michael@0	3919	if (k > 0) {
michael@0	3920	ds = tens[k&0xf];
michael@0	3921	j = k >> 4;
michael@0	3922	if (j & Bletch) {
michael@0	3923	/* prevent overflows */
michael@0	3924	j &= Bletch - 1;
michael@0	3925	dval(&u) /= bigtens[n_bigtens-1];
michael@0	3926	ieps++;
michael@0	3927	}
michael@0	3928	for(; j; j >>= 1, i++)
michael@0	3929	if (j & 1) {
michael@0	3930	ieps++;
michael@0	3931	ds *= bigtens[i];
michael@0	3932	}
michael@0	3933	dval(&u) /= ds;
michael@0	3934	}
michael@0	3935	else if ((j1 = -k)) {
michael@0	3936	dval(&u) *= tens[j1 & 0xf];
michael@0	3937	for(j = j1 >> 4; j; j >>= 1, i++)
michael@0	3938	if (j & 1) {
michael@0	3939	ieps++;
michael@0	3940	dval(&u) *= bigtens[i];
michael@0	3941	}
michael@0	3942	}
michael@0	3943	if (k_check && dval(&u) < 1. && ilim > 0) {
michael@0	3944	if (ilim1 <= 0)
michael@0	3945	goto fast_failed;
michael@0	3946	ilim = ilim1;
michael@0	3947	k--;
michael@0	3948	dval(&u) *= 10.;
michael@0	3949	ieps++;
michael@0	3950	}
michael@0	3951	dval(&eps) = ieps*dval(&u) + 7.;
michael@0	3952	word0(&eps) -= (P-1)*Exp_msk1;
michael@0	3953	if (ilim == 0) {
michael@0	3954	S = mhi = 0;
michael@0	3955	dval(&u) -= 5.;
michael@0	3956	if (dval(&u) > dval(&eps))
michael@0	3957	goto one_digit;
michael@0	3958	if (dval(&u) < -dval(&eps))
michael@0	3959	goto no_digits;
michael@0	3960	goto fast_failed;
michael@0	3961	}
michael@0	3962	#ifndef No_leftright
michael@0	3963	if (leftright) {
michael@0	3964	/* Use Steele & White method of only
michael@0	3965	* generating digits needed.
michael@0	3966	*/
michael@0	3967	dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
michael@0	3968	#ifdef IEEE_Arith
michael@0	3969	if (k0 < 0 && j1 >= 307) {
michael@0	3970	eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */
michael@0	3971	word0(&eps1) -= Exp_msk1 * (Bias+P-1);
michael@0	3972	dval(&eps1) *= tens[j1 & 0xf];
michael@0	3973	for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
michael@0	3974	if (j & 1)
michael@0	3975	dval(&eps1) *= bigtens[i];
michael@0	3976	if (eps.d < eps1.d)
michael@0	3977	eps.d = eps1.d;
michael@0	3978	}
michael@0	3979	#endif
michael@0	3980	for(i = 0;;) {
michael@0	3981	L = dval(&u);
michael@0	3982	dval(&u) -= L;
michael@0	3983	*s++ = '0' + (int)L;
michael@0	3984	if (1. - dval(&u) < dval(&eps))
michael@0	3985	goto bump_up;
michael@0	3986	if (dval(&u) < dval(&eps))
michael@0	3987	goto ret1;
michael@0	3988	if (++i >= ilim)
michael@0	3989	break;
michael@0	3990	dval(&eps) *= 10.;
michael@0	3991	dval(&u) *= 10.;
michael@0	3992	}
michael@0	3993	}
michael@0	3994	else {
michael@0	3995	#endif
michael@0	3996	/* Generate ilim digits, then fix them up. */
michael@0	3997	dval(&eps) *= tens[ilim-1];
michael@0	3998	for(i = 1;; i++, dval(&u) *= 10.) {
michael@0	3999	L = (Long)(dval(&u));
michael@0	4000	if (!(dval(&u) -= L))
michael@0	4001	ilim = i;
michael@0	4002	*s++ = '0' + (int)L;
michael@0	4003	if (i == ilim) {
michael@0	4004	if (dval(&u) > 0.5 + dval(&eps))
michael@0	4005	goto bump_up;
michael@0	4006	else if (dval(&u) < 0.5 - dval(&eps)) {
michael@0	4007	while(*--s == '0');
michael@0	4008	s++;
michael@0	4009	goto ret1;
michael@0	4010	}
michael@0	4011	break;
michael@0	4012	}
michael@0	4013	}
michael@0	4014	#ifndef No_leftright
michael@0	4015	}
michael@0	4016	#endif
michael@0	4017	fast_failed:
michael@0	4018	s = s0;
michael@0	4019	dval(&u) = dval(&d2);
michael@0	4020	k = k0;
michael@0	4021	ilim = ilim0;
michael@0	4022	}
michael@0	4023
michael@0	4024	/* Do we have a "small" integer? */
michael@0	4025
michael@0	4026	if (be >= 0 && k <= Int_max) {
michael@0	4027	/* Yes. */
michael@0	4028	ds = tens[k];
michael@0	4029	if (ndigits < 0 && ilim <= 0) {
michael@0	4030	S = mhi = 0;
michael@0	4031	if (ilim < 0 || dval(&u) <= 5*ds)
michael@0	4032	goto no_digits;
michael@0	4033	goto one_digit;
michael@0	4034	}
michael@0	4035	for(i = 1;; i++, dval(&u) *= 10.) {
michael@0	4036	L = (Long)(dval(&u) / ds);
michael@0	4037	dval(&u) -= L*ds;
michael@0	4038	#ifdef Check_FLT_ROUNDS
michael@0	4039	/* If FLT_ROUNDS == 2, L will usually be high by 1 */
michael@0	4040	if (dval(&u) < 0) {
michael@0	4041	L--;
michael@0	4042	dval(&u) += ds;
michael@0	4043	}
michael@0	4044	#endif
michael@0	4045	*s++ = '0' + (int)L;
michael@0	4046	if (!dval(&u)) {
michael@0	4047	#ifdef SET_INEXACT
michael@0	4048	inexact = 0;
michael@0	4049	#endif
michael@0	4050	break;
michael@0	4051	}
michael@0	4052	if (i == ilim) {
michael@0	4053	#ifdef Honor_FLT_ROUNDS
michael@0	4054	if (mode > 1)
michael@0	4055	switch(Rounding) {
michael@0	4056	case 0: goto ret1;
michael@0	4057	case 2: goto bump_up;
michael@0	4058	}
michael@0	4059	#endif
michael@0	4060	dval(&u) += dval(&u);
michael@0	4061	#ifdef ROUND_BIASED
michael@0	4062	if (dval(&u) >= ds)
michael@0	4063	#else
michael@0	4064	if (dval(&u) > ds || (dval(&u) == ds && L & 1))
michael@0	4065	#endif
michael@0	4066	{
michael@0	4067	bump_up:
michael@0	4068	while(*--s == '9')
michael@0	4069	if (s == s0) {
michael@0	4070	k++;
michael@0	4071	*s = '0';
michael@0	4072	break;
michael@0	4073	}
michael@0	4074	++*s++;
michael@0	4075	}
michael@0	4076	break;
michael@0	4077	}
michael@0	4078	}
michael@0	4079	goto ret1;
michael@0	4080	}
michael@0	4081
michael@0	4082	m2 = b2;
michael@0	4083	m5 = b5;
michael@0	4084	mhi = mlo = 0;
michael@0	4085	if (leftright) {
michael@0	4086	i =
michael@0	4087	#ifndef Sudden_Underflow
michael@0	4088	denorm ? be + (Bias + (P-1) - 1 + 1) :
michael@0	4089	#endif
michael@0	4090	#ifdef IBM
michael@0	4091	1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
michael@0	4092	#else
michael@0	4093	1 + P - bbits;
michael@0	4094	#endif
michael@0	4095	b2 += i;
michael@0	4096	s2 += i;
michael@0	4097	mhi = i2b(1);
michael@0	4098	}
michael@0	4099	if (m2 > 0 && s2 > 0) {
michael@0	4100	i = m2 < s2 ? m2 : s2;
michael@0	4101	b2 -= i;
michael@0	4102	m2 -= i;
michael@0	4103	s2 -= i;
michael@0	4104	}
michael@0	4105	if (b5 > 0) {
michael@0	4106	if (leftright) {
michael@0	4107	if (m5 > 0) {
michael@0	4108	mhi = pow5mult(mhi, m5);
michael@0	4109	b1 = mult(mhi, b);
michael@0	4110	Bfree(b);
michael@0	4111	b = b1;
michael@0	4112	}
michael@0	4113	if ((j = b5 - m5))
michael@0	4114	b = pow5mult(b, j);
michael@0	4115	}
michael@0	4116	else
michael@0	4117	b = pow5mult(b, b5);
michael@0	4118	}
michael@0	4119	S = i2b(1);
michael@0	4120	if (s5 > 0)
michael@0	4121	S = pow5mult(S, s5);
michael@0	4122
michael@0	4123	/* Check for special case that d is a normalized power of 2. */
michael@0	4124
michael@0	4125	spec_case = 0;
michael@0	4126	if ((mode < 2 || leftright)
michael@0	4127	#ifdef Honor_FLT_ROUNDS
michael@0	4128	&& Rounding == 1
michael@0	4129	#endif
michael@0	4130	) {
michael@0	4131	if (!word1(&u) && !(word0(&u) & Bndry_mask)
michael@0	4132	#ifndef Sudden_Underflow
michael@0	4133	&& word0(&u) & (Exp_mask & ~Exp_msk1)
michael@0	4134	#endif
michael@0	4135	) {
michael@0	4136	/* The special case */
michael@0	4137	b2 += Log2P;
michael@0	4138	s2 += Log2P;
michael@0	4139	spec_case = 1;
michael@0	4140	}
michael@0	4141	}
michael@0	4142
michael@0	4143	/* Arrange for convenient computation of quotients:
michael@0	4144	* shift left if necessary so divisor has 4 leading 0 bits.
michael@0	4145	*
michael@0	4146	* Perhaps we should just compute leading 28 bits of S once
michael@0	4147	* and for all and pass them and a shift to quorem, so it
michael@0	4148	* can do shifts and ors to compute the numerator for q.
michael@0	4149	*/
michael@0	4150	i = dshift(S, s2);
michael@0	4151	b2 += i;
michael@0	4152	m2 += i;
michael@0	4153	s2 += i;
michael@0	4154	if (b2 > 0)
michael@0	4155	b = lshift(b, b2);
michael@0	4156	if (s2 > 0)
michael@0	4157	S = lshift(S, s2);
michael@0	4158	if (k_check) {
michael@0	4159	if (cmp(b,S) < 0) {
michael@0	4160	k--;
michael@0	4161	b = multadd(b, 10, 0); /* we botched the k estimate */
michael@0	4162	if (leftright)
michael@0	4163	mhi = multadd(mhi, 10, 0);
michael@0	4164	ilim = ilim1;
michael@0	4165	}
michael@0	4166	}
michael@0	4167	if (ilim <= 0 && (mode == 3 || mode == 5)) {
michael@0	4168	if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
michael@0	4169	/* no digits, fcvt style */
michael@0	4170	no_digits:
michael@0	4171	k = -1 - ndigits;
michael@0	4172	goto ret;
michael@0	4173	}
michael@0	4174	one_digit:
michael@0	4175	*s++ = '1';
michael@0	4176	k++;
michael@0	4177	goto ret;
michael@0	4178	}
michael@0	4179	if (leftright) {
michael@0	4180	if (m2 > 0)
michael@0	4181	mhi = lshift(mhi, m2);
michael@0	4182
michael@0	4183	/* Compute mlo -- check for special case
michael@0	4184	* that d is a normalized power of 2.
michael@0	4185	*/
michael@0	4186
michael@0	4187	mlo = mhi;
michael@0	4188	if (spec_case) {
michael@0	4189	mhi = Balloc(mhi->k);
michael@0	4190	Bcopy(mhi, mlo);
michael@0	4191	mhi = lshift(mhi, Log2P);
michael@0	4192	}
michael@0	4193
michael@0	4194	for(i = 1;;i++) {
michael@0	4195	dig = quorem(b,S) + '0';
michael@0	4196	/* Do we yet have the shortest decimal string
michael@0	4197	* that will round to d?
michael@0	4198	*/
michael@0	4199	j = cmp(b, mlo);
michael@0	4200	delta = diff(S, mhi);
michael@0	4201	j1 = delta->sign ? 1 : cmp(b, delta);
michael@0	4202	Bfree(delta);
michael@0	4203	#ifndef ROUND_BIASED
michael@0	4204	if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
michael@0	4205	#ifdef Honor_FLT_ROUNDS
michael@0	4206	&& Rounding >= 1
michael@0	4207	#endif
michael@0	4208	) {
michael@0	4209	if (dig == '9')
michael@0	4210	goto round_9_up;
michael@0	4211	if (j > 0)
michael@0	4212	dig++;
michael@0	4213	#ifdef SET_INEXACT
michael@0	4214	else if (!b->x[0] && b->wds <= 1)
michael@0	4215	inexact = 0;
michael@0	4216	#endif
michael@0	4217	*s++ = dig;
michael@0	4218	goto ret;
michael@0	4219	}
michael@0	4220	#endif
michael@0	4221	if (j < 0 || (j == 0 && mode != 1
michael@0	4222	#ifndef ROUND_BIASED
michael@0	4223	&& !(word1(&u) & 1)
michael@0	4224	#endif
michael@0	4225	)) {
michael@0	4226	if (!b->x[0] && b->wds <= 1) {
michael@0	4227	#ifdef SET_INEXACT
michael@0	4228	inexact = 0;
michael@0	4229	#endif
michael@0	4230	goto accept_dig;
michael@0	4231	}
michael@0	4232	#ifdef Honor_FLT_ROUNDS
michael@0	4233	if (mode > 1)
michael@0	4234	switch(Rounding) {
michael@0	4235	case 0: goto accept_dig;
michael@0	4236	case 2: goto keep_dig;
michael@0	4237	}
michael@0	4238	#endif /*Honor_FLT_ROUNDS*/
michael@0	4239	if (j1 > 0) {
michael@0	4240	b = lshift(b, 1);
michael@0	4241	j1 = cmp(b, S);
michael@0	4242	#ifdef ROUND_BIASED
michael@0	4243	if (j1 >= 0 /*)*/
michael@0	4244	#else
michael@0	4245	if ((j1 > 0 || (j1 == 0 && dig & 1))
michael@0	4246	#endif
michael@0	4247	&& dig++ == '9')
michael@0	4248	goto round_9_up;
michael@0	4249	}
michael@0	4250	accept_dig:
michael@0	4251	*s++ = dig;
michael@0	4252	goto ret;
michael@0	4253	}
michael@0	4254	if (j1 > 0) {
michael@0	4255	#ifdef Honor_FLT_ROUNDS
michael@0	4256	if (!Rounding)
michael@0	4257	goto accept_dig;
michael@0	4258	#endif
michael@0	4259	if (dig == '9') { /* possible if i == 1 */
michael@0	4260	round_9_up:
michael@0	4261	*s++ = '9';
michael@0	4262	goto roundoff;
michael@0	4263	}
michael@0	4264	*s++ = dig + 1;
michael@0	4265	goto ret;
michael@0	4266	}
michael@0	4267	#ifdef Honor_FLT_ROUNDS
michael@0	4268	keep_dig:
michael@0	4269	#endif
michael@0	4270	*s++ = dig;
michael@0	4271	if (i == ilim)
michael@0	4272	break;
michael@0	4273	b = multadd(b, 10, 0);
michael@0	4274	if (mlo == mhi)
michael@0	4275	mlo = mhi = multadd(mhi, 10, 0);
michael@0	4276	else {
michael@0	4277	mlo = multadd(mlo, 10, 0);
michael@0	4278	mhi = multadd(mhi, 10, 0);
michael@0	4279	}
michael@0	4280	}
michael@0	4281	}
michael@0	4282	else
michael@0	4283	for(i = 1;; i++) {
michael@0	4284	*s++ = dig = quorem(b,S) + '0';
michael@0	4285	if (!b->x[0] && b->wds <= 1) {
michael@0	4286	#ifdef SET_INEXACT
michael@0	4287	inexact = 0;
michael@0	4288	#endif
michael@0	4289	goto ret;
michael@0	4290	}
michael@0	4291	if (i >= ilim)
michael@0	4292	break;
michael@0	4293	b = multadd(b, 10, 0);
michael@0	4294	}
michael@0	4295
michael@0	4296	/* Round off last digit */
michael@0	4297
michael@0	4298	#ifdef Honor_FLT_ROUNDS
michael@0	4299	switch(Rounding) {
michael@0	4300	case 0: goto trimzeros;
michael@0	4301	case 2: goto roundoff;
michael@0	4302	}
michael@0	4303	#endif
michael@0	4304	b = lshift(b, 1);
michael@0	4305	j = cmp(b, S);
michael@0	4306	#ifdef ROUND_BIASED
michael@0	4307	if (j >= 0)
michael@0	4308	#else
michael@0	4309	if (j > 0 || (j == 0 && dig & 1))
michael@0	4310	#endif
michael@0	4311	{
michael@0	4312	roundoff:
michael@0	4313	while(*--s == '9')
michael@0	4314	if (s == s0) {
michael@0	4315	k++;
michael@0	4316	*s++ = '1';
michael@0	4317	goto ret;
michael@0	4318	}
michael@0	4319	++*s++;
michael@0	4320	}
michael@0	4321	else {
michael@0	4322	#ifdef Honor_FLT_ROUNDS
michael@0	4323	trimzeros:
michael@0	4324	#endif
michael@0	4325	while(*--s == '0');
michael@0	4326	s++;
michael@0	4327	}
michael@0	4328	ret:
michael@0	4329	Bfree(S);
michael@0	4330	if (mhi) {
michael@0	4331	if (mlo && mlo != mhi)
michael@0	4332	Bfree(mlo);
michael@0	4333	Bfree(mhi);
michael@0	4334	}
michael@0	4335	ret1:
michael@0	4336	#ifdef SET_INEXACT
michael@0	4337	if (inexact) {
michael@0	4338	if (!oldinexact) {
michael@0	4339	word0(&u) = Exp_1 + (70 << Exp_shift);
michael@0	4340	word1(&u) = 0;
michael@0	4341	dval(&u) += 1.;
michael@0	4342	}
michael@0	4343	}
michael@0	4344	else if (!oldinexact)
michael@0	4345	clear_inexact();
michael@0	4346	#endif
michael@0	4347	Bfree(b);
michael@0	4348	*s = 0;
michael@0	4349	*decpt = k + 1;
michael@0	4350	if (rve)
michael@0	4351	*rve = s;
michael@0	4352	return s0;
michael@0	4353	}
michael@0	4354	#ifdef __cplusplus
michael@0	4355	}
michael@0	4356	#endif