# Python code coverage for Python/dtoa.c

# | count | content |
---|---|---|

1 | n/a | /**************************************************************** |

2 | n/a | * |

3 | n/a | * The author of this software is David M. Gay. |

4 | n/a | * |

5 | n/a | * Copyright (c) 1991, 2000, 2001 by Lucent Technologies. |

6 | n/a | * |

7 | n/a | * Permission to use, copy, modify, and distribute this software for any |

8 | n/a | * purpose without fee is hereby granted, provided that this entire notice |

9 | n/a | * is included in all copies of any software which is or includes a copy |

10 | n/a | * or modification of this software and in all copies of the supporting |

11 | n/a | * documentation for such software. |

12 | n/a | * |

13 | n/a | * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED |

14 | n/a | * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY |

15 | n/a | * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY |

16 | n/a | * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. |

17 | n/a | * |

18 | n/a | ***************************************************************/ |

19 | n/a | |

20 | n/a | /**************************************************************** |

21 | n/a | * This is dtoa.c by David M. Gay, downloaded from |

22 | n/a | * http://www.netlib.org/fp/dtoa.c on April 15, 2009 and modified for |

23 | n/a | * inclusion into the Python core by Mark E. T. Dickinson and Eric V. Smith. |

24 | n/a | * |

25 | n/a | * Please remember to check http://www.netlib.org/fp regularly (and especially |

26 | n/a | * before any Python release) for bugfixes and updates. |

27 | n/a | * |

28 | n/a | * The major modifications from Gay's original code are as follows: |

29 | n/a | * |

30 | n/a | * 0. The original code has been specialized to Python's needs by removing |

31 | n/a | * many of the #ifdef'd sections. In particular, code to support VAX and |

32 | n/a | * IBM floating-point formats, hex NaNs, hex floats, locale-aware |

33 | n/a | * treatment of the decimal point, and setting of the inexact flag have |

34 | n/a | * been removed. |

35 | n/a | * |

36 | n/a | * 1. We use PyMem_Malloc and PyMem_Free in place of malloc and free. |

37 | n/a | * |

38 | n/a | * 2. The public functions strtod, dtoa and freedtoa all now have |

39 | n/a | * a _Py_dg_ prefix. |

40 | n/a | * |

41 | n/a | * 3. Instead of assuming that PyMem_Malloc always succeeds, we thread |

42 | n/a | * PyMem_Malloc failures through the code. The functions |

43 | n/a | * |

44 | n/a | * Balloc, multadd, s2b, i2b, mult, pow5mult, lshift, diff, d2b |

45 | n/a | * |

46 | n/a | * of return type *Bigint all return NULL to indicate a malloc failure. |

47 | n/a | * Similarly, rv_alloc and nrv_alloc (return type char *) return NULL on |

48 | n/a | * failure. bigcomp now has return type int (it used to be void) and |

49 | n/a | * returns -1 on failure and 0 otherwise. _Py_dg_dtoa returns NULL |

50 | n/a | * on failure. _Py_dg_strtod indicates failure due to malloc failure |

51 | n/a | * by returning -1.0, setting errno=ENOMEM and *se to s00. |

52 | n/a | * |

53 | n/a | * 4. The static variable dtoa_result has been removed. Callers of |

54 | n/a | * _Py_dg_dtoa are expected to call _Py_dg_freedtoa to free |

55 | n/a | * the memory allocated by _Py_dg_dtoa. |

56 | n/a | * |

57 | n/a | * 5. The code has been reformatted to better fit with Python's |

58 | n/a | * C style guide (PEP 7). |

59 | n/a | * |

60 | n/a | * 6. A bug in the memory allocation has been fixed: to avoid FREEing memory |

61 | n/a | * that hasn't been MALLOC'ed, private_mem should only be used when k <= |

62 | n/a | * Kmax. |

63 | n/a | * |

64 | n/a | * 7. _Py_dg_strtod has been modified so that it doesn't accept strings with |

65 | n/a | * leading whitespace. |

66 | n/a | * |

67 | n/a | ***************************************************************/ |

68 | n/a | |

69 | n/a | /* Please send bug reports for the original dtoa.c code to David M. Gay (dmg |

70 | n/a | * at acm dot org, with " at " changed at "@" and " dot " changed to "."). |

71 | n/a | * Please report bugs for this modified version using the Python issue tracker |

72 | n/a | * (http://bugs.python.org). */ |

73 | n/a | |

74 | n/a | /* On a machine with IEEE extended-precision registers, it is |

75 | n/a | * necessary to specify double-precision (53-bit) rounding precision |

76 | n/a | * before invoking strtod or dtoa. If the machine uses (the equivalent |

77 | n/a | * of) Intel 80x87 arithmetic, the call |

78 | n/a | * _control87(PC_53, MCW_PC); |

79 | n/a | * does this with many compilers. Whether this or another call is |

80 | n/a | * appropriate depends on the compiler; for this to work, it may be |

81 | n/a | * necessary to #include "float.h" or another system-dependent header |

82 | n/a | * file. |

83 | n/a | */ |

84 | n/a | |

85 | n/a | /* strtod for IEEE-, VAX-, and IBM-arithmetic machines. |

86 | n/a | * |

87 | n/a | * This strtod returns a nearest machine number to the input decimal |

88 | n/a | * string (or sets errno to ERANGE). With IEEE arithmetic, ties are |

89 | n/a | * broken by the IEEE round-even rule. Otherwise ties are broken by |

90 | n/a | * biased rounding (add half and chop). |

91 | n/a | * |

92 | n/a | * Inspired loosely by William D. Clinger's paper "How to Read Floating |

93 | n/a | * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101]. |

94 | n/a | * |

95 | n/a | * Modifications: |

96 | n/a | * |

97 | n/a | * 1. We only require IEEE, IBM, or VAX double-precision |

98 | n/a | * arithmetic (not IEEE double-extended). |

99 | n/a | * 2. We get by with floating-point arithmetic in a case that |

100 | n/a | * Clinger missed -- when we're computing d * 10^n |

101 | n/a | * for a small integer d and the integer n is not too |

102 | n/a | * much larger than 22 (the maximum integer k for which |

103 | n/a | * we can represent 10^k exactly), we may be able to |

104 | n/a | * compute (d*10^k) * 10^(e-k) with just one roundoff. |

105 | n/a | * 3. Rather than a bit-at-a-time adjustment of the binary |

106 | n/a | * result in the hard case, we use floating-point |

107 | n/a | * arithmetic to determine the adjustment to within |

108 | n/a | * one bit; only in really hard cases do we need to |

109 | n/a | * compute a second residual. |

110 | n/a | * 4. Because of 3., we don't need a large table of powers of 10 |

111 | n/a | * for ten-to-e (just some small tables, e.g. of 10^k |

112 | n/a | * for 0 <= k <= 22). |

113 | n/a | */ |

114 | n/a | |

115 | n/a | /* Linking of Python's #defines to Gay's #defines starts here. */ |

116 | n/a | |

117 | n/a | #include "Python.h" |

118 | n/a | |

119 | n/a | /* if PY_NO_SHORT_FLOAT_REPR is defined, then don't even try to compile |

120 | n/a | the following code */ |

121 | n/a | #ifndef PY_NO_SHORT_FLOAT_REPR |

122 | n/a | |

123 | n/a | #include "float.h" |

124 | n/a | |

125 | n/a | #define MALLOC PyMem_Malloc |

126 | n/a | #define FREE PyMem_Free |

127 | n/a | |

128 | n/a | /* This code should also work for ARM mixed-endian format on little-endian |

129 | n/a | machines, where doubles have byte order 45670123 (in increasing address |

130 | n/a | order, 0 being the least significant byte). */ |

131 | n/a | #ifdef DOUBLE_IS_LITTLE_ENDIAN_IEEE754 |

132 | n/a | # define IEEE_8087 |

133 | n/a | #endif |

134 | n/a | #if defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) || \ |

135 | n/a | defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754) |

136 | n/a | # define IEEE_MC68k |

137 | n/a | #endif |

138 | n/a | #if defined(IEEE_8087) + defined(IEEE_MC68k) != 1 |

139 | n/a | #error "Exactly one of IEEE_8087 or IEEE_MC68k should be defined." |

140 | n/a | #endif |

141 | n/a | |

142 | n/a | /* The code below assumes that the endianness of integers matches the |

143 | n/a | endianness of the two 32-bit words of a double. Check this. */ |

144 | n/a | #if defined(WORDS_BIGENDIAN) && (defined(DOUBLE_IS_LITTLE_ENDIAN_IEEE754) || \ |

145 | n/a | defined(DOUBLE_IS_ARM_MIXED_ENDIAN_IEEE754)) |

146 | n/a | #error "doubles and ints have incompatible endianness" |

147 | n/a | #endif |

148 | n/a | |

149 | n/a | #if !defined(WORDS_BIGENDIAN) && defined(DOUBLE_IS_BIG_ENDIAN_IEEE754) |

150 | n/a | #error "doubles and ints have incompatible endianness" |

151 | n/a | #endif |

152 | n/a | |

153 | n/a | |

154 | n/a | typedef uint32_t ULong; |

155 | n/a | typedef int32_t Long; |

156 | n/a | typedef uint64_t ULLong; |

157 | n/a | |

158 | n/a | #undef DEBUG |

159 | n/a | #ifdef Py_DEBUG |

160 | n/a | #define DEBUG |

161 | n/a | #endif |

162 | n/a | |

163 | n/a | /* End Python #define linking */ |

164 | n/a | |

165 | n/a | #ifdef DEBUG |

166 | n/a | #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);} |

167 | n/a | #endif |

168 | n/a | |

169 | n/a | #ifndef PRIVATE_MEM |

170 | n/a | #define PRIVATE_MEM 2304 |

171 | n/a | #endif |

172 | n/a | #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) |

173 | n/a | static double private_mem[PRIVATE_mem], *pmem_next = private_mem; |

174 | n/a | |

175 | n/a | #ifdef __cplusplus |

176 | n/a | extern "C" { |

177 | n/a | #endif |

178 | n/a | |

179 | n/a | typedef union { double d; ULong L[2]; } U; |

180 | n/a | |

181 | n/a | #ifdef IEEE_8087 |

182 | n/a | #define word0(x) (x)->L[1] |

183 | n/a | #define word1(x) (x)->L[0] |

184 | n/a | #else |

185 | n/a | #define word0(x) (x)->L[0] |

186 | n/a | #define word1(x) (x)->L[1] |

187 | n/a | #endif |

188 | n/a | #define dval(x) (x)->d |

189 | n/a | |

190 | n/a | #ifndef STRTOD_DIGLIM |

191 | n/a | #define STRTOD_DIGLIM 40 |

192 | n/a | #endif |

193 | n/a | |

194 | n/a | /* maximum permitted exponent value for strtod; exponents larger than |

195 | n/a | MAX_ABS_EXP in absolute value get truncated to +-MAX_ABS_EXP. MAX_ABS_EXP |

196 | n/a | should fit into an int. */ |

197 | n/a | #ifndef MAX_ABS_EXP |

198 | n/a | #define MAX_ABS_EXP 1100000000U |

199 | n/a | #endif |

200 | n/a | /* Bound on length of pieces of input strings in _Py_dg_strtod; specifically, |

201 | n/a | this is used to bound the total number of digits ignoring leading zeros and |

202 | n/a | the number of digits that follow the decimal point. Ideally, MAX_DIGITS |

203 | n/a | should satisfy MAX_DIGITS + 400 < MAX_ABS_EXP; that ensures that the |

204 | n/a | exponent clipping in _Py_dg_strtod can't affect the value of the output. */ |

205 | n/a | #ifndef MAX_DIGITS |

206 | n/a | #define MAX_DIGITS 1000000000U |

207 | n/a | #endif |

208 | n/a | |

209 | n/a | /* Guard against trying to use the above values on unusual platforms with ints |

210 | n/a | * of width less than 32 bits. */ |

211 | n/a | #if MAX_ABS_EXP > INT_MAX |

212 | n/a | #error "MAX_ABS_EXP should fit in an int" |

213 | n/a | #endif |

214 | n/a | #if MAX_DIGITS > INT_MAX |

215 | n/a | #error "MAX_DIGITS should fit in an int" |

216 | n/a | #endif |

217 | n/a | |

218 | n/a | /* The following definition of Storeinc is appropriate for MIPS processors. |

219 | n/a | * An alternative that might be better on some machines is |

220 | n/a | * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff) |

221 | n/a | */ |

222 | n/a | #if defined(IEEE_8087) |

223 | n/a | #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \ |

224 | n/a | ((unsigned short *)a)[0] = (unsigned short)c, a++) |

225 | n/a | #else |

226 | n/a | #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \ |

227 | n/a | ((unsigned short *)a)[1] = (unsigned short)c, a++) |

228 | n/a | #endif |

229 | n/a | |

230 | n/a | /* #define P DBL_MANT_DIG */ |

231 | n/a | /* Ten_pmax = floor(P*log(2)/log(5)) */ |

232 | n/a | /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */ |

233 | n/a | /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */ |

234 | n/a | /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */ |

235 | n/a | |

236 | n/a | #define Exp_shift 20 |

237 | n/a | #define Exp_shift1 20 |

238 | n/a | #define Exp_msk1 0x100000 |

239 | n/a | #define Exp_msk11 0x100000 |

240 | n/a | #define Exp_mask 0x7ff00000 |

241 | n/a | #define P 53 |

242 | n/a | #define Nbits 53 |

243 | n/a | #define Bias 1023 |

244 | n/a | #define Emax 1023 |

245 | n/a | #define Emin (-1022) |

246 | n/a | #define Etiny (-1074) /* smallest denormal is 2**Etiny */ |

247 | n/a | #define Exp_1 0x3ff00000 |

248 | n/a | #define Exp_11 0x3ff00000 |

249 | n/a | #define Ebits 11 |

250 | n/a | #define Frac_mask 0xfffff |

251 | n/a | #define Frac_mask1 0xfffff |

252 | n/a | #define Ten_pmax 22 |

253 | n/a | #define Bletch 0x10 |

254 | n/a | #define Bndry_mask 0xfffff |

255 | n/a | #define Bndry_mask1 0xfffff |

256 | n/a | #define Sign_bit 0x80000000 |

257 | n/a | #define Log2P 1 |

258 | n/a | #define Tiny0 0 |

259 | n/a | #define Tiny1 1 |

260 | n/a | #define Quick_max 14 |

261 | n/a | #define Int_max 14 |

262 | n/a | |

263 | n/a | #ifndef Flt_Rounds |

264 | n/a | #ifdef FLT_ROUNDS |

265 | n/a | #define Flt_Rounds FLT_ROUNDS |

266 | n/a | #else |

267 | n/a | #define Flt_Rounds 1 |

268 | n/a | #endif |

269 | n/a | #endif /*Flt_Rounds*/ |

270 | n/a | |

271 | n/a | #define Rounding Flt_Rounds |

272 | n/a | |

273 | n/a | #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1)) |

274 | n/a | #define Big1 0xffffffff |

275 | n/a | |

276 | n/a | /* Standard NaN used by _Py_dg_stdnan. */ |

277 | n/a | |

278 | n/a | #define NAN_WORD0 0x7ff80000 |

279 | n/a | #define NAN_WORD1 0 |

280 | n/a | |

281 | n/a | /* Bits of the representation of positive infinity. */ |

282 | n/a | |

283 | n/a | #define POSINF_WORD0 0x7ff00000 |

284 | n/a | #define POSINF_WORD1 0 |

285 | n/a | |

286 | n/a | /* struct BCinfo is used to pass information from _Py_dg_strtod to bigcomp */ |

287 | n/a | |

288 | n/a | typedef struct BCinfo BCinfo; |

289 | n/a | struct |

290 | n/a | BCinfo { |

291 | n/a | int e0, nd, nd0, scale; |

292 | n/a | }; |

293 | n/a | |

294 | n/a | #define FFFFFFFF 0xffffffffUL |

295 | n/a | |

296 | n/a | #define Kmax 7 |

297 | n/a | |

298 | n/a | /* struct Bigint is used to represent arbitrary-precision integers. These |

299 | n/a | integers are stored in sign-magnitude format, with the magnitude stored as |

300 | n/a | an array of base 2**32 digits. Bigints are always normalized: if x is a |

301 | n/a | Bigint then x->wds >= 1, and either x->wds == 1 or x[wds-1] is nonzero. |

302 | n/a | |

303 | n/a | The Bigint fields are as follows: |

304 | n/a | |

305 | n/a | - next is a header used by Balloc and Bfree to keep track of lists |

306 | n/a | of freed Bigints; it's also used for the linked list of |

307 | n/a | powers of 5 of the form 5**2**i used by pow5mult. |

308 | n/a | - k indicates which pool this Bigint was allocated from |

309 | n/a | - maxwds is the maximum number of words space was allocated for |

310 | n/a | (usually maxwds == 2**k) |

311 | n/a | - sign is 1 for negative Bigints, 0 for positive. The sign is unused |

312 | n/a | (ignored on inputs, set to 0 on outputs) in almost all operations |

313 | n/a | involving Bigints: a notable exception is the diff function, which |

314 | n/a | ignores signs on inputs but sets the sign of the output correctly. |

315 | n/a | - wds is the actual number of significant words |

316 | n/a | - x contains the vector of words (digits) for this Bigint, from least |

317 | n/a | significant (x[0]) to most significant (x[wds-1]). |

318 | n/a | */ |

