1 | /* Convert string representing a number to float value, using given locale. |
---|---|

2 | Copyright (C) 1997-2012 Free Software Foundation, Inc. |

3 | This file is part of the GNU C Library. |

4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. |

5 | |

6 | The GNU C Library is free software; you can redistribute it and/or |

7 | modify it under the terms of the GNU Lesser General Public |

8 | License as published by the Free Software Foundation; either |

9 | version 2.1 of the License, or (at your option) any later version. |

10 | |

11 | The GNU C Library is distributed in the hope that it will be useful, |

12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

14 | Lesser General Public License for more details. |

15 | |

16 | You should have received a copy of the GNU Lesser General Public |

17 | License along with the GNU C Library; if not, see |

18 | <http://www.gnu.org/licenses/>. */ |

19 | |

20 | #include <config.h> |

21 | #include <stdarg.h> |

22 | #include <string.h> |

23 | #include <stdint.h> |

24 | #include <stdbool.h> |

25 | #include <float.h> |

26 | #include <math.h> |

27 | #define NDEBUG 1 |

28 | #include <assert.h> |

29 | #ifdef HAVE_ERRNO_H |

30 | #include <errno.h> |

31 | #endif |

32 | |

33 | #ifdef HAVE_FENV_H |

34 | #include <fenv.h> |

35 | #endif |

36 | |

37 | #ifdef HAVE_FENV_H |

38 | #include "quadmath-rounding-mode.h" |

39 | #endif |

40 | #include "../printf/quadmath-printf.h" |

41 | #include "../printf/fpioconst.h" |

42 | |

43 | #undef L_ |

44 | #ifdef USE_WIDE_CHAR |

45 | # define STRING_TYPE wchar_t |

46 | # define CHAR_TYPE wint_t |

47 | # define L_(Ch) L##Ch |

48 | # define ISSPACE(Ch) __iswspace_l ((Ch), loc) |

49 | # define ISDIGIT(Ch) __iswdigit_l ((Ch), loc) |

50 | # define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc) |

51 | # define TOLOWER(Ch) __towlower_l ((Ch), loc) |

52 | # define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr) |

53 | # define STRNCASECMP(S1, S2, N) \ |

54 | __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr) |

55 | # define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc) |

56 | #else |

57 | # define STRING_TYPE char |

58 | # define CHAR_TYPE char |

59 | # define L_(Ch) Ch |

60 | # define ISSPACE(Ch) isspace (Ch) |

61 | # define ISDIGIT(Ch) isdigit (Ch) |

62 | # define ISXDIGIT(Ch) isxdigit (Ch) |

63 | # define TOLOWER(Ch) tolower (Ch) |

64 | # define TOLOWER_C(Ch) \ |

65 | ({__typeof(Ch) __tlc = (Ch); \ |

66 | (__tlc >= 'A' && __tlc <= 'Z') ? __tlc - 'A' + 'a' : __tlc; }) |

67 | # define STRNCASECMP(S1, S2, N) \ |

68 | __quadmath_strncasecmp_c (S1, S2, N) |

69 | # ifdef HAVE_STRTOULL |

70 | # define STRTOULL(S, E, B) strtoull (S, E, B) |

71 | # else |

72 | # define STRTOULL(S, E, B) strtoul (S, E, B) |

73 | # endif |

74 | |

75 | static inline int |

76 | __quadmath_strncasecmp_c (const char *s1, const char *s2, size_t n) |

77 | { |

78 | const unsigned char *p1 = (const unsigned char *) s1; |

79 | const unsigned char *p2 = (const unsigned char *) s2; |

80 | int result; |

81 | if (p1 == p2 || n == 0) |

82 | return 0; |

83 | while ((result = TOLOWER_C (*p1) - TOLOWER_C (*p2++)) == 0) |

84 | if (*p1++ == '\0' || --n == 0) |

85 | break; |

86 | |

87 | return result; |

88 | } |

89 | #endif |

90 | |

91 | |

92 | /* Constants we need from float.h; select the set for the FLOAT precision. */ |

93 | #define MANT_DIG PASTE(FLT,_MANT_DIG) |

94 | #define DIG PASTE(FLT,_DIG) |

95 | #define MAX_EXP PASTE(FLT,_MAX_EXP) |

96 | #define MIN_EXP PASTE(FLT,_MIN_EXP) |

97 | #define MAX_10_EXP PASTE(FLT,_MAX_10_EXP) |

98 | #define MIN_10_EXP PASTE(FLT,_MIN_10_EXP) |

99 | #define MAX_VALUE PASTE(FLT,_MAX) |

100 | #define MIN_VALUE PASTE(FLT,_MIN) |

101 | |

102 | /* Extra macros required to get FLT expanded before the pasting. */ |

103 | #define PASTE(a,b) PASTE1(a,b) |

104 | #define PASTE1(a,b) a##b |

105 | |

106 | /* Function to construct a floating point number from an MP integer |

107 | containing the fraction bits, a base 2 exponent, and a sign flag. */ |

108 | extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative); |

109 | |

110 | /* Definitions according to limb size used. */ |

111 | #if BITS_PER_MP_LIMB == 32 |

112 | # define MAX_DIG_PER_LIMB 9 |

113 | # define MAX_FAC_PER_LIMB 1000000000UL |

114 | #elif BITS_PER_MP_LIMB == 64 |

115 | # define MAX_DIG_PER_LIMB 19 |

116 | # define MAX_FAC_PER_LIMB 10000000000000000000ULL |

117 | #else |

118 | # error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for" |

119 | #endif |

120 | |

121 | #define _tens_in_limb __quadmath_tens_in_limb |

122 | extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] attribute_hidden; |

123 | |

124 | #ifndef howmany |

125 | #define howmany(x,y) (((x)+((y)-1))/(y)) |

126 | #endif |

127 | #define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; }) |

128 | |

129 | #define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG) |

130 | #define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG) |

131 | #define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB) |

132 | |

133 | #define RETURN(val,end) \ |

134 | do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \ |

135 | return val; } while (0) |

136 | |

137 | /* Maximum size necessary for mpn integers to hold floating point |

138 | numbers. The largest number we need to hold is 10^n where 2^-n is |

139 | 1/4 ulp of the smallest representable value (that is, n = MANT_DIG |

140 | - MIN_EXP + 2). Approximate using 10^3 < 2^10. */ |

141 | #define MPNSIZE (howmany (1 + ((MANT_DIG - MIN_EXP + 2) * 10) / 3, \ |

142 | BITS_PER_MP_LIMB) + 2) |

143 | /* Declare an mpn integer variable that big. */ |

144 | #define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size |

145 | /* Copy an mpn integer value. */ |

146 | #define MPN_ASSIGN(dst, src) \ |

147 | memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t)) |

148 | |

149 | /* Set errno and return an overflowing value with sign specified by |

150 | NEGATIVE. */ |

151 | static FLOAT |

152 | overflow_value (int negative) |

153 | { |

154 | #if defined HAVE_ERRNO_H && defined ERANGE |

155 | errno = ERANGE; |

156 | #endif |

157 | FLOAT result = (negative ? -MAX_VALUE : MAX_VALUE) * MAX_VALUE; |

158 | return result; |

159 | } |

160 | |

161 | /* Set errno and return an underflowing value with sign specified by |

162 | NEGATIVE. */ |

163 | static FLOAT |

164 | underflow_value (int negative) |

165 | { |

166 | #if defined HAVE_ERRNO_H && defined ERANGE |

167 | errno = ERANGE; |

168 | #endif |

169 | FLOAT result = (negative ? -MIN_VALUE : MIN_VALUE) * MIN_VALUE; |

170 | return result; |

171 | } |

172 | |

173 | /* Return a floating point number of the needed type according to the given |

174 | multi-precision number after possible rounding. */ |

175 | static FLOAT |

176 | round_and_return (mp_limb_t *retval, intmax_t exponent, int negative, |

177 | mp_limb_t round_limb, mp_size_t round_bit, int more_bits) |

