| 1 | /****************************************************************
|
|---|
| 2 | *
|
|---|
| 3 | * The author of this software is David M. Gay.
|
|---|
| 4 | *
|
|---|
| 5 | * Copyright (c) 1991 by AT&T.
|
|---|
| 6 | *
|
|---|
| 7 | * Permission to use, copy, modify, and distribute this software for any
|
|---|
| 8 | * purpose without fee is hereby granted, provided that this entire notice
|
|---|
| 9 | * is included in all copies of any software which is or includes a copy
|
|---|
| 10 | * or modification of this software and in all copies of the supporting
|
|---|
| 11 | * documentation for such software.
|
|---|
| 12 | *
|
|---|
| 13 | * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
|---|
| 14 | * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
|
|---|
| 15 | * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
|---|
| 16 | * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
|---|
| 17 | *
|
|---|
| 18 | ***************************************************************/
|
|---|
| 19 |
|
|---|
| 20 | /* Please send bug reports to
|
|---|
| 21 | David M. Gay
|
|---|
| 22 | AT&T Bell Laboratories, Room 2C-463
|
|---|
| 23 | 600 Mountain Avenue
|
|---|
| 24 | Murray Hill, NJ 07974-2070
|
|---|
| 25 | U.S.A.
|
|---|
| 26 | dmg@research.att.com or research!dmg
|
|---|
| 27 | */
|
|---|
| 28 |
|
|---|
| 29 | #include "mprec.h"
|
|---|
| 30 | #include <string.h>
|
|---|
| 31 |
|
|---|
| 32 | static int
|
|---|
| 33 | _DEFUN (quorem,
|
|---|
| 34 | (b, S),
|
|---|
| 35 | _Jv_Bigint * b _AND _Jv_Bigint * S)
|
|---|
| 36 | {
|
|---|
| 37 | int n;
|
|---|
| 38 | long borrow, y;
|
|---|
| 39 | unsigned long carry, q, ys;
|
|---|
| 40 | unsigned long *bx, *bxe, *sx, *sxe;
|
|---|
| 41 | #ifdef Pack_32
|
|---|
| 42 | long z;
|
|---|
| 43 | unsigned long si, zs;
|
|---|
| 44 | #endif
|
|---|
| 45 |
|
|---|
| 46 | n = S->_wds;
|
|---|
| 47 | #ifdef DEBUG
|
|---|
| 48 | /*debug*/ if (b->_wds > n)
|
|---|
| 49 | /*debug*/ Bug ("oversize b in quorem");
|
|---|
| 50 | #endif
|
|---|
| 51 | if (b->_wds < n)
|
|---|
| 52 | return 0;
|
|---|
| 53 | sx = S->_x;
|
|---|
| 54 | sxe = sx + --n;
|
|---|
| 55 | bx = b->_x;
|
|---|
| 56 | bxe = bx + n;
|
|---|
| 57 | q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
|
|---|
| 58 | #ifdef DEBUG
|
|---|
| 59 | /*debug*/ if (q > 9)
|
|---|
| 60 | /*debug*/ Bug ("oversized quotient in quorem");
|
|---|
| 61 | #endif
|
|---|
| 62 | if (q)
|
|---|
| 63 | {
|
|---|
| 64 | borrow = 0;
|
|---|
| 65 | carry = 0;
|
|---|
| 66 | do
|
|---|
| 67 | {
|
|---|
| 68 | #ifdef Pack_32
|
|---|
| 69 | si = *sx++;
|
|---|
| 70 | ys = (si & 0xffff) * q + carry;
|
|---|
| 71 | zs = (si >> 16) * q + (ys >> 16);
|
|---|
| 72 | carry = zs >> 16;
|
|---|
| 73 | y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
|
|---|
| 74 | borrow = y >> 16;
|
|---|
| 75 | Sign_Extend (borrow, y);
|
|---|
| 76 | z = (*bx >> 16) - (zs & 0xffff) + borrow;
|
|---|
| 77 | borrow = z >> 16;
|
|---|
| 78 | Sign_Extend (borrow, z);
|
|---|
| 79 | Storeinc (bx, z, y);
|
|---|
| 80 | #else
|
|---|
| 81 | ys = *sx++ * q + carry;
|
|---|
| 82 | carry = ys >> 16;
|
|---|
| 83 | y = *bx - (ys & 0xffff) + borrow;
|
|---|
| 84 | borrow = y >> 16;
|
|---|
| 85 | Sign_Extend (borrow, y);
|
|---|
| 86 | *bx++ = y & 0xffff;
|
|---|
| 87 | #endif
|
|---|
| 88 | }
|
|---|
| 89 | while (sx <= sxe);
|
|---|
| 90 | if (!*bxe)
|
|---|
| 91 | {
|
|---|
| 92 | bx = b->_x;
|
|---|
| 93 | while (--bxe > bx && !*bxe)
|
|---|
| 94 | --n;
|
|---|
| 95 | b->_wds = n;
|
|---|
| 96 | }
|
|---|
| 97 | }
|
|---|
| 98 | if (cmp (b, S) >= 0)
|
|---|
| 99 | {
|
|---|
| 100 | q++;
|
|---|
| 101 | borrow = 0;
|
|---|
| 102 | carry = 0;
|
|---|
| 103 | bx = b->_x;
|
|---|
| 104 | sx = S->_x;
|
|---|
| 105 | do
|
|---|
| 106 | {
|
|---|
| 107 | #ifdef Pack_32
|
|---|
| 108 | si = *sx++;
|
|---|
| 109 | ys = (si & 0xffff) + carry;
|
|---|
| 110 | zs = (si >> 16) + (ys >> 16);
|
|---|
| 111 | carry = zs >> 16;
|
|---|
| 112 | y = (*bx & 0xffff) - (ys & 0xffff) + borrow;
|
|---|
| 113 | borrow = y >> 16;
|
|---|
| 114 | Sign_Extend (borrow, y);
|
|---|
| 115 | z = (*bx >> 16) - (zs & 0xffff) + borrow;
|
|---|
| 116 | borrow = z >> 16;
|
|---|
| 117 | Sign_Extend (borrow, z);
|
|---|
| 118 | Storeinc (bx, z, y);
|
|---|
| 119 | #else
|
|---|
| 120 | ys = *sx++ + carry;
|
|---|
| 121 | carry = ys >> 16;
|
|---|
| 122 | y = *bx - (ys & 0xffff) + borrow;
|
|---|
| 123 | borrow = y >> 16;
|
|---|
| 124 | Sign_Extend (borrow, y);
|
|---|
| 125 | *bx++ = y & 0xffff;
|
|---|
| 126 | #endif
|
|---|
| 127 | }
|
|---|
| 128 | while (sx <= sxe);
|
|---|
| 129 | bx = b->_x;
|
|---|
| 130 | bxe = bx + n;
|
|---|
| 131 | if (!*bxe)
|
|---|
| 132 | {
|
|---|
| 133 | while (--bxe > bx && !*bxe)
|
|---|
| 134 | --n;
|
|---|
| 135 | b->_wds = n;
|
|---|
| 136 | }
|
|---|
| 137 | }
|
|---|
| 138 | return q;
|
|---|
| 139 | }
|
|---|
| 140 |
|
|---|
| 141 | #ifdef DEBUG
|
|---|
| 142 | #include <stdio.h>
|
|---|
| 143 |
|
|---|
| 144 | void
|
|---|
| 145 | print (_Jv_Bigint * b)
|
|---|
| 146 | {
|
|---|
| 147 | int i, wds;
|
|---|
| 148 | unsigned long *x, y;
|
|---|
| 149 | wds = b->_wds;
|
|---|
| 150 | x = b->_x+wds;
|
|---|
| 151 | i = 0;
|
|---|
| 152 | do
|
|---|
| 153 | {
|
|---|
| 154 | x--;
|
|---|
| 155 | fprintf (stderr, "%08x", *x);
|
|---|
| 156 | }
|
|---|
| 157 | while (++i < wds);
|
|---|
| 158 | fprintf (stderr, "\n");
|
|---|
| 159 | }
|
|---|
| 160 | #endif
|
|---|
| 161 |
|
|---|
| 162 | /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
|
|---|
| 163 | *
|
|---|
| 164 | * Inspired by "How to Print Floating-Point Numbers Accurately" by
|
|---|
| 165 | * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].
|
|---|
| 166 | *
|
|---|
| 167 | * Modifications:
|
|---|
| 168 | * 1. Rather than iterating, we use a simple numeric overestimate
|
|---|
| 169 | * to determine k = floor(log10(d)). We scale relevant
|
|---|
| 170 | * quantities using O(log2(k)) rather than O(k) multiplications.
|
|---|
| 171 | * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
|
|---|
| 172 | * try to generate digits strictly left to right. Instead, we
|
|---|
| 173 | * compute with fewer bits and propagate the carry if necessary
|
|---|
| 174 | * when rounding the final digit up. This is often faster.
|
|---|
| 175 | * 3. Under the assumption that input will be rounded nearest,
|
|---|
| 176 | * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
|
|---|
| 177 | * That is, we allow equality in stopping tests when the
|
|---|
| 178 | * round-nearest rule will give the same floating-point value
|
|---|
| 179 | * as would satisfaction of the stopping test with strict
|
|---|
| 180 | * inequality.
|
|---|
| 181 | * 4. We remove common factors of powers of 2 from relevant
|
|---|
| 182 | * quantities.
|
|---|
| 183 | * 5. When converting floating-point integers less than 1e16,
|
|---|
| 184 | * we use floating-point arithmetic rather than resorting
|
|---|
| 185 | * to multiple-precision integers.
|
|---|
| 186 | * 6. When asked to produce fewer than 15 digits, we first try
|
|---|
| 187 | * to get by with floating-point arithmetic; we resort to
|
|---|
| 188 | * multiple-precision integer arithmetic only if we cannot
|
|---|
| 189 | * guarantee that the floating-point calculation has given
|
|---|
| 190 | * the correctly rounded result. For k requested digits and
|
|---|
| 191 | * "uniformly" distributed input, the probability is
|
|---|
| 192 | * something like 10^(k-15) that we must resort to the long
|
|---|
| 193 | * calculation.
|
|---|
| 194 | */
|
|---|
| 195 |
|
|---|
| 196 |
|
|---|
| 197 | char *
|
|---|
| 198 | _DEFUN (_dtoa_r,
|
|---|
| 199 | (ptr, _d, mode, ndigits, decpt, sign, rve, float_type),
|
|---|
| 200 | struct _Jv_reent *ptr _AND
|
|---|
| 201 | double _d _AND
|
|---|
| 202 | int mode _AND
|
|---|
| 203 | int ndigits _AND
|
|---|
| 204 | int *decpt _AND
|
|---|
| 205 | int *sign _AND
|
|---|
| 206 | char **rve _AND
|
|---|
| 207 | int float_type)
|
|---|
| 208 | {
|
|---|
| 209 | /*
|
|---|
| 210 | float_type == 0 for double precision, 1 for float.
|
|---|
| 211 |
|
|---|
| 212 | Arguments ndigits, decpt, sign are similar to those
|
|---|
| 213 | of ecvt and fcvt; trailing zeros are suppressed from
|
|---|
| 214 | the returned string. If not null, *rve is set to point
|
|---|
| 215 | to the end of the return value. If d is +-Infinity or NaN,
|
|---|
| 216 | then *decpt is set to 9999.
|
|---|
| 217 |
|
|---|
| 218 | mode:
|
|---|
| 219 | 0 ==> shortest string that yields d when read in
|
|---|
| 220 | and rounded to nearest.
|
|---|
| 221 | 1 ==> like 0, but with Steele & White stopping rule;
|
|---|
| 222 | e.g. with IEEE P754 arithmetic , mode 0 gives
|
|---|
| 223 | 1e23 whereas mode 1 gives 9.999999999999999e22.
|
|---|
| 224 | 2 ==> max(1,ndigits) significant digits. This gives a
|
|---|
| 225 | return value similar to that of ecvt, except
|
|---|
| 226 | that trailing zeros are suppressed.
|
|---|
| 227 | 3 ==> through ndigits past the decimal point. This
|
|---|
| 228 | gives a return value similar to that from fcvt,
|
|---|
| 229 | except that trailing zeros are suppressed, and
|
|---|
| 230 | ndigits can be negative.
|
|---|
| 231 | 4-9 should give the same return values as 2-3, i.e.,
|
|---|
| 232 | 4 <= mode <= 9 ==> same return as mode
|
|---|
| 233 | 2 + (mode & 1). These modes are mainly for
|
|---|
| 234 | debugging; often they run slower but sometimes
|
|---|
| 235 | faster than modes 2-3.
|
|---|
| 236 | 4,5,8,9 ==> left-to-right digit generation.
|
|---|
| 237 | 6-9 ==> don't try fast floating-point estimate
|
|---|
| 238 | (if applicable).
|
|---|
| 239 |
|
|---|
| 240 | > 16 ==> Floating-point arg is treated as single precision.
|
|---|
| 241 |
|
|---|
| 242 | Values of mode other than 0-9 are treated as mode 0.
|
|---|
| 243 |
|
|---|
| 244 | Sufficient space is allocated to the return value
|
|---|
| 245 | to hold the suppressed trailing zeros.
|
|---|
| 246 | */
|
|---|
| 247 |
|
|---|
| 248 | int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, j, j1, k, k0,
|
|---|
| 249 | k_check, leftright, m2, m5, s2, s5, spec_case, try_quick;
|
|---|
| 250 | union double_union d, d2, eps;
|
|---|
| 251 | long L;
|
|---|
| 252 | #ifndef Sudden_Underflow
|
|---|
| 253 | int denorm;
|
|---|
| 254 | unsigned long x;
|
|---|
| 255 | #endif
|
|---|
| 256 | _Jv_Bigint *b, *b1, *delta, *mlo, *mhi, *S;
|
|---|
| 257 | double ds;
|
|---|
| 258 | char *s, *s0;
|
|---|
| 259 |
|
|---|
| 260 | d.d = _d;
|
|---|
| 261 |
|
|---|
| 262 | if (ptr->_result)
|
|---|
| 263 | {
|
|---|
| 264 | ptr->_result->_k = ptr->_result_k;
|
|---|
| 265 | ptr->_result->_maxwds = 1 << ptr->_result_k;
|
|---|
| 266 | Bfree (ptr, ptr->_result);
|
|---|
| 267 | ptr->_result = 0;
|
|---|
| 268 | }
|
|---|
| 269 |
|
|---|
| 270 | if (word0 (d) & Sign_bit)
|
|---|
| 271 | {
|
|---|
| 272 | /* set sign for everything, including 0's and NaNs */
|
|---|
| 273 | *sign = 1;
|
|---|
| 274 | word0 (d) &= ~Sign_bit; /* clear sign bit */
|
|---|
| 275 | }
|
|---|
| 276 | else
|
|---|
| 277 | *sign = 0;
|
|---|
| 278 |
|
|---|
| 279 | #if defined(IEEE_Arith) + defined(VAX)
|
|---|
| 280 | #ifdef IEEE_Arith
|
|---|
| 281 | if ((word0 (d) & Exp_mask) == Exp_mask)
|
|---|
| 282 | #else
|
|---|
| 283 | if (word0 (d) == 0x8000)
|
|---|
| 284 | #endif
|
|---|
| 285 | {
|
|---|
| 286 | /* Infinity or NaN */
|
|---|
| 287 | *decpt = 9999;
|
|---|
| 288 | s =
|
|---|
| 289 | #ifdef IEEE_Arith
|
|---|
| 290 | !word1 (d) && !(word0 (d) & 0xfffff) ? "Infinity" :
|
|---|
| 291 | #endif
|
|---|
| 292 | "NaN";
|
|---|
| 293 | if (rve)
|
|---|
| 294 | *rve =
|
|---|
| 295 | #ifdef IEEE_Arith
|
|---|
| 296 | s[3] ? s + 8 :
|
|---|
| 297 | #endif
|
|---|
| 298 | s + 3;
|
|---|
| 299 | return s;
|
|---|
| 300 | }
|
|---|
| 301 | #endif
|
|---|
| 302 | #ifdef IBM
|
|---|
| 303 | d.d += 0; /* normalize */
|
|---|
| 304 | #endif
|
|---|
| 305 | if (!d.d)
|
|---|
| 306 | {
|
|---|
| 307 | *decpt = 1;
|
|---|
| 308 | s = "0";
|
|---|
| 309 | if (rve)
|
|---|
| 310 | *rve = s + 1;
|
|---|
| 311 | return s;
|
|---|
| 312 | }
|
|---|
| 313 |
|
|---|
| 314 | b = d2b (ptr, d.d, &be, &bbits);
|
|---|
| 315 | #ifdef Sudden_Underflow
|
|---|
| 316 | i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1));
|
|---|
| 317 | #else
|
|---|
| 318 | if ((i = (int) (word0 (d) >> Exp_shift1 & (Exp_mask >> Exp_shift1))))
|
|---|
| 319 | {
|
|---|
| 320 | #endif
|
|---|
| 321 | d2.d = d.d;
|
|---|
| 322 | word0 (d2) &= Frac_mask1;
|
|---|
| 323 | word0 (d2) |= Exp_11;
|
|---|
| 324 | #ifdef IBM
|
|---|
| 325 | if (j = 11 - hi0bits (word0 (d2) & Frac_mask))
|
|---|
| 326 | d2.d /= 1 << j;
|
|---|
| 327 | #endif
|
|---|
| 328 |
|
|---|
| 329 | /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
|
|---|
| 330 | * log10(x) = log(x) / log(10)
|
|---|
| 331 | * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
|
|---|
| 332 | * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
|
|---|
| 333 | *
|
|---|
| 334 | * This suggests computing an approximation k to log10(d) by
|
|---|
| 335 | *
|
|---|
| 336 | * k = (i - Bias)*0.301029995663981
|
|---|
| 337 | * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
|
|---|
| 338 | *
|
|---|
| 339 | * We want k to be too large rather than too small.
|
|---|
| 340 | * The error in the first-order Taylor series approximation
|
|---|
| 341 | * is in our favor, so we just round up the constant enough
|
|---|
| 342 | * to compensate for any error in the multiplication of
|
|---|
| 343 | * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
|
|---|
| 344 | * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
|
|---|
| 345 | * adding 1e-13 to the constant term more than suffices.
|
|---|
| 346 | * Hence we adjust the constant term to 0.1760912590558.
|
|---|
| 347 | * (We could get a more accurate k by invoking log10,
|
|---|
| 348 | * but this is probably not worthwhile.)
|
|---|
| 349 | */
|
|---|
| 350 |
|
|---|
| 351 | i -= Bias;
|
|---|
| 352 | #ifdef IBM
|
|---|
| 353 | i <<= 2;
|
|---|
| 354 | i += j;
|
|---|
| 355 | #endif
|
|---|
| 356 | #ifndef Sudden_Underflow
|
|---|
| 357 | denorm = 0;
|
|---|
| 358 | }
|
|---|
| 359 | else
|
|---|
| 360 | {
|
|---|
| 361 | /* d is denormalized */
|
|---|
| 362 |
|
|---|
| 363 | i = bbits + be + (Bias + (P - 1) - 1);
|
|---|
| 364 | x = i > 32 ? word0 (d) << (64 - i) | word1 (d) >> (i - 32)
|
|---|
| 365 | : word1 (d) << (32 - i);
|
|---|
| 366 | d2.d = x;
|
|---|
| 367 | word0 (d2) -= 31 * Exp_msk1; /* adjust exponent */
|
|---|
| 368 | i -= (Bias + (P - 1) - 1) + 1;
|
|---|
| 369 | denorm = 1;
|
|---|
| 370 | }
|
|---|
| 371 | #endif
|
|---|
| 372 | ds = (d2.d - 1.5) * 0.289529654602168 + 0.1760912590558 + i * 0.301029995663981;
|
|---|
| 373 | k = (int) ds;
|
|---|
| 374 | if (ds < 0. && ds != k)
|
|---|
| 375 | k--; /* want k = floor(ds) */
|
|---|
| 376 | k_check = 1;
|
|---|
| 377 | if (k >= 0 && k <= Ten_pmax)
|
|---|
| 378 | {
|
|---|
| 379 | if (d.d < tens[k])
|
|---|
| 380 | k--;
|
|---|
| 381 | k_check = 0;
|
|---|
| 382 | }
|
|---|
| 383 | j = bbits - i - 1;
|
|---|
| 384 | if (j >= 0)
|
|---|
| 385 | {
|
|---|
| 386 | b2 = 0;
|
|---|
| 387 | s2 = j;
|
|---|
| 388 | }
|
|---|
| 389 | else
|
|---|
| 390 | {
|
|---|
| 391 | b2 = -j;
|
|---|
| 392 | s2 = 0;
|
|---|
| 393 | }
|
|---|
| 394 | if (k >= 0)
|
|---|
| 395 | {
|
|---|
| 396 | b5 = 0;
|
|---|
| 397 | s5 = k;
|
|---|
| 398 | s2 += k;
|
|---|
| 399 | }
|
|---|
| 400 | else
|
|---|
| 401 | {
|
|---|
| 402 | b2 -= k;
|
|---|
| 403 | b5 = -k;
|
|---|
| 404 | s5 = 0;
|
|---|
| 405 | }
|
|---|
| 406 | if (mode < 0 || mode > 9)
|
|---|
| 407 | mode = 0;
|
|---|
| 408 | try_quick = 1;
|
|---|
| 409 | if (mode > 5)
|
|---|
| 410 | {
|
|---|
| 411 | mode -= 4;
|
|---|
| 412 | try_quick = 0;
|
|---|
| 413 | }
|
|---|
| 414 | leftright = 1;
|
|---|
| 415 | switch (mode)
|
|---|
| 416 | {
|
|---|
| 417 | case 0:
|
|---|
| 418 | case 1:
|
|---|
| 419 | ilim = ilim1 = -1;
|
|---|
| 420 | i = 18;
|
|---|
| 421 | ndigits = 0;
|
|---|
| 422 | break;
|
|---|
| 423 | case 2:
|
|---|
| 424 | leftright = 0;
|
|---|
| 425 | /* no break */
|
|---|
| 426 | case 4:
|
|---|
| 427 | if (ndigits <= 0)
|
|---|
| 428 | ndigits = 1;
|
|---|
| 429 | ilim = ilim1 = i = ndigits;
|
|---|
| 430 | break;
|
|---|
| 431 | case 3:
|
|---|
| 432 | leftright = 0;
|
|---|
| 433 | /* no break */
|
|---|
| 434 | case 5:
|
|---|
| 435 | i = ndigits + k + 1;
|
|---|
| 436 | ilim = i;
|
|---|
| 437 | ilim1 = i - 1;
|
|---|
| 438 | if (i <= 0)
|
|---|
| 439 | i = 1;
|
|---|
| 440 | }
|
|---|
| 441 | j = sizeof (unsigned long);
|
|---|
| 442 | for (ptr->_result_k = 0; (int) (sizeof (_Jv_Bigint) - sizeof (unsigned long)) + j <= i;
|
|---|
| 443 | j <<= 1)
|
|---|
| 444 | ptr->_result_k++;
|
|---|
| 445 | ptr->_result = Balloc (ptr, ptr->_result_k);
|
|---|
| 446 | s = s0 = (char *) ptr->_result;
|
|---|
| 447 |
|
|---|
| 448 | if (ilim >= 0 && ilim <= Quick_max && try_quick)
|
|---|
| 449 | {
|
|---|
| 450 | /* Try to get by with floating-point arithmetic. */
|
|---|
| 451 |
|
|---|
| 452 | i = 0;
|
|---|
| 453 | d2.d = d.d;
|
|---|
| 454 | k0 = k;
|
|---|
| 455 | ilim0 = ilim;
|
|---|
| 456 | ieps = 2; /* conservative */
|
|---|
| 457 | if (k > 0)
|
|---|
| 458 | {
|
|---|
| 459 | ds = tens[k & 0xf];
|
|---|
| 460 | j = k >> 4;
|
|---|
| 461 | if (j & Bletch)
|
|---|
| 462 | {
|
|---|
| 463 | /* prevent overflows */
|
|---|
| 464 | j &= Bletch - 1;
|
|---|
| 465 | d.d /= bigtens[n_bigtens - 1];
|
|---|
| 466 | ieps++;
|
|---|
| 467 | }
|
|---|
| 468 | for (; j; j >>= 1, i++)
|
|---|
| 469 | if (j & 1)
|
|---|
| 470 | {
|
|---|
| 471 | ieps++;
|
|---|
| 472 | ds *= bigtens[i];
|
|---|
| 473 | }
|
|---|
| 474 | d.d /= ds;
|
|---|
| 475 | }
|
|---|
| 476 | else if ((j1 = -k))
|
|---|
| 477 | {
|
|---|
| 478 | d.d *= tens[j1 & 0xf];
|
|---|
| 479 | for (j = j1 >> 4; j; j >>= 1, i++)
|
|---|
| 480 | if (j & 1)
|
|---|
| 481 | {
|
|---|
| 482 | ieps++;
|
|---|
| 483 | d.d *= bigtens[i];
|
|---|
| 484 | }
|
|---|
| 485 | }
|
|---|
| 486 | if (k_check && d.d < 1. && ilim > 0)
|
|---|
| 487 | {
|
|---|
| 488 | if (ilim1 <= 0)
|
|---|
| 489 | goto fast_failed;
|
|---|
| 490 | ilim = ilim1;
|
|---|
| 491 | k--;
|
|---|
| 492 | d.d *= 10.;
|
|---|
| 493 | ieps++;
|
|---|
| 494 | }
|
|---|
| 495 | eps.d = ieps * d.d + 7.;
|
|---|
| 496 | word0 (eps) -= (P - 1) * Exp_msk1;
|
|---|
| 497 | if (ilim == 0)
|
|---|
| 498 | {
|
|---|
| 499 | S = mhi = 0;
|
|---|
| 500 | d.d -= 5.;
|
|---|
| 501 | if (d.d > eps.d)
|
|---|
| 502 | goto one_digit;
|
|---|
| 503 | if (d.d < -eps.d)
|
|---|
| 504 | goto no_digits;
|
|---|
| 505 | goto fast_failed;
|
|---|
| 506 | }
|
|---|
| 507 | #ifndef No_leftright
|
|---|
| 508 | if (leftright)
|
|---|
| 509 | {
|
|---|
| 510 | /* Use Steele & White method of only
|
|---|
| 511 | * generating digits needed.
|
|---|
| 512 | */
|
|---|
| 513 | eps.d = 0.5 / tens[ilim - 1] - eps.d;
|
|---|
| 514 | for (i = 0;;)
|
|---|
| 515 | {
|
|---|
| 516 | L = d.d;
|
|---|
| 517 | d.d -= L;
|
|---|
| 518 | *s++ = '0' + (int) L;
|
|---|
| 519 | if (d.d < eps.d)
|
|---|
| 520 | goto ret1;
|
|---|
| 521 | if (1. - d.d < eps.d)
|
|---|
| 522 | goto bump_up;
|
|---|
| 523 | if (++i >= ilim)
|
|---|
| 524 | break;
|
|---|
| 525 | eps.d *= 10.;
|
|---|
| 526 | d.d *= 10.;
|
|---|
| 527 | }
|
|---|
| 528 | }
|
|---|
| 529 | else
|
|---|
| 530 | {
|
|---|
| 531 | #endif
|
|---|
| 532 | /* Generate ilim digits, then fix them up. */
|
|---|
| 533 | eps.d *= tens[ilim - 1];
|
|---|
| 534 | for (i = 1;; i++, d.d *= 10.)
|
|---|
| 535 | {
|
|---|
| 536 | L = d.d;
|
|---|
| 537 | d.d -= L;
|
|---|
| 538 | *s++ = '0' + (int) L;
|
|---|
| 539 | if (i == ilim)
|
|---|
| 540 | {
|
|---|
| 541 | if (d.d > 0.5 + eps.d)
|
|---|
| 542 | goto bump_up;
|
|---|
| 543 | else if (d.d < 0.5 - eps.d)
|
|---|
| 544 | {
|
|---|
| 545 | while (*--s == '0');
|
|---|
| 546 | s++;
|
|---|
| 547 | goto ret1;
|
|---|
| 548 | }
|
|---|
| 549 | break;
|
|---|
| 550 | }
|
|---|
| 551 | }
|
|---|
| 552 | #ifndef No_leftright
|
|---|
| 553 | }
|
|---|
| 554 | #endif
|
|---|
| 555 | fast_failed:
|
|---|
| 556 | s = s0;
|
|---|
| 557 | d.d = d2.d;
|
|---|
| 558 | k = k0;
|
|---|
| 559 | ilim = ilim0;
|
|---|
| 560 | }
|
|---|
| 561 |
|
|---|
| 562 | /* Do we have a "small" integer? */
|
|---|
| 563 |
|
|---|
| 564 | if (be >= 0 && k <= Int_max)
|
|---|
| 565 | {
|
|---|
| 566 | /* Yes. */
|
|---|
| 567 | ds = tens[k];
|
|---|
| 568 | if (ndigits < 0 && ilim <= 0)
|
|---|
| 569 | {
|
|---|
| 570 | S = mhi = 0;
|
|---|
| 571 | if (ilim < 0 || d.d <= 5 * ds)
|
|---|
| 572 | goto no_digits;
|
|---|
| 573 | goto one_digit;
|
|---|
| 574 | }
|
|---|
| 575 | for (i = 1;; i++)
|
|---|
| 576 | {
|
|---|
| 577 | L = d.d / ds;
|
|---|
| 578 | d.d -= L * ds;
|
|---|
| 579 | #ifdef Check_FLT_ROUNDS
|
|---|
| 580 | /* If FLT_ROUNDS == 2, L will usually be high by 1 */
|
|---|
| 581 | if (d.d < 0)
|
|---|
| 582 | {
|
|---|
| 583 | L--;
|
|---|
| 584 | d.d += ds;
|
|---|
| 585 | }
|
|---|
| 586 | #endif
|
|---|
| 587 | *s++ = '0' + (int) L;
|
|---|
| 588 | if (i == ilim)
|
|---|
| 589 | {
|
|---|
| 590 | d.d += d.d;
|
|---|
| 591 | if (d.d > ds || (d.d == ds && L & 1))
|
|---|
| 592 | {
|
|---|
| 593 | bump_up:
|
|---|
| 594 | while (*--s == '9')
|
|---|
| 595 | if (s == s0)
|
|---|
| 596 | {
|
|---|
| 597 | k++;
|
|---|
| 598 | *s = '0';
|
|---|
| 599 | break;
|
|---|
| 600 | }
|
|---|
| 601 | ++*s++;
|
|---|
| 602 | }
|
|---|
| 603 | break;
|
|---|
| 604 | }
|
|---|
| 605 | if (!(d.d *= 10.))
|
|---|
| 606 | break;
|
|---|
| 607 | }
|
|---|
| 608 | goto ret1;
|
|---|
| 609 | }
|
|---|
| 610 |
|
|---|
| 611 | m2 = b2;
|
|---|
| 612 | m5 = b5;
|
|---|
| 613 | mhi = mlo = 0;
|
|---|
| 614 | if (leftright)
|
|---|
| 615 | {
|
|---|
| 616 | if (mode < 2)
|
|---|
| 617 | {
|
|---|
| 618 | i =
|
|---|
| 619 | #ifndef Sudden_Underflow
|
|---|
| 620 | denorm ? be + (Bias + (P - 1) - 1 + 1) :
|
|---|
| 621 | #endif
|
|---|
| 622 | #ifdef IBM
|
|---|
| 623 | 1 + 4 * P - 3 - bbits + ((bbits + be - 1) & 3);
|
|---|
| 624 | #else
|
|---|
| 625 | 1 + P - bbits;
|
|---|
| 626 | #endif
|
|---|
| 627 | }
|
|---|
| 628 | else
|
|---|
| 629 | {
|
|---|
| 630 | j = ilim - 1;
|
|---|
| 631 | if (m5 >= j)
|
|---|
| 632 | m5 -= j;
|
|---|
| 633 | else
|
|---|
| 634 | {
|
|---|
| 635 | s5 += j -= m5;
|
|---|
| 636 | b5 += j;
|
|---|
| 637 | m5 = 0;
|
|---|
| 638 | }
|
|---|
| 639 | if ((i = ilim) < 0)
|
|---|
| 640 | {
|
|---|
| 641 | m2 -= i;
|
|---|
| 642 | i = 0;
|
|---|
| 643 | }
|
|---|
| 644 | }
|
|---|
| 645 | b2 += i;
|
|---|
| 646 | s2 += i;
|
|---|
| 647 | mhi = i2b (ptr, 1);
|
|---|
| 648 | }
|
|---|
| 649 | if (m2 > 0 && s2 > 0)
|
|---|
| 650 | {
|
|---|
| 651 | i = m2 < s2 ? m2 : s2;
|
|---|
| 652 | b2 -= i;
|
|---|
| 653 | m2 -= i;
|
|---|
| 654 | s2 -= i;
|
|---|
| 655 | }
|
|---|
| 656 | if (b5 > 0)
|
|---|
| 657 | {
|
|---|
| 658 | if (leftright)
|
|---|
| 659 | {
|
|---|
| 660 | if (m5 > 0)
|
|---|
| 661 | {
|
|---|
| 662 | mhi = pow5mult (ptr, mhi, m5);
|
|---|
| 663 | b1 = mult (ptr, mhi, b);
|
|---|
| 664 | Bfree (ptr, b);
|
|---|
| 665 | b = b1;
|
|---|
| 666 | }
|
|---|
| 667 | if ((j = b5 - m5))
|
|---|
| 668 | b = pow5mult (ptr, b, j);
|
|---|
| 669 | }
|
|---|
| 670 | else
|
|---|
| 671 | b = pow5mult (ptr, b, b5);
|
|---|
| 672 | }
|
|---|
| 673 | S = i2b (ptr, 1);
|
|---|
| 674 | if (s5 > 0)
|
|---|
| 675 | S = pow5mult (ptr, S, s5);
|
|---|
| 676 |
|
|---|
| 677 | /* Check for special case that d is a normalized power of 2. */
|
|---|
| 678 |
|
|---|
| 679 | if (mode < 2)
|
|---|
| 680 | {
|
|---|
| 681 | if (!word1 (d) && !(word0 (d) & Bndry_mask)
|
|---|
| 682 | #ifndef Sudden_Underflow
|
|---|
| 683 | && word0(d) & Exp_mask
|
|---|
| 684 | #endif
|
|---|
| 685 | )
|
|---|
| 686 | {
|
|---|
| 687 | /* The special case */
|
|---|
| 688 | b2 += Log2P;
|
|---|
| 689 | s2 += Log2P;
|
|---|
| 690 | spec_case = 1;
|
|---|
| 691 | }
|
|---|
| 692 | else
|
|---|
| 693 | spec_case = 0;
|
|---|
| 694 | }
|
|---|
| 695 |
|
|---|
| 696 | /* Arrange for convenient computation of quotients:
|
|---|
| 697 | * shift left if necessary so divisor has 4 leading 0 bits.
|
|---|
| 698 | *
|
|---|
| 699 | * Perhaps we should just compute leading 28 bits of S once
|
|---|
| 700 | * and for all and pass them and a shift to quorem, so it
|
|---|
| 701 | * can do shifts and ors to compute the numerator for q.
|
|---|
| 702 | */
|
|---|
| 703 |
|
|---|
| 704 | #ifdef Pack_32
|
|---|
| 705 | if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0x1f))
|
|---|
| 706 | i = 32 - i;
|
|---|
| 707 | #else
|
|---|
| 708 | if ((i = ((s5 ? 32 - hi0bits (S->_x[S->_wds - 1]) : 1) + s2) & 0xf))
|
|---|
| 709 | i = 16 - i;
|
|---|
| 710 | #endif
|
|---|
| 711 | if (i > 4)
|
|---|
| 712 | {
|
|---|
| 713 | i -= 4;
|
|---|
| 714 | b2 += i;
|
|---|
| 715 | m2 += i;
|
|---|
| 716 | s2 += i;
|
|---|
| 717 | }
|
|---|
| 718 | else if (i < 4)
|
|---|
| 719 | {
|
|---|
| 720 | i += 28;
|
|---|
| 721 | b2 += i;
|
|---|
| 722 | m2 += i;
|
|---|
| 723 | s2 += i;
|
|---|
| 724 | }
|
|---|
| 725 | if (b2 > 0)
|
|---|
| 726 | b = lshift (ptr, b, b2);
|
|---|
| 727 | if (s2 > 0)
|
|---|
| 728 | S = lshift (ptr, S, s2);
|
|---|
| 729 | if (k_check)
|
|---|
| 730 | {
|
|---|
| 731 | if (cmp (b, S) < 0)
|
|---|
| 732 | {
|
|---|
| 733 | k--;
|
|---|
| 734 | b = multadd (ptr, b, 10, 0); /* we botched the k estimate */
|
|---|
| 735 | if (leftright)
|
|---|
| 736 | mhi = multadd (ptr, mhi, 10, 0);
|
|---|
| 737 | ilim = ilim1;
|
|---|
| 738 | }
|
|---|
| 739 | }
|
|---|
| 740 | if (ilim <= 0 && mode > 2)
|
|---|
| 741 | {
|
|---|
| 742 | if (ilim < 0 || cmp (b, S = multadd (ptr, S, 5, 0)) <= 0)
|
|---|
| 743 | {
|
|---|
| 744 | /* no digits, fcvt style */
|
|---|
| 745 | no_digits:
|
|---|
| 746 | k = -1 - ndigits;
|
|---|
| 747 | goto ret;
|
|---|
| 748 | }
|
|---|
| 749 | one_digit:
|
|---|
| 750 | *s++ = '1';
|
|---|
| 751 | k++;
|
|---|
| 752 | goto ret;
|
|---|
| 753 | }
|
|---|
| 754 | if (leftright)
|
|---|
| 755 | {
|
|---|
| 756 | if (m2 > 0)
|
|---|
| 757 | mhi = lshift (ptr, mhi, m2);
|
|---|
| 758 |
|
|---|
| 759 | /* Single precision case, */
|
|---|
| 760 | if (float_type)
|
|---|
| 761 | mhi = lshift (ptr, mhi, 29);
|
|---|
| 762 |
|
|---|
| 763 | /* Compute mlo -- check for special case
|
|---|
| 764 | * that d is a normalized power of 2.
|
|---|
| 765 | */
|
|---|
| 766 |
|
|---|
| 767 | mlo = mhi;
|
|---|
| 768 | if (spec_case)
|
|---|
| 769 | {
|
|---|
| 770 | mhi = Balloc (ptr, mhi->_k);
|
|---|
| 771 | Bcopy (mhi, mlo);
|
|---|
| 772 | mhi = lshift (ptr, mhi, Log2P);
|
|---|
| 773 | }
|
|---|
| 774 |
|
|---|
| 775 | for (i = 1;; i++)
|
|---|
| 776 | {
|
|---|
| 777 | dig = quorem (b, S) + '0';
|
|---|
| 778 | /* Do we yet have the shortest decimal string
|
|---|
| 779 | * that will round to d?
|
|---|
| 780 | */
|
|---|
| 781 | j = cmp (b, mlo);
|
|---|
| 782 | delta = diff (ptr, S, mhi);
|
|---|
| 783 | j1 = delta->_sign ? 1 : cmp (b, delta);
|
|---|
| 784 | Bfree (ptr, delta);
|
|---|
| 785 | #ifndef ROUND_BIASED
|
|---|
| 786 | if (j1 == 0 && !mode && !(word1 (d) & 1))
|
|---|
| 787 | {
|
|---|
| 788 | if (dig == '9')
|
|---|
| 789 | goto round_9_up;
|
|---|
| 790 | if (j > 0)
|
|---|
| 791 | dig++;
|
|---|
| 792 | *s++ = dig;
|
|---|
| 793 | goto ret;
|
|---|
| 794 | }
|
|---|
| 795 | #endif
|
|---|
| 796 | if (j < 0 || (j == 0 && !mode
|
|---|
| 797 | #ifndef ROUND_BIASED
|
|---|
| 798 | && !(word1 (d) & 1)
|
|---|
| 799 | #endif
|
|---|
| 800 | ))
|
|---|
| 801 | {
|
|---|
| 802 | if (j1 > 0)
|
|---|
| 803 | {
|
|---|
| 804 | b = lshift (ptr, b, 1);
|
|---|
| 805 | j1 = cmp (b, S);
|
|---|
| 806 | if ((j1 > 0 || (j1 == 0 && dig & 1))
|
|---|
| 807 | && dig++ == '9')
|
|---|
| 808 | goto round_9_up;
|
|---|
| 809 | }
|
|---|
| 810 | *s++ = dig;
|
|---|
| 811 | goto ret;
|
|---|
| 812 | }
|
|---|
| 813 | if (j1 > 0)
|
|---|
| 814 | {
|
|---|
| 815 | if (dig == '9')
|
|---|
| 816 | { /* possible if i == 1 */
|
|---|
| 817 | round_9_up:
|
|---|
| 818 | *s++ = '9';
|
|---|
| 819 | goto roundoff;
|
|---|
| 820 | }
|
|---|
| 821 | *s++ = dig + 1;
|
|---|
| 822 | goto ret;
|
|---|
| 823 | }
|
|---|
| 824 | *s++ = dig;
|
|---|
| 825 | if (i == ilim)
|
|---|
| 826 | break;
|
|---|
| 827 | b = multadd (ptr, b, 10, 0);
|
|---|
| 828 | if (mlo == mhi)
|
|---|
| 829 | mlo = mhi = multadd (ptr, mhi, 10, 0);
|
|---|
| 830 | else
|
|---|
| 831 | {
|
|---|
| 832 | mlo = multadd (ptr, mlo, 10, 0);
|
|---|
| 833 | mhi = multadd (ptr, mhi, 10, 0);
|
|---|
| 834 | }
|
|---|
| 835 | }
|
|---|
| 836 | }
|
|---|
| 837 | else
|
|---|
| 838 | for (i = 1;; i++)
|
|---|
| 839 | {
|
|---|
| 840 | *s++ = dig = quorem (b, S) + '0';
|
|---|
| 841 | if (i >= ilim)
|
|---|
| 842 | break;
|
|---|
| 843 | b = multadd (ptr, b, 10, 0);
|
|---|
| 844 | }
|
|---|
| 845 |
|
|---|
| 846 | /* Round off last digit */
|
|---|
| 847 |
|
|---|
| 848 | b = lshift (ptr, b, 1);
|
|---|
| 849 | j = cmp (b, S);
|
|---|
| 850 | if (j > 0 || (j == 0 && dig & 1))
|
|---|
| 851 | {
|
|---|
| 852 | roundoff:
|
|---|
| 853 | while (*--s == '9')
|
|---|
| 854 | if (s == s0)
|
|---|
| 855 | {
|
|---|
| 856 | k++;
|
|---|
| 857 | *s++ = '1';
|
|---|
| 858 | goto ret;
|
|---|
| 859 | }
|
|---|
| 860 | ++*s++;
|
|---|
| 861 | }
|
|---|
| 862 | else
|
|---|
| 863 | {
|
|---|
| 864 | while (*--s == '0');
|
|---|
| 865 | s++;
|
|---|
| 866 | }
|
|---|
| 867 | ret:
|
|---|
| 868 | Bfree (ptr, S);
|
|---|
| 869 | if (mhi)
|
|---|
| 870 | {
|
|---|
| 871 | if (mlo && mlo != mhi)
|
|---|
| 872 | Bfree (ptr, mlo);
|
|---|
| 873 | Bfree (ptr, mhi);
|
|---|
| 874 | }
|
|---|
| 875 | ret1:
|
|---|
| 876 | Bfree (ptr, b);
|
|---|
| 877 | *s = 0;
|
|---|
| 878 | *decpt = k + 1;
|
|---|
| 879 | if (rve)
|
|---|
| 880 | *rve = s;
|
|---|
| 881 | return s0;
|
|---|
| 882 | }
|
|---|
| 883 |
|
|---|
| 884 |
|
|---|
| 885 | _VOID
|
|---|
| 886 | _DEFUN (_dtoa,
|
|---|
| 887 | (_d, mode, ndigits, decpt, sign, rve, buf, float_type),
|
|---|
| 888 | double _d _AND
|
|---|
| 889 | int mode _AND
|
|---|
| 890 | int ndigits _AND
|
|---|
| 891 | int *decpt _AND
|
|---|
| 892 | int *sign _AND
|
|---|
| 893 | char **rve _AND
|
|---|
| 894 | char *buf _AND
|
|---|
| 895 | int float_type)
|
|---|
| 896 | {
|
|---|
| 897 | struct _Jv_reent reent;
|
|---|
| 898 | char *p;
|
|---|
| 899 | memset (&reent, 0, sizeof reent);
|
|---|
| 900 |
|
|---|
| 901 | p = _dtoa_r (&reent, _d, mode, ndigits, decpt, sign, rve, float_type);
|
|---|
| 902 | strcpy (buf, p);
|
|---|
| 903 |
|
|---|
| 904 | return;
|
|---|
| 905 | }
|
|---|