319 | n/a | |

320 | n/a | struct |

321 | n/a | Bigint { |

322 | n/a | struct Bigint *next; |

323 | n/a | int k, maxwds, sign, wds; |

324 | n/a | ULong x[1]; |

325 | n/a | }; |

326 | n/a | |

327 | n/a | typedef struct Bigint Bigint; |

328 | n/a | |

329 | n/a | #ifndef Py_USING_MEMORY_DEBUGGER |

330 | n/a | |

331 | n/a | /* Memory management: memory is allocated from, and returned to, Kmax+1 pools |

332 | n/a | of memory, where pool k (0 <= k <= Kmax) is for Bigints b with b->maxwds == |

333 | n/a | 1 << k. These pools are maintained as linked lists, with freelist[k] |

334 | n/a | pointing to the head of the list for pool k. |

335 | n/a | |

336 | n/a | On allocation, if there's no free slot in the appropriate pool, MALLOC is |

337 | n/a | called to get more memory. This memory is not returned to the system until |

338 | n/a | Python quits. There's also a private memory pool that's allocated from |

339 | n/a | in preference to using MALLOC. |

340 | n/a | |

341 | n/a | For Bigints with more than (1 << Kmax) digits (which implies at least 1233 |

342 | n/a | decimal digits), memory is directly allocated using MALLOC, and freed using |

343 | n/a | FREE. |

344 | n/a | |

345 | n/a | XXX: it would be easy to bypass this memory-management system and |

346 | n/a | translate each call to Balloc into a call to PyMem_Malloc, and each |

347 | n/a | Bfree to PyMem_Free. Investigate whether this has any significant |

348 | n/a | performance on impact. */ |

349 | n/a | |

350 | n/a | static Bigint *freelist[Kmax+1]; |

351 | n/a | |

352 | n/a | /* Allocate space for a Bigint with up to 1<<k digits */ |

353 | n/a | |

354 | n/a | static Bigint * |

355 | n/a | Balloc(int k) |

356 | n/a | { |

357 | n/a | int x; |

358 | n/a | Bigint *rv; |

359 | n/a | unsigned int len; |

360 | n/a | |

361 | n/a | if (k <= Kmax && (rv = freelist[k])) |

362 | n/a | freelist[k] = rv->next; |

363 | n/a | else { |

364 | n/a | x = 1 << k; |

365 | n/a | len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) |

366 | n/a | /sizeof(double); |

367 | n/a | if (k <= Kmax && pmem_next - private_mem + len <= (Py_ssize_t)PRIVATE_mem) { |

368 | n/a | rv = (Bigint*)pmem_next; |

369 | n/a | pmem_next += len; |

370 | n/a | } |

371 | n/a | else { |

372 | n/a | rv = (Bigint*)MALLOC(len*sizeof(double)); |

373 | n/a | if (rv == NULL) |

374 | n/a | return NULL; |

375 | n/a | } |

376 | n/a | rv->k = k; |

377 | n/a | rv->maxwds = x; |

378 | n/a | } |

379 | n/a | rv->sign = rv->wds = 0; |

380 | n/a | return rv; |

381 | n/a | } |

382 | n/a | |

383 | n/a | /* Free a Bigint allocated with Balloc */ |

384 | n/a | |

385 | n/a | static void |

386 | n/a | Bfree(Bigint *v) |

387 | n/a | { |

388 | n/a | if (v) { |

389 | n/a | if (v->k > Kmax) |

390 | n/a | FREE((void*)v); |

391 | n/a | else { |

392 | n/a | v->next = freelist[v->k]; |

393 | n/a | freelist[v->k] = v; |

394 | n/a | } |

395 | n/a | } |

396 | n/a | } |

397 | n/a | |

398 | n/a | #else |

399 | n/a | |

400 | n/a | /* Alternative versions of Balloc and Bfree that use PyMem_Malloc and |

401 | n/a | PyMem_Free directly in place of the custom memory allocation scheme above. |

402 | n/a | These are provided for the benefit of memory debugging tools like |

403 | n/a | Valgrind. */ |

404 | n/a | |

405 | n/a | /* Allocate space for a Bigint with up to 1<<k digits */ |

406 | n/a | |

407 | n/a | static Bigint * |

408 | n/a | Balloc(int k) |

409 | n/a | { |

410 | n/a | int x; |

411 | n/a | Bigint *rv; |

412 | n/a | unsigned int len; |

413 | n/a | |

414 | n/a | x = 1 << k; |

415 | n/a | len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) |

416 | n/a | /sizeof(double); |

417 | n/a | |

418 | n/a | rv = (Bigint*)MALLOC(len*sizeof(double)); |

419 | n/a | if (rv == NULL) |

420 | n/a | return NULL; |

421 | n/a | |

422 | n/a | rv->k = k; |

423 | n/a | rv->maxwds = x; |

424 | n/a | rv->sign = rv->wds = 0; |

425 | n/a | return rv; |

426 | n/a | } |

427 | n/a | |

428 | n/a | /* Free a Bigint allocated with Balloc */ |

429 | n/a | |

430 | n/a | static void |

431 | n/a | Bfree(Bigint *v) |

432 | n/a | { |

433 | n/a | if (v) { |

434 | n/a | FREE((void*)v); |

435 | n/a | } |

436 | n/a | } |

437 | n/a | |

438 | n/a | #endif /* Py_USING_MEMORY_DEBUGGER */ |

439 | n/a | |

440 | n/a | #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \ |

441 | n/a | y->wds*sizeof(Long) + 2*sizeof(int)) |

442 | n/a | |

443 | n/a | /* Multiply a Bigint b by m and add a. Either modifies b in place and returns |

444 | n/a | a pointer to the modified b, or Bfrees b and returns a pointer to a copy. |

445 | n/a | On failure, return NULL. In this case, b will have been already freed. */ |

446 | n/a | |

447 | n/a | static Bigint * |

448 | n/a | multadd(Bigint *b, int m, int a) /* multiply by m and add a */ |

449 | n/a | { |

450 | n/a | int i, wds; |

451 | n/a | ULong *x; |

452 | n/a | ULLong carry, y; |

453 | n/a | Bigint *b1; |

454 | n/a | |

455 | n/a | wds = b->wds; |

456 | n/a | x = b->x; |

457 | n/a | i = 0; |

458 | n/a | carry = a; |

459 | n/a | do { |

460 | n/a | y = *x * (ULLong)m + carry; |

461 | n/a | carry = y >> 32; |

462 | n/a | *x++ = (ULong)(y & FFFFFFFF); |

463 | n/a | } |

464 | n/a | while(++i < wds); |

465 | n/a | if (carry) { |

466 | n/a | if (wds >= b->maxwds) { |

467 | n/a | b1 = Balloc(b->k+1); |

468 | n/a | if (b1 == NULL){ |

469 | n/a | Bfree(b); |

470 | n/a | return NULL; |

471 | n/a | } |

472 | n/a | Bcopy(b1, b); |

473 | n/a | Bfree(b); |

474 | n/a | b = b1; |

475 | n/a | } |

476 | n/a | b->x[wds++] = (ULong)carry; |

477 | n/a | b->wds = wds; |

478 | n/a | } |

479 | n/a | return b; |

480 | n/a | } |

481 | n/a | |

482 | n/a | /* convert a string s containing nd decimal digits (possibly containing a |

483 | n/a | decimal separator at position nd0, which is ignored) to a Bigint. This |

484 | n/a | function carries on where the parsing code in _Py_dg_strtod leaves off: on |

485 | n/a | entry, y9 contains the result of converting the first 9 digits. Returns |

486 | n/a | NULL on failure. */ |

487 | n/a | |

488 | n/a | static Bigint * |

489 | n/a | s2b(const char *s, int nd0, int nd, ULong y9) |

490 | n/a | { |

491 | n/a | Bigint *b; |

492 | n/a | int i, k; |

493 | n/a | Long x, y; |

494 | n/a | |

495 | n/a | x = (nd + 8) / 9; |

496 | n/a | for(k = 0, y = 1; x > y; y <<= 1, k++) ; |

497 | n/a | b = Balloc(k); |

498 | n/a | if (b == NULL) |

499 | n/a | return NULL; |

500 | n/a | b->x[0] = y9; |

501 | n/a | b->wds = 1; |

502 | n/a | |

503 | n/a | if (nd <= 9) |

504 | n/a | return b; |

505 | n/a | |

506 | n/a | s += 9; |

507 | n/a | for (i = 9; i < nd0; i++) { |

508 | n/a | b = multadd(b, 10, *s++ - '0'); |

509 | n/a | if (b == NULL) |

510 | n/a | return NULL; |

511 | n/a | } |

512 | n/a | s++; |

513 | n/a | for(; i < nd; i++) { |

514 | n/a | b = multadd(b, 10, *s++ - '0'); |

515 | n/a | if (b == NULL) |

516 | n/a | return NULL; |

517 | n/a | } |

518 | n/a | return b; |

519 | n/a | } |

520 | n/a | |

521 | n/a | /* count leading 0 bits in the 32-bit integer x. */ |

522 | n/a | |

523 | n/a | static int |

524 | n/a | hi0bits(ULong x) |

525 | n/a | { |

526 | n/a | int k = 0; |

527 | n/a | |

528 | n/a | if (!(x & 0xffff0000)) { |

529 | n/a | k = 16; |

530 | n/a | x <<= 16; |

531 | n/a | } |

532 | n/a | if (!(x & 0xff000000)) { |

533 | n/a | k += 8; |

534 | n/a | x <<= 8; |

535 | n/a | } |

536 | n/a | if (!(x & 0xf0000000)) { |

537 | n/a | k += 4; |

538 | n/a | x <<= 4; |

539 | n/a | } |

540 | n/a | if (!(x & 0xc0000000)) { |

541 | n/a | k += 2; |

542 | n/a | x <<= 2; |

543 | n/a | } |

544 | n/a | if (!(x & 0x80000000)) { |

545 | n/a | k++; |

546 | n/a | if (!(x & 0x40000000)) |

547 | n/a | return 32; |

548 | n/a | } |

549 | n/a | return k; |

550 | n/a | } |

551 | n/a | |

552 | n/a | /* count trailing 0 bits in the 32-bit integer y, and shift y right by that |

553 | n/a | number of bits. */ |

554 | n/a | |

555 | n/a | static int |

556 | n/a | lo0bits(ULong *y) |

557 | n/a | { |

558 | n/a | int k; |

559 | n/a | ULong x = *y; |

560 | n/a | |

561 | n/a | if (x & 7) { |

562 | n/a | if (x & 1) |

563 | n/a | return 0; |

564 | n/a | if (x & 2) { |

565 | n/a | *y = x >> 1; |

566 | n/a | return 1; |

567 | n/a | } |

568 | n/a | *y = x >> 2; |

569 | n/a | return 2; |

570 | n/a | } |

571 | n/a | k = 0; |

572 | n/a | if (!(x & 0xffff)) { |

573 | n/a | k = 16; |

574 | n/a | x >>= 16; |

575 | n/a | } |

576 | n/a | if (!(x & 0xff)) { |

577 | n/a | k += 8; |

578 | n/a | x >>= 8; |

579 | n/a | } |

580 | n/a | if (!(x & 0xf)) { |

581 | n/a | k += 4; |

582 | n/a | x >>= 4; |

583 | n/a | } |

584 | n/a | if (!(x & 0x3)) { |

585 | n/a | k += 2; |

586 | n/a | x >>= 2; |

587 | n/a | } |

588 | n/a | if (!(x & 1)) { |

589 | n/a | k++; |

590 | n/a | x >>= 1; |

591 | n/a | if (!x) |

592 | n/a | return 32; |

593 | n/a | } |

594 | n/a | *y = x; |

595 | n/a | return k; |

596 | n/a | } |

597 | n/a | |

598 | n/a | /* convert a small nonnegative integer to a Bigint */ |

599 | n/a | |

600 | n/a | static Bigint * |

601 | n/a | i2b(int i) |

602 | n/a | { |

603 | n/a | Bigint *b; |

604 | n/a | |

605 | n/a | b = Balloc(1); |

606 | n/a | if (b == NULL) |

607 | n/a | return NULL; |

608 | n/a | b->x[0] = i; |

609 | n/a | b->wds = 1; |

610 | n/a | return b; |

611 | n/a | } |

612 | n/a | |

613 | n/a | /* multiply two Bigints. Returns a new Bigint, or NULL on failure. Ignores |

614 | n/a | the signs of a and b. */ |

615 | n/a | |

616 | n/a | static Bigint * |

617 | n/a | mult(Bigint *a, Bigint *b) |

618 | n/a | { |

619 | n/a | Bigint *c; |

620 | n/a | int k, wa, wb, wc; |

621 | n/a | ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; |

622 | n/a | ULong y; |

623 | n/a | ULLong carry, z; |

624 | n/a | |

625 | n/a | if ((!a->x[0] && a->wds == 1) || (!b->x[0] && b->wds == 1)) { |

626 | n/a | c = Balloc(0); |

627 | n/a | if (c == NULL) |

628 | n/a | return NULL; |

629 | n/a | c->wds = 1; |

630 | n/a | c->x[0] = 0; |

631 | n/a | return c; |

632 | n/a | } |

633 | n/a | |

634 | n/a | if (a->wds < b->wds) { |

635 | n/a | c = a; |

636 | n/a | a = b; |

637 | n/a | b = c; |

638 | n/a | } |

639 | n/a | k = a->k; |

640 | n/a | wa = a->wds; |

641 | n/a | wb = b->wds; |

642 | n/a | wc = wa + wb; |

643 | n/a | if (wc > a->maxwds) |

644 | n/a | k++; |

645 | n/a | c = Balloc(k); |

646 | n/a | if (c == NULL) |

647 | n/a | return NULL; |

648 | n/a | for(x = c->x, xa = x + wc; x < xa; x++) |

649 | n/a | *x = 0; |

650 | n/a | xa = a->x; |

651 | n/a | xae = xa + wa; |

652 | n/a | xb = b->x; |

653 | n/a | xbe = xb + wb; |

654 | n/a | xc0 = c->x; |

655 | n/a | for(; xb < xbe; xc0++) { |

656 | n/a | if ((y = *xb++)) { |

657 | n/a | x = xa; |

658 | n/a | xc = xc0; |

659 | n/a | carry = 0; |

660 | n/a | do { |

661 | n/a | z = *x++ * (ULLong)y + *xc + carry; |

662 | n/a | carry = z >> 32; |

663 | n/a | *xc++ = (ULong)(z & FFFFFFFF); |

664 | n/a | } |

665 | n/a | while(x < xae); |

666 | n/a | *xc = (ULong)carry; |

667 | n/a | } |

668 | n/a | } |

669 | n/a | for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; |

670 | n/a | c->wds = wc; |

671 | n/a | return c; |

672 | n/a | } |

673 | n/a | |

674 | n/a | #ifndef Py_USING_MEMORY_DEBUGGER |

675 | n/a | |

676 | n/a | /* p5s is a linked list of powers of 5 of the form 5**(2**i), i >= 2 */ |

677 | n/a | |

678 | n/a | static Bigint *p5s; |

679 | n/a | |

680 | n/a | /* multiply the Bigint b by 5**k. Returns a pointer to the result, or NULL on |

681 | n/a | failure; if the returned pointer is distinct from b then the original |

682 | n/a | Bigint b will have been Bfree'd. Ignores the sign of b. */ |

683 | n/a | |

684 | n/a | static Bigint * |

685 | n/a | pow5mult(Bigint *b, int k) |

686 | n/a | { |

687 | n/a | Bigint *b1, *p5, *p51; |

688 | n/a | int i; |

689 | n/a | static const int p05[3] = { 5, 25, 125 }; |

690 | n/a | |

691 | n/a | if ((i = k & 3)) { |

692 | n/a | b = multadd(b, p05[i-1], 0); |

693 | n/a | if (b == NULL) |

694 | n/a | return NULL; |

695 | n/a | } |

696 | n/a | |

697 | n/a | if (!(k >>= 2)) |

698 | n/a | return b; |

699 | n/a | p5 = p5s; |

700 | n/a | if (!p5) { |

701 | n/a | /* first time */ |

702 | n/a | p5 = i2b(625); |

703 | n/a | if (p5 == NULL) { |

704 | n/a | Bfree(b); |

705 | n/a | return NULL; |

706 | n/a | } |

707 | n/a | p5s = p5; |

708 | n/a | p5->next = 0; |

709 | n/a | } |

710 | n/a | for(;;) { |

711 | n/a | if (k & 1) { |

712 | n/a | b1 = mult(b, p5); |

713 | n/a | Bfree(b); |

714 | n/a | b = b1; |

715 | n/a | if (b == NULL) |

716 | n/a | return NULL; |

717 | n/a | } |

718 | n/a | if (!(k >>= 1)) |

719 | n/a | break; |

720 | n/a | p51 = p5->next; |

721 | n/a | if (!p51) { |

722 | n/a | p51 = mult(p5,p5); |

723 | n/a | if (p51 == NULL) { |

724 | n/a | Bfree(b); |

725 | n/a | return NULL; |

726 | n/a | } |

727 | n/a | p51->next = 0; |

728 | n/a | p5->next = p51; |

729 | n/a | } |

730 | n/a | p5 = p51; |

731 | n/a | } |

732 | n/a | return b; |

733 | n/a | } |

