| ︙ | | |
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
-
+
|
* Derived quantities.
*/
#define TEN_PMAX 22 /* floor(FP_PRECISION*log(2)/log(5)) */
#define QUICK_MAX 14 /* floor((FP_PRECISION-1)*log(2)/log(10))-1 */
#define BLETCH 0x10 /* Highest power of two that is greater than
* DBL_MAX_10_EXP, divided by 16. */
#define DIGIT_GROUP 8 /* floor(DIGIT_BIT*log(2)/log(10)) */
#define DIGIT_GROUP 8 /* floor(MP_DIGIT_BIT*log(2)/log(10)) */
/*
* Union used to dismantle floating point numbers.
*/
typedef union Double {
struct {
|
| ︙ | | |
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
|
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
|
-
+
-
+
|
bignumRepPtr);
mp_add_d(bignumRepPtr, (mp_digit) digit, bignumRepPtr);
} else {
/*
* More than single digit multiplication. Multiply by the appropriate
* small powers of 5, and then shift. Large strings of zeroes are
* eaten 256 at a time; this is less efficient than it could be, but
* seems implausible. We presume that DIGIT_BIT is at least 27. The
* seems implausible. We presume that MP_DIGIT_BIT is at least 27. The
* first multiplication, by up to 10**7, is done with a one-DIGIT
* multiply (this presumes that DIGIT_BIT >= 24).
* multiply (this presumes that MP_DIGIT_BIT >= 24).
*/
n = numZeros + 1;
mp_mul_d(bignumRepPtr, (mp_digit) pow10_wide[n&0x7], bignumRepPtr);
for (i=3; i<=7; ++i) {
if (n & (1 << i)) {
mp_mul(bignumRepPtr, pow5+i, bignumRepPtr);
|
| ︙ | | |
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
|
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
|
-
+
-
+
|
/*
*----------------------------------------------------------------------
*
* ShouldBankerRoundUpPowD --
*
* Test whether bankers' rounding should round a digit up. Assumption is
* made that the denominator of the fraction being tested is a power of
* 2**DIGIT_BIT.
* 2**MP_DIGIT_BIT.
*
* Results:
* Returns 1 iff the fraction is more than 1/2, or if the fraction is
* exactly 1/2 and the digit is odd.
*
*----------------------------------------------------------------------
*/
static inline int
ShouldBankerRoundUpPowD(
mp_int *b, /* Numerator of the fraction. */
int sd, /* Denominator is 2**(sd*DIGIT_BIT). */
int sd, /* Denominator is 2**(sd*MP_DIGIT_BIT). */
int isodd) /* 1 if the digit is odd, 0 if even. */
{
int i;
static const mp_digit topbit = ((mp_digit)1) << (MP_DIGIT_BIT - 1);
if (b->used < sd || (b->dp[sd-1] & topbit) == 0) {
return 0;
|
| ︙ | | |
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219
|
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216
3217
3218
3219
|
-
+
-
+
|
*----------------------------------------------------------------------
*/
static inline int
ShouldBankerRoundUpToNextPowD(
mp_int *b, /* Numerator of the fraction. */
mp_int *m, /* Numerator of the rounding tolerance. */
int sd, /* Common denominator is 2**(sd*DIGIT_BIT). */
int sd, /* Common denominator is 2**(sd*MP_DIGIT_BIT). */
int convType, /* Conversion type: STEELE defeats
* round-to-even (not sure why one wants to do
* this; I copied it from Gay). FIXME */
int isodd, /* 1 if the integer significand is odd. */
mp_int *temp) /* Work area for the calculation. */
{
int i;
/*
* Compare B and S-m - which is the same as comparing B+m and S - which we
* do by computing b+m and doing a bitwhack compare against
* 2**(DIGIT_BIT*sd)
* 2**(MP_DIGIT_BIT*sd)
*/
mp_add(b, m, temp);
if (temp->used <= sd) { /* Too few digits to be > s */
return 0;
}
if (temp->used > sd+1 || temp->dp[sd] > 1) {
|
| ︙ | | |
3237
3238
3239
3240
3241
3242
3243
3244
3245
3246
3247
3248
3249
3250
3251
|
3237
3238
3239
3240
3241
3242
3243
3244
3245
3246
3247
3248
3249
3250
3251
|
-
+
|
*----------------------------------------------------------------------
*
* ShorteningBignumConversionPowD --
*
* Converts a double-precision number to the shortest string of digits
* that reconverts exactly to the given number, or to 'ilim' digits if
* that will yield a shorter result. The denominator in David Gay's
* conversion algorithm is known to be a power of 2**DIGIT_BIT, and hence
* conversion algorithm is known to be a power of 2**MP_DIGIT_BIT, and hence
* the division in the main loop may be replaced by a digit shift and
* mask.
*
* Results:
* Returns the string of significant decimal digits, in newly allocated
* memory
*
|
| ︙ | | |
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
3331
3332
3333
|
3319
3320
3321
3322
3323
3324
3325
3326
3327
3328
3329
3330
3331
3332
3333
|
-
+
|
if (m2plus > m2minus) {
mp_init_copy(&mplus, &mminus);
mp_mul_2d(&mplus, m2plus-m2minus, &mplus);
}
mp_init(&temp);
/*
* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
* Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
* by mp_digit extraction.
*/
i = 0;
for (;;) {
if (b.used <= sd) {
digit = 0;
|
| ︙ | | |
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
3443
3444
3445
|
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
3443
3444
3445
|
-
+
|
*----------------------------------------------------------------------
*
* StrictBignumConversionPowD --
*
* Converts a double-precision number to a fixed-lengt string of 'ilim'
* digits (or 'ilim1' if log10(d) has been overestimated). The
* denominator in David Gay's conversion algorithm is known to be a power
* of 2**DIGIT_BIT, and hence the division in the main loop may be
* of 2**MP_DIGIT_BIT, and hence the division in the main loop may be
* replaced by a digit shift and mask.
*
* Results:
* Returns the string of significant decimal digits, in newly allocated
* memory.
*
* Side effects:
|
| ︙ | | |
3492
3493
3494
3495
3496
3497
3498
3499
3500
3501
3502
3503
3504
3505
3506
|
3492
3493
3494
3495
3496
3497
3498
3499
3500
3501
3502
3503
3504
3505
3506
|
-
+
|
mp_mul_d(&b, 10, &b);
ilim = ilim1;
--k;
}
mp_init(&temp);
/*
* Loop through the digits. Do division and mod by s == 2**(sd*DIGIT_BIT)
* Loop through the digits. Do division and mod by s == 2**(sd*MP_DIGIT_BIT)
* by mp_digit extraction.
*/
i = 1;
for (;;) {
if (b.used <= sd) {
digit = 0;
|
| ︙ | | |
4258
4259
4260
4261
4262
4263
4264
4265
4266
4267
4268
4269
4270
4271
4272
|
4258
4259
4260
4261
4262
4263
4264
4265
4266
4267
4268
4269
4270
4271
4272
|
-
+
|
*/
return ShorteningInt64Conversion(&d, convType, bw, b2, b5, m2plus,
m2minus, m5, s2, s5, k, len, ilim, ilim1, decpt, endPtr);
} else if (s5 == 0) {
/*
* The denominator is a power of 2, so we can replace division by
* digit shifts. First we round up s2 to a multiple of DIGIT_BIT,
* digit shifts. First we round up s2 to a multiple of MP_DIGIT_BIT,
* and adjust m2 and b2 accordingly. Then we launch into a version
* of the comparison that's specialized for the 'power of mp_digit
* in the denominator' case.
*/
if (s2 % MP_DIGIT_BIT != 0) {
int delta = MP_DIGIT_BIT - (s2 % MP_DIGIT_BIT);
|
| ︙ | | |
4314
4315
4316
4317
4318
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
|
4314
4315
4316
4317
4318
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
|
-
+
|
*/
return StrictInt64Conversion(&d, convType, bw, b2, b5, s2, s5, k,
len, ilim, ilim1, decpt, endPtr);
} else if (s5 == 0) {
/*
* The denominator is a power of 2, so we can replace division by
* digit shifts. First we round up s2 to a multiple of DIGIT_BIT,
* digit shifts. First we round up s2 to a multiple of MP_DIGIT_BIT,
* and adjust m2 and b2 accordingly. Then we launch into a version
* of the comparison that's specialized for the 'power of mp_digit
* in the denominator' case.
*/
if (s2 % MP_DIGIT_BIT != 0) {
int delta = MP_DIGIT_BIT - (s2 % MP_DIGIT_BIT);
|
| ︙ | | |