178 | { |

179 | #ifdef HAVE_FENV_H |

180 | int mode = get_rounding_mode (); |

181 | #endif |

182 | |

183 | if (exponent < MIN_EXP - 1) |

184 | { |

185 | mp_size_t shift; |

186 | bool is_tiny; |

187 | |

188 | if (exponent < MIN_EXP - 1 - MANT_DIG) |

189 | return underflow_value (negative); |

190 | |

191 | shift = MIN_EXP - 1 - exponent; |

192 | is_tiny = true; |

193 | |

194 | more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0; |

195 | if (shift == MANT_DIG) |

196 | /* This is a special case to handle the very seldom case where |

197 | the mantissa will be empty after the shift. */ |

198 | { |

199 | int i; |

200 | |

201 | round_limb = retval[RETURN_LIMB_SIZE - 1]; |

202 | round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; |

203 | for (i = 0; i < RETURN_LIMB_SIZE; ++i) |

204 | more_bits |= retval[i] != 0; |

205 | MPN_ZERO (retval, RETURN_LIMB_SIZE); |

206 | } |

207 | else if (shift >= BITS_PER_MP_LIMB) |

208 | { |

209 | int i; |

210 | |

211 | round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB]; |

212 | round_bit = (shift - 1) % BITS_PER_MP_LIMB; |

213 | for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i) |

214 | more_bits |= retval[i] != 0; |

215 | more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) |

216 | != 0); |

217 | |

218 | (void) mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB], |

219 | RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB), |

220 | shift % BITS_PER_MP_LIMB); |

221 | MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)], |

222 | shift / BITS_PER_MP_LIMB); |

223 | } |

224 | else if (shift > 0) |

225 | { |

226 | #ifdef HAVE_FENV_H |

227 | if (TININESS_AFTER_ROUNDING && shift == 1) |

228 | { |

229 | /* Whether the result counts as tiny depends on whether, |

230 | after rounding to the normal precision, it still has |

231 | a subnormal exponent. */ |

232 | mp_limb_t retval_normal[RETURN_LIMB_SIZE]; |

233 | if (round_away (negative, |

234 | (retval[0] & 1) != 0, |

235 | (round_limb |

236 | & (((mp_limb_t) 1) << round_bit)) != 0, |

237 | (more_bits |

238 | || ((round_limb |

239 | & ((((mp_limb_t) 1) << round_bit) - 1)) |

240 | != 0)), |

241 | mode)) |

242 | { |

243 | mp_limb_t cy = mpn_add_1 (retval_normal, retval, |

244 | RETURN_LIMB_SIZE, 1); |

245 | |

246 | if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) || |

247 | ((MANT_DIG % BITS_PER_MP_LIMB) != 0 && |

248 | ((retval_normal[RETURN_LIMB_SIZE - 1] |

249 | & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) |

250 | != 0))) |

251 | is_tiny = false; |

252 | } |

253 | } |

254 | #endif |

255 | round_limb = retval[0]; |

256 | round_bit = shift - 1; |

257 | (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift); |

258 | } |

259 | /* This is a hook for the m68k long double format, where the |

260 | exponent bias is the same for normalized and denormalized |

261 | numbers. */ |

262 | #ifndef DENORM_EXP |

263 | # define DENORM_EXP (MIN_EXP - 2) |

264 | #endif |

265 | exponent = DENORM_EXP; |

266 | if (is_tiny |

267 | && ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0 |

268 | || more_bits |

269 | || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0)) |

270 | { |

271 | #if defined HAVE_ERRNO_H && defined ERANGE |

272 | errno = ERANGE; |

273 | #endif |

274 | volatile FLOAT force_underflow_exception = MIN_VALUE * MIN_VALUE; |

275 | (void) force_underflow_exception; |

276 | } |

277 | } |

278 | |

279 | if (exponent > MAX_EXP) |

280 | goto overflow; |

281 | |

282 | #ifdef HAVE_FENV_H |

283 | if (round_away (negative, |

284 | (retval[0] & 1) != 0, |

285 | (round_limb & (((mp_limb_t) 1) << round_bit)) != 0, |

286 | (more_bits |

287 | || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0), |

288 | mode)) |

289 | { |

290 | mp_limb_t cy = mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1); |

291 | |

292 | if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) || |

293 | ((MANT_DIG % BITS_PER_MP_LIMB) != 0 && |

294 | (retval[RETURN_LIMB_SIZE - 1] |

295 | & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0)) |

296 | { |

297 | ++exponent; |

298 | (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1); |

299 | retval[RETURN_LIMB_SIZE - 1] |

300 | |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB); |

301 | } |

302 | else if (exponent == DENORM_EXP |

303 | && (retval[RETURN_LIMB_SIZE - 1] |

304 | & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB))) |

305 | != 0) |

306 | /* The number was denormalized but now normalized. */ |

307 | exponent = MIN_EXP - 1; |

308 | } |

309 | #endif |

310 | |

311 | if (exponent > MAX_EXP) |

312 | overflow: |

313 | return overflow_value (negative); |

314 | |

315 | return MPN2FLOAT (retval, exponent, negative); |

316 | } |

317 | |

318 | |

319 | /* Read a multi-precision integer starting at STR with exactly DIGCNT digits |

320 | into N. Return the size of the number limbs in NSIZE at the first |

321 | character od the string that is not part of the integer as the function |

322 | value. If the EXPONENT is small enough to be taken as an additional |

323 | factor for the resulting number (see code) multiply by it. */ |

324 | static const STRING_TYPE * |

325 | str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize, |

326 | intmax_t *exponent |

327 | #ifndef USE_WIDE_CHAR |

328 | , const char *decimal, size_t decimal_len, const char *thousands |

329 | #endif |

330 | |

331 | ) |

332 | { |

333 | /* Number of digits for actual limb. */ |

334 | int cnt = 0; |

335 | mp_limb_t low = 0; |

336 | mp_limb_t start; |

337 | |

338 | *nsize = 0; |

339 | assert (digcnt > 0); |

340 | do |

341 | { |

342 | if (cnt == MAX_DIG_PER_LIMB) |

343 | { |

344 | if (*nsize == 0) |

345 | { |

346 | n[0] = low; |

347 | *nsize = 1; |

348 | } |

349 | else |

350 | { |

351 | mp_limb_t cy; |

352 | cy = mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB); |

353 | cy += mpn_add_1 (n, n, *nsize, low); |

354 | if (cy != 0) |

355 | { |

356 | assert (*nsize < MPNSIZE); |

357 | n[*nsize] = cy; |

358 | ++(*nsize); |

359 | } |

360 | } |

361 | cnt = 0; |

362 | low = 0; |

363 | } |

364 | |

365 | /* There might be thousands separators or radix characters in |

366 | the string. But these all can be ignored because we know the |

367 | format of the number is correct and we have an exact number |

368 | of characters to read. */ |

369 | #ifdef USE_WIDE_CHAR |

370 | if (*str < L'0' || *str > L'9') |

371 | ++str; |

372 | #else |

373 | if (*str < '0' || *str > '9') |

374 | { |

375 | int inner = 0; |

376 | if (thousands != NULL && *str == *thousands |

377 | && ({ for (inner = 1; thousands[inner] != '\0'; ++inner) |

378 | if (thousands[inner] != str[inner]) |

379 | break; |

380 | thousands[inner] == '\0'; })) |

381 | str += inner; |

382 | else |

383 | str += decimal_len; |

384 | } |

385 | #endif |

386 | low = low * 10 + *str++ - L_('0'); |

387 | ++cnt; |

388 | } |

389 | while (--digcnt > 0); |

390 | |

391 | if (*exponent > 0 && *exponent <= MAX_DIG_PER_LIMB - cnt) |

392 | { |

393 | low *= _tens_in_limb[*exponent]; |

394 | start = _tens_in_limb[cnt + *exponent]; |

395 | *exponent = 0; |

396 | } |

397 | else |

398 | start = _tens_in_limb[cnt]; |

399 | |

400 | if (*nsize == 0) |

401 | { |

402 | n[0] = low; |

403 | *nsize = 1; |

404 | } |

405 | else |

406 | { |

407 | mp_limb_t cy; |

408 | cy = mpn_mul_1 (n, n, *nsize, start); |

409 | cy += mpn_add_1 (n, n, *nsize, low); |

410 | if (cy != 0) |

411 | { |

412 | assert (*nsize < MPNSIZE); |

413 | n[(*nsize)++] = cy; |

414 | } |

415 | } |

416 | |

417 | return str; |

418 | } |

419 | |

420 | |

421 | /* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits |

422 | with the COUNT most significant bits of LIMB. |

423 | |

424 | Implemented as a macro, so that __builtin_constant_p works even at -O0. |

425 | |

426 | Tege doesn't like this macro so I have to write it here myself. :) |

427 | --drepper */ |

428 | #define mpn_lshift_1(ptr, size, count, limb) \ |

429 | do \ |

430 | { \ |

431 | mp_limb_t *__ptr = (ptr); \ |

432 | if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB) \ |

