Diff
Not logged in

Differences From Artifact [9ff45d8457]:

To Artifact [dc5c29ca69]:


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);