!68 Import freebsd strtod

Merge pull request !68 from yinchuang/import_freebsd_strtod

FreeBSD.gni

@@ -27,3 +27,4 @@
 # ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 FREEBSD_SYS_LIBKERN_SRC_FILES = [ "//third_party/FreeBSD/sys/libkern/crc32.c" ]
+FREEBSD_DIR = get_path_info(".", "abspath")

OAT.xml

@@ -83,6 +83,27 @@ please copy it to your project root dir and modify it refer to OpenHarmony/tools
                     uvwxyz license textB xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
 				 " desc=""/>
             </licensematcher-->
+            <licensematcher name="BSD-3-Clause" desc="">
+                <licensetext name="
+                    * Copyright (c) 1987 Regents of the University of California.
+                    * All rights reserved.
+                    *
+                    * Redistribution and use in source and binary forms are permitted
+                    * provided that the above copyright notice and this paragraph are
+                    * duplicated in all such forms and that any documentation,
+                    * advertising materials, and other materials related to such
+                    * distribution and use acknowledge that the software was developed
+                    * by the University of California, Berkeley.  The name of the
+                    * University may not be used to endorse or promote products derived
+                    * from this software without specific prior written permission.
+                    * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
+                    * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
+                    * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+				 " desc=""/>
+            </licensematcher>
+            <licensematcher name="MIT" desc="">
+                <licensetext name="The author of this software is David M. Gay." desc=""/>
+            </licensematcher>
         </licensematcherlist>
     </oatconfig>
 </configuration>

contrib/gdtoa/gdtoa.h

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#ifndef GDTOA_H_INCLUDED
+#define GDTOA_H_INCLUDED
+
+#include "arith.h"
+#include <stddef.h> /* for size_t */
+
+#ifndef Long
+#define Long int
+#endif
+#ifndef ULong
+typedef unsigned Long ULong;
+#endif
+#ifndef UShort
+typedef unsigned short UShort;
+#endif
+
+#ifndef ANSI
+#ifdef KR_headers
+#define ANSI(x) ()
+#define Void /*nothing*/
+#else
+#define ANSI(x) x
+#define Void void
+#endif
+#endif /* ANSI */
+
+#ifndef CONST
+#ifdef KR_headers
+#define CONST /* blank */
+#else
+#define CONST const
+#endif
+#endif /* CONST */
+
+ enum {	/* return values from strtodg */
+	STRTOG_Zero	= 0,
+	STRTOG_Normal	= 1,
+	STRTOG_Denormal	= 2,
+	STRTOG_Infinite	= 3,
+	STRTOG_NaN	= 4,
+	STRTOG_NaNbits	= 5,
+	STRTOG_NoNumber	= 6,
+	STRTOG_Retmask	= 7,
+
+	/* The following may be or-ed into one of the above values. */
+
+	STRTOG_Neg	= 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
+	STRTOG_Inexlo	= 0x10,	/* returned result rounded toward zero */
+	STRTOG_Inexhi	= 0x20, /* returned result rounded away from zero */
+	STRTOG_Inexact	= 0x30,
+	STRTOG_Underflow= 0x40,
+	STRTOG_Overflow	= 0x80
+	};
+
+ typedef struct
+FPI {
+	int nbits;
+	int emin;
+	int emax;
+	int rounding;
+	int sudden_underflow;
+	} FPI;
+
+enum {	/* FPI.rounding values: same as FLT_ROUNDS */
+	FPI_Round_zero = 0,
+	FPI_Round_near = 1,
+	FPI_Round_up = 2,
+	FPI_Round_down = 3
+	};
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern char* dtoa  ANSI((double d, int mode, int ndigits, int *decpt,
+			int *sign, char **rve));
+extern char* gdtoa ANSI((FPI *fpi, int be, ULong *bits, int *kindp,
+			int mode, int ndigits, int *decpt, char **rve));
+extern void freedtoa ANSI((char*));
+extern float  strtof ANSI((CONST char *, char **));
+extern double strtod ANSI((CONST char *, char **));
+extern int strtodg ANSI((CONST char*, char**, FPI*, Long*, ULong*));
+
+extern char*	g_ddfmt  ANSI((char*, double*, int, size_t));
+extern char*	g_dfmt   ANSI((char*, double*, int, size_t));
+extern char*	g_ffmt   ANSI((char*, float*,  int, size_t));
+extern char*	g_Qfmt   ANSI((char*, void*,   int, size_t));
+extern char*	g_xfmt   ANSI((char*, void*,   int, size_t));
+extern char*	g_xLfmt  ANSI((char*, void*,   int, size_t));
+
+extern int	strtoId  ANSI((CONST char*, char**, double*, double*));
+extern int	strtoIdd ANSI((CONST char*, char**, double*, double*));
+extern int	strtoIf  ANSI((CONST char*, char**, float*, float*));
+extern int	strtoIQ  ANSI((CONST char*, char**, void*, void*));
+extern int	strtoIx  ANSI((CONST char*, char**, void*, void*));
+extern int	strtoIxL ANSI((CONST char*, char**, void*, void*));
+extern int	strtord  ANSI((CONST char*, char**, int, double*));
+extern int	strtordd ANSI((CONST char*, char**, int, double*));
+extern int	strtorf  ANSI((CONST char*, char**, int, float*));
+extern int	strtorQ  ANSI((CONST char*, char**, int, void*));
+extern int	strtorx  ANSI((CONST char*, char**, int, void*));
+extern int	strtorxL ANSI((CONST char*, char**, int, void*));
+#if 1
+extern int	strtodI  ANSI((CONST char*, char**, double*));
+extern int	strtopd  ANSI((CONST char*, char**, double*));
+extern int	strtopdd ANSI((CONST char*, char**, double*));
+extern int	strtopf  ANSI((CONST char*, char**, float*));
+extern int	strtopQ  ANSI((CONST char*, char**, void*));
+extern int	strtopx  ANSI((CONST char*, char**, void*));
+extern int	strtopxL ANSI((CONST char*, char**, void*));
+#else
+#define strtopd(s,se,x) strtord(s,se,1,x)
+#define strtopdd(s,se,x) strtordd(s,se,1,x)
+#define strtopf(s,se,x) strtorf(s,se,1,x)
+#define strtopQ(s,se,x) strtorQ(s,se,1,x)
+#define strtopx(s,se,x) strtorx(s,se,1,x)
+#define strtopxL(s,se,x) strtorxL(s,se,1,x)
+#endif
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* GDTOA_H_INCLUDED */

contrib/gdtoa/gdtoa_fltrnds.h