433 | { \ |

434 | mp_size_t i; \ |

435 | for (i = (size) - 1; i > 0; --i) \ |

436 | __ptr[i] = __ptr[i - 1]; \ |

437 | __ptr[0] = (limb); \ |

438 | } \ |

439 | else \ |

440 | { \ |

441 | /* We assume count > 0 && count < BITS_PER_MP_LIMB here. */ \ |

442 | unsigned int __count = (count); \ |

443 | (void) mpn_lshift (__ptr, __ptr, size, __count); \ |

444 | __ptr[0] |= (limb) >> (BITS_PER_MP_LIMB - __count); \ |

445 | } \ |

446 | } \ |

447 | while (0) |

448 | |

449 | |

450 | #define INTERNAL(x) INTERNAL1(x) |

451 | #define INTERNAL1(x) __##x##_internal |

452 | #ifndef ____STRTOF_INTERNAL |

453 | # define ____STRTOF_INTERNAL INTERNAL (__STRTOF) |

454 | #endif |

455 | |

456 | /* This file defines a function to check for correct grouping. */ |

457 | #include "grouping.h" |

458 | |

459 | |

460 | /* Return a floating point number with the value of the given string NPTR. |

461 | Set *ENDPTR to the character after the last used one. If the number is |

462 | smaller than the smallest representable number, set `errno' to ERANGE and |

463 | return 0.0. If the number is too big to be represented, set `errno' to |

464 | ERANGE and return HUGE_VAL with the appropriate sign. */ |

465 | FLOAT |

466 | ____STRTOF_INTERNAL (nptr, endptr, group) |

467 | const STRING_TYPE *nptr; |

468 | STRING_TYPE **endptr; |

469 | int group; |

470 | { |

471 | int negative; /* The sign of the number. */ |

472 | MPN_VAR (num); /* MP representation of the number. */ |

473 | intmax_t exponent; /* Exponent of the number. */ |

474 | |

475 | /* Numbers starting `0X' or `0x' have to be processed with base 16. */ |

476 | int base = 10; |

477 | |

478 | /* When we have to compute fractional digits we form a fraction with a |

479 | second multi-precision number (and we sometimes need a second for |

480 | temporary results). */ |

481 | MPN_VAR (den); |

482 | |

483 | /* Representation for the return value. */ |

484 | mp_limb_t retval[RETURN_LIMB_SIZE]; |

485 | /* Number of bits currently in result value. */ |

486 | int bits; |

487 | |

488 | /* Running pointer after the last character processed in the string. */ |

489 | const STRING_TYPE *cp, *tp; |

490 | /* Start of significant part of the number. */ |

491 | const STRING_TYPE *startp, *start_of_digits; |

492 | /* Points at the character following the integer and fractional digits. */ |

493 | const STRING_TYPE *expp; |

494 | /* Total number of digit and number of digits in integer part. */ |

495 | size_t dig_no, int_no, lead_zero; |

496 | /* Contains the last character read. */ |

497 | CHAR_TYPE c; |

498 | |

499 | /* We should get wint_t from <stddef.h>, but not all GCC versions define it |

500 | there. So define it ourselves if it remains undefined. */ |

501 | #ifndef _WINT_T |

502 | typedef unsigned int wint_t; |

503 | #endif |

504 | /* The radix character of the current locale. */ |

505 | #ifdef USE_WIDE_CHAR |

506 | wchar_t decimal; |

507 | #else |

508 | const char *decimal; |

509 | size_t decimal_len; |

510 | #endif |

511 | /* The thousands character of the current locale. */ |

512 | #ifdef USE_WIDE_CHAR |

513 | wchar_t thousands = L'\0'; |

514 | #else |

515 | const char *thousands = NULL; |

516 | #endif |

517 | /* The numeric grouping specification of the current locale, |

518 | in the format described in <locale.h>. */ |

519 | const char *grouping; |

520 | /* Used in several places. */ |

521 | int cnt; |

522 | |

523 | #if defined USE_LOCALECONV && !defined USE_NL_LANGINFO |

524 | const struct lconv *lc = localeconv (); |

525 | #endif |

526 | |

527 | if (__builtin_expect (group, 0)) |

528 | { |

529 | #ifdef USE_NL_LANGINFO |

530 | grouping = nl_langinfo (GROUPING); |

531 | if (*grouping <= 0 || *grouping == CHAR_MAX) |

532 | grouping = NULL; |

533 | else |

534 | { |

535 | /* Figure out the thousands separator character. */ |

536 | #ifdef USE_WIDE_CHAR |

537 | thousands = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC); |

538 | if (thousands == L'\0') |

539 | grouping = NULL; |

540 | #else |

541 | thousands = nl_langinfo (THOUSANDS_SEP); |

542 | if (*thousands == '\0') |

543 | { |

544 | thousands = NULL; |

545 | grouping = NULL; |

546 | } |

547 | #endif |

548 | } |

549 | #elif defined USE_LOCALECONV |

550 | grouping = lc->grouping; |

551 | if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX) |

552 | grouping = NULL; |

553 | else |

554 | { |

555 | /* Figure out the thousands separator character. */ |

556 | thousands = lc->thousands_sep; |

557 | if (thousands == NULL || *thousands == '\0') |

558 | { |

559 | thousands = NULL; |

560 | grouping = NULL; |

561 | } |

562 | } |

563 | #else |

564 | grouping = NULL; |

565 | #endif |

566 | } |

567 | else |

568 | grouping = NULL; |

569 | |

570 | /* Find the locale's decimal point character. */ |

571 | #ifdef USE_WIDE_CHAR |

572 | decimal = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC); |

573 | assert (decimal != L'\0'); |

574 | # define decimal_len 1 |

575 | #else |

576 | #ifdef USE_NL_LANGINFO |

577 | decimal = nl_langinfo (DECIMAL_POINT); |

578 | decimal_len = strlen (decimal); |

579 | assert (decimal_len > 0); |

580 | #elif defined USE_LOCALECONV |

581 | decimal = lc->decimal_point; |

582 | if (decimal == NULL || *decimal == '\0') |

583 | decimal = "."; |

584 | decimal_len = strlen (decimal); |

585 | #else |

586 | decimal = "."; |

587 | decimal_len = 1; |

588 | #endif |

589 | #endif |

590 | |

591 | /* Prepare number representation. */ |

592 | exponent = 0; |

593 | negative = 0; |

594 | bits = 0; |

595 | |

596 | /* Parse string to get maximal legal prefix. We need the number of |

597 | characters of the integer part, the fractional part and the exponent. */ |

598 | cp = nptr - 1; |

599 | /* Ignore leading white space. */ |

600 | do |

601 | c = *++cp; |

602 | while (ISSPACE (c)); |

603 | |

604 | /* Get sign of the result. */ |

605 | if (c == L_('-')) |

606 | { |

607 | negative = 1; |

608 | c = *++cp; |

609 | } |

610 | else if (c == L_('+')) |

611 | c = *++cp; |

612 | |

613 | /* Return 0.0 if no legal string is found. |

614 | No character is used even if a sign was found. */ |

615 | #ifdef USE_WIDE_CHAR |

616 | if (c == (wint_t) decimal |

617 | && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9') |

618 | { |

619 | /* We accept it. This funny construct is here only to indent |

620 | the code correctly. */ |

621 | } |

622 | #else |

623 | for (cnt = 0; decimal[cnt] != '\0'; ++cnt) |

624 | if (cp[cnt] != decimal[cnt]) |

625 | break; |

626 | if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9') |

627 | { |

628 | /* We accept it. This funny construct is here only to indent |

629 | the code correctly. */ |

630 | } |

631 | #endif |

632 | else if (c < L_('0') || c > L_('9')) |

633 | { |

634 | /* Check for `INF' or `INFINITY'. */ |

635 | CHAR_TYPE lowc = TOLOWER_C (c); |

636 | |

637 | if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0) |

638 | { |

639 | /* Return +/- infinity. */ |

640 | if (endptr != NULL) |

641 | *endptr = (STRING_TYPE *) |

642 | (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0 |

643 | ? 8 : 3)); |

644 | |

645 | return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL; |

646 | } |

647 | |

648 | if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0) |