734 | n/a | |

735 | n/a | #else |

736 | n/a | |

737 | n/a | /* Version of pow5mult that doesn't cache powers of 5. Provided for |

738 | n/a | the benefit of memory debugging tools like Valgrind. */ |

739 | n/a | |

740 | n/a | static Bigint * |

741 | n/a | pow5mult(Bigint *b, int k) |

742 | n/a | { |

743 | n/a | Bigint *b1, *p5, *p51; |

744 | n/a | int i; |

745 | n/a | static const int p05[3] = { 5, 25, 125 }; |

746 | n/a | |

747 | n/a | if ((i = k & 3)) { |

748 | n/a | b = multadd(b, p05[i-1], 0); |

749 | n/a | if (b == NULL) |

750 | n/a | return NULL; |

751 | n/a | } |

752 | n/a | |

753 | n/a | if (!(k >>= 2)) |

754 | n/a | return b; |

755 | n/a | p5 = i2b(625); |

756 | n/a | if (p5 == NULL) { |

757 | n/a | Bfree(b); |

758 | n/a | return NULL; |

759 | n/a | } |

760 | n/a | |

761 | n/a | for(;;) { |

762 | n/a | if (k & 1) { |

763 | n/a | b1 = mult(b, p5); |

764 | n/a | Bfree(b); |

765 | n/a | b = b1; |

766 | n/a | if (b == NULL) { |

767 | n/a | Bfree(p5); |

768 | n/a | return NULL; |

769 | n/a | } |

770 | n/a | } |

771 | n/a | if (!(k >>= 1)) |

772 | n/a | break; |

773 | n/a | p51 = mult(p5, p5); |

774 | n/a | Bfree(p5); |

775 | n/a | p5 = p51; |

776 | n/a | if (p5 == NULL) { |

777 | n/a | Bfree(b); |

778 | n/a | return NULL; |

779 | n/a | } |

780 | n/a | } |

781 | n/a | Bfree(p5); |

782 | n/a | return b; |

783 | n/a | } |

784 | n/a | |

785 | n/a | #endif /* Py_USING_MEMORY_DEBUGGER */ |

786 | n/a | |

787 | n/a | /* shift a Bigint b left by k bits. Return a pointer to the shifted result, |

788 | n/a | or NULL on failure. If the returned pointer is distinct from b then the |

789 | n/a | original b will have been Bfree'd. Ignores the sign of b. */ |

790 | n/a | |

791 | n/a | static Bigint * |

792 | n/a | lshift(Bigint *b, int k) |

793 | n/a | { |

794 | n/a | int i, k1, n, n1; |

795 | n/a | Bigint *b1; |

796 | n/a | ULong *x, *x1, *xe, z; |

797 | n/a | |

798 | n/a | if (!k || (!b->x[0] && b->wds == 1)) |

799 | n/a | return b; |

800 | n/a | |

801 | n/a | n = k >> 5; |

802 | n/a | k1 = b->k; |

803 | n/a | n1 = n + b->wds + 1; |

804 | n/a | for(i = b->maxwds; n1 > i; i <<= 1) |

805 | n/a | k1++; |

806 | n/a | b1 = Balloc(k1); |

807 | n/a | if (b1 == NULL) { |

808 | n/a | Bfree(b); |

809 | n/a | return NULL; |

810 | n/a | } |

811 | n/a | x1 = b1->x; |

812 | n/a | for(i = 0; i < n; i++) |

813 | n/a | *x1++ = 0; |

814 | n/a | x = b->x; |

815 | n/a | xe = x + b->wds; |

816 | n/a | if (k &= 0x1f) { |

817 | n/a | k1 = 32 - k; |

818 | n/a | z = 0; |

819 | n/a | do { |

820 | n/a | *x1++ = *x << k | z; |

821 | n/a | z = *x++ >> k1; |

822 | n/a | } |

823 | n/a | while(x < xe); |

824 | n/a | if ((*x1 = z)) |

825 | n/a | ++n1; |

826 | n/a | } |

827 | n/a | else do |

828 | n/a | *x1++ = *x++; |

829 | n/a | while(x < xe); |

830 | n/a | b1->wds = n1 - 1; |

831 | n/a | Bfree(b); |

832 | n/a | return b1; |

833 | n/a | } |

834 | n/a | |

835 | n/a | /* Do a three-way compare of a and b, returning -1 if a < b, 0 if a == b and |

836 | n/a | 1 if a > b. Ignores signs of a and b. */ |

837 | n/a | |

838 | n/a | static int |

839 | n/a | cmp(Bigint *a, Bigint *b) |

840 | n/a | { |

841 | n/a | ULong *xa, *xa0, *xb, *xb0; |

842 | n/a | int i, j; |

843 | n/a | |

844 | n/a | i = a->wds; |

845 | n/a | j = b->wds; |

846 | n/a | #ifdef DEBUG |

847 | n/a | if (i > 1 && !a->x[i-1]) |

848 | n/a | Bug("cmp called with a->x[a->wds-1] == 0"); |

849 | n/a | if (j > 1 && !b->x[j-1]) |

850 | n/a | Bug("cmp called with b->x[b->wds-1] == 0"); |

851 | n/a | #endif |

852 | n/a | if (i -= j) |

853 | n/a | return i; |

854 | n/a | xa0 = a->x; |

855 | n/a | xa = xa0 + j; |

856 | n/a | xb0 = b->x; |

857 | n/a | xb = xb0 + j; |

858 | n/a | for(;;) { |

859 | n/a | if (*--xa != *--xb) |

860 | n/a | return *xa < *xb ? -1 : 1; |

861 | n/a | if (xa <= xa0) |

862 | n/a | break; |

863 | n/a | } |

864 | n/a | return 0; |

865 | n/a | } |

866 | n/a | |

867 | n/a | /* Take the difference of Bigints a and b, returning a new Bigint. Returns |

868 | n/a | NULL on failure. The signs of a and b are ignored, but the sign of the |

869 | n/a | result is set appropriately. */ |

870 | n/a | |

871 | n/a | static Bigint * |

872 | n/a | diff(Bigint *a, Bigint *b) |

873 | n/a | { |

874 | n/a | Bigint *c; |

875 | n/a | int i, wa, wb; |

876 | n/a | ULong *xa, *xae, *xb, *xbe, *xc; |

877 | n/a | ULLong borrow, y; |

878 | n/a | |

879 | n/a | i = cmp(a,b); |

880 | n/a | if (!i) { |

881 | n/a | c = Balloc(0); |

882 | n/a | if (c == NULL) |

883 | n/a | return NULL; |

884 | n/a | c->wds = 1; |

885 | n/a | c->x[0] = 0; |

886 | n/a | return c; |

887 | n/a | } |

888 | n/a | if (i < 0) { |

889 | n/a | c = a; |

890 | n/a | a = b; |

891 | n/a | b = c; |

892 | n/a | i = 1; |

893 | n/a | } |

894 | n/a | else |

895 | n/a | i = 0; |

896 | n/a | c = Balloc(a->k); |

897 | n/a | if (c == NULL) |

898 | n/a | return NULL; |

899 | n/a | c->sign = i; |

900 | n/a | wa = a->wds; |

901 | n/a | xa = a->x; |

902 | n/a | xae = xa + wa; |

903 | n/a | wb = b->wds; |

904 | n/a | xb = b->x; |

905 | n/a | xbe = xb + wb; |

906 | n/a | xc = c->x; |

907 | n/a | borrow = 0; |

908 | n/a | do { |

909 | n/a | y = (ULLong)*xa++ - *xb++ - borrow; |

910 | n/a | borrow = y >> 32 & (ULong)1; |

911 | n/a | *xc++ = (ULong)(y & FFFFFFFF); |

912 | n/a | } |

913 | n/a | while(xb < xbe); |

914 | n/a | while(xa < xae) { |

915 | n/a | y = *xa++ - borrow; |

916 | n/a | borrow = y >> 32 & (ULong)1; |

917 | n/a | *xc++ = (ULong)(y & FFFFFFFF); |

918 | n/a | } |

919 | n/a | while(!*--xc) |

920 | n/a | wa--; |

921 | n/a | c->wds = wa; |

922 | n/a | return c; |

923 | n/a | } |

924 | n/a | |

925 | n/a | /* Given a positive normal double x, return the difference between x and the |

926 | n/a | next double up. Doesn't give correct results for subnormals. */ |

927 | n/a | |

928 | n/a | static double |

929 | n/a | ulp(U *x) |

930 | n/a | { |

931 | n/a | Long L; |

932 | n/a | U u; |

933 | n/a | |

934 | n/a | L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1; |

935 | n/a | word0(&u) = L; |

936 | n/a | word1(&u) = 0; |

937 | n/a | return dval(&u); |

938 | n/a | } |

939 | n/a | |

940 | n/a | /* Convert a Bigint to a double plus an exponent */ |

941 | n/a | |

942 | n/a | static double |

943 | n/a | b2d(Bigint *a, int *e) |

944 | n/a | { |

945 | n/a | ULong *xa, *xa0, w, y, z; |

946 | n/a | int k; |

947 | n/a | U d; |

948 | n/a | |

949 | n/a | xa0 = a->x; |

950 | n/a | xa = xa0 + a->wds; |

951 | n/a | y = *--xa; |

952 | n/a | #ifdef DEBUG |

953 | n/a | if (!y) Bug("zero y in b2d"); |

954 | n/a | #endif |

955 | n/a | k = hi0bits(y); |

956 | n/a | *e = 32 - k; |

957 | n/a | if (k < Ebits) { |

958 | n/a | word0(&d) = Exp_1 | y >> (Ebits - k); |

959 | n/a | w = xa > xa0 ? *--xa : 0; |

960 | n/a | word1(&d) = y << ((32-Ebits) + k) | w >> (Ebits - k); |

961 | n/a | goto ret_d; |

962 | n/a | } |

963 | n/a | z = xa > xa0 ? *--xa : 0; |

964 | n/a | if (k -= Ebits) { |

965 | n/a | word0(&d) = Exp_1 | y << k | z >> (32 - k); |

966 | n/a | y = xa > xa0 ? *--xa : 0; |

967 | n/a | word1(&d) = z << k | y >> (32 - k); |

968 | n/a | } |

969 | n/a | else { |

970 | n/a | word0(&d) = Exp_1 | y; |

971 | n/a | word1(&d) = z; |

972 | n/a | } |

973 | n/a | ret_d: |

974 | n/a | return dval(&d); |

975 | n/a | } |

976 | n/a | |

977 | n/a | /* Convert a scaled double to a Bigint plus an exponent. Similar to d2b, |

978 | n/a | except that it accepts the scale parameter used in _Py_dg_strtod (which |

979 | n/a | should be either 0 or 2*P), and the normalization for the return value is |

980 | n/a | different (see below). On input, d should be finite and nonnegative, and d |

981 | n/a | / 2**scale should be exactly representable as an IEEE 754 double. |

982 | n/a | |

983 | n/a | Returns a Bigint b and an integer e such that |

984 | n/a | |

985 | n/a | dval(d) / 2**scale = b * 2**e. |

986 | n/a | |

987 | n/a | Unlike d2b, b is not necessarily odd: b and e are normalized so |

988 | n/a | that either 2**(P-1) <= b < 2**P and e >= Etiny, or b < 2**P |

989 | n/a | and e == Etiny. This applies equally to an input of 0.0: in that |

990 | n/a | case the return values are b = 0 and e = Etiny. |

991 | n/a | |

992 | n/a | The above normalization ensures that for all possible inputs d, |

993 | n/a | 2**e gives ulp(d/2**scale). |

994 | n/a | |

995 | n/a | Returns NULL on failure. |

996 | n/a | */ |

997 | n/a | |

998 | n/a | static Bigint * |

999 | n/a | sd2b(U *d, int scale, int *e) |

1000 | n/a | { |

1001 | n/a | Bigint *b; |

1002 | n/a | |

1003 | n/a | b = Balloc(1); |

1004 | n/a | if (b == NULL) |

1005 | n/a | return NULL; |

1006 | n/a | |

1007 | n/a | /* First construct b and e assuming that scale == 0. */ |

1008 | n/a | b->wds = 2; |

1009 | n/a | b->x[0] = word1(d); |

1010 | n/a | b->x[1] = word0(d) & Frac_mask; |

1011 | n/a | *e = Etiny - 1 + (int)((word0(d) & Exp_mask) >> Exp_shift); |

1012 | n/a | if (*e < Etiny) |

1013 | n/a | *e = Etiny; |

1014 | n/a | else |

1015 | n/a | b->x[1] |= Exp_msk1; |

1016 | n/a | |

1017 | n/a | /* Now adjust for scale, provided that b != 0. */ |

1018 | n/a | if (scale && (b->x[0] || b->x[1])) { |

1019 | n/a | *e -= scale; |

1020 | n/a | if (*e < Etiny) { |

1021 | n/a | scale = Etiny - *e; |

1022 | n/a | *e = Etiny; |

1023 | n/a | /* We can't shift more than P-1 bits without shifting out a 1. */ |

1024 | n/a | assert(0 < scale && scale <= P - 1); |

1025 | n/a | if (scale >= 32) { |

1026 | n/a | /* The bits shifted out should all be zero. */ |

1027 | n/a | assert(b->x[0] == 0); |

1028 | n/a | b->x[0] = b->x[1]; |

1029 | n/a | b->x[1] = 0; |

1030 | n/a | scale -= 32; |

1031 | n/a | } |

1032 | n/a | if (scale) { |

1033 | n/a | /* The bits shifted out should all be zero. */ |

1034 | n/a | assert(b->x[0] << (32 - scale) == 0); |

1035 | n/a | b->x[0] = (b->x[0] >> scale) | (b->x[1] << (32 - scale)); |

1036 | n/a | b->x[1] >>= scale; |

1037 | n/a | } |

1038 | n/a | } |

1039 | n/a | } |

1040 | n/a | /* Ensure b is normalized. */ |

1041 | n/a | if (!b->x[1]) |

1042 | n/a | b->wds = 1; |

1043 | n/a | |

1044 | n/a | return b; |

1045 | n/a | } |

1046 | n/a | |

1047 | n/a | /* Convert a double to a Bigint plus an exponent. Return NULL on failure. |

1048 | n/a | |

1049 | n/a | Given a finite nonzero double d, return an odd Bigint b and exponent *e |

1050 | n/a | such that fabs(d) = b * 2**e. On return, *bbits gives the number of |

1051 | n/a | significant bits of b; that is, 2**(*bbits-1) <= b < 2**(*bbits). |

1052 | n/a | |

1053 | n/a | If d is zero, then b == 0, *e == -1010, *bbits = 0. |

1054 | n/a | */ |

1055 | n/a | |

1056 | n/a | static Bigint * |

1057 | n/a | d2b(U *d, int *e, int *bits) |

1058 | n/a | { |

1059 | n/a | Bigint *b; |

1060 | n/a | int de, k; |

1061 | n/a | ULong *x, y, z; |

1062 | n/a | int i; |

1063 | n/a | |

1064 | n/a | b = Balloc(1); |

1065 | n/a | if (b == NULL) |

1066 | n/a | return NULL; |

1067 | n/a | x = b->x; |

1068 | n/a | |

1069 | n/a | z = word0(d) & Frac_mask; |

1070 | n/a | word0(d) &= 0x7fffffff; /* clear sign bit, which we ignore */ |

1071 | n/a | if ((de = (int)(word0(d) >> Exp_shift))) |

1072 | n/a | z |= Exp_msk1; |

1073 | n/a | if ((y = word1(d))) { |

1074 | n/a | if ((k = lo0bits(&y))) { |

1075 | n/a | x[0] = y | z << (32 - k); |

1076 | n/a | z >>= k; |

1077 | n/a | } |

1078 | n/a | else |

1079 | n/a | x[0] = y; |

1080 | n/a | i = |

1081 | n/a | b->wds = (x[1] = z) ? 2 : 1; |

1082 | n/a | } |

1083 | n/a | else { |

1084 | n/a | k = lo0bits(&z); |

1085 | n/a | x[0] = z; |

1086 | n/a | i = |

1087 | n/a | b->wds = 1; |

1088 | n/a | k += 32; |

1089 | n/a | } |

1090 | n/a | if (de) { |

1091 | n/a | *e = de - Bias - (P-1) + k; |

1092 | n/a | *bits = P - k; |

1093 | n/a | } |

1094 | n/a | else { |

1095 | n/a | *e = de - Bias - (P-1) + 1 + k; |

1096 | n/a | *bits = 32*i - hi0bits(x[i-1]); |

1097 | n/a | } |

1098 | n/a | return b; |

1099 | n/a | } |

1100 | n/a | |

1101 | n/a | /* Compute the ratio of two Bigints, as a double. The result may have an |

1102 | n/a | error of up to 2.5 ulps. */ |

1103 | n/a | |

1104 | n/a | static double |

1105 | n/a | ratio(Bigint *a, Bigint *b) |