+	FPI *fpi, fpi1;
+	int Rounding;
+#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
+	Rounding = Flt_Rounds;
+#else /*}{*/
+	Rounding = 1;
+	switch(fegetround()) {
+	  case FE_TOWARDZERO:	Rounding = 0; break;
+	  case FE_UPWARD:	Rounding = 2; break;
+	  case FE_DOWNWARD:	Rounding = 3;
+	  }
+#endif /*}}*/
+	fpi = &fpi0;
+	if (Rounding != 1) {
+		fpi1 = fpi0;
+		fpi = &fpi1;
+		fpi1.rounding = Rounding;
+		}

contrib/gdtoa/gethex.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+#ifdef USE_LOCALE
+#include "locale.h"
+#endif
+
+ int
+#ifdef KR_headers
+gethex(sp, fpi, exp, bp, sign)
+	CONST char **sp; FPI *fpi; Long *exp; Bigint **bp; int sign;
+#else
+gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
+#endif
+{
+	Bigint *b;
+	CONST unsigned char *decpt, *s0, *s, *s1;
+	int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
+	ULong L, lostbits, *x;
+	Long e, e1;
+#ifdef USE_LOCALE
+	int i;
+#ifdef NO_LOCALE_CACHE
+	const unsigned char *decimalpoint = (unsigned char*)localeconv()->decimal_point;
+#else
+	const unsigned char *decimalpoint;
+	static unsigned char *decimalpoint_cache;
+	if (!(s0 = decimalpoint_cache)) {
+		s0 = (unsigned char*)localeconv()->decimal_point;
+		if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
+			strcpy(decimalpoint_cache, s0);
+			s0 = decimalpoint_cache;
+			}
+		}
+	decimalpoint = s0;
+#endif
+#endif
+
+	if (!hexdig['0'])
+		hexdig_init_D2A();
+	*bp = 0;
+	havedig = 0;
+	s0 = *(CONST unsigned char **)sp + 2;
+	while(s0[havedig] == '0')
+		havedig++;
+	s0 += havedig;
+	s = s0;
+	decpt = 0;
+	zret = 0;
+	e = 0;
+	if (hexdig[*s])
+		havedig++;
+	else {
+		zret = 1;
+#ifdef USE_LOCALE
+		for(i = 0; decimalpoint[i]; ++i) {
+			if (s[i] != decimalpoint[i])
+				goto pcheck;
+			}
+		decpt = s += i;
+#else
+		if (*s != '.')
+			goto pcheck;
+		decpt = ++s;
+#endif
+		if (!hexdig[*s])
+			goto pcheck;
+		while(*s == '0')
+			s++;
+		if (hexdig[*s])
+			zret = 0;
+		havedig = 1;
+		s0 = s;
+		}
+	while(hexdig[*s])
+		s++;
+#ifdef USE_LOCALE
+	if (*s == *decimalpoint && !decpt) {
+		for(i = 1; decimalpoint[i]; ++i) {
+			if (s[i] != decimalpoint[i])
+				goto pcheck;
+			}
+		decpt = s += i;
+#else
+	if (*s == '.' && !decpt) {
+		decpt = ++s;
+#endif
+		while(hexdig[*s])
+			s++;
+		}/*}*/
+	if (decpt)
+		e = -(((Long)(s-decpt)) << 2);
+ pcheck:
+	s1 = s;
+	big = esign = 0;
+	switch(*s) {
+	  case 'p':
+	  case 'P':
+		switch(*++s) {
+		  case '-':
+			esign = 1;
+			/* no break */
+		  case '+':
+			s++;
+		  }
+		if ((n = hexdig[*s]) == 0 || n > 0x19) {
+			s = s1;
+			break;
+			}
+		e1 = n - 0x10;
+		while((n = hexdig[*++s]) !=0 && n <= 0x19) {
+			if (e1 & 0xf8000000)
+				big = 1;
+			e1 = 10*e1 + n - 0x10;
+			}
+		if (esign)
+			e1 = -e1;
+		e += e1;
+	  }
+	*sp = (char*)s;
+	if (!havedig)
+		*sp = (char*)s0 - 1;
+	if (zret)
+		return STRTOG_Zero;
+	if (big) {
+		if (esign) {
+			switch(fpi->rounding) {
+			  case FPI_Round_up:
+				if (sign)
+					break;
+				goto ret_tiny;
+			  case FPI_Round_down:
+				if (!sign)
+					break;
+				goto ret_tiny;
+			  }
+			goto retz;
+ ret_tiny:
+			b = Balloc(0);
+			b->wds = 1;
+			b->x[0] = 1;
+			goto dret;
+			}
+		switch(fpi->rounding) {
+		  case FPI_Round_near:
+			goto ovfl1;
+		  case FPI_Round_up:
+			if (!sign)
+				goto ovfl1;
+			goto ret_big;
+		  case FPI_Round_down:
+			if (sign)
+				goto ovfl1;
+			goto ret_big;
+		  }
+ ret_big:
+		nbits = fpi->nbits;
+		n0 = n = nbits >> kshift;
+		if (nbits & kmask)
+			++n;
+		for(j = n, k = 0; j >>= 1; ++k);
+		*bp = b = Balloc(k);
+		b->wds = n;
+		for(j = 0; j < n0; ++j)
+			b->x[j] = ALL_ON;
+		if (n > n0)
+			b->x[j] = ULbits >> (ULbits - (nbits & kmask));
+		*exp = fpi->emin;
+		return STRTOG_Normal | STRTOG_Inexlo;
+		}
+	n = s1 - s0 - 1;
+	for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
+		k++;
+	b = Balloc(k);
+	x = b->x;
+	n = 0;
+	L = 0;
+#ifdef USE_LOCALE
+	for(i = 0; decimalpoint[i+1]; ++i);
+#endif
+	while(s1 > s0) {
+#ifdef USE_LOCALE
+		if (*--s1 == decimalpoint[i]) {
+			s1 -= i;
+			continue;
+			}
+#else
+		if (*--s1 == '.')
+			continue;
+#endif
+		if (n == ULbits) {
+			*x++ = L;
+			L = 0;
+			n = 0;
+			}
+		L |= (hexdig[*s1] & 0x0f) << n;
+		n += 4;
+		}
+	*x++ = L;
+	b->wds = n = x - b->x;
+	n = ULbits*n - hi0bits(L);
+	nbits = fpi->nbits;
+	lostbits = 0;
+	x = b->x;
+	if (n > nbits) {
+		n -= nbits;
+		if (any_on(b,n)) {
+			lostbits = 1;
+			k = n - 1;
+			if (x[k>>kshift] & 1 << (k & kmask)) {
+				lostbits = 2;
+				if (k > 0 && any_on(b,k))
+					lostbits = 3;
+				}
+			}
+		rshift(b, n);
+		e += n;
+		}
+	else if (n < nbits) {
+		n = nbits - n;
+		b = lshift(b, n);
+		e -= n;
+		x = b->x;
+		}
+	if (e > fpi->emax) {
+ ovfl:
+		Bfree(b);
+ ovfl1:
+#ifndef NO_ERRNO
+		errno = ERANGE;
+#endif
+		return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
+		}
+	irv = STRTOG_Normal;
+	if (e < fpi->emin) {
+		irv = STRTOG_Denormal;
+		n = fpi->emin - e;
+		if (n >= nbits) {
+			switch (fpi->rounding) {
+			  case FPI_Round_near:
+				if (n == nbits && (n < 2 || any_on(b,n-1)))
+					goto one_bit;
+				break;
+			  case FPI_Round_up:
+				if (!sign)
+					goto one_bit;
+				break;
+			  case FPI_Round_down:
+				if (sign) {
+ one_bit:
+					x[0] = b->wds = 1;
+ dret:
+					*bp = b;
+					*exp = fpi->emin;
+#ifndef NO_ERRNO
+					errno = ERANGE;
+#endif
+					return STRTOG_Denormal | STRTOG_Inexhi
+						| STRTOG_Underflow;
+					}
+			  }
+			Bfree(b);
+ retz:
+#ifndef NO_ERRNO
+			errno = ERANGE;
+#endif
+			return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
+			}
+		k = n - 1;
+		if (lostbits)
+			lostbits = 1;
+		else if (k > 0)
+			lostbits = any_on(b,k);
+		if (x[k>>kshift] & 1 << (k & kmask))
+			lostbits |= 2;
+		nbits -= n;
+		rshift(b,n);
+		e = fpi->emin;
+		}
+	if (lostbits) {
+		up = 0;
+		switch(fpi->rounding) {
+		  case FPI_Round_zero:
+			break;
+		  case FPI_Round_near:
+			if (lostbits & 2
+			 && (lostbits | x[0]) & 1)
+				up = 1;
+			break;
+		  case FPI_Round_up:
+			up = 1 - sign;
+			break;
+		  case FPI_Round_down:
+			up = sign;
+		  }
+		if (up) {
+			k = b->wds;
+			b = increment(b);
+			x = b->x;
+			if (irv == STRTOG_Denormal) {
+				if (nbits == fpi->nbits - 1
+				 && x[nbits >> kshift] & 1 << (nbits & kmask))
+					irv =  STRTOG_Normal;
+				}
+			else if (b->wds > k
+			 || ((n = nbits & kmask) !=0
+			      && hi0bits(x[k-1]) < 32-n)) {
+				rshift(b,1);
+				if (++e > fpi->emax)
+					goto ovfl;
+				}
+			irv |= STRTOG_Inexhi;
+			}
+		else
+			irv |= STRTOG_Inexlo;
+		}
+	*bp = b;
+	*exp = e;
+	return irv;
+	}