649 | { |

650 | /* Return NaN. */ |

651 | FLOAT retval = NAN; |

652 | |

653 | cp += 3; |

654 | |

655 | /* Match `(n-char-sequence-digit)'. */ |

656 | if (*cp == L_('(')) |

657 | { |

658 | const STRING_TYPE *startp = cp; |

659 | do |

660 | ++cp; |

661 | while ((*cp >= L_('0') && *cp <= L_('9')) |

662 | || ({ CHAR_TYPE lo = TOLOWER (*cp); |

663 | lo >= L_('a') && lo <= L_('z'); }) |

664 | || *cp == L_('_')); |

665 | |

666 | if (*cp != L_(')')) |

667 | /* The closing brace is missing. Only match the NAN |

668 | part. */ |

669 | cp = startp; |

670 | else |

671 | { |

672 | /* This is a system-dependent way to specify the |

673 | bitmask used for the NaN. We expect it to be |

674 | a number which is put in the mantissa of the |

675 | number. */ |

676 | STRING_TYPE *endp; |

677 | unsigned long long int mant; |

678 | |

679 | mant = STRTOULL (startp + 1, &endp, 0); |

680 | if (endp == cp) |

681 | SET_MANTISSA (retval, mant); |

682 | |

683 | /* Consume the closing brace. */ |

684 | ++cp; |

685 | } |

686 | } |

687 | |

688 | if (endptr != NULL) |

689 | *endptr = (STRING_TYPE *) cp; |

690 | |

691 | return retval; |

692 | } |

693 | |

694 | /* It is really a text we do not recognize. */ |

695 | RETURN (0.0, nptr); |

696 | } |

697 | |

698 | /* First look whether we are faced with a hexadecimal number. */ |

699 | if (c == L_('0') && TOLOWER (cp[1]) == L_('x')) |

700 | { |

701 | /* Okay, it is a hexa-decimal number. Remember this and skip |

702 | the characters. BTW: hexadecimal numbers must not be |

703 | grouped. */ |

704 | base = 16; |

705 | cp += 2; |

706 | c = *cp; |

707 | grouping = NULL; |

708 | } |

709 | |

710 | /* Record the start of the digits, in case we will check their grouping. */ |

711 | start_of_digits = startp = cp; |

712 | |

713 | /* Ignore leading zeroes. This helps us to avoid useless computations. */ |

714 | #ifdef USE_WIDE_CHAR |

715 | while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands)) |

716 | c = *++cp; |

717 | #else |

718 | if (__builtin_expect (thousands == NULL, 1)) |

719 | while (c == '0') |

720 | c = *++cp; |

721 | else |

722 | { |

723 | /* We also have the multibyte thousands string. */ |

724 | while (1) |

725 | { |

726 | if (c != '0') |

727 | { |

728 | for (cnt = 0; thousands[cnt] != '\0'; ++cnt) |

729 | if (thousands[cnt] != cp[cnt]) |

730 | break; |

731 | if (thousands[cnt] != '\0') |

732 | break; |

733 | cp += cnt - 1; |

734 | } |

735 | c = *++cp; |

736 | } |

737 | } |

738 | #endif |

739 | |

740 | /* If no other digit but a '0' is found the result is 0.0. |

741 | Return current read pointer. */ |

742 | CHAR_TYPE lowc = TOLOWER (c); |

743 | if (!((c >= L_('0') && c <= L_('9')) |

744 | || (base == 16 && lowc >= L_('a') && lowc <= L_('f')) |

745 | || ( |

746 | #ifdef USE_WIDE_CHAR |

747 | c == (wint_t) decimal |

748 | #else |

749 | ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) |

750 | if (decimal[cnt] != cp[cnt]) |

751 | break; |

752 | decimal[cnt] == '\0'; }) |

753 | #endif |

754 | /* '0x.' alone is not a valid hexadecimal number. |

755 | '.' alone is not valid either, but that has been checked |

756 | already earlier. */ |

757 | && (base != 16 |

758 | || cp != start_of_digits |

759 | || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9')) |

760 | || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]); |

761 | lo >= L_('a') && lo <= L_('f'); }))) |

762 | || (base == 16 && (cp != start_of_digits |

763 | && lowc == L_('p'))) |

764 | || (base != 16 && lowc == L_('e')))) |

765 | { |

766 | #ifdef USE_WIDE_CHAR |

767 | tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands, |

768 | grouping); |

769 | #else |

770 | tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands, |

771 | grouping); |

772 | #endif |

773 | /* If TP is at the start of the digits, there was no correctly |

774 | grouped prefix of the string; so no number found. */ |

775 | RETURN (negative ? -0.0 : 0.0, |

776 | tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp); |

777 | } |

778 | |

779 | /* Remember first significant digit and read following characters until the |

780 | decimal point, exponent character or any non-FP number character. */ |

781 | startp = cp; |

782 | dig_no = 0; |

783 | while (1) |

784 | { |

785 | if ((c >= L_('0') && c <= L_('9')) |

786 | || (base == 16 |

787 | && ({ CHAR_TYPE lo = TOLOWER (c); |

788 | lo >= L_('a') && lo <= L_('f'); }))) |

789 | ++dig_no; |

790 | else |

791 | { |

792 | #ifdef USE_WIDE_CHAR |

793 | if (__builtin_expect ((wint_t) thousands == L'\0', 1) |

794 | || c != (wint_t) thousands) |

795 | /* Not a digit or separator: end of the integer part. */ |

796 | break; |

797 | #else |

798 | if (__builtin_expect (thousands == NULL, 1)) |

799 | break; |

800 | else |

801 | { |

802 | for (cnt = 0; thousands[cnt] != '\0'; ++cnt) |

803 | if (thousands[cnt] != cp[cnt]) |

804 | break; |

805 | if (thousands[cnt] != '\0') |

806 | break; |

807 | cp += cnt - 1; |

808 | } |

809 | #endif |

810 | } |

811 | c = *++cp; |

812 | } |

813 | |

814 | if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits) |

815 | { |

816 | /* Check the grouping of the digits. */ |

817 | #ifdef USE_WIDE_CHAR |

818 | tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands, |

819 | grouping); |

820 | #else |

821 | tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands, |

822 | grouping); |

823 | #endif |

824 | if (cp != tp) |

825 | { |

826 | /* Less than the entire string was correctly grouped. */ |

827 | |

828 | if (tp == start_of_digits) |

829 | /* No valid group of numbers at all: no valid number. */ |

830 | RETURN (0.0, nptr); |

831 | |

832 | if (tp < startp) |

833 | /* The number is validly grouped, but consists |

834 | only of zeroes. The whole value is zero. */ |

835 | RETURN (negative ? -0.0 : 0.0, tp); |

836 | |

837 | /* Recompute DIG_NO so we won't read more digits than |

838 | are properly grouped. */ |

839 | cp = tp; |

840 | dig_no = 0; |

841 | for (tp = startp; tp < cp; ++tp) |

842 | if (*tp >= L_('0') && *tp <= L_('9')) |

843 | ++dig_no; |

844 | |

845 | int_no = dig_no; |

846 | lead_zero = 0; |

847 | |

848 | goto number_parsed; |

849 | } |

850 | } |

851 | |

852 | /* We have the number of digits in the integer part. Whether these |

853 | are all or any is really a fractional digit will be decided |

854 | later. */ |

855 | int_no = dig_no; |

856 | lead_zero = int_no == 0 ? (size_t) -1 : 0; |

857 | |

858 | /* Read the fractional digits. A special case are the 'american |

859 | style' numbers like `16.' i.e. with decimal point but without |

860 | trailing digits. */ |

861 | if ( |

862 | #ifdef USE_WIDE_CHAR |

863 | c == (wint_t) decimal |

864 | #else |

865 | ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt) |

866 | if (decimal[cnt] != cp[cnt]) |

867 | break; |

868 | decimal[cnt] == '\0'; }) |

869 | #endif |

870 | ) |

871 | { |

872 | cp += decimal_len; |

873 | c = *cp; |

874 | while ((c >= L_('0') && c <= L_('9')) || |

875 | (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c); |

876 | lo >= L_('a') && lo <= L_('f'); }))) |

877 | { |

878 | if (c != L_('0') && lead_zero == (size_t) -1) |

879 | lead_zero = dig_no - int_no; |

880 | ++dig_no; |

881 | c = *++cp; |

882 | } |

883 | } |

884 | assert (dig_no <= (uintmax_t) INTMAX_MAX); |

885 | |

886 | /* Remember start of exponent (if any). */ |

887 | expp = cp; |

888 | |

889 | /* Read exponent. */ |

890 | lowc = TOLOWER (c); |

891 | if ((base == 16 && lowc == L_('p')) |

892 | || (base != 16 && lowc == L_('e'))) |