1106 | n/a | { |

1107 | n/a | U da, db; |

1108 | n/a | int k, ka, kb; |

1109 | n/a | |

1110 | n/a | dval(&da) = b2d(a, &ka); |

1111 | n/a | dval(&db) = b2d(b, &kb); |

1112 | n/a | k = ka - kb + 32*(a->wds - b->wds); |

1113 | n/a | if (k > 0) |

1114 | n/a | word0(&da) += k*Exp_msk1; |

1115 | n/a | else { |

1116 | n/a | k = -k; |

1117 | n/a | word0(&db) += k*Exp_msk1; |

1118 | n/a | } |

1119 | n/a | return dval(&da) / dval(&db); |

1120 | n/a | } |

1121 | n/a | |

1122 | n/a | static const double |

1123 | n/a | tens[] = { |

1124 | n/a | 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, |

1125 | n/a | 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, |

1126 | n/a | 1e20, 1e21, 1e22 |

1127 | n/a | }; |

1128 | n/a | |

1129 | n/a | static const double |

1130 | n/a | bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; |

1131 | n/a | static const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, |

1132 | n/a | 9007199254740992.*9007199254740992.e-256 |

1133 | n/a | /* = 2^106 * 1e-256 */ |

1134 | n/a | }; |

1135 | n/a | /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ |

1136 | n/a | /* flag unnecessarily. It leads to a song and dance at the end of strtod. */ |

1137 | n/a | #define Scale_Bit 0x10 |

1138 | n/a | #define n_bigtens 5 |

1139 | n/a | |

1140 | n/a | #define ULbits 32 |

1141 | n/a | #define kshift 5 |

1142 | n/a | #define kmask 31 |

1143 | n/a | |

1144 | n/a | |

1145 | n/a | static int |

1146 | n/a | dshift(Bigint *b, int p2) |

1147 | n/a | { |

1148 | n/a | int rv = hi0bits(b->x[b->wds-1]) - 4; |

1149 | n/a | if (p2 > 0) |

1150 | n/a | rv -= p2; |

1151 | n/a | return rv & kmask; |

1152 | n/a | } |

1153 | n/a | |

1154 | n/a | /* special case of Bigint division. The quotient is always in the range 0 <= |

1155 | n/a | quotient < 10, and on entry the divisor S is normalized so that its top 4 |

1156 | n/a | bits (28--31) are zero and bit 27 is set. */ |

1157 | n/a | |

1158 | n/a | static int |

1159 | n/a | quorem(Bigint *b, Bigint *S) |

1160 | n/a | { |

1161 | n/a | int n; |

1162 | n/a | ULong *bx, *bxe, q, *sx, *sxe; |

1163 | n/a | ULLong borrow, carry, y, ys; |

1164 | n/a | |

1165 | n/a | n = S->wds; |

1166 | n/a | #ifdef DEBUG |

1167 | n/a | /*debug*/ if (b->wds > n) |

1168 | n/a | /*debug*/ Bug("oversize b in quorem"); |

1169 | n/a | #endif |

1170 | n/a | if (b->wds < n) |

1171 | n/a | return 0; |

1172 | n/a | sx = S->x; |

1173 | n/a | sxe = sx + --n; |

1174 | n/a | bx = b->x; |

1175 | n/a | bxe = bx + n; |

1176 | n/a | q = *bxe / (*sxe + 1); /* ensure q <= true quotient */ |

1177 | n/a | #ifdef DEBUG |

1178 | n/a | /*debug*/ if (q > 9) |

1179 | n/a | /*debug*/ Bug("oversized quotient in quorem"); |

1180 | n/a | #endif |

1181 | n/a | if (q) { |

1182 | n/a | borrow = 0; |

1183 | n/a | carry = 0; |

1184 | n/a | do { |

1185 | n/a | ys = *sx++ * (ULLong)q + carry; |

1186 | n/a | carry = ys >> 32; |

1187 | n/a | y = *bx - (ys & FFFFFFFF) - borrow; |

1188 | n/a | borrow = y >> 32 & (ULong)1; |

1189 | n/a | *bx++ = (ULong)(y & FFFFFFFF); |

1190 | n/a | } |

1191 | n/a | while(sx <= sxe); |

1192 | n/a | if (!*bxe) { |

1193 | n/a | bx = b->x; |

1194 | n/a | while(--bxe > bx && !*bxe) |

1195 | n/a | --n; |

1196 | n/a | b->wds = n; |

1197 | n/a | } |

1198 | n/a | } |

1199 | n/a | if (cmp(b, S) >= 0) { |

1200 | n/a | q++; |

1201 | n/a | borrow = 0; |

1202 | n/a | carry = 0; |

1203 | n/a | bx = b->x; |

1204 | n/a | sx = S->x; |

1205 | n/a | do { |

1206 | n/a | ys = *sx++ + carry; |

1207 | n/a | carry = ys >> 32; |

1208 | n/a | y = *bx - (ys & FFFFFFFF) - borrow; |

1209 | n/a | borrow = y >> 32 & (ULong)1; |

1210 | n/a | *bx++ = (ULong)(y & FFFFFFFF); |

1211 | n/a | } |

1212 | n/a | while(sx <= sxe); |

1213 | n/a | bx = b->x; |

1214 | n/a | bxe = bx + n; |

1215 | n/a | if (!*bxe) { |

1216 | n/a | while(--bxe > bx && !*bxe) |

1217 | n/a | --n; |

1218 | n/a | b->wds = n; |

1219 | n/a | } |

1220 | n/a | } |

1221 | n/a | return q; |

1222 | n/a | } |

1223 | n/a | |

1224 | n/a | /* sulp(x) is a version of ulp(x) that takes bc.scale into account. |

1225 | n/a | |

1226 | n/a | Assuming that x is finite and nonnegative (positive zero is fine |

1227 | n/a | here) and x / 2^bc.scale is exactly representable as a double, |

1228 | n/a | sulp(x) is equivalent to 2^bc.scale * ulp(x / 2^bc.scale). */ |

1229 | n/a | |

1230 | n/a | static double |

1231 | n/a | sulp(U *x, BCinfo *bc) |

1232 | n/a | { |

1233 | n/a | U u; |

1234 | n/a | |

1235 | n/a | if (bc->scale && 2*P + 1 > (int)((word0(x) & Exp_mask) >> Exp_shift)) { |

1236 | n/a | /* rv/2^bc->scale is subnormal */ |

1237 | n/a | word0(&u) = (P+2)*Exp_msk1; |

1238 | n/a | word1(&u) = 0; |

1239 | n/a | return u.d; |

1240 | n/a | } |

1241 | n/a | else { |

1242 | n/a | assert(word0(x) || word1(x)); /* x != 0.0 */ |

1243 | n/a | return ulp(x); |

1244 | n/a | } |

1245 | n/a | } |

1246 | n/a | |

1247 | n/a | /* The bigcomp function handles some hard cases for strtod, for inputs |

1248 | n/a | with more than STRTOD_DIGLIM digits. It's called once an initial |

1249 | n/a | estimate for the double corresponding to the input string has |

1250 | n/a | already been obtained by the code in _Py_dg_strtod. |

1251 | n/a | |

1252 | n/a | The bigcomp function is only called after _Py_dg_strtod has found a |

1253 | n/a | double value rv such that either rv or rv + 1ulp represents the |

1254 | n/a | correctly rounded value corresponding to the original string. It |

1255 | n/a | determines which of these two values is the correct one by |

1256 | n/a | computing the decimal digits of rv + 0.5ulp and comparing them with |

1257 | n/a | the corresponding digits of s0. |

1258 | n/a | |

1259 | n/a | In the following, write dv for the absolute value of the number represented |

1260 | n/a | by the input string. |

1261 | n/a | |

1262 | n/a | Inputs: |

1263 | n/a | |

1264 | n/a | s0 points to the first significant digit of the input string. |

1265 | n/a | |

1266 | n/a | rv is a (possibly scaled) estimate for the closest double value to the |

1267 | n/a | value represented by the original input to _Py_dg_strtod. If |

1268 | n/a | bc->scale is nonzero, then rv/2^(bc->scale) is the approximation to |

1269 | n/a | the input value. |

1270 | n/a | |

1271 | n/a | bc is a struct containing information gathered during the parsing and |

1272 | n/a | estimation steps of _Py_dg_strtod. Description of fields follows: |

1273 | n/a | |

1274 | n/a | bc->e0 gives the exponent of the input value, such that dv = (integer |

1275 | n/a | given by the bd->nd digits of s0) * 10**e0 |

1276 | n/a | |

1277 | n/a | bc->nd gives the total number of significant digits of s0. It will |

1278 | n/a | be at least 1. |

1279 | n/a | |

1280 | n/a | bc->nd0 gives the number of significant digits of s0 before the |

1281 | n/a | decimal separator. If there's no decimal separator, bc->nd0 == |

1282 | n/a | bc->nd. |

1283 | n/a | |

1284 | n/a | bc->scale is the value used to scale rv to avoid doing arithmetic with |

1285 | n/a | subnormal values. It's either 0 or 2*P (=106). |

1286 | n/a | |

1287 | n/a | Outputs: |

1288 | n/a | |

1289 | n/a | On successful exit, rv/2^(bc->scale) is the closest double to dv. |

1290 | n/a | |

1291 | n/a | Returns 0 on success, -1 on failure (e.g., due to a failed malloc call). */ |

1292 | n/a | |

1293 | n/a | static int |

1294 | n/a | bigcomp(U *rv, const char *s0, BCinfo *bc) |

1295 | n/a | { |

1296 | n/a | Bigint *b, *d; |

1297 | n/a | int b2, d2, dd, i, nd, nd0, odd, p2, p5; |

1298 | n/a | |

1299 | n/a | nd = bc->nd; |

1300 | n/a | nd0 = bc->nd0; |

1301 | n/a | p5 = nd + bc->e0; |

1302 | n/a | b = sd2b(rv, bc->scale, &p2); |

1303 | n/a | if (b == NULL) |

1304 | n/a | return -1; |

1305 | n/a | |

1306 | n/a | /* record whether the lsb of rv/2^(bc->scale) is odd: in the exact halfway |

1307 | n/a | case, this is used for round to even. */ |

1308 | n/a | odd = b->x[0] & 1; |

1309 | n/a | |

1310 | n/a | /* left shift b by 1 bit and or a 1 into the least significant bit; |

1311 | n/a | this gives us b * 2**p2 = rv/2^(bc->scale) + 0.5 ulp. */ |

1312 | n/a | b = lshift(b, 1); |

1313 | n/a | if (b == NULL) |

1314 | n/a | return -1; |

1315 | n/a | b->x[0] |= 1; |

1316 | n/a | p2--; |

1317 | n/a | |

1318 | n/a | p2 -= p5; |

1319 | n/a | d = i2b(1); |

1320 | n/a | if (d == NULL) { |

1321 | n/a | Bfree(b); |

1322 | n/a | return -1; |

1323 | n/a | } |

1324 | n/a | /* Arrange for convenient computation of quotients: |

1325 | n/a | * shift left if necessary so divisor has 4 leading 0 bits. |

1326 | n/a | */ |

1327 | n/a | if (p5 > 0) { |

1328 | n/a | d = pow5mult(d, p5); |

1329 | n/a | if (d == NULL) { |

1330 | n/a | Bfree(b); |

1331 | n/a | return -1; |

1332 | n/a | } |

1333 | n/a | } |

1334 | n/a | else if (p5 < 0) { |

1335 | n/a | b = pow5mult(b, -p5); |

1336 | n/a | if (b == NULL) { |

1337 | n/a | Bfree(d); |

1338 | n/a | return -1; |

1339 | n/a | } |

1340 | n/a | } |

1341 | n/a | if (p2 > 0) { |

1342 | n/a | b2 = p2; |

1343 | n/a | d2 = 0; |

1344 | n/a | } |

1345 | n/a | else { |

1346 | n/a | b2 = 0; |

1347 | n/a | d2 = -p2; |

1348 | n/a | } |

1349 | n/a | i = dshift(d, d2); |

1350 | n/a | if ((b2 += i) > 0) { |

1351 | n/a | b = lshift(b, b2); |

1352 | n/a | if (b == NULL) { |

1353 | n/a | Bfree(d); |

1354 | n/a | return -1; |

1355 | n/a | } |

1356 | n/a | } |

1357 | n/a | if ((d2 += i) > 0) { |

1358 | n/a | d = lshift(d, d2); |

1359 | n/a | if (d == NULL) { |

1360 | n/a | Bfree(b); |

1361 | n/a | return -1; |

1362 | n/a | } |

1363 | n/a | } |

1364 | n/a | |

1365 | n/a | /* Compare s0 with b/d: set dd to -1, 0, or 1 according as s0 < b/d, s0 == |

1366 | n/a | * b/d, or s0 > b/d. Here the digits of s0 are thought of as representing |

1367 | n/a | * a number in the range [0.1, 1). */ |

1368 | n/a | if (cmp(b, d) >= 0) |

1369 | n/a | /* b/d >= 1 */ |

1370 | n/a | dd = -1; |

1371 | n/a | else { |

1372 | n/a | i = 0; |

1373 | n/a | for(;;) { |

1374 | n/a | b = multadd(b, 10, 0); |

1375 | n/a | if (b == NULL) { |

1376 | n/a | Bfree(d); |

1377 | n/a | return -1; |

1378 | n/a | } |

1379 | n/a | dd = s0[i < nd0 ? i : i+1] - '0' - quorem(b, d); |

1380 | n/a | i++; |

1381 | n/a | |

1382 | n/a | if (dd) |

1383 | n/a | break; |

1384 | n/a | if (!b->x[0] && b->wds == 1) { |

1385 | n/a | /* b/d == 0 */ |

1386 | n/a | dd = i < nd; |

1387 | n/a | break; |

1388 | n/a | } |

1389 | n/a | if (!(i < nd)) { |

1390 | n/a | /* b/d != 0, but digits of s0 exhausted */ |

1391 | n/a | dd = -1; |

1392 | n/a | break; |

1393 | n/a | } |

1394 | n/a | } |

1395 | n/a | } |

1396 | n/a | Bfree(b); |

1397 | n/a | Bfree(d); |

1398 | n/a | if (dd > 0 || (dd == 0 && odd)) |

1399 | n/a | dval(rv) += sulp(rv, bc); |

1400 | n/a | return 0; |

1401 | n/a | } |

1402 | n/a | |

1403 | n/a | /* Return a 'standard' NaN value. |

1404 | n/a | |

1405 | n/a | There are exactly two quiet NaNs that don't arise by 'quieting' signaling |

1406 | n/a | NaNs (see IEEE 754-2008, section 6.2.1). If sign == 0, return the one whose |

1407 | n/a | sign bit is cleared. Otherwise, return the one whose sign bit is set. |

1408 | n/a | */ |

1409 | n/a | |

1410 | n/a | double |

1411 | n/a | _Py_dg_stdnan(int sign) |

1412 | n/a | { |

1413 | n/a | U rv; |

1414 | n/a | word0(&rv) = NAN_WORD0; |

1415 | n/a | word1(&rv) = NAN_WORD1; |

1416 | n/a | if (sign) |

1417 | n/a | word0(&rv) |= Sign_bit; |

1418 | n/a | return dval(&rv); |

1419 | n/a | } |

1420 | n/a | |

1421 | n/a | /* Return positive or negative infinity, according to the given sign (0 for |

1422 | n/a | * positive infinity, 1 for negative infinity). */ |

1423 | n/a | |

1424 | n/a | double |

1425 | n/a | _Py_dg_infinity(int sign) |

1426 | n/a | { |

1427 | n/a | U rv; |

1428 | n/a | word0(&rv) = POSINF_WORD0; |

1429 | n/a | word1(&rv) = POSINF_WORD1; |

1430 | n/a | return sign ? -dval(&rv) : dval(&rv); |

1431 | n/a | } |

1432 | n/a | |

1433 | n/a | double |

1434 | n/a | _Py_dg_strtod(const char *s00, char **se) |

