Diff
Not logged in

Differences From Artifact [1d51246842]:

To Artifact [df6dd2545e]:


141
142
143
144
145
146
147
148

149
150
151
152
153
154
155
141
142
143
144
145
146
147

148
149
150
151
152
153
154
155







-
+







 * 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 {
1481
1482
1483
1484
1485
1486
1487
1488

1489
1490

1491
1492
1493
1494
1495
1496
1497
1481
1482
1483
1484
1485
1486
1487

1488
1489

1490
1491
1492
1493
1494
1495
1496
1497







-
+

-
+







		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);
3126
3127
3128
3129
3130
3131
3132
3133

3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144
3145

3146
3147
3148
3149
3150
3151
3152
3126
3127
3128
3129
3130
3131
3132

3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144

3145
3146
3147
3148
3149
3150
3151
3152







-
+











-
+







/*
 *----------------------------------------------------------------------
 *
 * 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;
3177
3178
3179
3180
3181
3182
3183
3184

3185
3186
3187
3188
3189
3190
3191
3192
3193

3194
3195
3196
3197
3198
3199
3200
3177
3178
3179
3180
3181
3182
3183

3184
3185
3186
3187
3188
3189
3190
3191
3192

3193
3194
3195
3196
3197
3198
3199
3200







-
+








-
+







 *----------------------------------------------------------------------
 */

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 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) {
3214
3215
3216
3217
3218
3219
3220
3221

3222
3223
3224
3225
3226
3227
3228
3214
3215
3216
3217
3218
3219
3220

3221
3222
3223
3224
3225
3226
3227
3228







-
+







 *----------------------------------------------------------------------
 *
 * 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
 *
3294
3295
3296
3297
3298
3299
3300
3301

3302
3303
3304
3305
3306
3307
3308
3294
3295
3296
3297
3298
3299
3300

3301
3302
3303
3304
3305
3306
3307
3308







-
+







    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;
3406
3407
3408
3409
3410
3411
3412
3413

3414
3415
3416
3417
3418
3419
3420
3406
3407
3408
3409
3410
3411
3412

3413
3414
3415
3416
3417
3418
3419
3420







-
+







 *----------------------------------------------------------------------
 *
 * 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:
3464
3465
3466
3467
3468
3469
3470
3471

3472
3473
3474
3475
3476
3477
3478
3464
3465
3466
3467
3468
3469
3470

3471
3472
3473
3474
3475
3476
3477
3478







-
+







    if (b.used <= sd) {
	mp_mul_d(&b, 10, &b);
	ilim = ilim1;
	--k;
    }

    /*
     * 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;
4207
4208
4209
4210
4211
4212
4213
4214

4215
4216
4217
4218
4219
4220
4221
4207
4208
4209
4210
4211
4212
4213

4214
4215
4216
4217
4218
4219
4220
4221







-
+







	     */

	    return ShorteningInt64Conversion(&d, 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);
4263
4264
4265
4266
4267
4268
4269
4270

4271
4272
4273
4274
4275
4276
4277
4263
4264
4265
4266
4267
4268
4269

4270
4271
4272
4273
4274
4275
4276
4277







-
+







	     */

	    return StrictInt64Conversion(&d, 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);