893 | { |

894 | int exp_negative = 0; |

895 | |

896 | c = *++cp; |

897 | if (c == L_('-')) |

898 | { |

899 | exp_negative = 1; |

900 | c = *++cp; |

901 | } |

902 | else if (c == L_('+')) |

903 | c = *++cp; |

904 | |

905 | if (c >= L_('0') && c <= L_('9')) |

906 | { |

907 | intmax_t exp_limit; |

908 | |

909 | /* Get the exponent limit. */ |

910 | if (base == 16) |

911 | { |

912 | if (exp_negative) |

913 | { |

914 | assert (int_no <= (uintmax_t) (INTMAX_MAX |

915 | + MIN_EXP - MANT_DIG) / 4); |

916 | exp_limit = -MIN_EXP + MANT_DIG + 4 * (intmax_t) int_no; |

917 | } |

918 | else |

919 | { |

920 | if (int_no) |

921 | { |

922 | assert (lead_zero == 0 |

923 | && int_no <= (uintmax_t) INTMAX_MAX / 4); |

924 | exp_limit = MAX_EXP - 4 * (intmax_t) int_no + 3; |

925 | } |

926 | else if (lead_zero == (size_t) -1) |

927 | { |

928 | /* The number is zero and this limit is |

929 | arbitrary. */ |

930 | exp_limit = MAX_EXP + 3; |

931 | } |

932 | else |

933 | { |

934 | assert (lead_zero |

935 | <= (uintmax_t) (INTMAX_MAX - MAX_EXP - 3) / 4); |

936 | exp_limit = (MAX_EXP |

937 | + 4 * (intmax_t) lead_zero |

938 | + 3); |

939 | } |

940 | } |

941 | } |

942 | else |

943 | { |

944 | if (exp_negative) |

945 | { |

946 | assert (int_no |

947 | <= (uintmax_t) (INTMAX_MAX + MIN_10_EXP - MANT_DIG)); |

948 | exp_limit = -MIN_10_EXP + MANT_DIG + (intmax_t) int_no; |

949 | } |

950 | else |

951 | { |

952 | if (int_no) |

953 | { |

954 | assert (lead_zero == 0 |

955 | && int_no <= (uintmax_t) INTMAX_MAX); |

956 | exp_limit = MAX_10_EXP - (intmax_t) int_no + 1; |

957 | } |

958 | else if (lead_zero == (size_t) -1) |

959 | { |

960 | /* The number is zero and this limit is |

961 | arbitrary. */ |

962 | exp_limit = MAX_10_EXP + 1; |

963 | } |

964 | else |

965 | { |

966 | assert (lead_zero |

967 | <= (uintmax_t) (INTMAX_MAX - MAX_10_EXP - 1)); |

968 | exp_limit = MAX_10_EXP + (intmax_t) lead_zero + 1; |

969 | } |

970 | } |

971 | } |

972 | |

973 | if (exp_limit < 0) |

974 | exp_limit = 0; |

975 | |

976 | do |

977 | { |

978 | if (__builtin_expect ((exponent > exp_limit / 10 |

979 | || (exponent == exp_limit / 10 |

980 | && c - L_('0') > exp_limit % 10)), 0)) |

981 | /* The exponent is too large/small to represent a valid |

982 | number. */ |

983 | { |

984 | FLOAT result; |

985 | |

986 | /* We have to take care for special situation: a joker |

987 | might have written "0.0e100000" which is in fact |

988 | zero. */ |

989 | if (lead_zero == (size_t) -1) |

990 | result = negative ? -0.0 : 0.0; |

991 | else |

992 | { |

993 | /* Overflow or underflow. */ |

994 | #if defined HAVE_ERRNO_H && defined ERANGE |

995 | errno = ERANGE; |

996 | #endif |

997 | result = (exp_negative ? (negative ? -0.0 : 0.0) : |

998 | negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL); |

999 | } |

1000 | |

1001 | /* Accept all following digits as part of the exponent. */ |

1002 | do |

1003 | ++cp; |

1004 | while (*cp >= L_('0') && *cp <= L_('9')); |

1005 | |

1006 | RETURN (result, cp); |

1007 | /* NOTREACHED */ |

1008 | } |

1009 | |

1010 | exponent *= 10; |

1011 | exponent += c - L_('0'); |

1012 | |

1013 | c = *++cp; |

1014 | } |

1015 | while (c >= L_('0') && c <= L_('9')); |

1016 | |

1017 | if (exp_negative) |

1018 | exponent = -exponent; |

1019 | } |

1020 | else |

1021 | cp = expp; |

1022 | } |

1023 | |

1024 | /* We don't want to have to work with trailing zeroes after the radix. */ |

1025 | if (dig_no > int_no) |

1026 | { |

1027 | while (expp[-1] == L_('0')) |

1028 | { |

1029 | --expp; |

1030 | --dig_no; |

1031 | } |

1032 | assert (dig_no >= int_no); |

1033 | } |

1034 | |

1035 | if (dig_no == int_no && dig_no > 0 && exponent < 0) |

1036 | do |

1037 | { |

1038 | while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1]))) |

1039 | --expp; |

1040 | |

1041 | if (expp[-1] != L_('0')) |

1042 | break; |

1043 | |

1044 | --expp; |

1045 | --dig_no; |

1046 | --int_no; |

1047 | exponent += base == 16 ? 4 : 1; |

1048 | } |

1049 | while (dig_no > 0 && exponent < 0); |

1050 | |

1051 | number_parsed: |

1052 | |

1053 | /* The whole string is parsed. Store the address of the next character. */ |

1054 | if (endptr) |

1055 | *endptr = (STRING_TYPE *) cp; |

1056 | |

1057 | if (dig_no == 0) |

1058 | return negative ? -0.0 : 0.0; |

1059 | |

1060 | if (lead_zero) |

1061 | { |

1062 | /* Find the decimal point */ |

1063 | #ifdef USE_WIDE_CHAR |

1064 | while (*startp != decimal) |

1065 | ++startp; |

1066 | #else |

1067 | while (1) |

1068 | { |

1069 | if (*startp == decimal[0]) |

1070 | { |

1071 | for (cnt = 1; decimal[cnt] != '\0'; ++cnt) |

1072 | if (decimal[cnt] != startp[cnt]) |

1073 | break; |

1074 | if (decimal[cnt] == '\0') |

1075 | break; |

1076 | } |

1077 | ++startp; |

1078 | } |

1079 | #endif |

1080 | startp += lead_zero + decimal_len; |

1081 | assert (lead_zero <= (base == 16 |

1082 | ? (uintmax_t) INTMAX_MAX / 4 |

1083 | : (uintmax_t) INTMAX_MAX)); |

1084 | assert (lead_zero <= (base == 16 |

1085 | ? ((uintmax_t) exponent |

1086 | - (uintmax_t) INTMAX_MIN) / 4 |

1087 | : ((uintmax_t) exponent - (uintmax_t) INTMAX_MIN))); |

1088 | exponent -= base == 16 ? 4 * (intmax_t) lead_zero : (intmax_t) lead_zero; |

1089 | dig_no -= lead_zero; |

1090 | } |

1091 | |

1092 | /* If the BASE is 16 we can use a simpler algorithm. */ |

1093 | if (base == 16) |

1094 | { |

1095 | static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3, |

1096 | 4, 4, 4, 4, 4, 4, 4, 4 }; |

1097 | int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB; |

1098 | int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB; |

1099 | mp_limb_t val; |

1100 | |

1101 | while (!ISXDIGIT (*startp)) |

1102 | ++startp; |

1103 | while (*startp == L_('0')) |

1104 | ++startp; |

1105 | if (ISDIGIT (*startp)) |

1106 | val = *startp++ - L_('0'); |

1107 | else |

1108 | val = 10 + TOLOWER (*startp++) - L_('a'); |

1109 | bits = nbits[val]; |

1110 | /* We cannot have a leading zero. */ |

1111 | assert (bits != 0); |

1112 | |

1113 | if (pos + 1 >= 4 || pos + 1 >= bits) |

1114 | { |

1115 | /* We don't have to care for wrapping. This is the normal |

1116 | case so we add the first clause in the `if' expression as |

1117 | an optimization. It is a compile-time constant and so does |

1118 | not cost anything. */ |

1119 | retval[idx] = val << (pos - bits + 1); |

1120 | pos -= bits; |

1121 | } |

1122 | else |

1123 | { |

1124 | retval[idx--] = val >> (bits - pos - 1); |

1125 | retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1)); |

1126 | pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1); |

1127 | } |

1128 | |

1129 | /* Adjust the exponent for the bits we are shifting in. */ |