contrib/gdtoa/gmisc.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+ void
+#ifdef KR_headers
+rshift(b, k) Bigint *b; int k;
+#else
+rshift(Bigint *b, int k)
+#endif
+{
+	ULong *x, *x1, *xe, y;
+	int n;
+
+	x = x1 = b->x;
+	n = k >> kshift;
+	if (n < b->wds) {
+		xe = x + b->wds;
+		x += n;
+		if (k &= kmask) {
+			n = ULbits - k;
+			y = *x++ >> k;
+			while(x < xe) {
+				*x1++ = (y | (*x << n)) & ALL_ON;
+				y = *x++ >> k;
+				}
+			if ((*x1 = y) !=0)
+				x1++;
+			}
+		else
+			while(x < xe)
+				*x1++ = *x++;
+		}
+	if ((b->wds = x1 - b->x) == 0)
+		b->x[0] = 0;
+	}
+
+ int
+#ifdef KR_headers
+trailz(b) Bigint *b;
+#else
+trailz(Bigint *b)
+#endif
+{
+	ULong L, *x, *xe;
+	int n = 0;
+
+	x = b->x;
+	xe = x + b->wds;
+	for(n = 0; x < xe && !*x; x++)
+		n += ULbits;
+	if (x < xe) {
+		L = *x;
+		n += lo0bits(&L);
+		}
+	return n;
+	}

contrib/gdtoa/hd_init.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 2000 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+ unsigned char hexdig[256];
+
+ static void
+#ifdef KR_headers
+htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
+#else
+htinit(unsigned char *h, unsigned char *s, int inc)
+#endif
+{
+	int i, j;
+	for(i = 0; (j = s[i]) !=0; i++)
+		h[j] = i + inc;
+	}
+
+ void
+hexdig_init_D2A(Void)
+{
+#define USC (unsigned char *)
+	htinit(hexdig, USC "0123456789", 0x10);
+	htinit(hexdig, USC "abcdef", 0x10 + 10);
+	htinit(hexdig, USC "ABCDEF", 0x10 + 10);
+	}

contrib/gdtoa/hexnan.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 2000 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+ static void
+#ifdef KR_headers
+L_shift(x, x1, i) ULong *x; ULong *x1; int i;
+#else
+L_shift(ULong *x, ULong *x1, int i)
+#endif
+{
+	int j;
+
+	i = 8 - i;
+	i <<= 2;
+	j = ULbits - i;
+	do {
+		*x |= x[1] << j;
+		x[1] >>= i;
+		} while(++x < x1);
+	}
+
+ int
+#ifdef KR_headers
+hexnan(sp, fpi, x0)
+	CONST char **sp; FPI *fpi; ULong *x0;
+#else
+hexnan( CONST char **sp, FPI *fpi, ULong *x0)
+#endif
+{
+	ULong c, h, *x, *x1, *xe;
+	CONST char *s;
+	int havedig, hd0, i, nbits;
+
+	if (!hexdig['0'])
+		hexdig_init_D2A();
+	nbits = fpi->nbits;
+	x = x0 + (nbits >> kshift);
+	if (nbits & kmask)
+		x++;
+	*--x = 0;
+	x1 = xe = x;
+	havedig = hd0 = i = 0;
+	s = *sp;
+	/* allow optional initial 0x or 0X */
+	while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
+		++s;
+	if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
+	 && *(CONST unsigned char*)(s+3) > ' ')
+		s += 2;
+	while((c = *(CONST unsigned char*)++s)) {
+		if (!(h = hexdig[c])) {
+			if (c <= ' ') {
+				if (hd0 < havedig) {
+					if (x < x1 && i < 8)
+						L_shift(x, x1, i);
+					if (x <= x0) {
+						i = 8;
+						continue;
+						}
+					hd0 = havedig;
+					*--x = 0;
+					x1 = x;
+					i = 0;
+					}
+				while(*(CONST unsigned char*)(s+1) <= ' ')
+					++s;
+				if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
+				 && *(CONST unsigned char*)(s+3) > ' ')
+					s += 2;
+				continue;
+				}
+			if (/*(*/ c == ')' && havedig) {
+				*sp = s + 1;
+				break;
+				}
+#ifndef GDTOA_NON_PEDANTIC_NANCHECK
+			do {
+				if (/*(*/ c == ')') {
+					*sp = s + 1;
+					break;
+					}
+				} while((c = *++s));
+#endif
+			return STRTOG_NaN;
+			}
+		havedig++;
+		if (++i > 8) {
+			if (x <= x0)
+				continue;
+			i = 1;
+			*--x = 0;
+			}
+		*x = (*x << 4) | (h & 0xf);
+		}
+	if (!havedig)
+		return STRTOG_NaN;
+	if (x < x1 && i < 8)
+		L_shift(x, x1, i);
+	if (x > x0) {
+		x1 = x0;
+		do *x1++ = *x++;
+			while(x <= xe);
+		do *x1++ = 0;
+			while(x1 <= xe);
+		}
+	else {
+		/* truncate high-order word if necessary */
+		if ( (i = nbits & (ULbits-1)) !=0)
+			*xe &= ((ULong)0xffffffff) >> (ULbits - i);
+		}
+	for(x1 = xe;; --x1) {
+		if (*x1 != 0)
+			break;
+		if (x1 == x0) {
+			*x1 = 1;
+			break;
+			}
+		}
+	return STRTOG_NaNbits;
+	}

