+ if (SCM_I_INUMP (n))
+ mpz_init_set_si (nn, SCM_I_INUM (n));
+ else
+ mpz_init_set (nn, SCM_I_BIG_MPZ (n));
+
+ neg = (mpz_sgn (nn) < 0) ^ (mpz_sgn (dd) < 0);
+ mpz_abs (nn, nn);
+ mpz_abs (dd, dd);
+
+ /* Now we need to find the value of e such that:
+
+ For e <= 0:
+ b^{p-1} - 1/2b <= b^-e n / d < b^p - 1/2 [1A]
+ (2 b^p - 1) <= 2 b b^-e n / d < (2 b^p - 1) b [2A]
+ (2 b^p - 1) d <= 2 b b^-e n < (2 b^p - 1) d b [3A]
+
+ For e >= 0:
+ b^{p-1} - 1/2b <= n / b^e d < b^p - 1/2 [1B]
+ (2 b^p - 1) <= 2 b n / b^e d < (2 b^p - 1) b [2B]
+ (2 b^p - 1) d b^e <= 2 b n < (2 b^p - 1) d b b^e [3B]
+
+ where: p = DBL_MANT_DIG
+ b = FLT_RADIX (here assumed to be 2)
+
+ After rounding, the mantissa must be an integer between b^{p-1} and
+ (b^p - 1), except for subnormal numbers. In the inequations [1A]
+ and [1B], the middle expression represents the mantissa *before*
+ rounding, and therefore is bounded by the range of values that will
+ round to a floating-point number with the exponent e. The upper
+ bound is (b^p - 1 + 1/2) = (b^p - 1/2), and is exclusive because
+ ties will round up to the next power of b. The lower bound is
+ (b^{p-1} - 1/2b), and is inclusive because ties will round toward
+ this power of b. Here we subtract 1/2b instead of 1/2 because it
+ is in the range of the next smaller exponent, where the
+ representable numbers are closer together by a factor of b.
+
+ Inequations [2A] and [2B] are derived from [1A] and [1B] by
+ multiplying by 2b, and in [3A] and [3B] we multiply by the
+ denominator of the middle value to obtain integer expressions.
+
+ In the code below, we refer to the three expressions in [3A] or
+ [3B] as lo, x, and hi. If the number is normalizable, we will
+ achieve the goal: lo <= x < hi */
+
+ /* Make an initial guess for e */
+ e = mpz_sizeinbase (nn, 2) - mpz_sizeinbase (dd, 2) - (DBL_MANT_DIG-1);
+ if (e < DBL_MIN_EXP - DBL_MANT_DIG)
+ e = DBL_MIN_EXP - DBL_MANT_DIG;
+
+ /* Compute the initial values of lo, x, and hi
+ based on the initial guess of e */
+ mpz_inits (lo, hi, x, NULL);
+ mpz_mul_2exp (x, nn, 2 + ((e < 0) ? -e : 0));
+ mpz_mul (lo, dd, scm_i_divide2double_lo2b);
+ if (e > 0)
+ mpz_mul_2exp (lo, lo, e);
+ mpz_mul_2exp (hi, lo, 1);
+
+ /* Adjust e as needed to satisfy the inequality lo <= x < hi,
+ (but without making e less then the minimum exponent) */
+ while (mpz_cmp (x, lo) < 0 && e > DBL_MIN_EXP - DBL_MANT_DIG)
+ {
+ mpz_mul_2exp (x, x, 1);
+ e--;
+ }
+ while (mpz_cmp (x, hi) >= 0)
+ {
+ /* If we ever used lo's value again,
+ we would need to double lo here. */
+ mpz_mul_2exp (hi, hi, 1);
+ e++;
+ }
+
+ /* Now compute the rounded mantissa:
+ n / b^e d (if e >= 0)
+ n b^-e / d (if e <= 0) */
HCoop / bpt / guile.git / blobdiff