1130 | assert (int_no <= (uintmax_t) (exponent < 0 |

1131 | ? (INTMAX_MAX - bits + 1) / 4 |

1132 | : (INTMAX_MAX - exponent - bits + 1) / 4)); |

1133 | exponent += bits - 1 + ((intmax_t) int_no - 1) * 4; |

1134 | |

1135 | while (--dig_no > 0 && idx >= 0) |

1136 | { |

1137 | if (!ISXDIGIT (*startp)) |

1138 | startp += decimal_len; |

1139 | if (ISDIGIT (*startp)) |

1140 | val = *startp++ - L_('0'); |

1141 | else |

1142 | val = 10 + TOLOWER (*startp++) - L_('a'); |

1143 | |

1144 | if (pos + 1 >= 4) |

1145 | { |

1146 | retval[idx] |= val << (pos - 4 + 1); |

1147 | pos -= 4; |

1148 | } |

1149 | else |

1150 | { |

1151 | retval[idx--] |= val >> (4 - pos - 1); |

1152 | val <<= BITS_PER_MP_LIMB - (4 - pos - 1); |

1153 | if (idx < 0) |

1154 | { |

1155 | int rest_nonzero = 0; |

1156 | while (--dig_no > 0) |

1157 | { |

1158 | if (*startp != L_('0')) |

1159 | { |

1160 | rest_nonzero = 1; |

1161 | break; |

1162 | } |

1163 | startp++; |

1164 | } |

1165 | return round_and_return (retval, exponent, negative, val, |

1166 | BITS_PER_MP_LIMB - 1, rest_nonzero); |

1167 | } |

1168 | |

1169 | retval[idx] = val; |

1170 | pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1); |

1171 | } |

1172 | } |

1173 | |

1174 | /* We ran out of digits. */ |

1175 | MPN_ZERO (retval, idx); |

1176 | |

1177 | return round_and_return (retval, exponent, negative, 0, 0, 0); |

1178 | } |

1179 | |

1180 | /* Now we have the number of digits in total and the integer digits as well |

1181 | as the exponent and its sign. We can decide whether the read digits are |

1182 | really integer digits or belong to the fractional part; i.e. we normalize |

1183 | 123e-2 to 1.23. */ |

1184 | { |

1185 | register intmax_t incr = (exponent < 0 |

1186 | ? MAX (-(intmax_t) int_no, exponent) |

1187 | : MIN ((intmax_t) dig_no - (intmax_t) int_no, |

1188 | exponent)); |

1189 | int_no += incr; |

1190 | exponent -= incr; |

1191 | } |

1192 | |

1193 | if (__builtin_expect (exponent > MAX_10_EXP + 1 - (intmax_t) int_no, 0)) |

1194 | return overflow_value (negative); |

1195 | |

1196 | if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0)) |

1197 | return underflow_value (negative); |

1198 | |

1199 | if (int_no > 0) |

1200 | { |

1201 | /* Read the integer part as a multi-precision number to NUM. */ |

1202 | startp = str_to_mpn (startp, int_no, num, &numsize, &exponent |

1203 | #ifndef USE_WIDE_CHAR |

1204 | , decimal, decimal_len, thousands |

1205 | #endif |

1206 | ); |

1207 | |

1208 | if (exponent > 0) |

1209 | { |

1210 | /* We now multiply the gained number by the given power of ten. */ |

1211 | mp_limb_t *psrc = num; |

1212 | mp_limb_t *pdest = den; |

1213 | int expbit = 1; |

1214 | const struct mp_power *ttab = &_fpioconst_pow10[0]; |

1215 | |

1216 | do |

1217 | { |

1218 | if ((exponent & expbit) != 0) |

1219 | { |

1220 | size_t size = ttab->arraysize - _FPIO_CONST_OFFSET; |

1221 | mp_limb_t cy; |

1222 | exponent ^= expbit; |

1223 | |

1224 | /* FIXME: not the whole multiplication has to be |

1225 | done. If we have the needed number of bits we |

1226 | only need the information whether more non-zero |

1227 | bits follow. */ |

1228 | if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET) |

1229 | cy = mpn_mul (pdest, psrc, numsize, |

1230 | &__tens[ttab->arrayoff |

1231 | + _FPIO_CONST_OFFSET], |

1232 | size); |

1233 | else |

1234 | cy = mpn_mul (pdest, &__tens[ttab->arrayoff |

1235 | + _FPIO_CONST_OFFSET], |

1236 | size, psrc, numsize); |

1237 | numsize += size; |

1238 | if (cy == 0) |

1239 | --numsize; |

1240 | (void) SWAP (psrc, pdest); |

1241 | } |

1242 | expbit <<= 1; |

1243 | ++ttab; |

1244 | } |

1245 | while (exponent != 0); |

1246 | |

1247 | if (psrc == den) |

1248 | memcpy (num, den, numsize * sizeof (mp_limb_t)); |

1249 | } |

1250 | |

1251 | /* Determine how many bits of the result we already have. */ |

1252 | count_leading_zeros (bits, num[numsize - 1]); |

1253 | bits = numsize * BITS_PER_MP_LIMB - bits; |

1254 | |

1255 | /* Now we know the exponent of the number in base two. |

1256 | Check it against the maximum possible exponent. */ |

1257 | if (__builtin_expect (bits > MAX_EXP, 0)) |

1258 | return overflow_value (negative); |

1259 | |

1260 | /* We have already the first BITS bits of the result. Together with |

1261 | the information whether more non-zero bits follow this is enough |

1262 | to determine the result. */ |

1263 | if (bits > MANT_DIG) |

1264 | { |

1265 | int i; |

1266 | const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB; |

1267 | const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB; |

1268 | const mp_size_t round_idx = least_bit == 0 ? least_idx - 1 |

1269 | : least_idx; |

1270 | const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1 |

1271 | : least_bit - 1; |

1272 | |

1273 | if (least_bit == 0) |

1274 | memcpy (retval, &num[least_idx], |

1275 | RETURN_LIMB_SIZE * sizeof (mp_limb_t)); |

1276 | else |

1277 | { |

1278 | for (i = least_idx; i < numsize - 1; ++i) |

1279 | retval[i - least_idx] = (num[i] >> least_bit) |

1280 | | (num[i + 1] |

1281 | << (BITS_PER_MP_LIMB - least_bit)); |

1282 | if (i - least_idx < RETURN_LIMB_SIZE) |

1283 | retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit; |

1284 | } |

1285 | |

1286 | /* Check whether any limb beside the ones in RETVAL are non-zero. */ |

1287 | for (i = 0; num[i] == 0; ++i) |

1288 | ; |

1289 | |

1290 | return round_and_return (retval, bits - 1, negative, |

1291 | num[round_idx], round_bit, |

1292 | int_no < dig_no || i < round_idx); |

1293 | /* NOTREACHED */ |

1294 | } |

1295 | else if (dig_no == int_no) |

1296 | { |

1297 | const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB; |

1298 | const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB; |

1299 | |

1300 | if (target_bit == is_bit) |

1301 | { |

1302 | memcpy (&retval[RETURN_LIMB_SIZE - numsize], num, |

1303 | numsize * sizeof (mp_limb_t)); |

1304 | /* FIXME: the following loop can be avoided if we assume a |

1305 | maximal MANT_DIG value. */ |

1306 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); |

1307 | } |

1308 | else if (target_bit > is_bit) |

1309 | { |

1310 | (void) mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize], |

1311 | num, numsize, target_bit - is_bit); |

1312 | /* FIXME: the following loop can be avoided if we assume a |

1313 | maximal MANT_DIG value. */ |

1314 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize); |

1315 | } |

1316 | else |

1317 | { |

1318 | mp_limb_t cy; |

1319 | assert (numsize < RETURN_LIMB_SIZE); |

1320 | |

1321 | cy = mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize], |

1322 | num, numsize, is_bit - target_bit); |

1323 | retval[RETURN_LIMB_SIZE - numsize - 1] = cy; |

1324 | /* FIXME: the following loop can be avoided if we assume a |

1325 | maximal MANT_DIG value. */ |

1326 | MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1); |

1327 | } |

1328 | |

1329 | return round_and_return (retval, bits - 1, negative, 0, 0, 0); |

1330 | /* NOTREACHED */ |

1331 | } |

1332 | |

1333 | /* Store the bits we already have. */ |

1334 | memcpy (retval, num, numsize * sizeof (mp_limb_t)); |

1335 | #if RETURN_LIMB_SIZE > 1 |

1336 | if (numsize < RETURN_LIMB_SIZE) |

1337 | # if RETURN_LIMB_SIZE == 2 |