contrib/gdtoa/misc.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998, 1999 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+ static Bigint *freelist[Kmax+1];
+#ifndef Omit_Private_Memory
+#ifndef PRIVATE_MEM
+#define PRIVATE_MEM 2304
+#endif
+#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
+static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
+#endif
+
+ Bigint *
+Balloc
+#ifdef KR_headers
+	(k) int k;
+#else
+	(int k)
+#endif
+{
+	int x;
+	Bigint *rv;
+#ifndef Omit_Private_Memory
+	unsigned int len;
+#endif
+
+	ACQUIRE_DTOA_LOCK(0);
+	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
+	/* but this case seems very unlikely. */
+	if (k <= Kmax && (rv = freelist[k]) !=0) {
+		freelist[k] = rv->next;
+		}
+	else {
+		x = 1 << k;
+#ifdef Omit_Private_Memory
+		rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
+#else
+		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
+			/sizeof(double);
+		if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
+			rv = (Bigint*)pmem_next;
+			pmem_next += len;
+			}
+		else
+			rv = (Bigint*)MALLOC(len*sizeof(double));
+#endif
+		rv->k = k;
+		rv->maxwds = x;
+		}
+	FREE_DTOA_LOCK(0);
+	rv->sign = rv->wds = 0;
+	return rv;
+	}
+
+ void
+Bfree
+#ifdef KR_headers
+	(v) Bigint *v;
+#else
+	(Bigint *v)
+#endif
+{
+	if (v) {
+		if (v->k > Kmax)
+#ifdef FREE
+			FREE((void*)v);
+#else
+			free((void*)v);
+#endif
+		else {
+			ACQUIRE_DTOA_LOCK(0);
+			v->next = freelist[v->k];
+			freelist[v->k] = v;
+			FREE_DTOA_LOCK(0);
+			}
+		}
+	}
+
+ int
+lo0bits
+#ifdef KR_headers
+	(y) ULong *y;
+#else
+	(ULong *y)
+#endif
+{
+	int k;
+	ULong x = *y;
+
+	if (x & 7) {
+		if (x & 1)
+			return 0;
+		if (x & 2) {
+			*y = x >> 1;
+			return 1;
+			}
+		*y = x >> 2;
+		return 2;
+		}
+	k = 0;
+	if (!(x & 0xffff)) {
+		k = 16;
+		x >>= 16;
+		}
+	if (!(x & 0xff)) {
+		k += 8;
+		x >>= 8;
+		}
+	if (!(x & 0xf)) {
+		k += 4;
+		x >>= 4;
+		}
+	if (!(x & 0x3)) {
+		k += 2;
+		x >>= 2;
+		}
+	if (!(x & 1)) {
+		k++;
+		x >>= 1;
+		if (!x)
+			return 32;
+		}
+	*y = x;
+	return k;
+	}
+
+ Bigint *
+multadd
+#ifdef KR_headers
+	(b, m, a) Bigint *b; int m, a;
+#else
+	(Bigint *b, int m, int a)	/* multiply by m and add a */
+#endif
+{
+	int i, wds;
+#ifdef ULLong
+	ULong *x;
+	ULLong carry, y;
+#else
+	ULong carry, *x, y;
+#ifdef Pack_32
+	ULong xi, z;
+#endif
+#endif
+	Bigint *b1;
+
+	wds = b->wds;
+	x = b->x;
+	i = 0;
+	carry = a;
+	do {
+#ifdef ULLong
+		y = *x * (ULLong)m + carry;
+		carry = y >> 32;
+		*x++ = y & 0xffffffffUL;
+#else
+#ifdef Pack_32
+		xi = *x;
+		y = (xi & 0xffff) * m + carry;
+		z = (xi >> 16) * m + (y >> 16);
+		carry = z >> 16;
+		*x++ = (z << 16) + (y & 0xffff);
+#else
+		y = *x * m + carry;
+		carry = y >> 16;
+		*x++ = y & 0xffff;
+#endif
+#endif
+		}
+		while(++i < wds);
+	if (carry) {
+		if (wds >= b->maxwds) {
+			b1 = Balloc(b->k+1);
+			Bcopy(b1, b);
+			Bfree(b);
+			b = b1;
+			}
+		b->x[wds++] = carry;
+		b->wds = wds;
+		}
+	return b;
+	}
+
+ int
+hi0bits_D2A
+#ifdef KR_headers
+	(x) ULong x;
+#else
+	(ULong x)
+#endif
+{
+	int k = 0;
+
+	if (!(x & 0xffff0000)) {
+		k = 16;
+		x <<= 16;
+		}
+	if (!(x & 0xff000000)) {
+		k += 8;
+		x <<= 8;
+		}
+	if (!(x & 0xf0000000)) {
+		k += 4;
+		x <<= 4;
+		}
+	if (!(x & 0xc0000000)) {
+		k += 2;
+		x <<= 2;
+		}
+	if (!(x & 0x80000000)) {
+		k++;
+		if (!(x & 0x40000000))
+			return 32;
+		}
+	return k;
+	}
+
+ Bigint *
+i2b
+#ifdef KR_headers
+	(i) int i;
+#else
+	(int i)
+#endif
+{
+	Bigint *b;
+
+	b = Balloc(1);
+	b->x[0] = i;
+	b->wds = 1;
+	return b;
+	}
+
+ Bigint *
+mult
+#ifdef KR_headers
+	(a, b) Bigint *a, *b;
+#else
+	(Bigint *a, Bigint *b)
+#endif
+{
+	Bigint *c;
+	int k, wa, wb, wc;
+	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
+	ULong y;
+#ifdef ULLong
+	ULLong carry, z;
+#else
+	ULong carry, z;
+#ifdef Pack_32
+	ULong z2;
+#endif
+#endif
+
+	if (a->wds < b->wds) {
+		c = a;
+		a = b;
+		b = c;
+		}
+	k = a->k;
+	wa = a->wds;
+	wb = b->wds;
+	wc = wa + wb;
+	if (wc > a->maxwds)
+		k++;
+	c = Balloc(k);
+	for(x = c->x, xa = x + wc; x < xa; x++)
+		*x = 0;
+	xa = a->x;
+	xae = xa + wa;
+	xb = b->x;
+	xbe = xb + wb;
+	xc0 = c->x;
+#ifdef ULLong
+	for(; xb < xbe; xc0++) {
+		if ( (y = *xb++) !=0) {
+			x = xa;
+			xc = xc0;
+			carry = 0;
+			do {
+				z = *x++ * (ULLong)y + *xc + carry;
+				carry = z >> 32;
+				*xc++ = z & 0xffffffffUL;
+				}
+				while(x < xae);
+			*xc = carry;
+			}
+		}
+#else
+#ifdef Pack_32
+	for(; xb < xbe; xb++, xc0++) {
+		if ( (y = *xb & 0xffff) !=0) {
+			x = xa;
+			xc = xc0;
+			carry = 0;
+			do {
+				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
+				carry = z >> 16;
+				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
+				carry = z2 >> 16;
+				Storeinc(xc, z2, z);
+				}
+				while(x < xae);
+			*xc = carry;
+			}
+		if ( (y = *xb >> 16) !=0) {
+			x = xa;
+			xc = xc0;
+			carry = 0;
+			z2 = *xc;
+			do {
+				z = (*x & 0xffff) * y + (*xc >> 16) + carry;
+				carry = z >> 16;
+				Storeinc(xc, z, z2);
+				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
+				carry = z2 >> 16;
+				}
+				while(x < xae);
+			*xc = z2;
+			}
+		}
+#else
+	for(; xb < xbe; xc0++) {
+		if ( (y = *xb++) !=0) {
+			x = xa;
+			xc = xc0;
+			carry = 0;
+			do {
+				z = *x++ * y + *xc + carry;
+				carry = z >> 16;
+				*xc++ = z & 0xffff;
+				}
+				while(x < xae);
+			*xc = carry;
+			}
+		}
+#endif
+#endif
+	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
+	c->wds = wc;
+	return c;
+	}
+
+ static Bigint *p5s;
+
+ Bigint *
+pow5mult
+#ifdef KR_headers
+	(b, k) Bigint *b; int k;
+#else
+	(Bigint *b, int k)
+#endif
+{
+	Bigint *b1, *p5, *p51;
+	int i;
+	static int p05[3] = { 5, 25, 125 };
+
+	if ( (i = k & 3) !=0)
+		b = multadd(b, p05[i-1], 0);
+
+	if (!(k >>= 2))
+		return b;
+	if ((p5 = p5s) == 0) {
+		/* first time */
+#ifdef MULTIPLE_THREADS
+		ACQUIRE_DTOA_LOCK(1);
+		if (!(p5 = p5s)) {
+			p5 = p5s = i2b(625);
+			p5->next = 0;
+			}
+		FREE_DTOA_LOCK(1);
+#else
+		p5 = p5s = i2b(625);
+		p5->next = 0;
+#endif
+		}
+	for(;;) {
+		if (k & 1) {
+			b1 = mult(b, p5);
+			Bfree(b);
+			b = b1;
+			}
+		if (!(k >>= 1))
+			break;
+		if ((p51 = p5->next) == 0) {
+#ifdef MULTIPLE_THREADS
+			ACQUIRE_DTOA_LOCK(1);
+			if (!(p51 = p5->next)) {
+				p51 = p5->next = mult(p5,p5);
+				p51->next = 0;
+				}
+			FREE_DTOA_LOCK(1);
+#else
+			p51 = p5->next = mult(p5,p5);
+			p51->next = 0;
+#endif
+			}
+		p5 = p51;
+		}
+	return b;
+	}
+
+ Bigint *
+lshift
+#ifdef KR_headers
+	(b, k) Bigint *b; int k;
+#else
+	(Bigint *b, int k)
+#endif
+{
+	int i, k1, n, n1;
+	Bigint *b1;
+	ULong *x, *x1, *xe, z;
+
+	n = k >> kshift;
+	k1 = b->k;
+	n1 = n + b->wds + 1;
+	for(i = b->maxwds; n1 > i; i <<= 1)
+		k1++;
+	b1 = Balloc(k1);
+	x1 = b1->x;
+	for(i = 0; i < n; i++)
+		*x1++ = 0;
+	x = b->x;
+	xe = x + b->wds;
+	if (k &= kmask) {
+#ifdef Pack_32
+		k1 = 32 - k;
+		z = 0;
+		do {
+			*x1++ = *x << k | z;
+			z = *x++ >> k1;
+			}
+			while(x < xe);
+		if ((*x1 = z) !=0)
+			++n1;
+#else
+		k1 = 16 - k;
+		z = 0;
+		do {
+			*x1++ = *x << k  & 0xffff | z;
+			z = *x++ >> k1;
+			}
+			while(x < xe);
+		if (*x1 = z)
+			++n1;
+#endif
+		}
+	else do
+		*x1++ = *x++;
+		while(x < xe);
+	b1->wds = n1 - 1;
+	Bfree(b);
+	return b1;
+	}
+
+ int
+cmp
+#ifdef KR_headers
+	(a, b) Bigint *a, *b;
+#else
+	(Bigint *a, Bigint *b)
+#endif
+{
+	ULong *xa, *xa0, *xb, *xb0;
+	int i, j;
+
+	i = a->wds;
+	j = b->wds;
+#ifdef DEBUG
+	if (i > 1 && !a->x[i-1])
+		Bug("cmp called with a->x[a->wds-1] == 0");
+	if (j > 1 && !b->x[j-1])
+		Bug("cmp called with b->x[b->wds-1] == 0");
+#endif
+	if (i -= j)
+		return i;
+	xa0 = a->x;
+	xa = xa0 + j;
+	xb0 = b->x;
+	xb = xb0 + j;
+	for(;;) {
+		if (*--xa != *--xb)
+			return *xa < *xb ? -1 : 1;
+		if (xa <= xa0)
+			break;
+		}
+	return 0;
+	}
+
+ Bigint *
+diff
+#ifdef KR_headers
+	(a, b) Bigint *a, *b;
+#else
+	(Bigint *a, Bigint *b)
+#endif
+{
+	Bigint *c;
+	int i, wa, wb;
+	ULong *xa, *xae, *xb, *xbe, *xc;
+#ifdef ULLong
+	ULLong borrow, y;
+#else
+	ULong borrow, y;
+#ifdef Pack_32
+	ULong z;
+#endif
+#endif
+
+	i = cmp(a,b);
+	if (!i) {
+		c = Balloc(0);
+		c->wds = 1;
+		c->x[0] = 0;
+		return c;
+		}
+	if (i < 0) {
+		c = a;
+		a = b;
+		b = c;
+		i = 1;
+		}
+	else
+		i = 0;
+	c = Balloc(a->k);
+	c->sign = i;
+	wa = a->wds;
+	xa = a->x;
+	xae = xa + wa;
+	wb = b->wds;
+	xb = b->x;
+	xbe = xb + wb;
+	xc = c->x;
+	borrow = 0;
+#ifdef ULLong
+	do {
+		y = (ULLong)*xa++ - *xb++ - borrow;
+		borrow = y >> 32 & 1UL;
+		*xc++ = y & 0xffffffffUL;
+		}
+		while(xb < xbe);
+	while(xa < xae) {
+		y = *xa++ - borrow;
+		borrow = y >> 32 & 1UL;
+		*xc++ = y & 0xffffffffUL;
+		}
+#else
+#ifdef Pack_32
+	do {
+		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
+		borrow = (y & 0x10000) >> 16;
+		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
+		borrow = (z & 0x10000) >> 16;
+		Storeinc(xc, z, y);
+		}
+		while(xb < xbe);
+	while(xa < xae) {
+		y = (*xa & 0xffff) - borrow;
+		borrow = (y & 0x10000) >> 16;
+		z = (*xa++ >> 16) - borrow;
+		borrow = (z & 0x10000) >> 16;
+		Storeinc(xc, z, y);
+		}
+#else
+	do {
+		y = *xa++ - *xb++ - borrow;
+		borrow = (y & 0x10000) >> 16;
+		*xc++ = y & 0xffff;
+		}
+		while(xb < xbe);
+	while(xa < xae) {
+		y = *xa++ - borrow;
+		borrow = (y & 0x10000) >> 16;
+		*xc++ = y & 0xffff;
+		}
+#endif
+#endif
+	while(!*--xc)
+		wa--;
+	c->wds = wa;
+	return c;
+	}
+
+ double
+b2d
+#ifdef KR_headers
+	(a, e) Bigint *a; int *e;
+#else
+	(Bigint *a, int *e)
+#endif
+{
+	ULong *xa, *xa0, w, y, z;
+	int k;
+	U d;
+#ifdef VAX
+	ULong d0, d1;
+#else
+#define d0 word0(&d)
+#define d1 word1(&d)
+#endif
+
+	xa0 = a->x;
+	xa = xa0 + a->wds;
+	y = *--xa;
+#ifdef DEBUG
+	if (!y) Bug("zero y in b2d");
+#endif
+	k = hi0bits(y);
+	*e = 32 - k;
+#ifdef Pack_32
+	if (k < Ebits) {
+		d0 = Exp_1 | y >> (Ebits - k);
+		w = xa > xa0 ? *--xa : 0;
+		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
+		goto ret_d;
+		}
+	z = xa > xa0 ? *--xa : 0;
+	if (k -= Ebits) {
+		d0 = Exp_1 | y << k | z >> (32 - k);
+		y = xa > xa0 ? *--xa : 0;
+		d1 = z << k | y >> (32 - k);
+		}
+	else {
+		d0 = Exp_1 | y;
+		d1 = z;
+		}
+#else
+	if (k < Ebits + 16) {
+		z = xa > xa0 ? *--xa : 0;
+		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
+		w = xa > xa0 ? *--xa : 0;
+		y = xa > xa0 ? *--xa : 0;
+		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
+		goto ret_d;
+		}
+	z = xa > xa0 ? *--xa : 0;
+	w = xa > xa0 ? *--xa : 0;
+	k -= Ebits + 16;
+	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
+	y = xa > xa0 ? *--xa : 0;
+	d1 = w << k + 16 | y << k;
+#endif
+ ret_d:
+#ifdef VAX
+	word0(&d) = d0 >> 16 | d0 << 16;
+	word1(&d) = d1 >> 16 | d1 << 16;
+#endif
+	return dval(&d);
+	}
+#undef d0
+#undef d1
+
+ Bigint *
+d2b
+#ifdef KR_headers
+	(dd, e, bits) double dd; int *e, *bits;
+#else
+	(double dd, int *e, int *bits)
+#endif
+{
+	Bigint *b;
+	U d;
+#ifndef Sudden_Underflow
+	int i;
+#endif
+	int de, k;
+	ULong *x, y, z;
+#ifdef VAX
+	ULong d0, d1;
+#else
+#define d0 word0(&d)
+#define d1 word1(&d)
+#endif
+	d.d = dd;
+#ifdef VAX
+	d0 = word0(&d) >> 16 | word0(&d) << 16;
+	d1 = word1(&d) >> 16 | word1(&d) << 16;
+#endif
+
+#ifdef Pack_32
+	b = Balloc(1);
+#else
+	b = Balloc(2);
+#endif
+	x = b->x;
+
+	z = d0 & Frac_mask;
+	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
+#ifdef Sudden_Underflow
+	de = (int)(d0 >> Exp_shift);
+#ifndef IBM
+	z |= Exp_msk11;
+#endif
+#else
+	if ( (de = (int)(d0 >> Exp_shift)) !=0)
+		z |= Exp_msk1;
+#endif
+#ifdef Pack_32
+	if ( (y = d1) !=0) {
+		if ( (k = lo0bits(&y)) !=0) {
+			x[0] = y | z << (32 - k);
+			z >>= k;
+			}
+		else
+			x[0] = y;
+#ifndef Sudden_Underflow
+		i =
+#endif
+		     b->wds = (x[1] = z) !=0 ? 2 : 1;
+		}
+	else {
+		k = lo0bits(&z);
+		x[0] = z;
+#ifndef Sudden_Underflow
+		i =
+#endif
+		    b->wds = 1;
+		k += 32;
+		}
+#else
+	if ( (y = d1) !=0) {
+		if ( (k = lo0bits(&y)) !=0)
+			if (k >= 16) {
+				x[0] = y | z << 32 - k & 0xffff;
+				x[1] = z >> k - 16 & 0xffff;
+				x[2] = z >> k;
+				i = 2;
+				}
+			else {
+				x[0] = y & 0xffff;
+				x[1] = y >> 16 | z << 16 - k & 0xffff;
+				x[2] = z >> k & 0xffff;
+				x[3] = z >> k+16;
+				i = 3;
+				}
+		else {
+			x[0] = y & 0xffff;
+			x[1] = y >> 16;
+			x[2] = z & 0xffff;
+			x[3] = z >> 16;
+			i = 3;
+			}
+		}
+	else {
+#ifdef DEBUG
+		if (!z)
+			Bug("Zero passed to d2b");
+#endif
+		k = lo0bits(&z);
+		if (k >= 16) {
+			x[0] = z;
+			i = 0;
+			}
+		else {
+			x[0] = z & 0xffff;
+			x[1] = z >> 16;
+			i = 1;
+			}
+		k += 32;
+		}
+	while(!x[i])
+		--i;
+	b->wds = i + 1;
+#endif
+#ifndef Sudden_Underflow
+	if (de) {
+#endif
+#ifdef IBM
+		*e = (de - Bias - (P-1) << 2) + k;
+		*bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
+#else
+		*e = de - Bias - (P-1) + k;
+		*bits = P - k;
+#endif
+#ifndef Sudden_Underflow
+		}
+	else {
+		*e = de - Bias - (P-1) + 1 + k;
+#ifdef Pack_32
+		*bits = 32*i - hi0bits(x[i-1]);
+#else
+		*bits = (i+2)*16 - hi0bits(x[i]);
+#endif
+		}
+#endif
+	return b;
+	}
+#undef d0
+#undef d1
+
+ CONST double
+#ifdef IEEE_Arith
+bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
+CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
+		};
+#else
+#ifdef IBM
+bigtens[] = { 1e16, 1e32, 1e64 };
+CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
+#else
+bigtens[] = { 1e16, 1e32 };
+CONST double tinytens[] = { 1e-16, 1e-32 };
+#endif
+#endif
+
+ CONST double
+tens[] = {
+		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
+		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
+		1e20, 1e21, 1e22
+#ifdef VAX
+		, 1e23, 1e24
+#endif
+		};
+
+ char *
+#ifdef KR_headers
+strcp_D2A(a, b) char *a; char *b;
+#else
+strcp_D2A(char *a, CONST char *b)
+#endif
+{
+	while((*a = *b++))
+		a++;
+	return a;
+	}
+
+#ifdef NO_STRING_H
+
+ Char *
+#ifdef KR_headers
+memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
+#else
+memcpy_D2A(void *a1, void *b1, size_t len)
+#endif
+{
+	char *a = (char*)a1, *ae = a + len;
+	char *b = (char*)b1, *a0 = a;
+	while(a < ae)
+		*a++ = *b++;
+	return a0;
+	}
+
+#endif /* NO_STRING_H */