1435 | n/a | { |

1436 | n/a | int bb2, bb5, bbe, bd2, bd5, bs2, c, dsign, e, e1, error; |

1437 | n/a | int esign, i, j, k, lz, nd, nd0, odd, sign; |

1438 | n/a | const char *s, *s0, *s1; |

1439 | n/a | double aadj, aadj1; |

1440 | n/a | U aadj2, adj, rv, rv0; |

1441 | n/a | ULong y, z, abs_exp; |

1442 | n/a | Long L; |

1443 | n/a | BCinfo bc; |

1444 | n/a | Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; |

1445 | n/a | size_t ndigits, fraclen; |

1446 | n/a | |

1447 | n/a | dval(&rv) = 0.; |

1448 | n/a | |

1449 | n/a | /* Start parsing. */ |

1450 | n/a | c = *(s = s00); |

1451 | n/a | |

1452 | n/a | /* Parse optional sign, if present. */ |

1453 | n/a | sign = 0; |

1454 | n/a | switch (c) { |

1455 | n/a | case '-': |

1456 | n/a | sign = 1; |

1457 | n/a | /* no break */ |

1458 | n/a | case '+': |

1459 | n/a | c = *++s; |

1460 | n/a | } |

1461 | n/a | |

1462 | n/a | /* Skip leading zeros: lz is true iff there were leading zeros. */ |

1463 | n/a | s1 = s; |

1464 | n/a | while (c == '0') |

1465 | n/a | c = *++s; |

1466 | n/a | lz = s != s1; |

1467 | n/a | |

1468 | n/a | /* Point s0 at the first nonzero digit (if any). fraclen will be the |

1469 | n/a | number of digits between the decimal point and the end of the |

1470 | n/a | digit string. ndigits will be the total number of digits ignoring |

1471 | n/a | leading zeros. */ |

1472 | n/a | s0 = s1 = s; |

1473 | n/a | while ('0' <= c && c <= '9') |

1474 | n/a | c = *++s; |

1475 | n/a | ndigits = s - s1; |

1476 | n/a | fraclen = 0; |

1477 | n/a | |

1478 | n/a | /* Parse decimal point and following digits. */ |

1479 | n/a | if (c == '.') { |

1480 | n/a | c = *++s; |

1481 | n/a | if (!ndigits) { |

1482 | n/a | s1 = s; |

1483 | n/a | while (c == '0') |

1484 | n/a | c = *++s; |

1485 | n/a | lz = lz || s != s1; |

1486 | n/a | fraclen += (s - s1); |

1487 | n/a | s0 = s; |

1488 | n/a | } |

1489 | n/a | s1 = s; |

1490 | n/a | while ('0' <= c && c <= '9') |

1491 | n/a | c = *++s; |

1492 | n/a | ndigits += s - s1; |

1493 | n/a | fraclen += s - s1; |

1494 | n/a | } |

1495 | n/a | |

1496 | n/a | /* Now lz is true if and only if there were leading zero digits, and |

1497 | n/a | ndigits gives the total number of digits ignoring leading zeros. A |

1498 | n/a | valid input must have at least one digit. */ |

1499 | n/a | if (!ndigits && !lz) { |

1500 | n/a | if (se) |

1501 | n/a | *se = (char *)s00; |

1502 | n/a | goto parse_error; |

1503 | n/a | } |

1504 | n/a | |

1505 | n/a | /* Range check ndigits and fraclen to make sure that they, and values |

1506 | n/a | computed with them, can safely fit in an int. */ |

1507 | n/a | if (ndigits > MAX_DIGITS || fraclen > MAX_DIGITS) { |

1508 | n/a | if (se) |

1509 | n/a | *se = (char *)s00; |

1510 | n/a | goto parse_error; |

1511 | n/a | } |

1512 | n/a | nd = (int)ndigits; |

1513 | n/a | nd0 = (int)ndigits - (int)fraclen; |

1514 | n/a | |

1515 | n/a | /* Parse exponent. */ |

1516 | n/a | e = 0; |

1517 | n/a | if (c == 'e' || c == 'E') { |

1518 | n/a | s00 = s; |

1519 | n/a | c = *++s; |

1520 | n/a | |

1521 | n/a | /* Exponent sign. */ |

1522 | n/a | esign = 0; |

1523 | n/a | switch (c) { |

1524 | n/a | case '-': |

1525 | n/a | esign = 1; |

1526 | n/a | /* no break */ |

1527 | n/a | case '+': |

1528 | n/a | c = *++s; |

1529 | n/a | } |

1530 | n/a | |

1531 | n/a | /* Skip zeros. lz is true iff there are leading zeros. */ |

1532 | n/a | s1 = s; |

1533 | n/a | while (c == '0') |

1534 | n/a | c = *++s; |

1535 | n/a | lz = s != s1; |

1536 | n/a | |

1537 | n/a | /* Get absolute value of the exponent. */ |

1538 | n/a | s1 = s; |

1539 | n/a | abs_exp = 0; |

1540 | n/a | while ('0' <= c && c <= '9') { |

1541 | n/a | abs_exp = 10*abs_exp + (c - '0'); |

1542 | n/a | c = *++s; |

1543 | n/a | } |

1544 | n/a | |

1545 | n/a | /* abs_exp will be correct modulo 2**32. But 10**9 < 2**32, so if |

1546 | n/a | there are at most 9 significant exponent digits then overflow is |

1547 | n/a | impossible. */ |

1548 | n/a | if (s - s1 > 9 || abs_exp > MAX_ABS_EXP) |

1549 | n/a | e = (int)MAX_ABS_EXP; |

1550 | n/a | else |

1551 | n/a | e = (int)abs_exp; |

1552 | n/a | if (esign) |

1553 | n/a | e = -e; |

1554 | n/a | |

1555 | n/a | /* A valid exponent must have at least one digit. */ |

1556 | n/a | if (s == s1 && !lz) |

1557 | n/a | s = s00; |

1558 | n/a | } |

1559 | n/a | |

1560 | n/a | /* Adjust exponent to take into account position of the point. */ |

1561 | n/a | e -= nd - nd0; |

1562 | n/a | if (nd0 <= 0) |

1563 | n/a | nd0 = nd; |

1564 | n/a | |

1565 | n/a | /* Finished parsing. Set se to indicate how far we parsed */ |

1566 | n/a | if (se) |

1567 | n/a | *se = (char *)s; |

1568 | n/a | |

1569 | n/a | /* If all digits were zero, exit with return value +-0.0. Otherwise, |

1570 | n/a | strip trailing zeros: scan back until we hit a nonzero digit. */ |

1571 | n/a | if (!nd) |

1572 | n/a | goto ret; |

1573 | n/a | for (i = nd; i > 0; ) { |

1574 | n/a | --i; |

1575 | n/a | if (s0[i < nd0 ? i : i+1] != '0') { |

1576 | n/a | ++i; |

1577 | n/a | break; |

1578 | n/a | } |

1579 | n/a | } |

1580 | n/a | e += nd - i; |

1581 | n/a | nd = i; |

1582 | n/a | if (nd0 > nd) |

1583 | n/a | nd0 = nd; |

1584 | n/a | |

1585 | n/a | /* Summary of parsing results. After parsing, and dealing with zero |

1586 | n/a | * inputs, we have values s0, nd0, nd, e, sign, where: |

1587 | n/a | * |

1588 | n/a | * - s0 points to the first significant digit of the input string |

1589 | n/a | * |

1590 | n/a | * - nd is the total number of significant digits (here, and |

1591 | n/a | * below, 'significant digits' means the set of digits of the |

1592 | n/a | * significand of the input that remain after ignoring leading |

1593 | n/a | * and trailing zeros). |

1594 | n/a | * |

1595 | n/a | * - nd0 indicates the position of the decimal point, if present; it |

1596 | n/a | * satisfies 1 <= nd0 <= nd. The nd significant digits are in |

1597 | n/a | * s0[0:nd0] and s0[nd0+1:nd+1] using the usual Python half-open slice |

1598 | n/a | * notation. (If nd0 < nd, then s0[nd0] contains a '.' character; if |

1599 | n/a | * nd0 == nd, then s0[nd0] could be any non-digit character.) |

1600 | n/a | * |

1601 | n/a | * - e is the adjusted exponent: the absolute value of the number |

1602 | n/a | * represented by the original input string is n * 10**e, where |

1603 | n/a | * n is the integer represented by the concatenation of |

1604 | n/a | * s0[0:nd0] and s0[nd0+1:nd+1] |

1605 | n/a | * |

1606 | n/a | * - sign gives the sign of the input: 1 for negative, 0 for positive |

1607 | n/a | * |

1608 | n/a | * - the first and last significant digits are nonzero |

1609 | n/a | */ |

1610 | n/a | |

1611 | n/a | /* put first DBL_DIG+1 digits into integer y and z. |

1612 | n/a | * |

1613 | n/a | * - y contains the value represented by the first min(9, nd) |

1614 | n/a | * significant digits |

1615 | n/a | * |

1616 | n/a | * - if nd > 9, z contains the value represented by significant digits |

1617 | n/a | * with indices in [9, min(16, nd)). So y * 10**(min(16, nd) - 9) + z |

1618 | n/a | * gives the value represented by the first min(16, nd) sig. digits. |

1619 | n/a | */ |

1620 | n/a | |

1621 | n/a | bc.e0 = e1 = e; |

1622 | n/a | y = z = 0; |

1623 | n/a | for (i = 0; i < nd; i++) { |

1624 | n/a | if (i < 9) |

1625 | n/a | y = 10*y + s0[i < nd0 ? i : i+1] - '0'; |

1626 | n/a | else if (i < DBL_DIG+1) |

1627 | n/a | z = 10*z + s0[i < nd0 ? i : i+1] - '0'; |

1628 | n/a | else |

1629 | n/a | break; |

1630 | n/a | } |

1631 | n/a | |

1632 | n/a | k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; |

1633 | n/a | dval(&rv) = y; |

1634 | n/a | if (k > 9) { |

1635 | n/a | dval(&rv) = tens[k - 9] * dval(&rv) + z; |

1636 | n/a | } |

1637 | n/a | bd0 = 0; |

1638 | n/a | if (nd <= DBL_DIG |

1639 | n/a | && Flt_Rounds == 1 |

1640 | n/a | ) { |

1641 | n/a | if (!e) |

1642 | n/a | goto ret; |

1643 | n/a | if (e > 0) { |

1644 | n/a | if (e <= Ten_pmax) { |

1645 | n/a | dval(&rv) *= tens[e]; |

1646 | n/a | goto ret; |

1647 | n/a | } |

1648 | n/a | i = DBL_DIG - nd; |

1649 | n/a | if (e <= Ten_pmax + i) { |

1650 | n/a | /* A fancier test would sometimes let us do |

1651 | n/a | * this for larger i values. |

1652 | n/a | */ |

1653 | n/a | e -= i; |

1654 | n/a | dval(&rv) *= tens[i]; |

1655 | n/a | dval(&rv) *= tens[e]; |

1656 | n/a | goto ret; |

1657 | n/a | } |

1658 | n/a | } |

1659 | n/a | else if (e >= -Ten_pmax) { |

1660 | n/a | dval(&rv) /= tens[-e]; |

1661 | n/a | goto ret; |

1662 | n/a | } |

1663 | n/a | } |

1664 | n/a | e1 += nd - k; |

1665 | n/a | |

1666 | n/a | bc.scale = 0; |

1667 | n/a | |

1668 | n/a | /* Get starting approximation = rv * 10**e1 */ |

1669 | n/a | |

1670 | n/a | if (e1 > 0) { |

1671 | n/a | if ((i = e1 & 15)) |

1672 | n/a | dval(&rv) *= tens[i]; |

1673 | n/a | if (e1 &= ~15) { |

1674 | n/a | if (e1 > DBL_MAX_10_EXP) |

1675 | n/a | goto ovfl; |

1676 | n/a | e1 >>= 4; |

1677 | n/a | for(j = 0; e1 > 1; j++, e1 >>= 1) |

1678 | n/a | if (e1 & 1) |

1679 | n/a | dval(&rv) *= bigtens[j]; |

1680 | n/a | /* The last multiplication could overflow. */ |

1681 | n/a | word0(&rv) -= P*Exp_msk1; |

1682 | n/a | dval(&rv) *= bigtens[j]; |

1683 | n/a | if ((z = word0(&rv) & Exp_mask) |

1684 | n/a | > Exp_msk1*(DBL_MAX_EXP+Bias-P)) |

1685 | n/a | goto ovfl; |

1686 | n/a | if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) { |

1687 | n/a | /* set to largest number */ |

1688 | n/a | /* (Can't trust DBL_MAX) */ |

1689 | n/a | word0(&rv) = Big0; |

1690 | n/a | word1(&rv) = Big1; |

1691 | n/a | } |

1692 | n/a | else |

1693 | n/a | word0(&rv) += P*Exp_msk1; |

1694 | n/a | } |

1695 | n/a | } |

1696 | n/a | else if (e1 < 0) { |

1697 | n/a | /* The input decimal value lies in [10**e1, 10**(e1+16)). |

1698 | n/a | |

1699 | n/a | If e1 <= -512, underflow immediately. |

1700 | n/a | If e1 <= -256, set bc.scale to 2*P. |

1701 | n/a | |

1702 | n/a | So for input value < 1e-256, bc.scale is always set; |

1703 | n/a | for input value >= 1e-240, bc.scale is never set. |

1704 | n/a | For input values in [1e-256, 1e-240), bc.scale may or may |

1705 | n/a | not be set. */ |

1706 | n/a | |

1707 | n/a | e1 = -e1; |

1708 | n/a | if ((i = e1 & 15)) |

1709 | n/a | dval(&rv) /= tens[i]; |

1710 | n/a | if (e1 >>= 4) { |

1711 | n/a | if (e1 >= 1 << n_bigtens) |

1712 | n/a | goto undfl; |

1713 | n/a | if (e1 & Scale_Bit) |

1714 | n/a | bc.scale = 2*P; |

1715 | n/a | for(j = 0; e1 > 0; j++, e1 >>= 1) |

1716 | n/a | if (e1 & 1) |

1717 | n/a | dval(&rv) *= tinytens[j]; |

1718 | n/a | if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask) |

1719 | n/a | >> Exp_shift)) > 0) { |

1720 | n/a | /* scaled rv is denormal; clear j low bits */ |

1721 | n/a | if (j >= 32) { |

1722 | n/a | word1(&rv) = 0; |

1723 | n/a | if (j >= 53) |

1724 | n/a | word0(&rv) = (P+2)*Exp_msk1; |

1725 | n/a | else |

1726 | n/a | word0(&rv) &= 0xffffffff << (j-32); |

1727 | n/a | } |

1728 | n/a | else |

1729 | n/a | word1(&rv) &= 0xffffffff << j; |

1730 | n/a | } |

1731 | n/a | if (!dval(&rv)) |

1732 | n/a | goto undfl; |

1733 | n/a | } |

1734 | n/a | } |

1735 | n/a | |

1736 | n/a | /* Now the hard part -- adjusting rv to the correct value.*/ |

1737 | n/a | |

1738 | n/a | /* Put digits into bd: true value = bd * 10^e */ |

1739 | n/a | |

1740 | n/a | bc.nd = nd; |

1741 | n/a | bc.nd0 = nd0; /* Only needed if nd > STRTOD_DIGLIM, but done here */ |

1742 | n/a | /* to silence an erroneous warning about bc.nd0 */ |

1743 | n/a | /* possibly not being initialized. */ |

1744 | n/a | if (nd > STRTOD_DIGLIM) { |

1745 | n/a | /* ASSERT(STRTOD_DIGLIM >= 18); 18 == one more than the */ |

1746 | n/a | /* minimum number of decimal digits to distinguish double values */ |

1747 | n/a | /* in IEEE arithmetic. */ |

1748 | n/a | |

1749 | n/a | /* Truncate input to 18 significant digits, then discard any trailing |

1750 | n/a | zeros on the result by updating nd, nd0, e and y suitably. (There's |

1751 | n/a | no need to update z; it's not reused beyond this point.) */ |

1752 | n/a | for (i = 18; i > 0; ) { |

1753 | n/a | /* scan back until we hit a nonzero digit. significant digit 'i' |

1754 | n/a | is s0[i] if i < nd0, s0[i+1] if i >= nd0. */ |

1755 | n/a | --i; |

1756 | n/a | if (s0[i < nd0 ? i : i+1] != '0') { |

1757 | n/a | ++i; |

1758 | n/a | break; |

1759 | n/a | } |

1760 | n/a | } |

1761 | n/a | e += nd - i; |

1762 | n/a | nd = i; |

1763 | n/a | if (nd0 > nd) |

1764 | n/a | nd0 = nd; |

1765 | n/a | if (nd < 9) { /* must recompute y */ |

1766 | n/a | y = 0; |

1767 | n/a | for(i = 0; i < nd0; ++i) |

1768 | n/a | y = 10*y + s0[i] - '0'; |

1769 | n/a | for(; i < nd; ++i) |

1770 | n/a | y = 10*y + s0[i+1] - '0'; |

1771 | n/a | } |

1772 | n/a | } |

1773 | n/a | bd0 = s2b(s0, nd0, nd, y); |

1774 | n/a | if (bd0 == NULL) |

1775 | n/a | goto failed_malloc; |

1776 | n/a | |

1777 | n/a | /* Notation for the comments below. Write: |

1778 | n/a | |

1779 | n/a | - dv for the absolute value of the number represented by the original |

1780 | n/a | decimal input string. |

1781 | n/a | |

1782 | n/a | - if we've truncated dv, write tdv for the truncated value. |

1783 | n/a | Otherwise, set tdv == dv. |

1784 | n/a | |

1785 | n/a | - srv for the quantity rv/2^bc.scale; so srv is the current binary |

1786 | n/a | approximation to tdv (and dv). It should be exactly representable |

1787 | n/a | in an IEEE 754 double. |

1788 | n/a | */ |

1789 | n/a | |