1338 | retval[numsize] = 0; |

1339 | # else |

1340 | MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize); |

1341 | # endif |

1342 | #endif |

1343 | } |

1344 | |

1345 | /* We have to compute at least some of the fractional digits. */ |

1346 | { |

1347 | /* We construct a fraction and the result of the division gives us |

1348 | the needed digits. The denominator is 1.0 multiplied by the |

1349 | exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and |

1350 | 123e-6 gives 123 / 1000000. */ |

1351 | |

1352 | int expbit; |

1353 | int neg_exp; |

1354 | int more_bits; |

1355 | int need_frac_digits; |

1356 | mp_limb_t cy; |

1357 | mp_limb_t *psrc = den; |

1358 | mp_limb_t *pdest = num; |

1359 | const struct mp_power *ttab = &_fpioconst_pow10[0]; |

1360 | |

1361 | assert (dig_no > int_no |

1362 | && exponent <= 0 |

1363 | && exponent >= MIN_10_EXP - (DIG + 1)); |

1364 | |

1365 | /* We need to compute MANT_DIG - BITS fractional bits that lie |

1366 | within the mantissa of the result, the following bit for |

1367 | rounding, and to know whether any subsequent bit is 0. |

1368 | Computing a bit with value 2^-n means looking at n digits after |

1369 | the decimal point. */ |

1370 | if (bits > 0) |

1371 | { |

1372 | /* The bits required are those immediately after the point. */ |

1373 | assert (int_no > 0 && exponent == 0); |

1374 | need_frac_digits = 1 + MANT_DIG - bits; |

1375 | } |

1376 | else |

1377 | { |

1378 | /* The number is in the form .123eEXPONENT. */ |

1379 | assert (int_no == 0 && *startp != L_('0')); |

1380 | /* The number is at least 10^(EXPONENT-1), and 10^3 < |

1381 | 2^10. */ |

1382 | int neg_exp_2 = ((1 - exponent) * 10) / 3 + 1; |

1383 | /* The number is at least 2^-NEG_EXP_2. We need up to |

1384 | MANT_DIG bits following that bit. */ |

1385 | need_frac_digits = neg_exp_2 + MANT_DIG; |

1386 | /* However, we never need bits beyond 1/4 ulp of the smallest |

1387 | representable value. (That 1/4 ulp bit is only needed to |

1388 | determine tinyness on machines where tinyness is determined |

1389 | after rounding.) */ |

1390 | if (need_frac_digits > MANT_DIG - MIN_EXP + 2) |

1391 | need_frac_digits = MANT_DIG - MIN_EXP + 2; |

1392 | /* At this point, NEED_FRAC_DIGITS is the total number of |

1393 | digits needed after the point, but some of those may be |

1394 | leading 0s. */ |

1395 | need_frac_digits += exponent; |

1396 | /* Any cases underflowing enough that none of the fractional |

1397 | digits are needed should have been caught earlier (such |

1398 | cases are on the order of 10^-n or smaller where 2^-n is |

1399 | the least subnormal). */ |

1400 | assert (need_frac_digits > 0); |

1401 | } |

1402 | |

1403 | if (need_frac_digits > (intmax_t) dig_no - (intmax_t) int_no) |

1404 | need_frac_digits = (intmax_t) dig_no - (intmax_t) int_no; |

1405 | |

1406 | if ((intmax_t) dig_no > (intmax_t) int_no + need_frac_digits) |

1407 | { |

1408 | dig_no = int_no + need_frac_digits; |

1409 | more_bits = 1; |

1410 | } |

1411 | else |

1412 | more_bits = 0; |

1413 | |

1414 | neg_exp = (intmax_t) dig_no - (intmax_t) int_no - exponent; |

1415 | |

1416 | /* Construct the denominator. */ |

1417 | densize = 0; |

1418 | expbit = 1; |

1419 | do |

1420 | { |

1421 | if ((neg_exp & expbit) != 0) |

1422 | { |

1423 | mp_limb_t cy; |

1424 | neg_exp ^= expbit; |

1425 | |

1426 | if (densize == 0) |

1427 | { |

1428 | densize = ttab->arraysize - _FPIO_CONST_OFFSET; |

1429 | memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET], |

1430 | densize * sizeof (mp_limb_t)); |

1431 | } |

1432 | else |

1433 | { |

1434 | cy = mpn_mul (pdest, &__tens[ttab->arrayoff |

1435 | + _FPIO_CONST_OFFSET], |

1436 | ttab->arraysize - _FPIO_CONST_OFFSET, |

1437 | psrc, densize); |

1438 | densize += ttab->arraysize - _FPIO_CONST_OFFSET; |

1439 | if (cy == 0) |

1440 | --densize; |

1441 | (void) SWAP (psrc, pdest); |

1442 | } |

1443 | } |

1444 | expbit <<= 1; |

1445 | ++ttab; |

1446 | } |

1447 | while (neg_exp != 0); |

1448 | |

1449 | if (psrc == num) |

1450 | memcpy (den, num, densize * sizeof (mp_limb_t)); |

1451 | |

1452 | /* Read the fractional digits from the string. */ |

1453 | (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent |

1454 | #ifndef USE_WIDE_CHAR |

1455 | , decimal, decimal_len, thousands |

1456 | #endif |

1457 | ); |

1458 | |

1459 | /* We now have to shift both numbers so that the highest bit in the |

1460 | denominator is set. In the same process we copy the numerator to |

1461 | a high place in the array so that the division constructs the wanted |

1462 | digits. This is done by a "quasi fix point" number representation. |

1463 | |

1464 | num: ddddddddddd . 0000000000000000000000 |

1465 | |--- m ---| |

1466 | den: ddddddddddd n >= m |

1467 | |--- n ---| |

1468 | */ |

1469 | |

1470 | count_leading_zeros (cnt, den[densize - 1]); |

1471 | |

1472 | if (cnt > 0) |

1473 | { |

1474 | /* Don't call `mpn_shift' with a count of zero since the specification |

1475 | does not allow this. */ |

1476 | (void) mpn_lshift (den, den, densize, cnt); |

1477 | cy = mpn_lshift (num, num, numsize, cnt); |

1478 | if (cy != 0) |

1479 | num[numsize++] = cy; |

1480 | } |

1481 | |

1482 | /* Now we are ready for the division. But it is not necessary to |

1483 | do a full multi-precision division because we only need a small |

1484 | number of bits for the result. So we do not use mpn_divmod |

1485 | here but instead do the division here by hand and stop whenever |

1486 | the needed number of bits is reached. The code itself comes |

1487 | from the GNU MP Library by Torbj\"orn Granlund. */ |

1488 | |

1489 | exponent = bits; |

1490 | |

1491 | switch (densize) |

1492 | { |

1493 | case 1: |

1494 | { |

1495 | mp_limb_t d, n, quot; |

1496 | int used = 0; |

1497 | |

1498 | n = num[0]; |

1499 | d = den[0]; |

1500 | assert (numsize == 1 && n < d); |

1501 | |

1502 | do |

1503 | { |

1504 | udiv_qrnnd (quot, n, n, 0, d); |

1505 | |

1506 | #define got_limb \ |

1507 | if (bits == 0) \ |

1508 | { \ |

1509 | register int cnt; \ |

1510 | if (quot == 0) \ |

1511 | cnt = BITS_PER_MP_LIMB; \ |

1512 | else \ |

1513 | count_leading_zeros (cnt, quot); \ |

1514 | exponent -= cnt; \ |

1515 | if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \ |

1516 | { \ |

1517 | used = MANT_DIG + cnt; \ |

1518 | retval[0] = quot >> (BITS_PER_MP_LIMB - used); \ |

1519 | bits = MANT_DIG + 1; \ |

1520 | } \ |

1521 | else \ |

1522 | { \ |

1523 | /* Note that we only clear the second element. */ \ |

1524 | /* The conditional is determined at compile time. */ \ |

1525 | if (RETURN_LIMB_SIZE > 1) \ |

1526 | retval[1] = 0; \ |

1527 | retval[0] = quot; \ |

1528 | bits = -cnt; \ |

1529 | } \ |

1530 | } \ |

1531 | else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \ |

1532 | mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \ |

1533 | quot); \ |

1534 | else \ |

1535 | { \ |

1536 | used = MANT_DIG - bits; \ |

1537 | if (used > 0) \ |

1538 | mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \ |

1539 | } \ |

1540 | bits += BITS_PER_MP_LIMB |

1541 | |

1542 | got_limb; |

1543 | } |

1544 | while (bits <= MANT_DIG); |

1545 | |

1546 | return round_and_return (retval, exponent - 1, negative, |

1547 | quot, BITS_PER_MP_LIMB - 1 - used, |

1548 | more_bits || n != 0); |

1549 | } |