contrib/gdtoa/smisc.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998, 1999 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+ Bigint *
+s2b
+#ifdef KR_headers
+	(s, nd0, nd, y9, dplen) CONST char *s; int dplen, nd0, nd; ULong y9;
+#else
+	(CONST char *s, int nd0, int nd, ULong y9, int dplen)
+#endif
+{
+	Bigint *b;
+	int i, k;
+	Long x, y;
+
+	x = (nd + 8) / 9;
+	for(k = 0, y = 1; x > y; y <<= 1, k++) ;
+#ifdef Pack_32
+	b = Balloc(k);
+	b->x[0] = y9;
+	b->wds = 1;
+#else
+	b = Balloc(k+1);
+	b->x[0] = y9 & 0xffff;
+	b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
+#endif
+
+	i = 9;
+	if (9 < nd0) {
+		s += 9;
+		do b = multadd(b, 10, *s++ - '0');
+			while(++i < nd0);
+		s += dplen;
+		}
+	else
+		s += dplen + 9;
+	for(; i < nd; i++)
+		b = multadd(b, 10, *s++ - '0');
+	return b;
+	}
+
+ double
+ratio
+#ifdef KR_headers
+	(a, b) Bigint *a, *b;
+#else
+	(Bigint *a, Bigint *b)
+#endif
+{
+	U da, db;
+	int k, ka, kb;
+
+	dval(&da) = b2d(a, &ka);
+	dval(&db) = b2d(b, &kb);
+	k = ka - kb + ULbits*(a->wds - b->wds);
+#ifdef IBM
+	if (k > 0) {
+		word0(&da) += (k >> 2)*Exp_msk1;
+		if (k &= 3)
+			dval(&da) *= 1 << k;
+		}
+	else {
+		k = -k;
+		word0(&db) += (k >> 2)*Exp_msk1;
+		if (k &= 3)
+			dval(&db) *= 1 << k;
+		}
+#else
+	if (k > 0)
+		word0(&da) += k*Exp_msk1;
+	else {
+		k = -k;
+		word0(&db) += k*Exp_msk1;
+		}
+#endif
+	return dval(&da) / dval(&db);
+	}
+
+#ifdef INFNAN_CHECK
+
+ int
+match
+#ifdef KR_headers
+	(sp, t) char **sp, *t;
+#else
+	(CONST char **sp, char *t)
+#endif
+{
+	int c, d;
+	CONST char *s = *sp;
+
+	while( (d = *t++) !=0) {
+		if ((c = *++s) >= 'A' && c <= 'Z')
+			c += 'a' - 'A';
+		if (c != d)
+			return 0;
+		}
+	*sp = s + 1;
+	return 1;
+	}
+#endif /* INFNAN_CHECK */
+
+ void
+#ifdef KR_headers
+copybits(c, n, b) ULong *c; int n; Bigint *b;
+#else
+copybits(ULong *c, int n, Bigint *b)
+#endif
+{
+	ULong *ce, *x, *xe;
+#ifdef Pack_16
+	int nw, nw1;
+#endif
+
+	ce = c + ((n-1) >> kshift) + 1;
+	x = b->x;
+#ifdef Pack_32
+	xe = x + b->wds;
+	while(x < xe)
+		*c++ = *x++;
+#else
+	nw = b->wds;
+	nw1 = nw & 1;
+	for(xe = x + (nw - nw1); x < xe; x += 2)
+		Storeinc(c, x[1], x[0]);
+	if (nw1)
+		*c++ = *x;
+#endif
+	while(c < ce)
+		*c++ = 0;
+	}
+
+ ULong
+#ifdef KR_headers
+any_on(b, k) Bigint *b; int k;
+#else
+any_on(Bigint *b, int k)
+#endif
+{
+	int n, nwds;
+	ULong *x, *x0, x1, x2;
+
+	x = b->x;
+	nwds = b->wds;
+	n = k >> kshift;
+	if (n > nwds)
+		n = nwds;
+	else if (n < nwds && (k &= kmask)) {
+		x1 = x2 = x[n];
+		x1 >>= k;
+		x1 <<= k;
+		if (x1 != x2)
+			return 1;
+		}
+	x0 = x;
+	x += n;
+	while(x > x0)
+		if (*--x)
+			return 1;
+	return 0;
+	}