1790 | n/a | for(;;) { |

1791 | n/a | |

1792 | n/a | /* This is the main correction loop for _Py_dg_strtod. |

1793 | n/a | |

1794 | n/a | We've got a decimal value tdv, and a floating-point approximation |

1795 | n/a | srv=rv/2^bc.scale to tdv. The aim is to determine whether srv is |

1796 | n/a | close enough (i.e., within 0.5 ulps) to tdv, and to compute a new |

1797 | n/a | approximation if not. |

1798 | n/a | |

1799 | n/a | To determine whether srv is close enough to tdv, compute integers |

1800 | n/a | bd, bb and bs proportional to tdv, srv and 0.5 ulp(srv) |

1801 | n/a | respectively, and then use integer arithmetic to determine whether |

1802 | n/a | |tdv - srv| is less than, equal to, or greater than 0.5 ulp(srv). |

1803 | n/a | */ |

1804 | n/a | |

1805 | n/a | bd = Balloc(bd0->k); |

1806 | n/a | if (bd == NULL) { |

1807 | n/a | Bfree(bd0); |

1808 | n/a | goto failed_malloc; |

1809 | n/a | } |

1810 | n/a | Bcopy(bd, bd0); |

1811 | n/a | bb = sd2b(&rv, bc.scale, &bbe); /* srv = bb * 2^bbe */ |

1812 | n/a | if (bb == NULL) { |

1813 | n/a | Bfree(bd); |

1814 | n/a | Bfree(bd0); |

1815 | n/a | goto failed_malloc; |

1816 | n/a | } |

1817 | n/a | /* Record whether lsb of bb is odd, in case we need this |

1818 | n/a | for the round-to-even step later. */ |

1819 | n/a | odd = bb->x[0] & 1; |

1820 | n/a | |

1821 | n/a | /* tdv = bd * 10**e; srv = bb * 2**bbe */ |

1822 | n/a | bs = i2b(1); |

1823 | n/a | if (bs == NULL) { |

1824 | n/a | Bfree(bb); |

1825 | n/a | Bfree(bd); |

1826 | n/a | Bfree(bd0); |

1827 | n/a | goto failed_malloc; |

1828 | n/a | } |

1829 | n/a | |

1830 | n/a | if (e >= 0) { |

1831 | n/a | bb2 = bb5 = 0; |

1832 | n/a | bd2 = bd5 = e; |

1833 | n/a | } |

1834 | n/a | else { |

1835 | n/a | bb2 = bb5 = -e; |

1836 | n/a | bd2 = bd5 = 0; |

1837 | n/a | } |

1838 | n/a | if (bbe >= 0) |

1839 | n/a | bb2 += bbe; |

1840 | n/a | else |

1841 | n/a | bd2 -= bbe; |

1842 | n/a | bs2 = bb2; |

1843 | n/a | bb2++; |

1844 | n/a | bd2++; |

1845 | n/a | |

1846 | n/a | /* At this stage bd5 - bb5 == e == bd2 - bb2 + bbe, bb2 - bs2 == 1, |

1847 | n/a | and bs == 1, so: |

1848 | n/a | |

1849 | n/a | tdv == bd * 10**e = bd * 2**(bbe - bb2 + bd2) * 5**(bd5 - bb5) |

1850 | n/a | srv == bb * 2**bbe = bb * 2**(bbe - bb2 + bb2) |

1851 | n/a | 0.5 ulp(srv) == 2**(bbe-1) = bs * 2**(bbe - bb2 + bs2) |

1852 | n/a | |

1853 | n/a | It follows that: |

1854 | n/a | |

1855 | n/a | M * tdv = bd * 2**bd2 * 5**bd5 |

1856 | n/a | M * srv = bb * 2**bb2 * 5**bb5 |

1857 | n/a | M * 0.5 ulp(srv) = bs * 2**bs2 * 5**bb5 |

1858 | n/a | |

1859 | n/a | for some constant M. (Actually, M == 2**(bb2 - bbe) * 5**bb5, but |

1860 | n/a | this fact is not needed below.) |

1861 | n/a | */ |

1862 | n/a | |

1863 | n/a | /* Remove factor of 2**i, where i = min(bb2, bd2, bs2). */ |

1864 | n/a | i = bb2 < bd2 ? bb2 : bd2; |

1865 | n/a | if (i > bs2) |

1866 | n/a | i = bs2; |

1867 | n/a | if (i > 0) { |

1868 | n/a | bb2 -= i; |

1869 | n/a | bd2 -= i; |

1870 | n/a | bs2 -= i; |

1871 | n/a | } |

1872 | n/a | |

1873 | n/a | /* Scale bb, bd, bs by the appropriate powers of 2 and 5. */ |

1874 | n/a | if (bb5 > 0) { |

1875 | n/a | bs = pow5mult(bs, bb5); |

1876 | n/a | if (bs == NULL) { |

1877 | n/a | Bfree(bb); |

1878 | n/a | Bfree(bd); |

1879 | n/a | Bfree(bd0); |

1880 | n/a | goto failed_malloc; |

1881 | n/a | } |

1882 | n/a | bb1 = mult(bs, bb); |

1883 | n/a | Bfree(bb); |

1884 | n/a | bb = bb1; |

1885 | n/a | if (bb == NULL) { |

1886 | n/a | Bfree(bs); |

1887 | n/a | Bfree(bd); |

1888 | n/a | Bfree(bd0); |

1889 | n/a | goto failed_malloc; |

1890 | n/a | } |

1891 | n/a | } |

1892 | n/a | if (bb2 > 0) { |

1893 | n/a | bb = lshift(bb, bb2); |

1894 | n/a | if (bb == NULL) { |

1895 | n/a | Bfree(bs); |

1896 | n/a | Bfree(bd); |

1897 | n/a | Bfree(bd0); |

1898 | n/a | goto failed_malloc; |

1899 | n/a | } |

1900 | n/a | } |

1901 | n/a | if (bd5 > 0) { |

1902 | n/a | bd = pow5mult(bd, bd5); |

1903 | n/a | if (bd == NULL) { |

1904 | n/a | Bfree(bb); |

1905 | n/a | Bfree(bs); |

1906 | n/a | Bfree(bd0); |

1907 | n/a | goto failed_malloc; |

1908 | n/a | } |

1909 | n/a | } |

1910 | n/a | if (bd2 > 0) { |

1911 | n/a | bd = lshift(bd, bd2); |

1912 | n/a | if (bd == NULL) { |

1913 | n/a | Bfree(bb); |

1914 | n/a | Bfree(bs); |

1915 | n/a | Bfree(bd0); |

1916 | n/a | goto failed_malloc; |

1917 | n/a | } |

1918 | n/a | } |

1919 | n/a | if (bs2 > 0) { |

1920 | n/a | bs = lshift(bs, bs2); |

1921 | n/a | if (bs == NULL) { |

1922 | n/a | Bfree(bb); |

1923 | n/a | Bfree(bd); |

1924 | n/a | Bfree(bd0); |

1925 | n/a | goto failed_malloc; |

1926 | n/a | } |

1927 | n/a | } |

1928 | n/a | |

1929 | n/a | /* Now bd, bb and bs are scaled versions of tdv, srv and 0.5 ulp(srv), |

1930 | n/a | respectively. Compute the difference |tdv - srv|, and compare |

1931 | n/a | with 0.5 ulp(srv). */ |

1932 | n/a | |

1933 | n/a | delta = diff(bb, bd); |

1934 | n/a | if (delta == NULL) { |

1935 | n/a | Bfree(bb); |

1936 | n/a | Bfree(bs); |

1937 | n/a | Bfree(bd); |

1938 | n/a | Bfree(bd0); |

1939 | n/a | goto failed_malloc; |

1940 | n/a | } |

1941 | n/a | dsign = delta->sign; |

1942 | n/a | delta->sign = 0; |

1943 | n/a | i = cmp(delta, bs); |

1944 | n/a | if (bc.nd > nd && i <= 0) { |

1945 | n/a | if (dsign) |

1946 | n/a | break; /* Must use bigcomp(). */ |

1947 | n/a | |

1948 | n/a | /* Here rv overestimates the truncated decimal value by at most |

1949 | n/a | 0.5 ulp(rv). Hence rv either overestimates the true decimal |

1950 | n/a | value by <= 0.5 ulp(rv), or underestimates it by some small |

1951 | n/a | amount (< 0.1 ulp(rv)); either way, rv is within 0.5 ulps of |

1952 | n/a | the true decimal value, so it's possible to exit. |

1953 | n/a | |

1954 | n/a | Exception: if scaled rv is a normal exact power of 2, but not |

1955 | n/a | DBL_MIN, then rv - 0.5 ulp(rv) takes us all the way down to the |

1956 | n/a | next double, so the correctly rounded result is either rv - 0.5 |

1957 | n/a | ulp(rv) or rv; in this case, use bigcomp to distinguish. */ |

1958 | n/a | |

1959 | n/a | if (!word1(&rv) && !(word0(&rv) & Bndry_mask)) { |

1960 | n/a | /* rv can't be 0, since it's an overestimate for some |

1961 | n/a | nonzero value. So rv is a normal power of 2. */ |

1962 | n/a | j = (int)(word0(&rv) & Exp_mask) >> Exp_shift; |

1963 | n/a | /* rv / 2^bc.scale = 2^(j - 1023 - bc.scale); use bigcomp if |

1964 | n/a | rv / 2^bc.scale >= 2^-1021. */ |

1965 | n/a | if (j - bc.scale >= 2) { |

1966 | n/a | dval(&rv) -= 0.5 * sulp(&rv, &bc); |

1967 | n/a | break; /* Use bigcomp. */ |

1968 | n/a | } |

1969 | n/a | } |

1970 | n/a | |

1971 | n/a | { |

1972 | n/a | bc.nd = nd; |

1973 | n/a | i = -1; /* Discarded digits make delta smaller. */ |

1974 | n/a | } |

1975 | n/a | } |

1976 | n/a | |

1977 | n/a | if (i < 0) { |

1978 | n/a | /* Error is less than half an ulp -- check for |

1979 | n/a | * special case of mantissa a power of two. |

1980 | n/a | */ |

1981 | n/a | if (dsign || word1(&rv) || word0(&rv) & Bndry_mask |

1982 | n/a | || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1 |

1983 | n/a | ) { |

1984 | n/a | break; |

1985 | n/a | } |

1986 | n/a | if (!delta->x[0] && delta->wds <= 1) { |

1987 | n/a | /* exact result */ |

1988 | n/a | break; |

1989 | n/a | } |

1990 | n/a | delta = lshift(delta,Log2P); |

1991 | n/a | if (delta == NULL) { |

1992 | n/a | Bfree(bb); |

1993 | n/a | Bfree(bs); |

1994 | n/a | Bfree(bd); |

1995 | n/a | Bfree(bd0); |

1996 | n/a | goto failed_malloc; |

1997 | n/a | } |

1998 | n/a | if (cmp(delta, bs) > 0) |

1999 | n/a | goto drop_down; |

2000 | n/a | break; |

2001 | n/a | } |

2002 | n/a | if (i == 0) { |

2003 | n/a | /* exactly half-way between */ |

2004 | n/a | if (dsign) { |

2005 | n/a | if ((word0(&rv) & Bndry_mask1) == Bndry_mask1 |

2006 | n/a | && word1(&rv) == ( |

2007 | n/a | (bc.scale && |

2008 | n/a | (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1) ? |

2009 | n/a | (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : |

2010 | n/a | 0xffffffff)) { |

2011 | n/a | /*boundary case -- increment exponent*/ |

2012 | n/a | word0(&rv) = (word0(&rv) & Exp_mask) |

2013 | n/a | + Exp_msk1 |

2014 | n/a | ; |

2015 | n/a | word1(&rv) = 0; |

2016 | n/a | /* dsign = 0; */ |

2017 | n/a | break; |

2018 | n/a | } |

2019 | n/a | } |

2020 | n/a | else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) { |

2021 | n/a | drop_down: |

2022 | n/a | /* boundary case -- decrement exponent */ |

2023 | n/a | if (bc.scale) { |

2024 | n/a | L = word0(&rv) & Exp_mask; |

2025 | n/a | if (L <= (2*P+1)*Exp_msk1) { |

2026 | n/a | if (L > (P+2)*Exp_msk1) |

2027 | n/a | /* round even ==> */ |

2028 | n/a | /* accept rv */ |

2029 | n/a | break; |

2030 | n/a | /* rv = smallest denormal */ |

2031 | n/a | if (bc.nd > nd) |

2032 | n/a | break; |

2033 | n/a | goto undfl; |

2034 | n/a | } |

2035 | n/a | } |

2036 | n/a | L = (word0(&rv) & Exp_mask) - Exp_msk1; |

2037 | n/a | word0(&rv) = L | Bndry_mask1; |

2038 | n/a | word1(&rv) = 0xffffffff; |

2039 | n/a | break; |

2040 | n/a | } |

2041 | n/a | if (!odd) |

2042 | n/a | break; |

2043 | n/a | if (dsign) |

2044 | n/a | dval(&rv) += sulp(&rv, &bc); |

2045 | n/a | else { |

2046 | n/a | dval(&rv) -= sulp(&rv, &bc); |

2047 | n/a | if (!dval(&rv)) { |

2048 | n/a | if (bc.nd >nd) |

2049 | n/a | break; |

2050 | n/a | goto undfl; |

2051 | n/a | } |

2052 | n/a | } |

2053 | n/a | /* dsign = 1 - dsign; */ |

2054 | n/a | break; |

2055 | n/a | } |

2056 | n/a | if ((aadj = ratio(delta, bs)) <= 2.) { |

2057 | n/a | if (dsign) |

2058 | n/a | aadj = aadj1 = 1.; |

2059 | n/a | else if (word1(&rv) || word0(&rv) & Bndry_mask) { |

2060 | n/a | if (word1(&rv) == Tiny1 && !word0(&rv)) { |

2061 | n/a | if (bc.nd >nd) |

2062 | n/a | break; |

2063 | n/a | goto undfl; |

2064 | n/a | } |

2065 | n/a | aadj = 1.; |

2066 | n/a | aadj1 = -1.; |

2067 | n/a | } |

2068 | n/a | else { |

2069 | n/a | /* special case -- power of FLT_RADIX to be */ |

2070 | n/a | /* rounded down... */ |

2071 | n/a | |

2072 | n/a | if (aadj < 2./FLT_RADIX) |

2073 | n/a | aadj = 1./FLT_RADIX; |

2074 | n/a | else |

2075 | n/a | aadj *= 0.5; |

2076 | n/a | aadj1 = -aadj; |

2077 | n/a | } |

2078 | n/a | } |

2079 | n/a | else { |

2080 | n/a | aadj *= 0.5; |

2081 | n/a | aadj1 = dsign ? aadj : -aadj; |

2082 | n/a | if (Flt_Rounds == 0) |

2083 | n/a | aadj1 += 0.5; |

2084 | n/a | } |

2085 | n/a | y = word0(&rv) & Exp_mask; |

2086 | n/a | |

2087 | n/a | /* Check for overflow */ |

2088 | n/a | |

2089 | n/a | if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { |

2090 | n/a | dval(&rv0) = dval(&rv); |

2091 | n/a | word0(&rv) -= P*Exp_msk1; |

2092 | n/a | adj.d = aadj1 * ulp(&rv); |

2093 | n/a | dval(&rv) += adj.d; |

2094 | n/a | if ((word0(&rv) & Exp_mask) >= |

2095 | n/a | Exp_msk1*(DBL_MAX_EXP+Bias-P)) { |

2096 | n/a | if (word0(&rv0) == Big0 && word1(&rv0) == Big1) { |

2097 | n/a | Bfree(bb); |

2098 | n/a | Bfree(bd); |

2099 | n/a | Bfree(bs); |

2100 | n/a | Bfree(bd0); |

2101 | n/a | Bfree(delta); |

2102 | n/a | goto ovfl; |

2103 | n/a | } |

2104 | n/a | word0(&rv) = Big0; |

2105 | n/a | word1(&rv) = Big1; |

2106 | n/a | goto cont; |

2107 | n/a | } |

2108 | n/a | else |

2109 | n/a | word0(&rv) += P*Exp_msk1; |

2110 | n/a | } |

2111 | n/a | else { |

2112 | n/a | if (bc.scale && y <= 2*P*Exp_msk1) { |

2113 | n/a | if (aadj <= 0x7fffffff) { |

2114 | n/a | if ((z = (ULong)aadj) <= 0) |

2115 | n/a | z = 1; |

2116 | n/a | aadj = z; |

2117 | n/a | aadj1 = dsign ? aadj : -aadj; |

2118 | n/a | } |

2119 | n/a | dval(&aadj2) = aadj1; |

2120 | n/a | word0(&aadj2) += (2*P+1)*Exp_msk1 - y; |

2121 | n/a | aadj1 = dval(&aadj2); |

2122 | n/a | } |

2123 | n/a | adj.d = aadj1 * ulp(&rv); |

2124 | n/a | dval(&rv) += adj.d; |

2125 | n/a | } |

2126 | n/a | z = word0(&rv) & Exp_mask; |

2127 | n/a | if (bc.nd == nd) { |

2128 | n/a | if (!bc.scale) |

2129 | n/a | if (y == z) { |

2130 | n/a | /* Can we stop now? */ |

2131 | n/a | L = (Long)aadj; |

2132 | n/a | aadj -= L; |

2133 | n/a | /* The tolerances below are conservative. */ |

2134 | n/a | if (dsign || word1(&rv) || word0(&rv) & Bndry_mask) { |

2135 | n/a | if (aadj < .4999999 || aadj > .5000001) |

2136 | n/a | break; |

2137 | n/a | } |

2138 | n/a | else if (aadj < .4999999/FLT_RADIX) |

2139 | n/a | break; |

2140 | n/a | } |

2141 | n/a | } |

2142 | n/a | cont: |

2143 | n/a | Bfree(bb); |

2144 | n/a | Bfree(bd); |

2145 | n/a | Bfree(bs); |

2146 | n/a | Bfree(delta); |

2147 | n/a | } |

2148 | n/a | Bfree(bb); |

2149 | n/a | Bfree(bd); |

2150 | n/a | Bfree(bs); |

2151 | n/a | Bfree(bd0); |

2152 | n/a | Bfree(delta); |

2153 | n/a | if (bc.nd > nd) { |

2154 | n/a | error = bigcomp(&rv, s0, &bc); |

2155 | n/a | if (error) |

2156 | n/a | goto failed_malloc; |

2157 | n/a | } |

2158 | n/a | |

2159 | n/a | if (bc.scale) { |

2160 | n/a | word0(&rv0) = Exp_1 - 2*P*Exp_msk1; |

2161 | n/a | word1(&rv0) = 0; |

2162 | n/a | dval(&rv) *= dval(&rv0); |

2163 | n/a | } |

2164 | n/a | |

2165 | n/a | ret: |

2166 | n/a | return sign ? -dval(&rv) : dval(&rv); |

2167 | n/a | |

2168 | n/a | parse_error: |

2169 | n/a | return 0.0; |

2170 | n/a | |

2171 | n/a | failed_malloc: |

2172 | n/a | errno = ENOMEM; |

2173 | n/a | return -1.0; |

2174 | n/a | |

2175 | n/a | undfl: |

2176 | n/a | return sign ? -0.0 : 0.0; |

2177 | n/a | |

2178 | n/a | ovfl: |

2179 | n/a | errno = ERANGE; |

2180 | n/a | /* Can't trust HUGE_VAL */ |

2181 | n/a | word0(&rv) = Exp_mask; |

2182 | n/a | word1(&rv) = 0; |

2183 | n/a | return sign ? -dval(&rv) : dval(&rv); |

2184 | n/a | |

2185 | n/a | } |