1550 | case 2: |

1551 | { |

1552 | mp_limb_t d0, d1, n0, n1; |

1553 | mp_limb_t quot = 0; |

1554 | int used = 0; |

1555 | |

1556 | d0 = den[0]; |

1557 | d1 = den[1]; |

1558 | |

1559 | if (numsize < densize) |

1560 | { |

1561 | if (num[0] >= d1) |

1562 | { |

1563 | /* The numerator of the number occupies fewer bits than |

1564 | the denominator but the one limb is bigger than the |

1565 | high limb of the numerator. */ |

1566 | n1 = 0; |

1567 | n0 = num[0]; |

1568 | } |

1569 | else |

1570 | { |

1571 | if (bits <= 0) |

1572 | exponent -= BITS_PER_MP_LIMB; |

1573 | else |

1574 | { |

1575 | if (bits + BITS_PER_MP_LIMB <= MANT_DIG) |

1576 | mpn_lshift_1 (retval, RETURN_LIMB_SIZE, |

1577 | BITS_PER_MP_LIMB, 0); |

1578 | else |

1579 | { |

1580 | used = MANT_DIG - bits; |

1581 | if (used > 0) |

1582 | mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); |

1583 | } |

1584 | bits += BITS_PER_MP_LIMB; |

1585 | } |

1586 | n1 = num[0]; |

1587 | n0 = 0; |

1588 | } |

1589 | } |

1590 | else |

1591 | { |

1592 | n1 = num[1]; |

1593 | n0 = num[0]; |

1594 | } |

1595 | |

1596 | while (bits <= MANT_DIG) |

1597 | { |

1598 | mp_limb_t r; |

1599 | |

1600 | if (n1 == d1) |

1601 | { |

1602 | /* QUOT should be either 111..111 or 111..110. We need |

1603 | special treatment of this rare case as normal division |

1604 | would give overflow. */ |

1605 | quot = ~(mp_limb_t) 0; |

1606 | |

1607 | r = n0 + d1; |

1608 | if (r < d1) /* Carry in the addition? */ |

1609 | { |

1610 | add_ssaaaa (n1, n0, r - d0, 0, 0, d0); |

1611 | goto have_quot; |

1612 | } |

1613 | n1 = d0 - (d0 != 0); |

1614 | n0 = -d0; |

1615 | } |

1616 | else |

1617 | { |

1618 | udiv_qrnnd (quot, r, n1, n0, d1); |

1619 | umul_ppmm (n1, n0, d0, quot); |

1620 | } |

1621 | |

1622 | q_test: |

1623 | if (n1 > r || (n1 == r && n0 > 0)) |

1624 | { |

1625 | /* The estimated QUOT was too large. */ |

1626 | --quot; |

1627 | |

1628 | sub_ddmmss (n1, n0, n1, n0, 0, d0); |

1629 | r += d1; |

1630 | if (r >= d1) /* If not carry, test QUOT again. */ |

1631 | goto q_test; |

1632 | } |

1633 | sub_ddmmss (n1, n0, r, 0, n1, n0); |

1634 | |

1635 | have_quot: |

1636 | got_limb; |

1637 | } |

1638 | |

1639 | return round_and_return (retval, exponent - 1, negative, |

1640 | quot, BITS_PER_MP_LIMB - 1 - used, |

1641 | more_bits || n1 != 0 || n0 != 0); |

1642 | } |

1643 | default: |

1644 | { |

1645 | int i; |

1646 | mp_limb_t cy, dX, d1, n0, n1; |

1647 | mp_limb_t quot = 0; |

1648 | int used = 0; |

1649 | |

1650 | dX = den[densize - 1]; |

1651 | d1 = den[densize - 2]; |

1652 | |

1653 | /* The division does not work if the upper limb of the two-limb |

1654 | numerator is greater than the denominator. */ |

1655 | if (mpn_cmp (num, &den[densize - numsize], numsize) > 0) |

1656 | num[numsize++] = 0; |

1657 | |

1658 | if (numsize < densize) |

1659 | { |

1660 | mp_size_t empty = densize - numsize; |

1661 | register int i; |

1662 | |

1663 | if (bits <= 0) |

1664 | exponent -= empty * BITS_PER_MP_LIMB; |

1665 | else |

1666 | { |

1667 | if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG) |

1668 | { |

1669 | /* We make a difference here because the compiler |

1670 | cannot optimize the `else' case that good and |

1671 | this reflects all currently used FLOAT types |

1672 | and GMP implementations. */ |

1673 | #if RETURN_LIMB_SIZE <= 2 |

1674 | assert (empty == 1); |

1675 | mpn_lshift_1 (retval, RETURN_LIMB_SIZE, |

1676 | BITS_PER_MP_LIMB, 0); |

1677 | #else |

1678 | for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i) |

1679 | retval[i] = retval[i - empty]; |

1680 | while (i >= 0) |

1681 | retval[i--] = 0; |

1682 | #endif |

1683 | } |

1684 | else |

1685 | { |

1686 | used = MANT_DIG - bits; |

1687 | if (used >= BITS_PER_MP_LIMB) |

1688 | { |

1689 | register int i; |

1690 | (void) mpn_lshift (&retval[used |

1691 | / BITS_PER_MP_LIMB], |

1692 | retval, |

1693 | (RETURN_LIMB_SIZE |

1694 | - used / BITS_PER_MP_LIMB), |

1695 | used % BITS_PER_MP_LIMB); |

1696 | for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i) |

1697 | retval[i] = 0; |

1698 | } |

1699 | else if (used > 0) |

1700 | mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0); |

1701 | } |

1702 | bits += empty * BITS_PER_MP_LIMB; |

1703 | } |

1704 | for (i = numsize; i > 0; --i) |

1705 | num[i + empty] = num[i - 1]; |

1706 | MPN_ZERO (num, empty + 1); |

1707 | } |

1708 | else |

1709 | { |

1710 | int i; |

1711 | assert (numsize == densize); |

1712 | for (i = numsize; i > 0; --i) |

1713 | num[i] = num[i - 1]; |

1714 | num[0] = 0; |

1715 | } |

1716 | |

1717 | den[densize] = 0; |

1718 | n0 = num[densize]; |

1719 | |

1720 | while (bits <= MANT_DIG) |

1721 | { |

1722 | if (n0 == dX) |

1723 | /* This might over-estimate QUOT, but it's probably not |

1724 | worth the extra code here to find out. */ |

1725 | quot = ~(mp_limb_t) 0; |

1726 | else |

1727 | { |

1728 | mp_limb_t r; |

1729 | |

1730 | udiv_qrnnd (quot, r, n0, num[densize - 1], dX); |

1731 | umul_ppmm (n1, n0, d1, quot); |

1732 | |

1733 | while (n1 > r || (n1 == r && n0 > num[densize - 2])) |

1734 | { |

1735 | --quot; |

1736 | r += dX; |

1737 | if (r < dX) /* I.e. "carry in previous addition?" */ |

1738 | break; |

1739 | n1 -= n0 < d1; |

1740 | n0 -= d1; |

1741 | } |

1742 | } |

1743 | |

1744 | /* Possible optimization: We already have (q * n0) and (1 * n1) |

1745 | after the calculation of QUOT. Taking advantage of this, we |

1746 | could make this loop make two iterations less. */ |

1747 | |

1748 | cy = mpn_submul_1 (num, den, densize + 1, quot); |

1749 | |

1750 | if (num[densize] != cy) |

1751 | { |

1752 | cy = mpn_add_n (num, num, den, densize); |

1753 | assert (cy != 0); |

1754 | --quot; |

1755 | } |

1756 | n0 = num[densize] = num[densize - 1]; |

1757 | for (i = densize - 1; i > 0; --i) |

1758 | num[i] = num[i - 1]; |

1759 | num[0] = 0; |

1760 | |

1761 | got_limb; |

1762 | } |

1763 | |

1764 | for (i = densize; num[i] == 0 && i >= 0; --i) |

1765 | ; |

1766 | return round_and_return (retval, exponent - 1, negative, |

1767 | quot, BITS_PER_MP_LIMB - 1 - used, |

1768 | more_bits || i >= 0); |

1769 | } |

1770 | } |

1771 | } |

1772 | |

1773 | /* NOTREACHED */ |

1774 | } |

1775 |