contrib/gdtoa/strtof.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998, 2000 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+/* $FreeBSD$ */
+
+#include "gdtoaimp.h"
+
+static float
+#ifdef KR_headers
+strtof_l(s, sp, loc) CONST char *s; char **sp; locale_t loc;
+#else
+strtof_l(CONST char *s, char **sp, locale_t loc)
+#endif
+{
+	static FPI fpi0 = { 24, 1-127-24+1,  254-127-24+1, 1, SI };
+	ULong bits[1];
+	Long exp;
+	int k;
+	union { ULong L[1]; float f; } u;
+#ifdef Honor_FLT_ROUNDS
+#include "gdtoa_fltrnds.h"
+#else
+#define fpi &fpi0
+#endif
+
+	k = strtodg_l(s, sp, fpi, &exp, bits, loc);
+	switch(k & STRTOG_Retmask) {
+	  case STRTOG_NoNumber:
+	  case STRTOG_Zero:
+		u.L[0] = 0;
+		break;
+
+	  case STRTOG_Normal:
+		u.L[0] = (bits[0] & 0x7fffff) | ((exp + 0x7f + 23) << 23);
+		break;
+
+	  case STRTOG_NaNbits:
+		/* FreeBSD local: always return a quiet NaN */
+		u.L[0] = bits[0] | 0x7fc00000;
+		break;	      
+
+	  case STRTOG_Denormal:
+		u.L[0] = bits[0];
+		break;
+
+	  case STRTOG_Infinite:
+		u.L[0] = 0x7f800000;
+		break;
+
+	  case STRTOG_NaN:
+		u.L[0] = f_QNAN;
+	  }
+	if (k & STRTOG_Neg)
+		u.L[0] |= 0x80000000L;
+	return u.f;
+	}
+ float
+#ifdef KR_headers
+strtof(s, sp) CONST char *s; char **sp;
+#else
+strtof(CONST char *s, char **sp)
+#endif
+{
+	// return strtof_l(s, sp, __get_locale());
+	return strtof_l(s, sp, NULL);
+}
+

