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