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