contrib/gdtoa/strtord.c

+/****************************************************************
+
+The author of this software is David M. Gay.
+
+Copyright (C) 1998, 2000 by Lucent Technologies
+All Rights Reserved
+
+Permission to use, copy, modify, and distribute this software and
+its documentation for any purpose and without fee is hereby
+granted, provided that the above copyright notice appear in all
+copies and that both that the copyright notice and this
+permission notice and warranty disclaimer appear in supporting
+documentation, and that the name of Lucent or any of its entities
+not be used in advertising or publicity pertaining to
+distribution of the software without specific, written prior
+permission.
+
+LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
+INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
+IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
+SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
+IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
+ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
+THIS SOFTWARE.
+
+****************************************************************/
+
+/* Please send bug reports to David M. Gay (dmg at acm dot org,
+ * with " at " changed at "@" and " dot " changed to ".").	*/
+
+#include "gdtoaimp.h"
+
+ void
+#ifdef KR_headers
+ULtod(L, bits, exp, k) ULong *L; ULong *bits; Long exp; int k;
+#else
+ULtod(ULong *L, ULong *bits, Long exp, int k)
+#endif
+{
+	switch(k & STRTOG_Retmask) {
+	  case STRTOG_NoNumber:
+	  case STRTOG_Zero:
+		L[0] = L[1] = 0;
+		break;
+
+	  case STRTOG_Denormal:
+		L[_1] = bits[0];
+		L[_0] = bits[1];
+		break;
+
+	  case STRTOG_Normal:
+	  case STRTOG_NaNbits:
+		L[_1] = bits[0];
+		L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20);
+		break;
+
+	  case STRTOG_Infinite:
+		L[_0] = 0x7ff00000;
+		L[_1] = 0;
+		break;
+
+	  case STRTOG_NaN:
+		L[0] = d_QNAN0;
+		L[1] = d_QNAN1;
+	  }
+	if (k & STRTOG_Neg)
+		L[_0] |= 0x80000000L;
+	}
+
+ int
+#ifdef KR_headers
+strtord_l(s, sp, rounding, d, locale) CONST char *s; char **sp; int rounding;
+double *d; locale_t locale;
+#else
+strtord_l(CONST char *s, char **sp, int rounding, double *d, locale_t locale)
+#endif
+{
+	static FPI fpi0 = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
+	FPI *fpi, fpi1;
+	ULong bits[2];
+	Long exp;
+	int k;
+
+	fpi = &fpi0;
+	if (rounding != FPI_Round_near) {
+		fpi1 = fpi0;
+		fpi1.rounding = rounding;
+		fpi = &fpi1;
+		}
+	k = strtodg_l(s, sp, fpi, &exp, bits, locale);
+	ULtod((ULong*)d, bits, exp, k);
+	return k;
+	}
+