2186 | n/a | |

2187 | n/a | static char * |

2188 | n/a | rv_alloc(int i) |

2189 | n/a | { |

2190 | n/a | int j, k, *r; |

2191 | n/a | |

2192 | n/a | j = sizeof(ULong); |

2193 | n/a | for(k = 0; |

2194 | n/a | sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= (unsigned)i; |

2195 | n/a | j <<= 1) |

2196 | n/a | k++; |

2197 | n/a | r = (int*)Balloc(k); |

2198 | n/a | if (r == NULL) |

2199 | n/a | return NULL; |

2200 | n/a | *r = k; |

2201 | n/a | return (char *)(r+1); |

2202 | n/a | } |

2203 | n/a | |

2204 | n/a | static char * |

2205 | n/a | nrv_alloc(const char *s, char **rve, int n) |

2206 | n/a | { |

2207 | n/a | char *rv, *t; |

2208 | n/a | |

2209 | n/a | rv = rv_alloc(n); |

2210 | n/a | if (rv == NULL) |

2211 | n/a | return NULL; |

2212 | n/a | t = rv; |

2213 | n/a | while((*t = *s++)) t++; |

2214 | n/a | if (rve) |

2215 | n/a | *rve = t; |

2216 | n/a | return rv; |

2217 | n/a | } |

2218 | n/a | |

2219 | n/a | /* freedtoa(s) must be used to free values s returned by dtoa |

2220 | n/a | * when MULTIPLE_THREADS is #defined. It should be used in all cases, |

2221 | n/a | * but for consistency with earlier versions of dtoa, it is optional |

2222 | n/a | * when MULTIPLE_THREADS is not defined. |

2223 | n/a | */ |

2224 | n/a | |

2225 | n/a | void |

2226 | n/a | _Py_dg_freedtoa(char *s) |

2227 | n/a | { |

2228 | n/a | Bigint *b = (Bigint *)((int *)s - 1); |

2229 | n/a | b->maxwds = 1 << (b->k = *(int*)b); |

2230 | n/a | Bfree(b); |

2231 | n/a | } |

2232 | n/a | |

2233 | n/a | /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string. |

2234 | n/a | * |

2235 | n/a | * Inspired by "How to Print Floating-Point Numbers Accurately" by |

2236 | n/a | * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126]. |

2237 | n/a | * |

2238 | n/a | * Modifications: |

2239 | n/a | * 1. Rather than iterating, we use a simple numeric overestimate |

2240 | n/a | * to determine k = floor(log10(d)). We scale relevant |

2241 | n/a | * quantities using O(log2(k)) rather than O(k) multiplications. |

2242 | n/a | * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't |

2243 | n/a | * try to generate digits strictly left to right. Instead, we |

2244 | n/a | * compute with fewer bits and propagate the carry if necessary |

2245 | n/a | * when rounding the final digit up. This is often faster. |

2246 | n/a | * 3. Under the assumption that input will be rounded nearest, |

2247 | n/a | * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22. |

2248 | n/a | * That is, we allow equality in stopping tests when the |

2249 | n/a | * round-nearest rule will give the same floating-point value |

2250 | n/a | * as would satisfaction of the stopping test with strict |

2251 | n/a | * inequality. |

2252 | n/a | * 4. We remove common factors of powers of 2 from relevant |

2253 | n/a | * quantities. |

2254 | n/a | * 5. When converting floating-point integers less than 1e16, |

2255 | n/a | * we use floating-point arithmetic rather than resorting |

2256 | n/a | * to multiple-precision integers. |

2257 | n/a | * 6. When asked to produce fewer than 15 digits, we first try |

2258 | n/a | * to get by with floating-point arithmetic; we resort to |

2259 | n/a | * multiple-precision integer arithmetic only if we cannot |

2260 | n/a | * guarantee that the floating-point calculation has given |

2261 | n/a | * the correctly rounded result. For k requested digits and |

2262 | n/a | * "uniformly" distributed input, the probability is |

2263 | n/a | * something like 10^(k-15) that we must resort to the Long |

2264 | n/a | * calculation. |

2265 | n/a | */ |

2266 | n/a | |

2267 | n/a | /* Additional notes (METD): (1) returns NULL on failure. (2) to avoid memory |

2268 | n/a | leakage, a successful call to _Py_dg_dtoa should always be matched by a |

2269 | n/a | call to _Py_dg_freedtoa. */ |

2270 | n/a | |

2271 | n/a | char * |

2272 | n/a | _Py_dg_dtoa(double dd, int mode, int ndigits, |

2273 | n/a | int *decpt, int *sign, char **rve) |

2274 | n/a | { |

2275 | n/a | /* Arguments ndigits, decpt, sign are similar to those |

2276 | n/a | of ecvt and fcvt; trailing zeros are suppressed from |

2277 | n/a | the returned string. If not null, *rve is set to point |

2278 | n/a | to the end of the return value. If d is +-Infinity or NaN, |

2279 | n/a | then *decpt is set to 9999. |

2280 | n/a | |

2281 | n/a | mode: |

2282 | n/a | 0 ==> shortest string that yields d when read in |

2283 | n/a | and rounded to nearest. |

2284 | n/a | 1 ==> like 0, but with Steele & White stopping rule; |

2285 | n/a | e.g. with IEEE P754 arithmetic , mode 0 gives |

2286 | n/a | 1e23 whereas mode 1 gives 9.999999999999999e22. |

2287 | n/a | 2 ==> max(1,ndigits) significant digits. This gives a |

2288 | n/a | return value similar to that of ecvt, except |

2289 | n/a | that trailing zeros are suppressed. |

2290 | n/a | 3 ==> through ndigits past the decimal point. This |

2291 | n/a | gives a return value similar to that from fcvt, |

2292 | n/a | except that trailing zeros are suppressed, and |

2293 | n/a | ndigits can be negative. |

2294 | n/a | 4,5 ==> similar to 2 and 3, respectively, but (in |

2295 | n/a | round-nearest mode) with the tests of mode 0 to |

2296 | n/a | possibly return a shorter string that rounds to d. |

2297 | n/a | With IEEE arithmetic and compilation with |

2298 | n/a | -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same |

2299 | n/a | as modes 2 and 3 when FLT_ROUNDS != 1. |

2300 | n/a | 6-9 ==> Debugging modes similar to mode - 4: don't try |

2301 | n/a | fast floating-point estimate (if applicable). |

2302 | n/a | |

2303 | n/a | Values of mode other than 0-9 are treated as mode 0. |

2304 | n/a | |

2305 | n/a | Sufficient space is allocated to the return value |

2306 | n/a | to hold the suppressed trailing zeros. |

2307 | n/a | */ |

2308 | n/a | |

2309 | n/a | int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1, |

2310 | n/a | j, j1, k, k0, k_check, leftright, m2, m5, s2, s5, |

2311 | n/a | spec_case, try_quick; |

2312 | n/a | Long L; |

2313 | n/a | int denorm; |

2314 | n/a | ULong x; |

2315 | n/a | Bigint *b, *b1, *delta, *mlo, *mhi, *S; |

2316 | n/a | U d2, eps, u; |

2317 | n/a | double ds; |

2318 | n/a | char *s, *s0; |

2319 | n/a | |

2320 | n/a | /* set pointers to NULL, to silence gcc compiler warnings and make |

2321 | n/a | cleanup easier on error */ |

2322 | n/a | mlo = mhi = S = 0; |

2323 | n/a | s0 = 0; |

2324 | n/a | |

2325 | n/a | u.d = dd; |

2326 | n/a | if (word0(&u) & Sign_bit) { |

2327 | n/a | /* set sign for everything, including 0's and NaNs */ |

2328 | n/a | *sign = 1; |

2329 | n/a | word0(&u) &= ~Sign_bit; /* clear sign bit */ |

2330 | n/a | } |

2331 | n/a | else |

2332 | n/a | *sign = 0; |

2333 | n/a | |

2334 | n/a | /* quick return for Infinities, NaNs and zeros */ |

2335 | n/a | if ((word0(&u) & Exp_mask) == Exp_mask) |

2336 | n/a | { |

2337 | n/a | /* Infinity or NaN */ |

2338 | n/a | *decpt = 9999; |

2339 | n/a | if (!word1(&u) && !(word0(&u) & 0xfffff)) |

2340 | n/a | return nrv_alloc("Infinity", rve, 8); |

2341 | n/a | return nrv_alloc("NaN", rve, 3); |

2342 | n/a | } |

2343 | n/a | if (!dval(&u)) { |

2344 | n/a | *decpt = 1; |

2345 | n/a | return nrv_alloc("0", rve, 1); |

2346 | n/a | } |

2347 | n/a | |

2348 | n/a | /* compute k = floor(log10(d)). The computation may leave k |

2349 | n/a | one too large, but should never leave k too small. */ |

2350 | n/a | b = d2b(&u, &be, &bbits); |

2351 | n/a | if (b == NULL) |

2352 | n/a | goto failed_malloc; |

2353 | n/a | if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) { |

2354 | n/a | dval(&d2) = dval(&u); |

2355 | n/a | word0(&d2) &= Frac_mask1; |

2356 | n/a | word0(&d2) |= Exp_11; |

2357 | n/a | |

2358 | n/a | /* log(x) ~=~ log(1.5) + (x-1.5)/1.5 |

2359 | n/a | * log10(x) = log(x) / log(10) |

2360 | n/a | * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10)) |

2361 | n/a | * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2) |

2362 | n/a | * |

2363 | n/a | * This suggests computing an approximation k to log10(d) by |

2364 | n/a | * |

2365 | n/a | * k = (i - Bias)*0.301029995663981 |

2366 | n/a | * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 ); |

2367 | n/a | * |

2368 | n/a | * We want k to be too large rather than too small. |

2369 | n/a | * The error in the first-order Taylor series approximation |

2370 | n/a | * is in our favor, so we just round up the constant enough |

2371 | n/a | * to compensate for any error in the multiplication of |

2372 | n/a | * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077, |

2373 | n/a | * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14, |

2374 | n/a | * adding 1e-13 to the constant term more than suffices. |

2375 | n/a | * Hence we adjust the constant term to 0.1760912590558. |

2376 | n/a | * (We could get a more accurate k by invoking log10, |

2377 | n/a | * but this is probably not worthwhile.) |

2378 | n/a | */ |

2379 | n/a | |

2380 | n/a | i -= Bias; |

2381 | n/a | denorm = 0; |

2382 | n/a | } |

2383 | n/a | else { |

2384 | n/a | /* d is denormalized */ |

2385 | n/a | |

2386 | n/a | i = bbits + be + (Bias + (P-1) - 1); |

2387 | n/a | x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32) |

2388 | n/a | : word1(&u) << (32 - i); |

2389 | n/a | dval(&d2) = x; |

2390 | n/a | word0(&d2) -= 31*Exp_msk1; /* adjust exponent */ |

2391 | n/a | i -= (Bias + (P-1) - 1) + 1; |

2392 | n/a | denorm = 1; |

2393 | n/a | } |

2394 | n/a | ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + |

2395 | n/a | i*0.301029995663981; |

2396 | n/a | k = (int)ds; |

2397 | n/a | if (ds < 0. && ds != k) |

2398 | n/a | k--; /* want k = floor(ds) */ |

2399 | n/a | k_check = 1; |

2400 | n/a | if (k >= 0 && k <= Ten_pmax) { |

2401 | n/a | if (dval(&u) < tens[k]) |

2402 | n/a | k--; |

2403 | n/a | k_check = 0; |

2404 | n/a | } |

2405 | n/a | j = bbits - i - 1; |

2406 | n/a | if (j >= 0) { |

2407 | n/a | b2 = 0; |

2408 | n/a | s2 = j; |

2409 | n/a | } |

2410 | n/a | else { |

2411 | n/a | b2 = -j; |

2412 | n/a | s2 = 0; |

2413 | n/a | } |

2414 | n/a | if (k >= 0) { |

2415 | n/a | b5 = 0; |

2416 | n/a | s5 = k; |

2417 | n/a | s2 += k; |

2418 | n/a | } |

2419 | n/a | else { |

2420 | n/a | b2 -= k; |

2421 | n/a | b5 = -k; |

2422 | n/a | s5 = 0; |

2423 | n/a | } |

2424 | n/a | if (mode < 0 || mode > 9) |

2425 | n/a | mode = 0; |

2426 | n/a | |

2427 | n/a | try_quick = 1; |

2428 | n/a | |

2429 | n/a | if (mode > 5) { |

2430 | n/a | mode -= 4; |

2431 | n/a | try_quick = 0; |

2432 | n/a | } |

2433 | n/a | leftright = 1; |

2434 | n/a | ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */ |

2435 | n/a | /* silence erroneous "gcc -Wall" warning. */ |

2436 | n/a | switch(mode) { |

2437 | n/a | case 0: |

2438 | n/a | case 1: |

2439 | n/a | i = 18; |

2440 | n/a | ndigits = 0; |

2441 | n/a | break; |

2442 | n/a | case 2: |

2443 | n/a | leftright = 0; |

2444 | n/a | /* no break */ |

2445 | n/a | case 4: |

2446 | n/a | if (ndigits <= 0) |

2447 | n/a | ndigits = 1; |

2448 | n/a | ilim = ilim1 = i = ndigits; |

2449 | n/a | break; |

2450 | n/a | case 3: |

2451 | n/a | leftright = 0; |

2452 | n/a | /* no break */ |

2453 | n/a | case 5: |

2454 | n/a | i = ndigits + k + 1; |

2455 | n/a | ilim = i; |

2456 | n/a | ilim1 = i - 1; |

2457 | n/a | if (i <= 0) |

2458 | n/a | i = 1; |

2459 | n/a | } |

2460 | n/a | s0 = rv_alloc(i); |

2461 | n/a | if (s0 == NULL) |

2462 | n/a | goto failed_malloc; |

2463 | n/a | s = s0; |

2464 | n/a | |

2465 | n/a | |

