~ chicken-core (chicken-5) 15fd900137c5007d60e6300973df7e73b2934110
commit 15fd900137c5007d60e6300973df7e73b2934110
Author: Peter Bex <peter@more-magic.net>
AuthorDate: Sun Mar 15 17:56:56 2015 +0100
Commit: Peter Bex <peter@more-magic.net>
CommitDate: Sun May 31 14:55:23 2015 +0200
Convert exact->inexact to a horribly long and ugly C function, making it inlineable (needed for converting + and - to be inlineable)
diff --git a/chicken.h b/chicken.h
index f0387662..2970689e 100644
--- a/chicken.h
+++ b/chicken.h
@@ -2183,7 +2183,6 @@ C_fctexport double C_fcall C_bignum_to_double(C_word bignum) C_regparm;
C_fctexport C_word C_fcall C_a_i_cpu_time(C_word **a, int c, C_word buf) C_regparm;
/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
C_fctexport C_word C_fcall C_a_i_string_to_number(C_word **a, int c, C_word str, C_word radix) C_regparm;
-/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
C_fctexport C_word C_fcall C_a_i_exact_to_inexact(C_word **a, int c, C_word n) C_regparm;
C_fctexport C_word C_fcall C_i_file_exists_p(C_word name, C_word file, C_word dir) C_regparm;
diff --git a/library.scm b/library.scm
index 3c6f756f..e9b0190c 100644
--- a/library.scm
+++ b/library.scm
@@ -36,7 +36,7 @@
##sys#format-here-doc-warning
exit-in-progress
maximal-string-length find-ratio-between find-ratio
- make-complex flonum->ratnum ratnum rat+/- minimum-denorm-flonum-expt
+ make-complex flonum->ratnum ratnum rat+/-
+maximum-allowed-exponent+ mantexp->dbl ldexp round-quotient
##sys#string->compnum ##sys#bignum-extract-digits ##sys#internal-gcd)
(not inline ##sys#user-read-hook ##sys#error-hook ##sys#signal-hook ##sys#schedule
@@ -1130,58 +1130,8 @@ EOF
(inexact->exact (%cplxnum-imag x))))
(else (##sys#error-bad-number x 'inexact->exact))))
-;; Exponent of the lowest allowed flonum; if we get any lower we get zero.
-;; In other words, this is the first (smallest) flonum after 0.
-;; Equal to (expt 2.0 (- flonum-minimum-exponent flonum-precision))
-(define minimum-denorm-flonum-expt (fx- flonum-minimum-exponent flonum-precision))
-
(define (exact->inexact x)
- (cond ((##core#inline "C_fixnump" x)
- (##core#inline_allocate ("C_a_i_fix_to_flo" 4) x))
- ((##core#inline "C_i_flonump" x) x)
- ((##core#inline "C_i_bignump" x)
- (##core#inline_allocate ("C_a_u_i_big_to_flo" 4) x))
- ((ratnum? x)
- ;; This tries to keep the numbers within representable ranges
- ;; and tries to drop as few significant digits as possible by
- ;; bringing the two numbers to within the same powers of two.
- ;; See algorithms M & N in Knuth, 4.2.1
- (let* ((n1 (%ratnum-numerator x))
- (an (##sys#integer-abs n1))
- (d1 (%ratnum-denominator x))
- ;; Approximate distance between the numbers in powers
- ;; of 2 ie, 2^e-1 < n/d < 2^e+1 (e is the *un*biased
- ;; value of e_w in M2)
- ;; XXX: What if b != 2 (ie, flonum-radix is not 2)?
- (e (fx- (integer-length an) (integer-length d1)))
- (rnd (lambda (n d e) ; Here, 1 <= n/d < 2 (normalized) [N5]
- ;; Cannot shift above the available precision,
- ;; and can't have an exponent that's below the
- ;; minimum flonum exponent.
- (let* ((s (min (fx- flonum-precision 1)
- (fx- e minimum-denorm-flonum-expt)))
- (norm (##sys#/-2 (##sys#integer-shift n s) d))
- (r (round norm))
- (fraction (exact->inexact r))
- (exp (fx- e s)))
- (let ((res (fp* fraction (expt 2.0 exp))))
- (if (negative? n1) (##sys#--2 0 res) res)))))
- (scale (lambda (n d) ; Here, 1/2 <= n/d < 2 [N3]
- (if (##sys#<-2 n d) ; n/d < 1?
- ;; Scale left [N3]; only needed once (see note in M3)
- (rnd (##sys#integer-shift n 1) d (fx- e 1))
- ;; Already normalized
- (rnd n d e)))))
- ;; After this step, which shifts the smaller number to
- ;; align with the larger, "f" in algorithm N is represented
- ;; in the procedures above by n/d.
- (if (negative? e)
- (scale (##sys#integer-shift an (##sys#--2 0 e)) d1)
- (scale an (##sys#integer-shift d1 e)))))
- ((cplxnum? x)
- (make-complex (exact->inexact (%cplxnum-real x))
- (exact->inexact (%cplxnum-imag x))))
- (else (##sys#error-bad-number x 'exact->inexact))))
+ (##core#inline_allocate ("C_a_i_exact_to_inexact" 12) x))
(define ##sys#exact->inexact exact->inexact)
(define ##sys#inexact->exact inexact->exact)
diff --git a/runtime.c b/runtime.c
index 5de43743..b6db3da1 100644
--- a/runtime.c
+++ b/runtime.c
@@ -9777,17 +9777,171 @@ void C_ccall C_string_to_symbol(C_word c, C_word closure, C_word k, C_word strin
C_kontinue(k, s);
}
+C_inline int integer_length_abs(C_word x)
+{
+ if (x & C_FIXNUM_BIT) {
+ return C_ilen(labs(C_unfix(x)));
+ } else {
+ C_uword result = (C_bignum_size(x) - 1) * C_BIGNUM_DIGIT_LENGTH,
+ *last_digit = C_bignum_digits(x) + C_bignum_size(x) - 1,
+ last_digit_length = C_ilen(*last_digit);
+ return result + last_digit_length;
+ }
+}
-/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
+/* This will usually return a flonum, but it may also return a cplxnum
+ * consisting of two flonums, making for a total of 12 words.
+ */
C_regparm C_word C_fcall
-C_a_i_exact_to_inexact(C_word **a, int c, C_word n)
-{
- if(n & C_FIXNUM_BIT)
- return C_flonum(a, (double)C_unfix(n));
- else if(C_immediatep(n) || C_block_header(n) != C_FLONUM_TAG)
- barf(C_BAD_ARGUMENT_TYPE_ERROR, "exact->inexact", n);
+C_a_i_exact_to_inexact(C_word **ptr, int c, C_word n)
+{
+ if (n & C_FIXNUM_BIT) {
+ return C_flonum(ptr, (double)C_unfix(n));
+ } else if (C_immediatep(n)) {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "exact->inexact", n);
+ } else if (C_block_header(n) == C_FLONUM_TAG) {
+ return n;
+ } else if (C_truep(C_bignump(n))) {
+ return C_a_u_i_big_to_flo(ptr, c, n);
+ } else if (C_block_header(n) == C_STRUCTURE3_TAG &&
+ (C_block_item(n, 0) == C_cplxnum_type_tag)) {
+ return C_cplxnum(ptr, C_a_i_exact_to_inexact(ptr, 1, C_block_item(n, 1)),
+ C_a_i_exact_to_inexact(ptr, 1, C_block_item(n, 2)));
+ /* The horribly painful case: ratnums */
+ } else if (C_block_header(n) == C_STRUCTURE3_TAG &&
+ (C_block_item(n, 0) == C_ratnum_type_tag)) {
+ /* This tries to keep the numbers within representable ranges and
+ * tries to drop as few significant digits as possible by bringing
+ * the two numbers to within the same powers of two. See
+ * algorithms M & N in Knuth, 4.2.1.
+ */
+ C_word num = C_block_item(n, 1), denom = C_block_item(n, 2),
+ /* e = approx. distance between the numbers in powers of 2.
+ * ie, 2^e-1 < n/d < 2^e+1 (e is the *un*biased value of
+ * e_w in M2. TODO: What if b!=2 (ie, flonum-radix isn't 2)?
+ */
+ e = integer_length_abs(num) - integer_length_abs(denom),
+ ab[C_SIZEOF_FIX_BIGNUM*4], *a = ab, tmp1 = 0, tmp2 = 0, tmp3 = 0,
+ shift_amount, negp = C_i_integer_negativep(num), q, r, len;
+ C_uword *d;
+ double res, fraction;
+
+ /* Simplify logic by ensuring bignums */
+ if (num & C_FIXNUM_BIT) num = C_a_u_i_fix_to_big(&a, num);
+ if (denom & C_FIXNUM_BIT) denom = C_a_u_i_fix_to_big(&a, denom);
+
+ /* Align numbers by shifting the smaller to the same size of the
+ * larger. After this, "f" in alg. N is represented by num/denom.
+ */
+ if (e < 0) {
+ tmp1 = allocate_tmp_bignum(C_fix(C_bignum_size(denom)),
+ C_SCHEME_FALSE, C_SCHEME_TRUE);
+ d = C_bignum_digits(tmp1) - e / C_BIGNUM_DIGIT_LENGTH;
+ C_memcpy(d, C_bignum_digits(num), C_wordstobytes(C_bignum_size(num)));
+ shift_amount = -e % C_BIGNUM_DIGIT_LENGTH;
+ if(shift_amount > 0) {
+ bignum_digits_destructive_shift_left(
+ d, C_bignum_digits(tmp1) + C_bignum_size(tmp1), shift_amount);
+ }
+ num = tmp1;
+ } else if (e > 0) {
+ tmp1 = allocate_tmp_bignum(C_fix(C_bignum_size(num)),
+ C_SCHEME_FALSE, C_SCHEME_TRUE);
+ d = C_bignum_digits(tmp1) + e / C_BIGNUM_DIGIT_LENGTH;
+ C_memcpy(d, C_bignum_digits(denom), C_wordstobytes(C_bignum_size(denom)));
+ shift_amount = e % C_BIGNUM_DIGIT_LENGTH;
+ if(shift_amount > 0) {
+ bignum_digits_destructive_shift_left(
+ d, C_bignum_digits(tmp1) + C_bignum_size(tmp1), shift_amount);
+ }
+ denom = tmp1;
+ }
+ /* From here on, 1/2 <= n/d < 2 [N3] */
+ if (C_truep(C_i_integer_lessp(num, denom))) { /* n/d < 1? */
+ len = C_bignum_size(num) + 1;
+ tmp2 = allocate_tmp_bignum(C_fix(len), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_digits_destructive_copy(tmp2, num);
+ d = C_bignum_digits(tmp2);
+ d[len-1] = 0; /* Init most significant digit */
+ bignum_digits_destructive_shift_left(d, d + len, 1);
+ num = tmp2;
+ e -= 1;
+ }
- return n;
+ /* Here, 1 <= n/d < 2 (normalized) [N5] */
+ shift_amount = nmin(DBL_MANT_DIG-1, e - (DBL_MIN_EXP - DBL_MANT_DIG));
+
+ len = C_bignum_size(num) + shift_amount / C_BIGNUM_DIGIT_LENGTH + 1;
+ tmp3 = allocate_tmp_bignum(C_fix(len), C_SCHEME_FALSE, C_SCHEME_TRUE);
+ d = C_bignum_digits(tmp3) + shift_amount / C_BIGNUM_DIGIT_LENGTH;
+ C_memcpy(d, C_bignum_digits(num), C_wordstobytes(C_bignum_size(num)));
+ shift_amount = shift_amount % C_BIGNUM_DIGIT_LENGTH;
+ if (shift_amount > 0) {
+ bignum_digits_destructive_shift_left(
+ d, C_bignum_digits(tmp3) + len, shift_amount);
+ }
+ num = tmp3;
+
+ /* Now, calculate round(num/denom). We start with a quotient&remainder */
+ switch(bignum_cmp_unsigned(num, denom)) {
+ case 0: /* q = 1, r = 0 */
+ q = allocate_tmp_bignum(C_fix(1), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ *(C_bignum_digits(q)) = 1;
+ r = allocate_tmp_bignum(C_fix(0), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ break;
+
+ case -1: /* q = 0, r = num */
+ q = allocate_tmp_bignum(C_fix(0), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ len = C_bignum_size(num) + 1; /* Ensure we can shift left by one */
+ r = allocate_tmp_bignum(C_fix(len), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_digits_destructive_copy(r, num);
+ d = C_bignum_digits(r);
+ d[len-1] = 0; /* Initialize most significant digit */
+ bignum_digits_destructive_shift_left(d, d + len, 1);
+ break;
+
+ case 1:
+ default:
+ len = C_bignum_size(num) + 1 - C_bignum_size(denom);
+ q = allocate_tmp_bignum(C_fix(len), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ len = C_bignum_size(num) + 1; /* LEN */
+ r = allocate_tmp_bignum(C_fix(len), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(num, denom, q, r, C_SCHEME_TRUE);
+ d = C_bignum_digits(r);
+ /* There should always be room to shift left by 1 because of LEN */
+ assert(C_ilen(d[len-1]) < C_BIGNUM_DIGIT_LENGTH);
+ bignum_digits_destructive_shift_left(d, d + len, 1);
+ break;
+ }
+
+ /* Now q is the quotient, but to "round" result we need to
+ * adjust. This follows the semantics of the "round" procedure:
+ * Round away from zero on positive numbers (this is never
+ * negative). In case of exactly halfway, we round up if odd.
+ */
+ fraction = C_bignum_to_double(q);
+ switch (basic_cmp(C_bignum_simplify(r), C_bignum_simplify(denom), "", 0)) {
+ case C_fix(0):
+ if (*(C_bignum_digits(q)) & 1) fraction += 1.0;
+ break;
+ case C_fix(1):
+ fraction += 1.0;
+ break;
+ default: /* if r <= denom, we're done */ break;
+ }
+
+ free_tmp_bignum(q);
+ free_tmp_bignum(r);
+ if (tmp1) free_tmp_bignum(tmp1);
+ if (tmp2) free_tmp_bignum(tmp2);
+ if (tmp3) free_tmp_bignum(tmp3);
+
+ shift_amount = nmin(DBL_MANT_DIG-1, e - (DBL_MIN_EXP - DBL_MANT_DIG));
+ res = ldexp(fraction, e - shift_amount);
+ return C_flonum(ptr, C_truep(negp) ? -res : res);
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "exact->inexact", n);
+ }
}
diff --git a/types.db b/types.db
index dbe39e66..58e938bd 100644
--- a/types.db
+++ b/types.db
@@ -532,7 +532,8 @@
(exact->inexact (#(procedure #:clean #:enforce #:foldable) exact->inexact (number) (or float cplxnum))
((float) (float) #(1))
- ((fixnum) (float) (##core#inline_allocate ("C_a_i_fix_to_flo" 4) #(1))))
+ ((fixnum) (float) (##core#inline_allocate ("C_a_i_fix_to_flo" 4) #(1)))
+ ((number) (##core#inline_allocate ("C_a_i_exact_to_inexact" 12) #(1))))
(inexact->exact (#(procedure #:clean #:enforce #:foldable) inexact->exact (number) (or integer ratnum))
((fixnum) (fixnum) #(1))
@@ -811,8 +812,10 @@
(##core#inline_allocate ("C_a_i_fix_to_flo" 4) #(1))))
((cplxnum) (##core#inline_allocate
("C_a_i_flonum_atan2" 4)
- (##sys#exact->inexact (##sys#slot #(1) '2))
- (##sys#exact->inexact (##sys#slot #(1) '1)))))
+ (##core#inline_allocate ("C_a_i_exact_to_inexact" 12)
+ (##sys#slot #(1) '2))
+ (##core#inline_allocate ("C_a_i_exact_to_inexact" 12)
+ (##sys#slot #(1) '1)))))
(numerator (#(procedure #:clean #:enforce #:foldable) numerator ((or float integer ratnum)) (or float integer))
((fixnum) (fixnum) #(1))
Trap