lib/libc/aarch64/arith.h

+/*
+ * MD header for contrib/gdtoa
+ *
+ * $FreeBSD$
+ */
+
+/*
+ * NOTE: The definitions in this file must be correct or strtod(3) and
+ * floating point formats in printf(3) will break!  The file can be
+ * generated by running contrib/gdtoa/arithchk.c on the target
+ * architecture.  See contrib/gdtoa/gdtoaimp.h for details.
+ */
+
+#define IEEE_8087
+#define Arith_Kind_ASL 1
+#define Long int
+#define Intcast (int)(long)
+#define Double_Align
+#define X64_bit_pointers

lib/libc/aarch64/gd_qnan.h

+/*
+ * MD header for contrib/gdtoa
+ *
+ * This file can be generated by compiling and running contrib/gdtoa/qnan.c
+ * on the target architecture after arith.h has been generated.
+ *
+ * $FreeBSD$
+ */
+
+#define f_QNAN 0x7fc00000
+#define d_QNAN0 0x0
+#define d_QNAN1 0x7ff80000
+#define ld_QNAN0 0x0
+#define ld_QNAN1 0x0
+#define ld_QNAN2 0x0
+#define ld_QNAN3 0x7fff8000
+#define ldus_QNAN0 0x0
+#define ldus_QNAN1 0x0
+#define ldus_QNAN2 0x0
+#define ldus_QNAN3 0x0
+#define ldus_QNAN4 0x0