dtoa.c revision 5a0b399aa9d75d4d208394d80b4997c164435ecb
1bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/**************************************************************** 2bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 3bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * The author of this software is David M. Gay. 4bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 5bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. 6bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 7bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Permission to use, copy, modify, and distribute this software for any 8bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * purpose without fee is hereby granted, provided that this entire notice 9bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * is included in all copies of any software which is or includes a copy 10bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * or modification of this software and in all copies of the supporting 11bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * documentation for such software. 12bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 13bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 14bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY 15bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 16bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 17bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 18bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ***************************************************************/ 19bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 20bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/**************************************************************** 21bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * This is dtoa.c by David M. Gay, downloaded from 22bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for 23bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith. 24bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 25bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Please remember to check http://www.netlib.org/fp regularly (and especially 26bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * before any Python release) for bugfixes and updates. 27bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 28bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * The major modifications from Gay's original code are as follows: 29bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 30bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 0. The original code has been specialized to Python's needs by removing 31bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * many of the #ifdef'd sections. In particular, code to support VAX and 32bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * IBM floating-point formats, hex NaNs, hex floats, locale-aware 33bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * treatment of the decimal point, and setting of the inexact flag have 34bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * been removed. 35bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 36bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free. 37bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 38bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2. The public functions strtod, dtoa and freedtoa all now have 39bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * a _Py_dg_ prefix. 40bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 41bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread 42bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * PyMem_Malloc failures through the code. The functions 43bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 44bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b 45bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 46bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * of return type *Bigint all return NULL to indicate a malloc failure. 47bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on 48bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * failure. bigcomp now has return type int (it used to be void) and 49bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * returns -1 on failure and 0 otherwise. _Py_dg_dtoa returns NULL 50bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * on failure. _Py_dg_strtod indicates failure due to malloc failure 51bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * by returning -1.0, setting errno=ENOMEM and *se to s00. 52bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 53bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 4. The static variable dtoa_result has been removed. Callers of 54bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * _Py_dg_dtoa are expected to call _Py_dg_freedtoa to free 55bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * the memory allocated by _Py_dg_dtoa. 56bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 57bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 5. The code has been reformatted to better fit with Python's 58bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * C style guide (PEP 7). 59bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 60bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory 61bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * that hasn't been MALLOC'ed, private_mem should only be used when k <= 62bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Kmax. 63bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 64bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with 65bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * leading whitespace. 66bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 67bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ***************************************************************/ 68bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 69bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Please send bug reports for the original dtoa.c code to David M. Gay (dmg 70bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * at acm dot org, with " at " changed at "@" and " dot " changed to "."). 71bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Please report bugs for this modified version using the Python issue tracker 72bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * (http://bugs.python.org). */ 73bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 74bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* On a machine with IEEE extended-precision registers, it is 75bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * necessary to specify double-precision (53-bit) rounding precision 76bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * before invoking strtod or dtoa. If the machine uses (the equivalent 77bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * of) Intel 80x87 arithmetic, the call 78bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * _control87(PC_53, MCW_PC); 79bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * does this with many compilers. Whether this or another call is 80bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * appropriate depends on the compiler; for this to work, it may be 81bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * necessary to #include "float.h" or another system-dependent header 82bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * file. 83bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 84bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 85bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* strtod for IEEE-, VAX-, and IBM-arithmetic machines. 86bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 87bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * This strtod returns a nearest machine number to the input decimal 88bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * string (or sets errno to ERANGE). With IEEE arithmetic, ties are 89bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * broken by the IEEE round-even rule. Otherwise ties are broken by 90bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * biased rounding (add half and chop). 91bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 92bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Inspired loosely by William D. Clinger's paper "How to Read Floating 93bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. 94bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 95bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Modifications: 96bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 97bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 1. We only require IEEE, IBM, or VAX double-precision 98bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * arithmetic (not IEEE double-extended). 99bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2. We get by with floating-point arithmetic in a case that 100bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Clinger missed -- when we're computing d * 10^n 101bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * for a small integer d and the integer n is not too 102bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * much larger than 22 (the maximum integer k for which 103bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * we can represent 10^k exactly), we may be able to 104bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * compute (d*10^k) * 10^(e-k) with just one roundoff. 105bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 3. Rather than a bit-at-a-time adjustment of the binary 106bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * result in the hard case, we use floating-point 107bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * arithmetic to determine the adjustment to within 108bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * one bit; only in really hard cases do we need to 109bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * compute a second residual. 110bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 4. Because of 3., we don't need a large table of powers of 10 111bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * for ten-to-e (just some small tables, e.g. of 10^k 112bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * for 0 <= k <= 22). 113bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 114bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 115bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Linking of Python's #defines to Gay's #defines starts here. */ 116bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 117bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#include "Python.h" 118bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 119bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile 120bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson the following code */ 121bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifndef PY_NO_SHORT_FLOAT_REPR 122bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 123bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#include "float.h" 124bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 125bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define MALLOC PyMem_Malloc 126bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define FREE PyMem_Free 127bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 128bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* This code should also work for ARM mixed-endian format on little-endian 129bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson machines, where doubles have byte order 45670123 (in increasing address 130bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson order, 0 being the least significant byte). */ 131bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754 132bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson# define IEEE_8087 133bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 134bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \ 135bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754) 136bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson# define IEEE_MC68k 137bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 138bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if defined(IEEE_8087) + defined(IEEE_MC68k) != 1 139bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined." 140bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 141bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 142bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* The code below assumes that the endianness of integers matches the 143bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson endianness of the two 32-bit words of a double. Check this. */ 144bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \ 145bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)) 146bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#error "doubles and ints have incompatible endianness" 147bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 148bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 149bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) 150bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#error "doubles and ints have incompatible endianness" 151bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 152bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 153bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 154bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if defined(HAVE_UINT32_T) && defined(HAVE_INT32_T) 155bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsontypedef PY_UINT32_T ULong; 156bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsontypedef PY_INT32_T Long; 157bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 158bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#error "Failed to find an exact-width 32-bit integer type" 159bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 160bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 161bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if defined(HAVE_UINT64_T) 162bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define ULLong PY_UINT64_T 163bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 164bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#undef ULLong 165bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 166bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 167bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#undef DEBUG 168bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef Py_DEBUG 169bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define DEBUG 170bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 171bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 172bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* End Python #define linking */ 173bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 174bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef DEBUG 175bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} 176bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 177bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 178bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifndef PRIVATE_MEM 179bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define PRIVATE_MEM 2304 180bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 181bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 182bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic double private_mem[PRIVATE_mem], *pmem_next = private_mem; 183bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 184bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef __cplusplus 185bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonextern "C" { 186bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 187bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 188bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsontypedef union { double d; ULong L[2]; } U; 189bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 190bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef IEEE_8087 191bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define word0(x) (x)->L[1] 192bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define word1(x) (x)->L[0] 193bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 194bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define word0(x) (x)->L[0] 195bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define word1(x) (x)->L[1] 196bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 197bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define dval(x) (x)->d 198bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 199bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifndef STRTOD_DIGLIM 200bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define STRTOD_DIGLIM 40 201bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 202bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 203bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* The following definition of Storeinc is appropriate for MIPS processors. 204bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * An alternative that might be better on some machines is 205bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) 206bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 207bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#if defined(IEEE_8087) 208bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ 209bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ((unsigned short *)a)[0] = (unsigned short)c, a++) 210bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 211bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ 212bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ((unsigned short *)a)[1] = (unsigned short)c, a++) 213bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 214bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 215bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* #define P DBL_MANT_DIG */ 216bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Ten_pmax = floor(P*log(2)/log(5)) */ 217bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ 218bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ 219bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ 220bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 221bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_shift 20 222bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_shift1 20 223bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_msk1 0x100000 224bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_msk11 0x100000 225bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_mask 0x7ff00000 226bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define P 53 227bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Nbits 53 228bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Bias 1023 229bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Emax 1023 230bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Emin (-1022) 231bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_1 0x3ff00000 232bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Exp_11 0x3ff00000 233bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Ebits 11 234bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Frac_mask 0xfffff 235bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Frac_mask1 0xfffff 236bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Ten_pmax 22 237bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Bletch 0x10 238bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Bndry_mask 0xfffff 239bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Bndry_mask1 0xfffff 240bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define LSB 1 241bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Sign_bit 0x80000000 242bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Log2P 1 243bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Tiny0 0 244bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Tiny1 1 245bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Quick_max 14 246bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Int_max 14 247bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 248bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifndef Flt_Rounds 249bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef FLT_ROUNDS 250bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Flt_Rounds FLT_ROUNDS 251bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 252bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Flt_Rounds 1 253bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 254bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif /*Flt_Rounds*/ 255bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 256bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Rounding Flt_Rounds 257bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 258bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) 259bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Big1 0xffffffff 260bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 261bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */ 262bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 263bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsontypedef struct BCinfo BCinfo; 264bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstruct 265bb28285ea2f01e97a26bc595d49da43fbee62913Mark DickinsonBCinfo { 2665a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson int dp0, dp1, dplen, dsign, e0, nd, nd0, scale; 267bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson}; 268bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 269bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define FFFFFFFF 0xffffffffUL 270bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 271bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Kmax 7 272bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 273bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* struct Bigint is used to represent arbitrary-precision integers. These 274bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson integers are stored in sign-magnitude format, with the magnitude stored as 275bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson an array of base 2**32 digits. Bigints are always normalized: if x is a 276bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero. 277bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 278bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson The Bigint fields are as follows: 279bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 280bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson - next is a header used by Balloc and Bfree to keep track of lists 281bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson of freed Bigints; it's also used for the linked list of 282bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson powers of 5 of the form 5**2**i used by pow5mult. 283bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson - k indicates which pool this Bigint was allocated from 284bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson - maxwds is the maximum number of words space was allocated for 285bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson (usually maxwds == 2**k) 286bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson - sign is 1 for negative Bigints, 0 for positive. The sign is unused 287bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson (ignored on inputs, set to 0 on outputs) in almost all operations 288bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson involving Bigints: a notable exception is the diff function, which 289bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ignores signs on inputs but sets the sign of the output correctly. 290bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson - wds is the actual number of significant words 291bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson - x contains the vector of words (digits) for this Bigint, from least 292bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson significant (x[0]) to most significant (x[wds-1]). 293bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson*/ 294bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 295bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstruct 296bb28285ea2f01e97a26bc595d49da43fbee62913Mark DickinsonBigint { 297bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson struct Bigint *next; 298bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int k, maxwds, sign, wds; 299bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong x[1]; 300bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson}; 301bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 302bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsontypedef struct Bigint Bigint; 303bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 304bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Memory management: memory is allocated from, and returned to, Kmax+1 pools 305bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds == 306bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1 << k. These pools are maintained as linked lists, with freelist[k] 307bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson pointing to the head of the list for pool k. 308bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 309bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson On allocation, if there's no free slot in the appropriate pool, MALLOC is 310bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson called to get more memory. This memory is not returned to the system until 311bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Python quits. There's also a private memory pool that's allocated from 312bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson in preference to using MALLOC. 313bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 314bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson For Bigints with more than (1 << Kmax) digits (which implies at least 1233 315bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson decimal digits), memory is directly allocated using MALLOC, and freed using 316bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson FREE. 317bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 318bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson XXX: it would be easy to bypass this memory-management system and 319bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson translate each call to Balloc into a call to PyMem_Malloc, and each 320bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree to PyMem_Free. Investigate whether this has any significant 321bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson performance on impact. */ 322bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 323bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint *freelist[Kmax+1]; 324bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 325bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Allocate space for a Bigint with up to 1<<k digits */ 326bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 327bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 328bb28285ea2f01e97a26bc595d49da43fbee62913Mark DickinsonBalloc(int k) 329bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 330bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int x; 331bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *rv; 332bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson unsigned int len; 333bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 334bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k <= Kmax && (rv = freelist[k])) 335bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson freelist[k] = rv->next; 336bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 337bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = 1 << k; 338bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 339bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /sizeof(double); 340bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 341bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv = (Bigint*)pmem_next; 342bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson pmem_next += len; 343bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 344bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 345bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv = (Bigint*)MALLOC(len*sizeof(double)); 346bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (rv == NULL) 347bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 348bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 349bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv->k = k; 350bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv->maxwds = x; 351bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 352bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv->sign = rv->wds = 0; 353bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return rv; 354bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 355bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 356bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Free a Bigint allocated with Balloc */ 357bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 358bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic void 359bb28285ea2f01e97a26bc595d49da43fbee62913Mark DickinsonBfree(Bigint *v) 360bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 361bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (v) { 362bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (v->k > Kmax) 363bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson FREE((void*)v); 364bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 365bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson v->next = freelist[v->k]; 366bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson freelist[v->k] = v; 367bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 368bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 369bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 370bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 371bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ 372bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y->wds*sizeof(Long) + 2*sizeof(int)) 373bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 374bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Multiply a Bigint b by m and add a. Either modifies b in place and returns 375bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson a pointer to the modified b, or Bfrees b and returns a pointer to a copy. 376bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson On failure, return NULL. In this case, b will have been already freed. */ 377bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 378bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 379bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonmultadd(Bigint *b, int m, int a) /* multiply by m and add a */ 380bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 381bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i, wds; 382bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 383bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *x; 384bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULLong carry, y; 385bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 386bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong carry, *x, y; 387bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong xi, z; 388bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 389bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b1; 390bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 391bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wds = b->wds; 392bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = b->x; 393bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 0; 394bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = a; 395bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 396bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 397bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = *x * (ULLong)m + carry; 398bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = y >> 32; 399bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *x++ = (ULong)(y & FFFFFFFF); 400bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 401bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xi = *x; 402bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = (xi & 0xffff) * m + carry; 403bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (xi >> 16) * m + (y >> 16); 404bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = z >> 16; 405bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *x++ = (z << 16) + (y & 0xffff); 406bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 407bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 408bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(++i < wds); 409bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (carry) { 410bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (wds >= b->maxwds) { 411bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b1 = Balloc(b->k+1); 412bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b1 == NULL){ 413bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 414bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 415bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 416bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bcopy(b1, b); 417bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 418bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = b1; 419bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 420bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->x[wds++] = (ULong)carry; 421bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = wds; 422bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 423bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b; 424bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 425bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 426bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* convert a string s containing nd decimal digits (possibly containing a 427bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson decimal separator at position nd0, which is ignored) to a Bigint. This 428bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson function carries on where the parsing code in _Py_dg_strtod leaves off: on 429bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson entry, y9 contains the result of converting the first 9 digits. Returns 430bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson NULL on failure. */ 431bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 432bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 433bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsons2b(const char *s, int nd0, int nd, ULong y9, int dplen) 434bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 435bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b; 436bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i, k; 437bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Long x, y; 438bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 439bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = (nd + 8) / 9; 440bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(k = 0, y = 1; x > y; y <<= 1, k++) ; 441bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = Balloc(k); 442bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 443bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 444bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->x[0] = y9; 445bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = 1; 446bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 447bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 9; 448bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (9 < nd0) { 449bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s += 9; 450bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 451bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, *s++ - '0'); 452bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 453bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 454bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } while(++i < nd0); 455bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s += dplen; 456bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 457bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 458bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s += dplen + 9; 459bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(; i < nd; i++) { 460bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, *s++ - '0'); 461bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 462bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 463bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 464bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b; 465bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 466bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 467bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* count leading 0 bits in the 32-bit integer x. */ 468bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 469bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic int 470bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonhi0bits(ULong x) 471bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 472bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int k = 0; 473bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 474bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xffff0000)) { 475bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = 16; 476bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x <<= 16; 477bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 478bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xff000000)) { 479bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 8; 480bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x <<= 8; 481bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 482bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xf0000000)) { 483bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 4; 484bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x <<= 4; 485bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 486bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xc0000000)) { 487bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 2; 488bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x <<= 2; 489bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 490bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0x80000000)) { 491bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 492bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0x40000000)) 493bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 32; 494bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 495bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return k; 496bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 497bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 498bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* count trailing 0 bits in the 32-bit integer y, and shift y right by that 499bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson number of bits. */ 500bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 501bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic int 502bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonlo0bits(ULong *y) 503bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 504bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int k; 505bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong x = *y; 506bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 507bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (x & 7) { 508bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (x & 1) 509bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 0; 510bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (x & 2) { 511bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *y = x >> 1; 512bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 1; 513bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 514bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *y = x >> 2; 515bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 2; 516bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 517bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = 0; 518bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xffff)) { 519bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = 16; 520bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x >>= 16; 521bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 522bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xff)) { 523bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 8; 524bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x >>= 8; 525bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 526bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0xf)) { 527bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 4; 528bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x >>= 4; 529bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 530bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 0x3)) { 531bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 2; 532bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x >>= 2; 533bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 534bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(x & 1)) { 535bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 536bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x >>= 1; 537bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!x) 538bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 32; 539bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 540bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *y = x; 541bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return k; 542bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 543bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 544bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* convert a small nonnegative integer to a Bigint */ 545bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 546bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 547bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsoni2b(int i) 548bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 549bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b; 550bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 551bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = Balloc(1); 552bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 553bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 554bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->x[0] = i; 555bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = 1; 556bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b; 557bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 558bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 559bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores 560bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson the signs of a and b. */ 561bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 562bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 563bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonmult(Bigint *a, Bigint *b) 564bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 565bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *c; 566bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int k, wa, wb, wc; 567bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 568bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong y; 569bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 570bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULLong carry, z; 571bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 572bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong carry, z; 573bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong z2; 574bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 575bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 576bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (a->wds < b->wds) { 577bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = a; 578bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson a = b; 579bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = c; 580bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 581bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = a->k; 582bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wa = a->wds; 583bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wb = b->wds; 584bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wc = wa + wb; 585bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (wc > a->maxwds) 586bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 587bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = Balloc(k); 588bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c == NULL) 589bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 590bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(x = c->x, xa = x + wc; x < xa; x++) 591bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *x = 0; 592bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xa = a->x; 593bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xae = xa + wa; 594bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xb = b->x; 595bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xbe = xb + wb; 596bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xc0 = c->x; 597bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 598bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(; xb < xbe; xc0++) { 599bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((y = *xb++)) { 600bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = xa; 601bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xc = xc0; 602bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = 0; 603bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 604bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = *x++ * (ULLong)y + *xc + carry; 605bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = z >> 32; 606bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *xc++ = (ULong)(z & FFFFFFFF); 607bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 608bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(x < xae); 609bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *xc = (ULong)carry; 610bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 611bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 612bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 613bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(; xb < xbe; xb++, xc0++) { 614bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (y = *xb & 0xffff) { 615bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = xa; 616bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xc = xc0; 617bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = 0; 618bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 619bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 620bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = z >> 16; 621bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 622bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = z2 >> 16; 623bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Storeinc(xc, z2, z); 624bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 625bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(x < xae); 626bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *xc = carry; 627bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 628bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (y = *xb >> 16) { 629bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = xa; 630bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xc = xc0; 631bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = 0; 632bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z2 = *xc; 633bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 634bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (*x & 0xffff) * y + (*xc >> 16) + carry; 635bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = z >> 16; 636bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Storeinc(xc, z, z2); 637bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 638bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = z2 >> 16; 639bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 640bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(x < xae); 641bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *xc = z2; 642bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 643bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 644bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 645bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 646bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c->wds = wc; 647bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return c; 648bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 649bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 650bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */ 651bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 652bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint *p5s; 653bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 654bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on 655bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson failure; if the returned pointer is distinct from b then the original 656bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint b will have been Bfree'd. Ignores the sign of b. */ 657bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 658bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 659bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonpow5mult(Bigint *b, int k) 660bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 661bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b1, *p5, *p51; 662bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i; 663bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson static int p05[3] = { 5, 25, 125 }; 664bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 665bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((i = k & 3)) { 666bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, p05[i-1], 0); 667bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 668bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 669bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 670bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 671bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(k >>= 2)) 672bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b; 673bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5 = p5s; 674bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!p5) { 675bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* first time */ 676bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5 = i2b(625); 677bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (p5 == NULL) { 678bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 679bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 680bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 681bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5s = p5; 682bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5->next = 0; 683bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 684bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(;;) { 685bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k & 1) { 686bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b1 = mult(b, p5); 687bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 688bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = b1; 689bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 690bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 691bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 692bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(k >>= 1)) 693bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 694bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p51 = p5->next; 695bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!p51) { 696bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p51 = mult(p5,p5); 697bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (p51 == NULL) { 698bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 699bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 700bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 701bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p51->next = 0; 702bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5->next = p51; 703bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 704bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5 = p51; 705bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 706bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b; 707bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 708bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 709bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* shift a Bigint b left by k bits. Return a pointer to the shifted result, 710bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson or NULL on failure. If the returned pointer is distinct from b then the 711bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson original b will have been Bfree'd. Ignores the sign of b. */ 712bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 713bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 714bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonlshift(Bigint *b, int k) 715bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 716bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i, k1, n, n1; 717bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b1; 718bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *x, *x1, *xe, z; 719bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 720bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson n = k >> 5; 721bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k1 = b->k; 722bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson n1 = n + b->wds + 1; 723bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = b->maxwds; n1 > i; i <<= 1) 724bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k1++; 725bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b1 = Balloc(k1); 726bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b1 == NULL) { 727bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 728bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 729bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 730bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x1 = b1->x; 731bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 0; i < n; i++) 732bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *x1++ = 0; 733bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = b->x; 734bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xe = x + b->wds; 735bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k &= 0x1f) { 736bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k1 = 32 - k; 737bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = 0; 738bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 739bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *x1++ = *x << k | z; 740bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = *x++ >> k1; 741bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 742bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(x < xe); 743bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((*x1 = z)) 744bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ++n1; 745bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 746bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else do 747bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *x1++ = *x++; 748bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(x < xe); 749bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b1->wds = n1 - 1; 750bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 751bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b1; 752bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 753bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 754bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and 755bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1 if a > b. Ignores signs of a and b. */ 756bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 757bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic int 758bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsoncmp(Bigint *a, Bigint *b) 759bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 760bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *xa, *xa0, *xb, *xb0; 761bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i, j; 762bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 763bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = a->wds; 764bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = b->wds; 765bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef DEBUG 766bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i > 1 && !a->x[i-1]) 767bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bug("cmp called with a->x[a->wds-1] == 0"); 768bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j > 1 && !b->x[j-1]) 769bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bug("cmp called with b->x[b->wds-1] == 0"); 770bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 771bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i -= j) 772bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return i; 773bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xa0 = a->x; 774bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xa = xa0 + j; 775bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xb0 = b->x; 776bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xb = xb0 + j; 777bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(;;) { 778bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (*--xa != *--xb) 779bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return *xa < *xb ? -1 : 1; 780bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (xa <= xa0) 781bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 782bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 783bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 0; 784bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 785bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 786bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Take the difference of Bigints a and b, returning a new Bigint. Returns 787bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson NULL on failure. The signs of a and b are ignored, but the sign of the 788bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson result is set appropriately. */ 789bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 790bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 791bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsondiff(Bigint *a, Bigint *b) 792bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 793bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *c; 794bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i, wa, wb; 795bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *xa, *xae, *xb, *xbe, *xc; 796bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 797bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULLong borrow, y; 798bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 799bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong borrow, y; 800bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong z; 801bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 802bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 803bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = cmp(a,b); 804bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!i) { 805bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = Balloc(0); 806bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c == NULL) 807bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 808bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c->wds = 1; 809bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c->x[0] = 0; 810bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return c; 811bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 812bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i < 0) { 813bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = a; 814bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson a = b; 815bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = c; 816bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 1; 817bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 818bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 819bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 0; 820bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = Balloc(a->k); 821bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c == NULL) 822bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 823bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c->sign = i; 824bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wa = a->wds; 825bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xa = a->x; 826bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xae = xa + wa; 827bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wb = b->wds; 828bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xb = b->x; 829bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xbe = xb + wb; 830bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xc = c->x; 831bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = 0; 832bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 833bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 834bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = (ULLong)*xa++ - *xb++ - borrow; 835bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = y >> 32 & (ULong)1; 836bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *xc++ = (ULong)(y & FFFFFFFF); 837bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 838bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(xb < xbe); 839bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(xa < xae) { 840bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = *xa++ - borrow; 841bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = y >> 32 & (ULong)1; 842bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *xc++ = (ULong)(y & FFFFFFFF); 843bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 844bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 845bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 846bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 847bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (y & 0x10000) >> 16; 848bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 849bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (z & 0x10000) >> 16; 850bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Storeinc(xc, z, y); 851bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 852bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(xb < xbe); 853bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(xa < xae) { 854bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = (*xa & 0xffff) - borrow; 855bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (y & 0x10000) >> 16; 856bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (*xa++ >> 16) - borrow; 857bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (z & 0x10000) >> 16; 858bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Storeinc(xc, z, y); 859bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 860bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 861bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(!*--xc) 862bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson wa--; 863bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c->wds = wa; 864bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return c; 865bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 866bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 867bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Given a positive normal double x, return the difference between x and the next 868bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson double up. Doesn't give correct results for subnormals. */ 869bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 870bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic double 871bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonulp(U *x) 872bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 873bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Long L; 874bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson U u; 875bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 876bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; 877bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&u) = L; 878bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&u) = 0; 879bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return dval(&u); 880bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 881bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 882bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Convert a Bigint to a double plus an exponent */ 883bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 884bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic double 885bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonb2d(Bigint *a, int *e) 886bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 887bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *xa, *xa0, w, y, z; 888bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int k; 889bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson U d; 890bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 891bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xa0 = a->x; 892bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson xa = xa0 + a->wds; 893bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = *--xa; 894bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef DEBUG 895bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!y) Bug("zero y in b2d"); 896bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 897bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = hi0bits(y); 898bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *e = 32 - k; 899bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k < Ebits) { 900bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&d) = Exp_1 | y >> (Ebits - k); 901bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson w = xa > xa0 ? *--xa : 0; 902bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k); 903bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret_d; 904bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 905bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = xa > xa0 ? *--xa : 0; 906bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k -= Ebits) { 907bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&d) = Exp_1 | y << k | z >> (32 - k); 908bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = xa > xa0 ? *--xa : 0; 909bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&d) = z << k | y >> (32 - k); 910bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 911bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 912bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&d) = Exp_1 | y; 913bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&d) = z; 914bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 915bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ret_d: 916bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return dval(&d); 917bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 918bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 919bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Convert a double to a Bigint plus an exponent. Return NULL on failure. 920bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 921bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Given a finite nonzero double d, return an odd Bigint b and exponent *e 922bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson such that fabs(d) = b * 2**e. On return, *bbits gives the number of 9232bcd17727027b6c923340209f8c6b92e34c69556Mark Dickinson significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits). 924bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 925bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson If d is zero, then b == 0, *e == -1010, *bbits = 0. 926bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 927bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 928bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 929bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic Bigint * 930bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsond2b(U *d, int *e, int *bits) 931bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 932bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b; 933bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int de, k; 934bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *x, y, z; 935bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int i; 936bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 937bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = Balloc(1); 938bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 939bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 940bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = b->x; 941bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 942bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = word0(d) & Frac_mask; 943bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */ 944bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((de = (int)(word0(d) >> Exp_shift))) 945bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z |= Exp_msk1; 946bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((y = word1(d))) { 947bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((k = lo0bits(&y))) { 948bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x[0] = y | z << (32 - k); 949bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z >>= k; 950bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 951bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 952bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x[0] = y; 953bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 954bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = (x[1] = z) ? 2 : 1; 955bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 956bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 957bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = lo0bits(&z); 958bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x[0] = z; 959bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 960bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = 1; 961bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k += 32; 962bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 963bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (de) { 964bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *e = de - Bias - (P-1) + k; 965bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *bits = P - k; 966bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 967bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 968bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *e = de - Bias - (P-1) + 1 + k; 969bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *bits = 32*i - hi0bits(x[i-1]); 970bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 971bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return b; 972bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 973bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 974bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Compute the ratio of two Bigints, as a double. The result may have an 975bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson error of up to 2.5 ulps. */ 976bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 977bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic double 978bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonratio(Bigint *a, Bigint *b) 979bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 980bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson U da, db; 981bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int k, ka, kb; 982bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 983bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&da) = b2d(a, &ka); 984bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&db) = b2d(b, &kb); 985bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = ka - kb + 32*(a->wds - b->wds); 986bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k > 0) 987bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&da) += k*Exp_msk1; 988bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 989bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = -k; 990bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&db) += k*Exp_msk1; 991bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 992bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return dval(&da) / dval(&db); 993bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 994bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 995bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic const double 996bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsontens[] = { 997bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 998bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 999bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1e20, 1e21, 1e22 1000bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson}; 1001bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1002bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic const double 1003bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonbigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 1004bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1005bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 9007199254740992.*9007199254740992.e-256 1006bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* = 2^106 * 1e-256 */ 1007bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson}; 1008bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ 1009bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ 1010bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define Scale_Bit 0x10 1011bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define n_bigtens 5 1012bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1013bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define ULbits 32 1014bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define kshift 5 1015bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define kmask 31 1016bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1017bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1018bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic int 1019bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsondshift(Bigint *b, int p2) 1020bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1021bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int rv = hi0bits(b->x[b->wds-1]) - 4; 1022bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (p2 > 0) 1023bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv -= p2; 1024bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return rv & kmask; 1025bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1026bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1027bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* special case of Bigint division. The quotient is always in the range 0 <= 1028bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson quotient < 10, and on entry the divisor S is normalized so that its top 4 1029bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bits (28--31) are zero and bit 27 is set. */ 1030bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1031bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic int 1032bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonquorem(Bigint *b, Bigint *S) 1033bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1034bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int n; 1035bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong *bx, *bxe, q, *sx, *sxe; 1036bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 1037bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULLong borrow, carry, y, ys; 1038bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 1039bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong borrow, carry, y, ys; 1040bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong si, z, zs; 1041bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 1042bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1043bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson n = S->wds; 1044bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef DEBUG 1045bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /*debug*/ if (b->wds > n) 1046bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /*debug*/ Bug("oversize b in quorem"); 1047bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 1048bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b->wds < n) 1049bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 0; 1050bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson sx = S->x; 1051bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson sxe = sx + --n; 1052bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bx = b->x; 1053bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bxe = bx + n; 1054bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ 1055bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef DEBUG 1056bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /*debug*/ if (q > 9) 1057bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /*debug*/ Bug("oversized quotient in quorem"); 1058bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 1059bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (q) { 1060bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = 0; 1061bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = 0; 1062bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 1063bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 1064bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ys = *sx++ * (ULLong)q + carry; 1065bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = ys >> 32; 1066bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = *bx - (ys & FFFFFFFF) - borrow; 1067bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = y >> 32 & (ULong)1; 1068bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *bx++ = (ULong)(y & FFFFFFFF); 1069bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 1070bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson si = *sx++; 1071bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ys = (si & 0xffff) * q + carry; 1072bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson zs = (si >> 16) * q + (ys >> 16); 1073bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = zs >> 16; 1074bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 1075bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (y & 0x10000) >> 16; 1076bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (*bx >> 16) - (zs & 0xffff) - borrow; 1077bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (z & 0x10000) >> 16; 1078bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Storeinc(bx, z, y); 1079bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 1080bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1081bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(sx <= sxe); 1082bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!*bxe) { 1083bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bx = b->x; 1084bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(--bxe > bx && !*bxe) 1085bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson --n; 1086bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = n; 1087bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1088bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1089bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (cmp(b, S) >= 0) { 1090bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson q++; 1091bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = 0; 1092bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = 0; 1093bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bx = b->x; 1094bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson sx = S->x; 1095bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson do { 1096bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef ULLong 1097bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ys = *sx++ + carry; 1098bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = ys >> 32; 1099bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = *bx - (ys & FFFFFFFF) - borrow; 1100bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = y >> 32 & (ULong)1; 1101bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *bx++ = (ULong)(y & FFFFFFFF); 1102bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#else 1103bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson si = *sx++; 1104bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ys = (si & 0xffff) + carry; 1105bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson zs = (si >> 16) + (ys >> 16); 1106bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson carry = zs >> 16; 1107bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = (*bx & 0xffff) - (ys & 0xffff) - borrow; 1108bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (y & 0x10000) >> 16; 1109bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = (*bx >> 16) - (zs & 0xffff) - borrow; 1110bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson borrow = (z & 0x10000) >> 16; 1111bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Storeinc(bx, z, y); 1112bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 1113bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1114bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(sx <= sxe); 1115bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bx = b->x; 1116bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bxe = bx + n; 1117bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!*bxe) { 1118bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(--bxe > bx && !*bxe) 1119bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson --n; 1120bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->wds = n; 1121bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1122bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1123bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return q; 1124bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1125bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1126bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1127bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* return 0 on success, -1 on failure */ 1128bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1129bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic int 1130bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonbigcomp(U *rv, const char *s0, BCinfo *bc) 1131bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1132bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b, *d; 1133bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase; 1134bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1135bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dsign = bc->dsign; 1136bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nd = bc->nd; 1137bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nd0 = bc->nd0; 1138bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p5 = nd + bc->e0 - 1; 1139bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson speccase = 0; 1140bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (rv->d == 0.) { /* special case: value near underflow-to-zero */ 1141bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* threshold was rounded to zero */ 1142bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = i2b(1); 1143bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 1144bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1145bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p2 = Emin - P + 1; 1146bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bbits = 1; 1147bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(rv) = (P+2) << Exp_shift; 1148bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 0; 1149bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson { 1150bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson speccase = 1; 1151bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson --p2; 1152bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dsign = 0; 1153bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto have_i; 1154bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1155bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1156bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1157bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson { 1158bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = d2b(rv, &p2, &bbits); 1159bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 1160bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1161bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1162bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p2 -= bc->scale; 1163bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* floor(log2(rv)) == bbits - 1 + p2 */ 1164bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Check for denormal case. */ 1165bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = P - bbits; 1166bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i > (j = P - Emin - 1 + p2)) { 1167bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = j; 1168bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1169bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson { 1170bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = lshift(b, ++i); 1171bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 1172bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1173bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->x[0] |= 1; 1174bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1175bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson have_i: 1176bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson p2 -= p5 + i; 1177bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson d = i2b(1); 1178bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (d == NULL) { 1179bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 1180bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1181bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1182bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Arrange for convenient computation of quotients: 1183bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * shift left if necessary so divisor has 4 leading 0 bits. 1184bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1185bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (p5 > 0) { 1186bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson d = pow5mult(d, p5); 1187bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (d == NULL) { 1188bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 1189bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1190bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1191bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1192bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (p5 < 0) { 1193bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = pow5mult(b, -p5); 1194bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) { 1195bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(d); 1196bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1197bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1198bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1199bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (p2 > 0) { 1200bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 = p2; 1201bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson d2 = 0; 1202bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1203bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1204bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 = 0; 1205bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson d2 = -p2; 1206bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1207bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = dshift(d, d2); 1208bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((b2 += i) > 0) { 1209bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = lshift(b, b2); 1210bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) { 1211bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(d); 1212bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1213bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1214bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1215bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((d2 += i) > 0) { 1216bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson d = lshift(d, d2); 1217bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (d == NULL) { 1218bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 1219bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1220bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1221bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1222bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1223bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Now b/d = exactly half-way between the two floating-point values */ 1224bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* on either side of the input string. Compute first digit of b/d. */ 1225bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1226bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(dig = quorem(b,d))) { 1227bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, 0); /* very unlikely */ 1228bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) { 1229bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(d); 1230bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1231bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1232bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dig = quorem(b,d); 1233bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1234bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1235bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Compare b/d with s0 */ 1236bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1237bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson assert(nd > 0); 1238bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dd = 9999; /* silence gcc compiler warning */ 1239bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 0; i < nd0; ) { 1240bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((dd = s0[i++] - '0' - dig)) 1241bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1242bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!b->x[0] && b->wds == 1) { 1243bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i < nd) 1244bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dd = 1; 1245bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1246bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1247bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, 0); 1248bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) { 1249bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(d); 1250bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1251bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1252bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dig = quorem(b,d); 1253bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1254bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(j = bc->dp1; i++ < nd;) { 1255bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((dd = s0[j++] - '0' - dig)) 1256bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1257bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!b->x[0] && b->wds == 1) { 1258bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i < nd) 1259bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dd = 1; 1260bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1261bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1262bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, 0); 1263bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) { 1264bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(d); 1265bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1; 1266bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1267bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dig = quorem(b,d); 1268bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1269bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b->x[0] || b->wds > 1) 1270bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dd = -1; 1271bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ret: 1272bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 1273bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(d); 1274bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (speccase) { 1275bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dd <= 0) 1276bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv->d = 0.; 1277bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1278bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (dd < 0) { 1279bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!dsign) /* does not happen for round-near */ 1280bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson retlow1: 1281bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(rv) -= ulp(rv); 1282bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1283bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (dd > 0) { 1284bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dsign) { 1285bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rethi1: 1286bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(rv) += ulp(rv); 1287bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1288bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1289bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1290bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Exact half-way case: apply round-even rule. */ 1291bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (word1(rv) & 1) { 1292bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dsign) 1293bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto rethi1; 1294bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto retlow1; 1295bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1296bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1297bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1298bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return 0; 1299bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1300bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1301bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsondouble 1302bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson_Py_dg_strtod(const char *s00, char **se) 1303bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1304bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1, error; 1305bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int esign, i, j, k, nd, nd0, nf, nz, nz0, sign; 1306bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson const char *s, *s0, *s1; 1307bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson double aadj, aadj1; 1308bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Long L; 1309bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson U aadj2, adj, rv, rv0; 1310bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong y, z; 1311bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson BCinfo bc; 1312bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 1313bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 13145a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson sign = nz0 = nz = bc.dplen = 0; 1315bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) = 0.; 1316bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(s = s00;;s++) switch(*s) { 1317bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '-': 1318bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson sign = 1; 1319bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* no break */ 1320bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '+': 1321bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (*++s) 1322bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto break2; 1323bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* no break */ 1324bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 0: 1325bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret0; 1326bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* modify original dtoa.c so that it doesn't accept leading whitespace 1327bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '\t': 1328bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '\n': 1329bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '\v': 1330bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '\f': 1331bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '\r': 1332bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case ' ': 1333bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson continue; 1334bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1335bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson default: 1336bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto break2; 1337bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1338bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break2: 1339bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (*s == '0') { 1340bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nz0 = 1; 1341bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(*++s == '0') ; 1342bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!*s) 1343bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1344bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1345bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s0 = s; 1346bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = z = 0; 1347bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 1348bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (nd < 9) 1349bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = 10*y + c - '0'; 1350bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (nd < 16) 1351bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = 10*z + c - '0'; 1352bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nd0 = nd; 1353bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.dp0 = bc.dp1 = s - s0; 1354bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c == '.') { 1355bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = *++s; 1356bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.dp1 = s - s0; 1357bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.dplen = bc.dp1 - bc.dp0; 1358bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!nd) { 1359bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(; c == '0'; c = *++s) 1360bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nz++; 1361bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c > '0' && c <= '9') { 1362bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s0 = s; 1363bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nf += nz; 1364bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nz = 0; 1365bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto have_dig; 1366bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1367bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto dig_done; 1368bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1369bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(; c >= '0' && c <= '9'; c = *++s) { 1370bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson have_dig: 1371bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nz++; 1372bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c -= '0') { 1373bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nf += nz; 1374bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 1; i < nz; i++) 1375bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (nd++ < 9) 1376bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y *= 10; 1377bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (nd <= DBL_DIG + 1) 1378bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z *= 10; 1379bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (nd++ < 9) 1380bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = 10*y + c; 1381bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (nd <= DBL_DIG + 1) 1382bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = 10*z + c; 1383bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nz = 0; 1384bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1385bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1386bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1387bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dig_done: 1388bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e = 0; 1389bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c == 'e' || c == 'E') { 1390bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!nd && !nz && !nz0) { 1391bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret0; 1392bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1393bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s00 = s; 1394bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson esign = 0; 1395bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson switch(c = *++s) { 1396bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '-': 1397bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson esign = 1; 1398bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case '+': 1399bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = *++s; 1400bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1401bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c >= '0' && c <= '9') { 1402bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(c == '0') 1403bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson c = *++s; 1404bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (c > '0' && c <= '9') { 1405bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = c - '0'; 1406bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s1 = s; 1407bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while((c = *++s) >= '0' && c <= '9') 1408bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = 10*L + c - '0'; 1409bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s - s1 > 8 || L > 19999) 1410bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Avoid confusion from exponents 1411bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * so large that e might overflow. 1412bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1413bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e = 19999; /* safe for 16 bit ints */ 1414bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1415bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e = (int)L; 1416bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (esign) 1417bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e = -e; 1418bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1419bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1420bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e = 0; 1421bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1422bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1423bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s = s00; 1424bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1425bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!nd) { 1426bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!nz && !nz0) { 1427bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ret0: 1428bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s = s00; 1429bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson sign = 0; 1430bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1431bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1432bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1433bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.e0 = e1 = e -= nf; 1434bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1435bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Now we have nd0 digits, starting at s0, followed by a 1436bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * decimal point, followed by nd-nd0 digits. The number we're 1437bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * after is the integer represented by those digits times 1438bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 10**e */ 1439bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1440bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!nd0) 1441bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nd0 = nd; 1442bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 1443bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) = y; 1444bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k > 9) { 1445bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) = tens[k - 9] * dval(&rv) + z; 1446bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1447bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd0 = 0; 1448bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (nd <= DBL_DIG 1449bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson && Flt_Rounds == 1 1450bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ) { 1451bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!e) 1452bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1453bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e > 0) { 1454bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e <= Ten_pmax) { 1455bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= tens[e]; 1456bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1457bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1458bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = DBL_DIG - nd; 1459bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e <= Ten_pmax + i) { 1460bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* A fancier test would sometimes let us do 1461bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * this for larger i values. 1462bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1463bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e -= i; 1464bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= tens[i]; 1465bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= tens[e]; 1466bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1467bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1468bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1469bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (e >= -Ten_pmax) { 1470bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) /= tens[-e]; 1471bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1472bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1473bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1474bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e1 += nd - k; 1475bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1476bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.scale = 0; 1477bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1478bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Get starting approximation = rv * 10**e1 */ 1479bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1480bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 > 0) { 1481bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((i = e1 & 15)) 1482bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= tens[i]; 1483bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 &= ~15) { 1484bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 > DBL_MAX_10_EXP) { 1485bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ovfl: 1486bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson errno = ERANGE; 1487bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Can't trust HUGE_VAL */ 1488bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) = Exp_mask; 1489bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) = 0; 1490bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1491bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1492bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e1 >>= 4; 1493bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(j = 0; e1 > 1; j++, e1 >>= 1) 1494bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 & 1) 1495bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= bigtens[j]; 1496bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* The last multiplication could overflow. */ 1497bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) -= P*Exp_msk1; 1498bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= bigtens[j]; 1499bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((z = word0(&rv) & Exp_mask) 1500bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson > Exp_msk1*(DBL_MAX_EXP+Bias-P)) 1501bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ovfl; 1502bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { 1503bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* set to largest number */ 1504bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* (Can't trust DBL_MAX) */ 1505bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) = Big0; 1506bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) = Big1; 1507bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1508bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1509bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) += P*Exp_msk1; 1510bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1511bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1512bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (e1 < 0) { 1513bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e1 = -e1; 1514bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((i = e1 & 15)) 1515bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) /= tens[i]; 1516bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 >>= 4) { 1517bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 >= 1 << n_bigtens) 1518bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto undfl; 1519bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 & Scale_Bit) 1520bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.scale = 2*P; 1521bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(j = 0; e1 > 0; j++, e1 >>= 1) 1522bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e1 & 1) 1523bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= tinytens[j]; 1524bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) 1525bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson >> Exp_shift)) > 0) { 1526bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* scaled rv is denormal; clear j low bits */ 1527bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j >= 32) { 1528bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) = 0; 1529bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j >= 53) 1530bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) = (P+2)*Exp_msk1; 1531bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1532bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) &= 0xffffffff << (j-32); 1533bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1534bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1535bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) &= 0xffffffff << j; 1536bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1537bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!dval(&rv)) { 1538bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson undfl: 1539bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) = 0.; 1540bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson errno = ERANGE; 1541bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 1542bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1543bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1544bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1545bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1546bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Now the hard part -- adjusting rv to the correct value.*/ 1547bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1548bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Put digits into bd: true value = bd * 10^e */ 1549bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1550bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.nd = nd; 15515a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */ 1552bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* to silence an erroneous warning about bc.nd0 */ 1553bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* possibly not being initialized. */ 15545a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson if (nd > STRTOD_DIGLIM) { 15555a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */ 1556bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* minimum number of decimal digits to distinguish double values */ 1557bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* in IEEE arithmetic. */ 1558bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = j = 18; 1559bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i > nd0) 1560bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j += bc.dplen; 1561bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(;;) { 1562bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (--j <= bc.dp1 && j >= bc.dp0) 1563bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = bc.dp0 - 1; 1564bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s0[j] != '0') 1565bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1566bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson --i; 1567bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1568bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e += nd - i; 1569bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nd = i; 1570bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (nd0 > nd) 1571bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson nd0 = nd; 1572bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (nd < 9) { /* must recompute y */ 1573bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = 0; 1574bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 0; i < nd0; ++i) 1575bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = 10*y + s0[i] - '0'; 1576bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(j = bc.dp1; i < nd; ++i) 1577bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = 10*y + s0[j++] - '0'; 1578bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1579bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1580bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd0 = s2b(s0, nd0, nd, y, bc.dplen); 1581bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bd0 == NULL) 1582bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1583bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1584bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(;;) { 1585bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd = Balloc(bd0->k); 1586bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bd == NULL) { 1587bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1588bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1589bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1590bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bcopy(bd, bd0); 1591bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */ 1592bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bb == NULL) { 1593bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1594bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1595bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1596bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1597bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bs = i2b(1); 1598bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bs == NULL) { 1599bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1600bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1601bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1602bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1603bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1604bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1605bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (e >= 0) { 1606bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb2 = bb5 = 0; 1607bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd2 = bd5 = e; 1608bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1609bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1610bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb2 = bb5 = -e; 1611bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd2 = bd5 = 0; 1612bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1613bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bbe >= 0) 1614bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb2 += bbe; 1615bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1616bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd2 -= bbe; 1617bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bs2 = bb2; 1618bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = bbe - bc.scale; 1619bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = j + bbbits - 1; /* logb(rv) */ 1620bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i < Emin) /* denormal */ 1621bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j += P - Emin; 1622bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1623bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = P + 1 - bbbits; 1624bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb2 += j; 1625bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd2 += j; 1626bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd2 += bc.scale; 1627bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = bb2 < bd2 ? bb2 : bd2; 1628bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i > bs2) 1629bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = bs2; 1630bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i > 0) { 1631bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb2 -= i; 1632bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd2 -= i; 1633bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bs2 -= i; 1634bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1635bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bb5 > 0) { 1636bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bs = pow5mult(bs, bb5); 1637bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bs == NULL) { 1638bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1639bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1640bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1641bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1642bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1643bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb1 = mult(bs, bb); 1644bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1645bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb = bb1; 1646bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bb == NULL) { 1647bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1648bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1649bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1650bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1651bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1652bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1653bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bb2 > 0) { 1654bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bb = lshift(bb, bb2); 1655bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bb == NULL) { 1656bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1657bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1658bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1659bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1660bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1661bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1662bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bd5 > 0) { 1663bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd = pow5mult(bd, bd5); 1664bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bd == NULL) { 1665bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1666bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1667bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1668bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1669bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1670bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1671bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bd2 > 0) { 1672bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bd = lshift(bd, bd2); 1673bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bd == NULL) { 1674bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1675bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1676bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1677bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1678bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1679bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1680bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bs2 > 0) { 1681bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bs = lshift(bs, bs2); 1682bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bs == NULL) { 1683bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1684bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1685bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1686bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1687bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1688bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1689bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson delta = diff(bb, bd); 1690bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (delta == NULL) { 1691bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1692bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1693bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1694bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1695bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1696bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1697bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.dsign = delta->sign; 1698bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson delta->sign = 0; 1699bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = cmp(delta, bs); 1700bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.nd > nd && i <= 0) { 1701bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.dsign) 1702bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; /* Must use bigcomp(). */ 1703bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson { 1704bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.nd = nd; 1705bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = -1; /* Discarded digits make delta smaller. */ 1706bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1707bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1708bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1709bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i < 0) { 1710bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Error is less than half an ulp -- check for 1711bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * special case of mantissa a power of two. 1712bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1713bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask 1714bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 1715bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ) { 1716bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1717bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1718bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!delta->x[0] && delta->wds <= 1) { 1719bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* exact result */ 1720bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1721bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1722bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson delta = lshift(delta,Log2P); 1723bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (delta == NULL) { 1724bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1725bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1726bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1727bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1728bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1729bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1730bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (cmp(delta, bs) > 0) 1731bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto drop_down; 1732bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1733bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1734bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i == 0) { 1735bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* exactly half-way between */ 1736bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.dsign) { 1737bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 1738bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson && word1(&rv) == ( 1739bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson (bc.scale && 1740bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ? 1741bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : 1742bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 0xffffffff)) { 1743bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /*boundary case -- increment exponent*/ 1744bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) = (word0(&rv) & Exp_mask) 1745bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson + Exp_msk1 1746bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ; 1747bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) = 0; 1748bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.dsign = 0; 1749bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1750bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1751bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1752bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { 1753bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson drop_down: 1754bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* boundary case -- decrement exponent */ 1755bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.scale) { 1756bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = word0(&rv) & Exp_mask; 1757bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (L <= (2*P+1)*Exp_msk1) { 1758bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (L > (P+2)*Exp_msk1) 1759bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* round even ==> */ 1760bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* accept rv */ 1761bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1762bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* rv = smallest denormal */ 17635a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson if (bc.nd >nd) 1764bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1765bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto undfl; 1766bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1767bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1768bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = (word0(&rv) & Exp_mask) - Exp_msk1; 1769bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) = L | Bndry_mask1; 1770bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) = 0xffffffff; 1771bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1772bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1773bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(word1(&rv) & LSB)) 1774bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1775bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.dsign) 1776bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) += ulp(&rv); 1777bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1778bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) -= ulp(&rv); 1779bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!dval(&rv)) { 17805a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson if (bc.nd >nd) 1781bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1782bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto undfl; 1783bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1784bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1785bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bc.dsign = 1 - bc.dsign; 1786bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1787bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1788bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((aadj = ratio(delta, bs)) <= 2.) { 1789bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.dsign) 1790bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj = aadj1 = 1.; 1791bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (word1(&rv) || word0(&rv) & Bndry_mask) { 1792bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (word1(&rv) == Tiny1 && !word0(&rv)) { 17935a0b399aa9d75d4d208394d80b4997c164435ecbMark Dickinson if (bc.nd >nd) 1794bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1795bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto undfl; 1796bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1797bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj = 1.; 1798bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj1 = -1.; 1799bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1800bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1801bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* special case -- power of FLT_RADIX to be */ 1802bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* rounded down... */ 1803bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1804bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (aadj < 2./FLT_RADIX) 1805bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj = 1./FLT_RADIX; 1806bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1807bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj *= 0.5; 1808bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj1 = -aadj; 1809bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1810bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1811bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1812bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj *= 0.5; 1813bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj1 = bc.dsign ? aadj : -aadj; 1814bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (Flt_Rounds == 0) 1815bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj1 += 0.5; 1816bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1817bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson y = word0(&rv) & Exp_mask; 1818bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1819bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Check for overflow */ 1820bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1821bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { 1822bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv0) = dval(&rv); 1823bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) -= P*Exp_msk1; 1824bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson adj.d = aadj1 * ulp(&rv); 1825bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) += adj.d; 1826bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((word0(&rv) & Exp_mask) >= 1827bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Exp_msk1*(DBL_MAX_EXP+Bias-P)) { 1828bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (word0(&rv0) == Big0 && word1(&rv0) == Big1) 1829bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ovfl; 1830bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) = Big0; 1831bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv) = Big1; 1832bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto cont; 1833bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1834bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 1835bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv) += P*Exp_msk1; 1836bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1837bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 1838bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.scale && y <= 2*P*Exp_msk1) { 1839bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (aadj <= 0x7fffffff) { 1840bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((z = (ULong)aadj) <= 0) 1841bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = 1; 1842bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj = z; 1843bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj1 = bc.dsign ? aadj : -aadj; 1844bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1845bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&aadj2) = aadj1; 1846bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&aadj2) += (2*P+1)*Exp_msk1 - y; 1847bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj1 = dval(&aadj2); 1848bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1849bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson adj.d = aadj1 * ulp(&rv); 1850bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) += adj.d; 1851bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1852bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson z = word0(&rv) & Exp_mask; 1853bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.nd == nd) { 1854bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!bc.scale) 1855bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (y == z) { 1856bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Can we stop now? */ 1857bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = (Long)aadj; 1858bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson aadj -= L; 1859bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* The tolerances below are conservative. */ 1860bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) { 1861bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (aadj < .4999999 || aadj > .5000001) 1862bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1863bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1864bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (aadj < .4999999/FLT_RADIX) 1865bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 1866bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1867bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1868bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson cont: 1869bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1870bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1871bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1872bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(delta); 1873bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1874bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bb); 1875bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd); 1876bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bs); 1877bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(bd0); 1878bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(delta); 1879bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.nd > nd) { 1880bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson error = bigcomp(&rv, s0, &bc); 1881bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (error) 1882bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 1883bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1884bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1885bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (bc.scale) { 1886bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&rv0) = Exp_1 - 2*P*Exp_msk1; 1887bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word1(&rv0) = 0; 1888bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&rv) *= dval(&rv0); 1889bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* try to avoid the bug of testing an 8087 register value */ 1890bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(word0(&rv) & Exp_mask)) 1891bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson errno = ERANGE; 1892bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 1893bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ret: 1894bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (se) 1895bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *se = (char *)s; 1896bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return sign ? -dval(&rv) : dval(&rv); 1897bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1898bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson failed_malloc: 1899bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (se) 1900bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *se = (char *)s00; 1901bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson errno = ENOMEM; 1902bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return -1.0; 1903bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1904bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1905bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic char * 1906bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonrv_alloc(int i) 1907bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1908bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int j, k, *r; 1909bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1910bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = sizeof(ULong); 1911bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(k = 0; 1912bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; 1913bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j <<= 1) 1914bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 1915bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson r = (int*)Balloc(k); 1916bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (r == NULL) 1917bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 1918bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *r = k; 1919bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return (char *)(r+1); 1920bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1921bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1922bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonstatic char * 1923bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonnrv_alloc(char *s, char **rve, int n) 1924bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1925bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson char *rv, *t; 1926bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1927bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson rv = rv_alloc(n); 1928bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (rv == NULL) 1929bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 1930bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson t = rv; 1931bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while((*t = *s++)) t++; 1932bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (rve) 1933bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *rve = t; 1934bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return rv; 1935bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1936bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1937bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* freedtoa(s) must be used to free values s returned by dtoa 1938bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * when MULTIPLE_THREADS is #defined. It should be used in all cases, 1939bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * but for consistency with earlier versions of dtoa, it is optional 1940bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * when MULTIPLE_THREADS is not defined. 1941bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1942bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1943bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonvoid 1944bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson_Py_dg_freedtoa(char *s) 1945bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1946bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b = (Bigint *)((int *)s - 1); 1947bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b->maxwds = 1 << (b->k = *(int*)b); 1948bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 1949bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 1950bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1951bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. 1952bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 1953bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Inspired by "How to Print Floating-Point Numbers Accurately" by 1954bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. 1955bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 1956bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Modifications: 1957bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 1. Rather than iterating, we use a simple numeric overestimate 1958bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * to determine k = floor(log10(d)). We scale relevant 1959bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * quantities using O(log2(k)) rather than O(k) multiplications. 1960bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't 1961bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * try to generate digits strictly left to right. Instead, we 1962bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * compute with fewer bits and propagate the carry if necessary 1963bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * when rounding the final digit up. This is often faster. 1964bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 3. Under the assumption that input will be rounded nearest, 1965bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. 1966bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * That is, we allow equality in stopping tests when the 1967bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * round-nearest rule will give the same floating-point value 1968bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * as would satisfaction of the stopping test with strict 1969bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * inequality. 1970bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 4. We remove common factors of powers of 2 from relevant 1971bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * quantities. 1972bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 5. When converting floating-point integers less than 1e16, 1973bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * we use floating-point arithmetic rather than resorting 1974bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * to multiple-precision integers. 1975bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 6. When asked to produce fewer than 15 digits, we first try 1976bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * to get by with floating-point arithmetic; we resort to 1977bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * multiple-precision integer arithmetic only if we cannot 1978bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * guarantee that the floating-point calculation has given 1979bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * the correctly rounded result. For k requested digits and 1980bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * "uniformly" distributed input, the probability is 1981bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * something like 10^(k-15) that we must resort to the Long 1982bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * calculation. 1983bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 1984bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1985bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson/* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory 1986bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson leakage, a successful call to _Py_dg_dtoa should always be matched by a 1987bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson call to _Py_dg_freedtoa. */ 1988bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1989bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinsonchar * 1990bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson_Py_dg_dtoa(double dd, int mode, int ndigits, 1991bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int *decpt, int *sign, char **rve) 1992bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson{ 1993bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Arguments ndigits, decpt, sign are similar to those 1994bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson of ecvt and fcvt; trailing zeros are suppressed from 1995bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson the returned string. If not null, *rve is set to point 1996bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson to the end of the return value. If d is +-Infinity or NaN, 1997bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson then *decpt is set to 9999. 1998bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1999bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mode: 2000bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 0 ==> shortest string that yields d when read in 2001bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson and rounded to nearest. 2002bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1 ==> like 0, but with Steele & White stopping rule; 2003bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson e.g. with IEEE P754 arithmetic , mode 0 gives 2004bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1e23 whereas mode 1 gives 9.999999999999999e22. 2005bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2 ==> max(1,ndigits) significant digits. This gives a 2006bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return value similar to that of ecvt, except 2007bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson that trailing zeros are suppressed. 2008bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 3 ==> through ndigits past the decimal point. This 2009bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson gives a return value similar to that from fcvt, 2010bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson except that trailing zeros are suppressed, and 2011bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ndigits can be negative. 2012bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 4,5 ==> similar to 2 and 3, respectively, but (in 2013bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson round-nearest mode) with the tests of mode 0 to 2014bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson possibly return a shorter string that rounds to d. 2015bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson With IEEE arithmetic and compilation with 2016bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same 2017bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson as modes 2 and 3 when FLT_ROUNDS != 1. 2018bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 6-9 ==> Debugging modes similar to mode - 4: don't try 2019bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson fast floating-point estimate (if applicable). 2020bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2021bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Values of mode other than 0-9 are treated as mode 0. 2022bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2023bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Sufficient space is allocated to the return value 2024bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson to hold the suppressed trailing zeros. 2025bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 2026bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2027bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, 2028bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, 2029bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson spec_case, try_quick; 2030bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Long L; 2031bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson int denorm; 2032bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ULong x; 2033bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bigint *b, *b1, *delta, *mlo, *mhi, *S; 2034bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson U d2, eps, u; 2035bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson double ds; 2036bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson char *s, *s0; 2037bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2038bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* set pointers to NULL, to silence gcc compiler warnings and make 2039bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson cleanup easier on error */ 2040bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mlo = mhi = b = S = 0; 2041bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s0 = 0; 2042bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2043bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson u.d = dd; 2044bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (word0(&u) & Sign_bit) { 2045bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* set sign for everything, including 0's and NaNs */ 2046bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *sign = 1; 2047bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&u) &= ~Sign_bit; /* clear sign bit */ 2048bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2049bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 2050bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *sign = 0; 2051bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2052bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* quick return for Infinities, NaNs and zeros */ 2053bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((word0(&u) & Exp_mask) == Exp_mask) 2054bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson { 2055bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Infinity or NaN */ 2056bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *decpt = 9999; 2057bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!word1(&u) && !(word0(&u) & 0xfffff)) 2058bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return nrv_alloc("Infinity", rve, 8); 2059bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return nrv_alloc("NaN", rve, 3); 2060bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2061bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!dval(&u)) { 2062bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *decpt = 1; 2063bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return nrv_alloc("0", rve, 1); 2064bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2065bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2066bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* compute k = floor(log10(d)). The computation may leave k 2067bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson one too large, but should never leave k too small. */ 2068bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = d2b(&u, &be, &bbits); 2069bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2070bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2071bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { 2072bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&d2) = dval(&u); 2073bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&d2) &= Frac_mask1; 2074bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&d2) |= Exp_11; 2075bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2076bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 2077bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * log10(x) = log(x) / log(10) 2078bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) 2079bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) 2080bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2081bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * This suggests computing an approximation k to log10(d) by 2082bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2083bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * k = (i - Bias)*0.301029995663981 2084bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); 2085bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2086bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * We want k to be too large rather than too small. 2087bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * The error in the first-order Taylor series approximation 2088bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * is in our favor, so we just round up the constant enough 2089bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * to compensate for any error in the multiplication of 2090bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, 2091bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, 2092bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * adding 1e-13 to the constant term more than suffices. 2093bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Hence we adjust the constant term to 0.1760912590558. 2094bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * (We could get a more accurate k by invoking log10, 2095bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * but this is probably not worthwhile.) 2096bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 2097bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2098bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i -= Bias; 2099bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson denorm = 0; 2100bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2101bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2102bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* d is denormalized */ 2103bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2104bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = bbits + be + (Bias + (P-1) - 1); 2105bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) 2106bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson : word1(&u) << (32 - i); 2107bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&d2) = x; 2108bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ 2109bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i -= (Bias + (P-1) - 1) + 1; 2110bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson denorm = 1; 2111bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2112bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + 2113bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i*0.301029995663981; 2114bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = (int)ds; 2115bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ds < 0. && ds != k) 2116bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k--; /* want k = floor(ds) */ 2117bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k_check = 1; 2118bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k >= 0 && k <= Ten_pmax) { 2119bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dval(&u) < tens[k]) 2120bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k--; 2121bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k_check = 0; 2122bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2123bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = bbits - i - 1; 2124bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j >= 0) { 2125bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 = 0; 2126bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 = j; 2127bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2128bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2129bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 = -j; 2130bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 = 0; 2131bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2132bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k >= 0) { 2133bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b5 = 0; 2134bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s5 = k; 2135bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 += k; 2136bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2137bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2138bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 -= k; 2139bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b5 = -k; 2140bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s5 = 0; 2141bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2142bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mode < 0 || mode > 9) 2143bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mode = 0; 2144bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2145bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson try_quick = 1; 2146bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2147bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mode > 5) { 2148bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mode -= 4; 2149bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson try_quick = 0; 2150bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2151bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson leftright = 1; 2152bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ 2153bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* silence erroneous "gcc -Wall" warning. */ 2154bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson switch(mode) { 2155bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 0: 2156bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 1: 2157bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 18; 2158bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ndigits = 0; 2159bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2160bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 2: 2161bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson leftright = 0; 2162bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* no break */ 2163bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 4: 2164bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ndigits <= 0) 2165bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ndigits = 1; 2166bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = ilim1 = i = ndigits; 2167bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2168bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 3: 2169bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson leftright = 0; 2170bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* no break */ 2171bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson case 5: 2172bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = ndigits + k + 1; 2173bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = i; 2174bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim1 = i - 1; 2175bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i <= 0) 2176bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 1; 2177bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2178bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s0 = rv_alloc(i); 2179bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s0 == NULL) 2180bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2181bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s = s0; 2182bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2183bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2184bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ilim >= 0 && ilim <= Quick_max && try_quick) { 2185bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2186bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Try to get by with floating-point arithmetic. */ 2187bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2188bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 0; 2189bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&d2) = dval(&u); 2190bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k0 = k; 2191bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim0 = ilim; 2192bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ieps = 2; /* conservative */ 2193bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k > 0) { 2194bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ds = tens[k&0xf]; 2195bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = k >> 4; 2196bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j & Bletch) { 2197bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* prevent overflows */ 2198bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j &= Bletch - 1; 2199bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) /= bigtens[n_bigtens-1]; 2200bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ieps++; 2201bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2202bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(; j; j >>= 1, i++) 2203bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j & 1) { 2204bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ieps++; 2205bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ds *= bigtens[i]; 2206bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2207bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) /= ds; 2208bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2209bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if ((j1 = -k)) { 2210bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) *= tens[j1 & 0xf]; 2211bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(j = j1 >> 4; j; j >>= 1, i++) 2212bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j & 1) { 2213bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ieps++; 2214bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) *= bigtens[i]; 2215bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2216bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2217bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k_check && dval(&u) < 1. && ilim > 0) { 2218bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ilim1 <= 0) 2219bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto fast_failed; 2220bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = ilim1; 2221bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k--; 2222bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) *= 10.; 2223bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ieps++; 2224bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2225bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&eps) = ieps*dval(&u) + 7.; 2226bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson word0(&eps) -= (P-1)*Exp_msk1; 2227bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ilim == 0) { 2228bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson S = mhi = 0; 2229bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) -= 5.; 2230bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dval(&u) > dval(&eps)) 2231bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto one_digit; 2232bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dval(&u) < -dval(&eps)) 2233bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto no_digits; 2234bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto fast_failed; 2235bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2236bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (leftright) { 2237bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Use Steele & White method of only 2238bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * generating digits needed. 2239bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 2240bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); 2241bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 0;;) { 2242bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = (Long)dval(&u); 2243bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) -= L; 2244bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = '0' + (int)L; 2245bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dval(&u) < dval(&eps)) 2246bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret1; 2247bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (1. - dval(&u) < dval(&eps)) 2248bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto bump_up; 2249bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (++i >= ilim) 2250bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2251bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&eps) *= 10.; 2252bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) *= 10.; 2253bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2254bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2255bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2256bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Generate ilim digits, then fix them up. */ 2257bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&eps) *= tens[ilim-1]; 2258bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 1;; i++, dval(&u) *= 10.) { 2259bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = (Long)(dval(&u)); 2260bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!(dval(&u) -= L)) 2261bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = i; 2262bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = '0' + (int)L; 2263bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i == ilim) { 2264bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dval(&u) > 0.5 + dval(&eps)) 2265bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto bump_up; 2266bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else if (dval(&u) < 0.5 - dval(&eps)) { 2267bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(*--s == '0'); 2268bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s++; 2269bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret1; 2270bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2271bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2272bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2273bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2274bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2275bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson fast_failed: 2276bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s = s0; 2277bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) = dval(&d2); 2278bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = k0; 2279bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = ilim0; 2280bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2281bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2282bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Do we have a "small" integer? */ 2283bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2284bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (be >= 0 && k <= Int_max) { 2285bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Yes. */ 2286bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ds = tens[k]; 2287bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ndigits < 0 && ilim <= 0) { 2288bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson S = mhi = 0; 2289bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ilim < 0 || dval(&u) <= 5*ds) 2290bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto no_digits; 2291bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto one_digit; 2292bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2293bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 1;; i++, dval(&u) *= 10.) { 2294bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson L = (Long)(dval(&u) / ds); 2295bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) -= L*ds; 2296bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = '0' + (int)L; 2297bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!dval(&u)) { 2298bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2299bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2300bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i == ilim) { 2301bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dval(&u) += dval(&u); 2302bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dval(&u) > ds || (dval(&u) == ds && L & 1)) { 2303bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson bump_up: 2304bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(*--s == '9') 2305bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s == s0) { 2306bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 2307bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s = '0'; 2308bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2309bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2310bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ++*s++; 2311bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2312bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2313bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2314bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2315bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret1; 2316bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2317bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2318bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson m2 = b2; 2319bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson m5 = b5; 2320bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (leftright) { 2321bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 2322bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson denorm ? be + (Bias + (P-1) - 1 + 1) : 2323bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 1 + P - bbits; 2324bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 += i; 2325bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 += i; 2326bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = i2b(1); 2327bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2328bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2329bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2330bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (m2 > 0 && s2 > 0) { 2331bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = m2 < s2 ? m2 : s2; 2332bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 -= i; 2333bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson m2 -= i; 2334bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 -= i; 2335bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2336bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b5 > 0) { 2337bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (leftright) { 2338bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (m5 > 0) { 2339bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = pow5mult(mhi, m5); 2340bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2341bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2342bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b1 = mult(mhi, b); 2343bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 2344bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = b1; 2345bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2346bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2347bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2348bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((j = b5 - m5)) { 2349bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = pow5mult(b, j); 2350bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2351bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2352bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2353bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2354bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2355bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = pow5mult(b, b5); 2356bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2357bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2358bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2359bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2360bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson S = i2b(1); 2361bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (S == NULL) 2362bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2363bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s5 > 0) { 2364bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson S = pow5mult(S, s5); 2365bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (S == NULL) 2366bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2367bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2368bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2369bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Check for special case that d is a normalized power of 2. */ 2370bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2371bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson spec_case = 0; 2372bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((mode < 2 || leftright) 2373bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ) { 2374bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!word1(&u) && !(word0(&u) & Bndry_mask) 2375bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson && word0(&u) & (Exp_mask & ~Exp_msk1) 2376bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ) { 2377bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* The special case */ 2378bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 += Log2P; 2379bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 += Log2P; 2380bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson spec_case = 1; 2381bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2382bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2383bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2384bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Arrange for convenient computation of quotients: 2385bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * shift left if necessary so divisor has 4 leading 0 bits. 2386bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * 2387bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * Perhaps we should just compute leading 28 bits of S once 2388bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * and for all and pass them and a shift to quorem, so it 2389bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * can do shifts and ors to compute the numerator for q. 2390bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 2391bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((i = ((s5 ? 32 - hi0bits(S->x[S->wds-1]) : 1) + s2) & 0x1f)) 2392bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = 32 - i; 2393bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#define iInc 28 2394bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson i = dshift(S, s2); 2395bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b2 += i; 2396bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson m2 += i; 2397bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s2 += i; 2398bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b2 > 0) { 2399bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = lshift(b, b2); 2400bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2401bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2402bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2403bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s2 > 0) { 2404bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson S = lshift(S, s2); 2405bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (S == NULL) 2406bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2407bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2408bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (k_check) { 2409bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (cmp(b,S) < 0) { 2410bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k--; 2411bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, 0); /* we botched the k estimate */ 2412bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2413bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2414bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (leftright) { 2415bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = multadd(mhi, 10, 0); 2416bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2417bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2418bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2419bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ilim = ilim1; 2420bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2421bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2422bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ilim <= 0 && (mode == 3 || mode == 5)) { 2423bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (ilim < 0) { 2424bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* no digits, fcvt style */ 2425bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson no_digits: 2426bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k = -1 - ndigits; 2427bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2428bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2429bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2430bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson S = multadd(S, 5, 0); 2431bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (S == NULL) 2432bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2433bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (cmp(b, S) <= 0) 2434bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto no_digits; 2435bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2436bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson one_digit: 2437bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = '1'; 2438bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 2439bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2440bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2441bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (leftright) { 2442bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (m2 > 0) { 2443bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = lshift(mhi, m2); 2444bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2445bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2446bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2447bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2448bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Compute mlo -- check for special case 2449bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * that d is a normalized power of 2. 2450bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 2451bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2452bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mlo = mhi; 2453bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (spec_case) { 2454bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = Balloc(mhi->k); 2455bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2456bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2457bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bcopy(mhi, mlo); 2458bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = lshift(mhi, Log2P); 2459bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2460bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2461bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2462bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2463bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 1;;i++) { 2464bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dig = quorem(b,S) + '0'; 2465bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Do we yet have the shortest decimal string 2466bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson * that will round to d? 2467bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson */ 2468bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = cmp(b, mlo); 2469bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson delta = diff(S, mhi); 2470bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (delta == NULL) 2471bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2472bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j1 = delta->sign ? 1 : cmp(b, delta); 2473bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(delta); 2474bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j1 == 0 && mode != 1 && !(word1(&u) & 1) 2475bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ) { 2476bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dig == '9') 2477bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto round_9_up; 2478bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j > 0) 2479bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson dig++; 2480bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = dig; 2481bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2482bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2483bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j < 0 || (j == 0 && mode != 1 2484bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson && !(word1(&u) & 1) 2485bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson )) { 2486bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!b->x[0] && b->wds <= 1) { 2487bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto accept_dig; 2488bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2489bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j1 > 0) { 2490bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = lshift(b, 1); 2491bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2492bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2493bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j1 = cmp(b, S); 2494bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if ((j1 > 0 || (j1 == 0 && dig & 1)) 2495bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson && dig++ == '9') 2496bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto round_9_up; 2497bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2498bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson accept_dig: 2499bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = dig; 2500bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2501bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2502bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j1 > 0) { 2503bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (dig == '9') { /* possible if i == 1 */ 2504bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson round_9_up: 2505bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = '9'; 2506bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto roundoff; 2507bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2508bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = dig + 1; 2509bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2510bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2511bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = dig; 2512bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i == ilim) 2513bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2514bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, 0); 2515bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2516bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2517bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mlo == mhi) { 2518bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mlo = mhi = multadd(mhi, 10, 0); 2519bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mlo == NULL) 2520bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2521bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2522bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2523bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mlo = multadd(mlo, 10, 0); 2524bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mlo == NULL) 2525bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2526bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson mhi = multadd(mhi, 10, 0); 2527bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi == NULL) 2528bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2529bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2530bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2531bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2532bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else 2533bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson for(i = 1;; i++) { 2534bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = dig = quorem(b,S) + '0'; 2535bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (!b->x[0] && b->wds <= 1) { 2536bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2537bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2538bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (i >= ilim) 2539bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson break; 2540bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = multadd(b, 10, 0); 2541bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2542bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2543bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2544bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2545bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson /* Round off last digit */ 2546bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2547bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson b = lshift(b, 1); 2548bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b == NULL) 2549bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto failed_malloc; 2550bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson j = cmp(b, S); 2551bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (j > 0 || (j == 0 && dig & 1)) { 2552bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson roundoff: 2553bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(*--s == '9') 2554bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s == s0) { 2555bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson k++; 2556bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s++ = '1'; 2557bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson goto ret; 2558bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2559bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ++*s++; 2560bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2561bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson else { 2562bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson while(*--s == '0'); 2563bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson s++; 2564bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2565bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ret: 2566bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(S); 2567bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi) { 2568bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mlo && mlo != mhi) 2569bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(mlo); 2570bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(mhi); 2571bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson } 2572bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson ret1: 2573bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 2574bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *s = 0; 2575bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *decpt = k + 1; 2576bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (rve) 2577bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson *rve = s; 2578bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return s0; 2579bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson failed_malloc: 2580bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (S) 2581bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(S); 2582bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mlo && mlo != mhi) 2583bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(mlo); 2584bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (mhi) 2585bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(mhi); 2586bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (b) 2587bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson Bfree(b); 2588bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson if (s0) 2589bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson _Py_dg_freedtoa(s0); 2590bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson return NULL; 2591bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 2592bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#ifdef __cplusplus 2593bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson} 2594bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif 2595bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson 2596bb28285ea2f01e97a26bc595d49da43fbee62913Mark Dickinson#endif /* PY_NO_SHORT_FLOAT_REPR */ 2597