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