2466 | n/a | if (ilim >= 0 && ilim <= Quick_max && try_quick) { |

2467 | n/a | |

2468 | n/a | /* Try to get by with floating-point arithmetic. */ |

2469 | n/a | |

2470 | n/a | i = 0; |

2471 | n/a | dval(&d2) = dval(&u); |

2472 | n/a | k0 = k; |

2473 | n/a | ilim0 = ilim; |

2474 | n/a | ieps = 2; /* conservative */ |

2475 | n/a | if (k > 0) { |

2476 | n/a | ds = tens[k&0xf]; |

2477 | n/a | j = k >> 4; |

2478 | n/a | if (j & Bletch) { |

2479 | n/a | /* prevent overflows */ |

2480 | n/a | j &= Bletch - 1; |

2481 | n/a | dval(&u) /= bigtens[n_bigtens-1]; |

2482 | n/a | ieps++; |

2483 | n/a | } |

2484 | n/a | for(; j; j >>= 1, i++) |

2485 | n/a | if (j & 1) { |

2486 | n/a | ieps++; |

2487 | n/a | ds *= bigtens[i]; |

2488 | n/a | } |

2489 | n/a | dval(&u) /= ds; |

2490 | n/a | } |

2491 | n/a | else if ((j1 = -k)) { |

2492 | n/a | dval(&u) *= tens[j1 & 0xf]; |

2493 | n/a | for(j = j1 >> 4; j; j >>= 1, i++) |

2494 | n/a | if (j & 1) { |

2495 | n/a | ieps++; |

2496 | n/a | dval(&u) *= bigtens[i]; |

2497 | n/a | } |

2498 | n/a | } |

2499 | n/a | if (k_check && dval(&u) < 1. && ilim > 0) { |

2500 | n/a | if (ilim1 <= 0) |

2501 | n/a | goto fast_failed; |

2502 | n/a | ilim = ilim1; |

2503 | n/a | k--; |

2504 | n/a | dval(&u) *= 10.; |

2505 | n/a | ieps++; |

2506 | n/a | } |

2507 | n/a | dval(&eps) = ieps*dval(&u) + 7.; |

2508 | n/a | word0(&eps) -= (P-1)*Exp_msk1; |

2509 | n/a | if (ilim == 0) { |

2510 | n/a | S = mhi = 0; |

2511 | n/a | dval(&u) -= 5.; |

2512 | n/a | if (dval(&u) > dval(&eps)) |

2513 | n/a | goto one_digit; |

2514 | n/a | if (dval(&u) < -dval(&eps)) |

2515 | n/a | goto no_digits; |

2516 | n/a | goto fast_failed; |

2517 | n/a | } |

2518 | n/a | if (leftright) { |

2519 | n/a | /* Use Steele & White method of only |

2520 | n/a | * generating digits needed. |

2521 | n/a | */ |

2522 | n/a | dval(&eps) = 0.5/tens[ilim-1] - dval(&eps); |

2523 | n/a | for(i = 0;;) { |

2524 | n/a | L = (Long)dval(&u); |

2525 | n/a | dval(&u) -= L; |

2526 | n/a | *s++ = '0' + (int)L; |

2527 | n/a | if (dval(&u) < dval(&eps)) |

2528 | n/a | goto ret1; |

2529 | n/a | if (1. - dval(&u) < dval(&eps)) |

2530 | n/a | goto bump_up; |

2531 | n/a | if (++i >= ilim) |

2532 | n/a | break; |

2533 | n/a | dval(&eps) *= 10.; |

2534 | n/a | dval(&u) *= 10.; |

2535 | n/a | } |

2536 | n/a | } |

2537 | n/a | else { |

2538 | n/a | /* Generate ilim digits, then fix them up. */ |

2539 | n/a | dval(&eps) *= tens[ilim-1]; |

2540 | n/a | for(i = 1;; i++, dval(&u) *= 10.) { |

2541 | n/a | L = (Long)(dval(&u)); |

2542 | n/a | if (!(dval(&u) -= L)) |

2543 | n/a | ilim = i; |

2544 | n/a | *s++ = '0' + (int)L; |

2545 | n/a | if (i == ilim) { |

2546 | n/a | if (dval(&u) > 0.5 + dval(&eps)) |

2547 | n/a | goto bump_up; |

2548 | n/a | else if (dval(&u) < 0.5 - dval(&eps)) { |

2549 | n/a | while(*--s == '0'); |

2550 | n/a | s++; |

2551 | n/a | goto ret1; |

2552 | n/a | } |

2553 | n/a | break; |

2554 | n/a | } |

2555 | n/a | } |

2556 | n/a | } |

2557 | n/a | fast_failed: |

2558 | n/a | s = s0; |

2559 | n/a | dval(&u) = dval(&d2); |

2560 | n/a | k = k0; |

2561 | n/a | ilim = ilim0; |

2562 | n/a | } |

2563 | n/a | |

2564 | n/a | /* Do we have a "small" integer? */ |

2565 | n/a | |

2566 | n/a | if (be >= 0 && k <= Int_max) { |

2567 | n/a | /* Yes. */ |

2568 | n/a | ds = tens[k]; |

2569 | n/a | if (ndigits < 0 && ilim <= 0) { |

2570 | n/a | S = mhi = 0; |

2571 | n/a | if (ilim < 0 || dval(&u) <= 5*ds) |

2572 | n/a | goto no_digits; |

2573 | n/a | goto one_digit; |

2574 | n/a | } |

2575 | n/a | for(i = 1;; i++, dval(&u) *= 10.) { |

2576 | n/a | L = (Long)(dval(&u) / ds); |

2577 | n/a | dval(&u) -= L*ds; |

2578 | n/a | *s++ = '0' + (int)L; |

2579 | n/a | if (!dval(&u)) { |

2580 | n/a | break; |

2581 | n/a | } |

2582 | n/a | if (i == ilim) { |

2583 | n/a | dval(&u) += dval(&u); |

2584 | n/a | if (dval(&u) > ds || (dval(&u) == ds && L & 1)) { |

2585 | n/a | bump_up: |

2586 | n/a | while(*--s == '9') |

2587 | n/a | if (s == s0) { |

2588 | n/a | k++; |

2589 | n/a | *s = '0'; |

2590 | n/a | break; |

2591 | n/a | } |

2592 | n/a | ++*s++; |

2593 | n/a | } |

2594 | n/a | break; |

2595 | n/a | } |

2596 | n/a | } |

2597 | n/a | goto ret1; |

2598 | n/a | } |

2599 | n/a | |

2600 | n/a | m2 = b2; |

2601 | n/a | m5 = b5; |

2602 | n/a | if (leftright) { |

2603 | n/a | i = |

2604 | n/a | denorm ? be + (Bias + (P-1) - 1 + 1) : |

2605 | n/a | 1 + P - bbits; |

2606 | n/a | b2 += i; |

2607 | n/a | s2 += i; |

2608 | n/a | mhi = i2b(1); |

2609 | n/a | if (mhi == NULL) |

2610 | n/a | goto failed_malloc; |

2611 | n/a | } |

2612 | n/a | if (m2 > 0 && s2 > 0) { |

2613 | n/a | i = m2 < s2 ? m2 : s2; |

2614 | n/a | b2 -= i; |

2615 | n/a | m2 -= i; |

2616 | n/a | s2 -= i; |

2617 | n/a | } |

2618 | n/a | if (b5 > 0) { |

2619 | n/a | if (leftright) { |

2620 | n/a | if (m5 > 0) { |

2621 | n/a | mhi = pow5mult(mhi, m5); |

2622 | n/a | if (mhi == NULL) |

2623 | n/a | goto failed_malloc; |

2624 | n/a | b1 = mult(mhi, b); |

2625 | n/a | Bfree(b); |

2626 | n/a | b = b1; |

2627 | n/a | if (b == NULL) |

2628 | n/a | goto failed_malloc; |

2629 | n/a | } |

2630 | n/a | if ((j = b5 - m5)) { |

2631 | n/a | b = pow5mult(b, j); |

2632 | n/a | if (b == NULL) |

2633 | n/a | goto failed_malloc; |

2634 | n/a | } |

2635 | n/a | } |

2636 | n/a | else { |

2637 | n/a | b = pow5mult(b, b5); |

2638 | n/a | if (b == NULL) |

2639 | n/a | goto failed_malloc; |

2640 | n/a | } |

2641 | n/a | } |

2642 | n/a | S = i2b(1); |

2643 | n/a | if (S == NULL) |

2644 | n/a | goto failed_malloc; |

2645 | n/a | if (s5 > 0) { |

2646 | n/a | S = pow5mult(S, s5); |

2647 | n/a | if (S == NULL) |

2648 | n/a | goto failed_malloc; |

2649 | n/a | } |

2650 | n/a | |

2651 | n/a | /* Check for special case that d is a normalized power of 2. */ |

2652 | n/a | |

2653 | n/a | spec_case = 0; |

2654 | n/a | if ((mode < 2 || leftright) |

2655 | n/a | ) { |

2656 | n/a | if (!word1(&u) && !(word0(&u) & Bndry_mask) |

2657 | n/a | && word0(&u) & (Exp_mask & ~Exp_msk1) |

2658 | n/a | ) { |

2659 | n/a | /* The special case */ |

2660 | n/a | b2 += Log2P; |

2661 | n/a | s2 += Log2P; |

2662 | n/a | spec_case = 1; |

2663 | n/a | } |

2664 | n/a | } |

2665 | n/a | |

2666 | n/a | /* Arrange for convenient computation of quotients: |

2667 | n/a | * shift left if necessary so divisor has 4 leading 0 bits. |

2668 | n/a | * |

2669 | n/a | * Perhaps we should just compute leading 28 bits of S once |

2670 | n/a | * and for all and pass them and a shift to quorem, so it |

2671 | n/a | * can do shifts and ors to compute the numerator for q. |

2672 | n/a | */ |

2673 | n/a | #define iInc 28 |

2674 | n/a | i = dshift(S, s2); |

2675 | n/a | b2 += i; |

2676 | n/a | m2 += i; |

2677 | n/a | s2 += i; |

2678 | n/a | if (b2 > 0) { |

2679 | n/a | b = lshift(b, b2); |

2680 | n/a | if (b == NULL) |

2681 | n/a | goto failed_malloc; |

2682 | n/a | } |

2683 | n/a | if (s2 > 0) { |

2684 | n/a | S = lshift(S, s2); |

2685 | n/a | if (S == NULL) |

2686 | n/a | goto failed_malloc; |

2687 | n/a | } |

2688 | n/a | if (k_check) { |

2689 | n/a | if (cmp(b,S) < 0) { |

2690 | n/a | k--; |

2691 | n/a | b = multadd(b, 10, 0); /* we botched the k estimate */ |

2692 | n/a | if (b == NULL) |

2693 | n/a | goto failed_malloc; |

2694 | n/a | if (leftright) { |

2695 | n/a | mhi = multadd(mhi, 10, 0); |

2696 | n/a | if (mhi == NULL) |

2697 | n/a | goto failed_malloc; |

2698 | n/a | } |

2699 | n/a | ilim = ilim1; |

2700 | n/a | } |

2701 | n/a | } |

2702 | n/a | if (ilim <= 0 && (mode == 3 || mode == 5)) { |

2703 | n/a | if (ilim < 0) { |

2704 | n/a | /* no digits, fcvt style */ |

2705 | n/a | no_digits: |

2706 | n/a | k = -1 - ndigits; |

2707 | n/a | goto ret; |

2708 | n/a | } |

2709 | n/a | else { |

2710 | n/a | S = multadd(S, 5, 0); |

2711 | n/a | if (S == NULL) |

2712 | n/a | goto failed_malloc; |

2713 | n/a | if (cmp(b, S) <= 0) |

2714 | n/a | goto no_digits; |

2715 | n/a | } |

2716 | n/a | one_digit: |

2717 | n/a | *s++ = '1'; |

2718 | n/a | k++; |

2719 | n/a | goto ret; |

2720 | n/a | } |

2721 | n/a | if (leftright) { |

2722 | n/a | if (m2 > 0) { |

2723 | n/a | mhi = lshift(mhi, m2); |

2724 | n/a | if (mhi == NULL) |

2725 | n/a | goto failed_malloc; |

2726 | n/a | } |

2727 | n/a | |

2728 | n/a | /* Compute mlo -- check for special case |

2729 | n/a | * that d is a normalized power of 2. |

2730 | n/a | */ |

2731 | n/a | |

2732 | n/a | mlo = mhi; |

2733 | n/a | if (spec_case) { |

2734 | n/a | mhi = Balloc(mhi->k); |

2735 | n/a | if (mhi == NULL) |

2736 | n/a | goto failed_malloc; |

2737 | n/a | Bcopy(mhi, mlo); |

2738 | n/a | mhi = lshift(mhi, Log2P); |

2739 | n/a | if (mhi == NULL) |

2740 | n/a | goto failed_malloc; |

2741 | n/a | } |

2742 | n/a | |

2743 | n/a | for(i = 1;;i++) { |

2744 | n/a | dig = quorem(b,S) + '0'; |

2745 | n/a | /* Do we yet have the shortest decimal string |

2746 | n/a | * that will round to d? |

2747 | n/a | */ |

2748 | n/a | j = cmp(b, mlo); |

2749 | n/a | delta = diff(S, mhi); |

2750 | n/a | if (delta == NULL) |

2751 | n/a | goto failed_malloc; |

2752 | n/a | j1 = delta->sign ? 1 : cmp(b, delta); |

2753 | n/a | Bfree(delta); |

2754 | n/a | if (j1 == 0 && mode != 1 && !(word1(&u) & 1) |

2755 | n/a | ) { |

2756 | n/a | if (dig == '9') |

2757 | n/a | goto round_9_up; |

2758 | n/a | if (j > 0) |

2759 | n/a | dig++; |

2760 | n/a | *s++ = dig; |

2761 | n/a | goto ret; |

2762 | n/a | } |

2763 | n/a | if (j < 0 || (j == 0 && mode != 1 |

2764 | n/a | && !(word1(&u) & 1) |

2765 | n/a | )) { |

2766 | n/a | if (!b->x[0] && b->wds <= 1) { |

2767 | n/a | goto accept_dig; |

2768 | n/a | } |

2769 | n/a | if (j1 > 0) { |

2770 | n/a | b = lshift(b, 1); |

2771 | n/a | if (b == NULL) |

2772 | n/a | goto failed_malloc; |

2773 | n/a | j1 = cmp(b, S); |

2774 | n/a | if ((j1 > 0 || (j1 == 0 && dig & 1)) |

2775 | n/a | && dig++ == '9') |

2776 | n/a | goto round_9_up; |

2777 | n/a | } |

2778 | n/a | accept_dig: |

2779 | n/a | *s++ = dig; |

2780 | n/a | goto ret; |

2781 | n/a | } |

2782 | n/a | if (j1 > 0) { |

2783 | n/a | if (dig == '9') { /* possible if i == 1 */ |

2784 | n/a | round_9_up: |

2785 | n/a | *s++ = '9'; |

2786 | n/a | goto roundoff; |

2787 | n/a | } |

2788 | n/a | *s++ = dig + 1; |

2789 | n/a | goto ret; |

2790 | n/a | } |

2791 | n/a | *s++ = dig; |

2792 | n/a | if (i == ilim) |

2793 | n/a | break; |

2794 | n/a | b = multadd(b, 10, 0); |

2795 | n/a | if (b == NULL) |

2796 | n/a | goto failed_malloc; |

2797 | n/a | if (mlo == mhi) { |

2798 | n/a | mlo = mhi = multadd(mhi, 10, 0); |

2799 | n/a | if (mlo == NULL) |

2800 | n/a | goto failed_malloc; |

2801 | n/a | } |

2802 | n/a | else { |

2803 | n/a | mlo = multadd(mlo, 10, 0); |

2804 | n/a | if (mlo == NULL) |

2805 | n/a | goto failed_malloc; |

2806 | n/a | mhi = multadd(mhi, 10, 0); |

2807 | n/a | if (mhi == NULL) |

2808 | n/a | goto failed_malloc; |

2809 | n/a | } |

2810 | n/a | } |

2811 | n/a | } |

2812 | n/a | else |

2813 | n/a | for(i = 1;; i++) { |

2814 | n/a | *s++ = dig = quorem(b,S) + '0'; |

2815 | n/a | if (!b->x[0] && b->wds <= 1) { |

2816 | n/a | goto ret; |

2817 | n/a | } |

2818 | n/a | if (i >= ilim) |

2819 | n/a | break; |

2820 | n/a | b = multadd(b, 10, 0); |

2821 | n/a | if (b == NULL) |

2822 | n/a | goto failed_malloc; |

2823 | n/a | } |

2824 | n/a | |

2825 | n/a | /* Round off last digit */ |

2826 | n/a | |

2827 | n/a | b = lshift(b, 1); |

2828 | n/a | if (b == NULL) |

2829 | n/a | goto failed_malloc; |

2830 | n/a | j = cmp(b, S); |

2831 | n/a | if (j > 0 || (j == 0 && dig & 1)) { |

2832 | n/a | roundoff: |

2833 | n/a | while(*--s == '9') |

2834 | n/a | if (s == s0) { |

2835 | n/a | k++; |

2836 | n/a | *s++ = '1'; |

2837 | n/a | goto ret; |

2838 | n/a | } |

2839 | n/a | ++*s++; |

2840 | n/a | } |

2841 | n/a | else { |

2842 | n/a | while(*--s == '0'); |

2843 | n/a | s++; |

2844 | n/a | } |

2845 | n/a | ret: |

2846 | n/a | Bfree(S); |

2847 | n/a | if (mhi) { |

2848 | n/a | if (mlo && mlo != mhi) |

2849 | n/a | Bfree(mlo); |

2850 | n/a | Bfree(mhi); |

2851 | n/a | } |

2852 | n/a | ret1: |

2853 | n/a | Bfree(b); |

2854 | n/a | *s = 0; |

2855 | n/a | *decpt = k + 1; |

2856 | n/a | if (rve) |

2857 | n/a | *rve = s; |

2858 | n/a | return s0; |

2859 | n/a | failed_malloc: |

2860 | n/a | if (S) |

2861 | n/a | Bfree(S); |

2862 | n/a | if (mlo && mlo != mhi) |

2863 | n/a | Bfree(mlo); |

2864 | n/a | if (mhi) |

2865 | n/a | Bfree(mhi); |

2866 | n/a | if (b) |

2867 | n/a | Bfree(b); |

2868 | n/a | if (s0) |

2869 | n/a | _Py_dg_freedtoa(s0); |

2870 | n/a | return NULL; |

2871 | n/a | } |

2872 | n/a | #ifdef __cplusplus |

2873 | n/a | } |

2874 | n/a | #endif |

2875 | n/a | |

2876 | n/a | #endif /* PY_NO_SHORT_FLOAT_REPR */ |