~ chicken-core (chicken-5) 404627f287754ff03ae16a2aa6cb9b0cfa203632
commit 404627f287754ff03ae16a2aa6cb9b0cfa203632
Author: Peter Bex <peter@more-magic.net>
AuthorDate: Sat Mar 28 19:28:56 2015 +0100
Commit: Peter Bex <peter@more-magic.net>
CommitDate: Sun May 31 14:55:24 2015 +0200
Make generic dyadic * inlineable! Restore old-style compiler specialization rewrites for dyadic *.
diff --git a/c-platform.scm b/c-platform.scm
index b6091cbe..b9b5b05a 100644
--- a/c-platform.scm
+++ b/c-platform.scm
@@ -640,6 +640,7 @@
(rewrite '+ 16 2 "C_s_a_i_plus" #t 36)
(rewrite '- 16 2 "C_s_a_i_minus" #t 36)
+(rewrite '* 16 2 "C_s_a_i_times" #t 40)
(rewrite '= 17 2 "C_i_nequalp")
(rewrite '> 17 2 "C_i_greaterp")
diff --git a/chicken.h b/chicken.h
index 241cda78..52e802f1 100644
--- a/chicken.h
+++ b/chicken.h
@@ -1958,9 +1958,7 @@ C_fctexport void C_ccall C_values(C_word c, C_word closure, C_word k, ...) C_nor
C_fctexport void C_ccall C_apply_values(C_word c, C_word closure, C_word k, C_word lst) C_noret;
C_fctexport void C_ccall C_call_with_values(C_word c, C_word closure, C_word k, C_word thunk, C_word kont) C_noret;
C_fctexport void C_ccall C_u_call_with_values(C_word c, C_word closure, C_word k, C_word thunk, C_word kont) C_noret;
-/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
C_fctexport void C_ccall C_times(C_word c, C_word closure, C_word k, ...) C_noret;
-C_fctexport void C_ccall C_2_basic_times(C_word c, C_word self, C_word k, C_word x, C_word y) C_noret;
C_fctexport void C_ccall C_plus(C_word c, C_word closure, C_word k, ...) C_noret;
C_fctexport void C_ccall C_minus(C_word c, C_word closure, C_word k, C_word n1, ...) C_noret;
/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
@@ -2186,6 +2184,7 @@ C_fctexport C_word C_fcall C_s_a_u_i_integer_negate(C_word **ptr, C_word n, C_wo
C_fctexport C_word C_fcall C_s_a_u_i_integer_minus(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
C_fctexport C_word C_fcall C_s_a_i_plus(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
C_fctexport C_word C_fcall C_s_a_u_i_integer_plus(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
+C_fctexport C_word C_fcall C_s_a_i_times(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
C_fctexport C_word C_fcall C_s_a_u_i_integer_times(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
C_fctexport C_word C_fcall C_s_a_i_arithmetic_shift(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
C_fctexport C_word C_fcall C_s_a_u_i_integer_gcd(C_word **ptr, C_word n, C_word x, C_word y) C_regparm;
diff --git a/library.scm b/library.scm
index 4fd9045a..27ddfdd6 100644
--- a/library.scm
+++ b/library.scm
@@ -995,8 +995,8 @@ EOF
(##sys#check-real phi 'make-polar)
(let ((fphi (exact->inexact phi)))
(make-complex
- (##sys#*-2 r (##core#inline_allocate ("C_a_i_cos" 4) fphi))
- (##sys#*-2 r (##core#inline_allocate ("C_a_i_sin" 4) fphi)))))
+ (* r (##core#inline_allocate ("C_a_i_cos" 4) fphi))
+ (* r (##core#inline_allocate ("C_a_i_sin" 4) fphi)))))
(define (real-part x)
(cond ((cplxnum? x) (%cplxnum-real x))
@@ -1019,7 +1019,7 @@ EOF
(cond ((cplxnum? x)
(let ((r (%cplxnum-real x))
(i (%cplxnum-imag x)) )
- (sqrt (+ (##sys#*-2 r r) (##sys#*-2 i i))) ))
+ (sqrt (+ (* r r) (* i i))) ))
((number? x) (abs x))
(else (##sys#error-bad-number x 'magnitude))))
@@ -1076,14 +1076,14 @@ EOF
(define (deliver y d)
(let* ((q (##sys#integer-power 2 (float-fraction-length y)))
- (scaled-y (##sys#*-2 y (exact->inexact q))))
+ (scaled-y (* y (exact->inexact q))))
(if (finite? scaled-y) ; Shouldn't this always be true?
(##sys#/-2 (##sys#/-2 ((##core#primitive "C_u_flo_to_int") scaled-y) q) d)
(##sys#error-bad-inexact x 'inexact->exact))))
(if (and (fp< x 1.0) ; Watch out for denormalized numbers
(fp> x -1.0)) ; XXX: Needs a test, it seems pointless
- (deliver (##sys#*-2 x (expt 2.0 flonum-precision))
+ (deliver (* x (expt 2.0 flonum-precision))
;; Can be bignum (is on 32-bit), so must wait until after init.
;; We shouldn't need to calculate this every single time, tho..
(##sys#integer-power 2 flonum-precision))
@@ -1112,8 +1112,6 @@ EOF
(define (abs x) (##core#inline_allocate ("C_s_a_i_abs" 10) x))
-(define ##sys#*-2 (##core#primitive "C_2_basic_times"))
-
(define (* . args)
(if (null? args)
1
@@ -1122,45 +1120,14 @@ EOF
(if (null? args)
(if (number? x) x (##sys#error-bad-number x '*))
(let loop ((args (##sys#slot args 1))
- (x (##sys#*-2 x (##sys#slot args 0))))
+ (x (##core#inline_allocate
+ ("C_s_a_i_times" 40) x (##sys#slot args 0))))
(if (null? args)
x
(loop (##sys#slot args 1)
- (##sys#*-2 x (##sys#slot args 0))) ) ) ) ) ) )
-
-(define (##sys#extended-times x y)
- (define (nonrat*rat x y)
- ;; a/b * c/d = a*c / b*d [with b = 1]
- ;; = ((a / g) * c) / (d / g)
- ;; With g = gcd(a, d) and a = x [Knuth, 4.5.1]
- (let* ((d (%ratnum-denominator y))
- (g (##sys#internal-gcd '* x d)))
- (ratnum (##sys#*-2 (quotient x g) (%ratnum-numerator y))
- (quotient d g))))
-
- (cond ((or (cplxnum? x) (cplxnum? y))
- (let* ((a (real-part x)) (b (imag-part x))
- (c (real-part y)) (d (imag-part y))
- (r (- (##sys#*-2 a c) (##sys#*-2 b d)))
- (i (+ (##sys#*-2 a d) (##sys#*-2 b c))) )
- (make-complex r i) ) )
- ((or (##core#inline "C_i_flonump" x) (##core#inline "C_i_flonump" y))
- ;; This may be incorrect when one is a ratnum consisting of bignums
- (fp* (exact->inexact y) (exact->inexact x))) ; loc?
- ((ratnum? x)
- (if (ratnum? y)
- ;; a/b * c/d = a*c / b*d [generic]
- ;; = ((a / g1) * (c / g2)) / ((b / g2) * (d / g1))
- ;; With g1 = gcd(a, d) and g2 = gcd(b, c) [Knuth, 4.5.1]
- (let* ((a (%ratnum-numerator x)) (b (%ratnum-denominator x))
- (c (%ratnum-numerator y)) (d (%ratnum-denominator y))
- (g1 (##sys#integer-gcd a d))
- (g2 (##sys#integer-gcd b c)))
- (ratnum (##sys#*-2 (quotient a g1) (quotient c g2))
- (##sys#*-2 (quotient b g2) (quotient d g1))))
- (nonrat*rat y x)))
- ((ratnum? y) (nonrat*rat x y))
- (else (##sys#error-bad-number x '*))))
+ (##core#inline_allocate
+ ("C_s_a_i_times" 40)
+ x (##sys#slot args 0))) ) ) ) ) ) )
(define (/ arg1 . args)
(if (null? args)
@@ -1182,9 +1149,9 @@ EOF
((or (cplxnum? x) (cplxnum? y))
(let* ((a (real-part x)) (b (imag-part x))
(c (real-part y)) (d (imag-part y))
- (r (+ (##sys#*-2 c c) (##sys#*-2 d d)))
- (x (##sys#/-2 (+ (##sys#*-2 a c) (##sys#*-2 b d)) r))
- (y (##sys#/-2 (- (##sys#*-2 b c) (##sys#*-2 a d)) r)) )
+ (r (+ (* c c) (* d d)))
+ (x (##sys#/-2 (+ (* a c) (* b d)) r))
+ (y (##sys#/-2 (- (* b c) (* a d)) r)) )
(make-complex x y) ))
((or (##core#inline "C_i_flonump" x) (##core#inline "C_i_flonump" y))
;; This may be incorrect when one is a ratnum consisting of bignums
@@ -1198,15 +1165,15 @@ EOF
(c (%ratnum-numerator y)) (d (%ratnum-denominator y))
(g1 (##sys#integer-gcd a c))
(g2 (##sys#integer-gcd b d)))
- (ratnum (##sys#*-2 (quotient a g1) (quotient d g2))
- (##sys#*-2 (quotient b g2) (quotient c g1))))
+ (ratnum (* (quotient a g1) (quotient d g2))
+ (* (quotient b g2) (quotient c g1))))
;; a/b / c/d = a*d / b*c [with d = 1]
;; = ((a / g) * sign(a)) / abs(b * (c / g))
;; With g = gcd(a, c) and c = y [Knuth, 4.5.1 ex. 4]
(let* ((a (%ratnum-numerator x))
(g (##sys#internal-gcd '/ a y))
(num (quotient a g))
- (denom (##sys#*-2 (%ratnum-denominator x) (quotient y g))))
+ (denom (* (%ratnum-denominator x) (quotient y g))))
(if (##core#inline "C_i_flonump" denom)
(##sys#/-2 num denom)
(ratnum num denom)))))
@@ -1216,7 +1183,7 @@ EOF
;; With g1 = gcd(a, c) and a = x [Knuth, 4.5.1 ex. 4]
(let* ((c (%ratnum-numerator y))
(g (##sys#internal-gcd '/ x c))
- (num (##sys#*-2 (quotient x g) (%ratnum-denominator y)))
+ (num (* (quotient x g) (%ratnum-denominator y)))
(denom (quotient c g)))
(if (##core#inline "C_i_flonump" denom)
(##sys#/-2 num denom)
@@ -1258,7 +1225,7 @@ EOF
(values (- (arithmetic-shift 1 base*n) 1) ; B^n-1
(+ (- a12 (arithmetic-shift b1 base*n)) b1))))
(let ((r1a3 (+ (arithmetic-shift r1 (digit-bits n)) a3)))
- (let lp ((r^ (- r1a3 (##sys#*-2 q^ b2)))
+ (let lp ((r^ (- r1a3 (* q^ b2)))
(q^ q^))
(if (negative? r^)
(lp (+ r^ b) (- q^ 1))
@@ -1365,7 +1332,7 @@ EOF
((= fx fy)
(let ((rat (sr (##sys#/-2 1 (- y fy))
(##sys#/-2 1 (- x fx)))))
- (list (+ (cadr rat) (##sys#*-2 fx (car rat)))
+ (list (+ (cadr rat) (* fx (car rat)))
(car rat))))
(else (list (+ 1 fx) 1)))))
(cond ((< y x) (find-ratio-between y x))
@@ -1434,14 +1401,12 @@ EOF
(define (exp n)
(##sys#check-number n 'exp)
(if (cplxnum? n)
- (##sys#*-2 (##core#inline_allocate ("C_a_i_exp" 4)
- (exact->inexact (%cplxnum-real n)))
- (let ((p (%cplxnum-imag n)))
- (make-complex
- (##core#inline_allocate
- ("C_a_i_cos" 4) (exact->inexact p))
- (##core#inline_allocate
- ("C_a_i_sin" 4) (exact->inexact p)) ) ) )
+ (* (##core#inline_allocate ("C_a_i_exp" 4)
+ (exact->inexact (%cplxnum-real n)))
+ (let ((p (%cplxnum-imag n)))
+ (make-complex
+ (##core#inline_allocate ("C_a_i_cos" 4) (exact->inexact p))
+ (##core#inline_allocate ("C_a_i_sin" 4) (exact->inexact p)) ) ) )
(##core#inline_allocate ("C_a_i_flonum_exp" 4) (exact->inexact n))))
(define (##sys#log-1 x) ; log_e(x)
@@ -1451,7 +1416,7 @@ EOF
;; avoid calling inexact->exact on X here (to avoid overflow?)
((or (cplxnum? x) (negative? x)) ; General case
(+ (##sys#log-1 (magnitude x))
- (##sys#*-2 (make-complex 0 1) (angle x))))
+ (* (make-complex 0 1) (angle x))))
(else ; Real number case (< already ensured the argument type is a number)
(##core#inline_allocate ("C_a_i_log" 4) (exact->inexact x)))))
@@ -1467,14 +1432,14 @@ EOF
(define (sin n)
(##sys#check-number n 'sin)
(if (cplxnum? n)
- (let ((in (##sys#*-2 %i n)))
+ (let ((in (* %i n)))
(##sys#/-2 (- (exp in) (exp (- in))) %i2))
(##core#inline_allocate ("C_a_i_sin" 4) (exact->inexact n))))
(define (cos n)
(##sys#check-number n 'cos)
(if (cplxnum? n)
- (let ((in (##sys#*-2 %i n)))
+ (let ((in (* %i n)))
(##sys#/-2 (+ (exp in) (exp (- in))) 2) )
(##core#inline_allocate ("C_a_i_cos" 4) (exact->inexact n))))
@@ -1494,10 +1459,9 @@ EOF
(##core#inline_allocate
("C_a_i_fix_to_flo" 4) n)))
;; General definition can return compnums
- (else (##sys#*-2 %-i
- (##sys#log-1
- (+ (##sys#*-2 %i n)
- (##sys#sqrt/loc 'asin (- 1 (##sys#*-2 n n)))))))))
+ (else (* %-i (##sys#log-1
+ (+ (* %i n)
+ (##sys#sqrt/loc 'asin (- 1 (* n n)))))))))
;; General case:
;; cos^{-1}(z) = 1/2\pi + i\ln(iz + \sqrt{1-z^2}) = 1/2\pi - sin^{-1}(z) = sin(1) - sin(z)
@@ -1519,7 +1483,7 @@ EOF
(cond ((cplxnum? n)
(if b
(##sys#error-bad-real n 'atan)
- (let ((in (##sys#*-2 %i n)))
+ (let ((in (* %i n)))
(##sys#/-2 (- (##sys#log-1 (+ 1 in))
(##sys#log-1 (- 1 in))) %i2))))
(b
@@ -1553,7 +1517,7 @@ EOF
(+ (arithmetic-shift r^ len/4) a1)
(arithmetic-shift s^ 1)))
((s) (+ (arithmetic-shift s^ len/4) q))
- ((r) (+ (arithmetic-shift u len/4) (- a0 (##sys#*-2 q q)))))
+ ((r) (+ (arithmetic-shift u len/4) (- a0 (* q q)))))
(if (negative? r)
(values (- s 1)
(- (+ r (arithmetic-shift s 1)) 1))
@@ -1569,7 +1533,7 @@ EOF
(cond ((cplxnum? n) ; Must be checked before we call "negative?"
(let ((p (##sys#/-2 (angle n) 2))
(m (##core#inline_allocate ("C_a_i_sqrt" 4) (magnitude n))) )
- (make-complex (##sys#*-2 m (cos p)) (##sys#*-2 m (sin p)) ) ))
+ (make-complex (* m (cos p)) (* m (sin p)) ) ))
((negative? n)
(make-complex .0 (##core#inline_allocate
("C_a_i_sqrt" 4) (exact->inexact (- n)))))
@@ -1610,44 +1574,44 @@ EOF
(n-1 (- n 1)))
(let lp ((g0 g0)
(g1 (quotient
- (+ (##sys#*-2 n-1 g0)
+ (+ (* n-1 g0)
(quotient k (##sys#integer-power g0 n-1)))
n)))
(if (< g1 g0)
(lp g1 (quotient
- (+ (##sys#*-2 n-1 g1)
+ (+ (* n-1 g1)
(quotient k (##sys#integer-power g1 n-1)))
n))
(values g0 (- k (##sys#integer-power g0 n))))))))))
(define (##sys#integer-power base e)
- (define (square x) (##sys#*-2 x x))
+ (define (square x) (* x x))
(if (negative? e)
(##sys#/-2 1 (##sys#integer-power base (integer-negate e)))
(let lp ((res 1) (e2 e))
(cond
((eq? e2 0) res)
((even? e2) ; recursion is faster than iteration here
- (##sys#*-2 res (square (lp 1 (arithmetic-shift e2 -1)))))
+ (* res (square (lp 1 (arithmetic-shift e2 -1)))))
(else
- (lp (##sys#*-2 res base) (- e2 1)))))))
+ (lp (* res base) (- e2 1)))))))
(define (expt a b)
(define (log-expt a b)
- (exp (##sys#*-2 b (##sys#log-1 a))))
+ (exp (* b (##sys#log-1 a))))
(define (slow-expt a b)
(if (eq? 0 a)
(##sys#signal-hook
#:arithmetic-error 'expt
"exponent of exact 0 with complex argument is undefined" a b)
- (exp (##sys#*-2 b (##sys#log-1 a)))))
+ (exp (* b (##sys#log-1 a)))))
(cond ((not (number? a)) (##sys#error-bad-number a 'expt))
((not (number? b)) (##sys#error-bad-number b 'expt))
((and (ratnum? a) (not (inexact? b)))
;; (n*d)^b = n^b * d^b = n^b * x^{-b} | x = 1/b
;; Hopefully faster than integer-power
- (##sys#*-2 (expt (%ratnum-numerator a) b)
- (expt (%ratnum-denominator a) (- b))))
+ (* (expt (%ratnum-numerator a) b)
+ (expt (%ratnum-denominator a) (- b))))
((ratnum? b)
;; x^{a/b} = (x^{1/b})^a
(cond
@@ -1728,7 +1692,7 @@ EOF
(if (integer? head) (abs head) (##sys#error-bad-integer head 'lcm))
(let* ((n2 (##sys#slot next 0))
(gcd (##sys#internal-gcd 'lcm head n2)))
- (loop (quotient (##sys#*-2 head n2) gcd)
+ (loop (quotient (* head n2) gcd)
(##sys#slot next 1)) ) ) ) ) )
;; This simple enough idea is from
@@ -1774,7 +1738,7 @@ EOF
;; by Aubrey Jaffer.
(define (mantexp->dbl mant point)
(if (not (negative? point))
- (exact->inexact (##sys#*-2 mant (##sys#integer-power 10 point)))
+ (exact->inexact (* mant (##sys#integer-power 10 point)))
(let* ((scl (##sys#integer-power 10 (abs point)))
(bex (fx- (fx- (integer-length mant) (integer-length scl))
flonum-precision)))
@@ -1784,17 +1748,17 @@ EOF
(cond ((> (integer-length quo) flonum-precision)
;; Too many bits of quotient; readjust
(set! bex (fx+ 1 bex))
- (set! quo (round-quotient num (##sys#*-2 scl 2)))))
+ (set! quo (round-quotient num (* scl 2)))))
(ldexp (exact->inexact quo) bex))
;; Fall back to exact calculation in extreme cases
- (##sys#*-2 mant (##sys#integer-power 10 point))))))
+ (* mant (##sys#integer-power 10 point))))))
(define ldexp (foreign-lambda double "ldexp" double int))
;; Should we export this?
(define (round-quotient n d)
(let ((q (##sys#integer-quotient n d)))
- (if ((if (even? q) > >=) (##sys#*-2 (abs (remainder n d)) 2) (abs d))
+ (if ((if (even? q) > >=) (* (abs (remainder n d)) 2) (abs d))
(+ q (if (eqv? (negative? n) (negative? d)) 1 -1))
q)))
@@ -1815,7 +1779,7 @@ EOF
((< e (fxneg +maximum-allowed-exponent+))
(and (eq? exactness 'i) +0.0))
((eq? exactness 'i) (mantexp->dbl value e))
- (else (##sys#*-2 value (##sys#integer-power 10 e))))))
+ (else (* value (##sys#integer-power 10 e))))))
(define (make-nan)
;; Return fresh NaNs, so eqv? returns #f on two read NaNs. This
;; is not mandated by the standard, but compatible with earlier
diff --git a/runtime.c b/runtime.c
index e01d26c4..65b92b0e 100644
--- a/runtime.c
+++ b/runtime.c
@@ -518,6 +518,9 @@ static C_word bignum_plus_unsigned(C_word **ptr, C_word x, C_word y, C_word negp
static C_word rat_plusmin_integer(C_word **ptr, C_word rat, C_word i, integer_plusmin_op plusmin_op);
static C_word integer_minus_rat(C_word **ptr, C_word i, C_word rat);
static C_word rat_plusmin_rat(C_word **ptr, C_word x, C_word y, integer_plusmin_op plusmin_op);
+static C_word rat_times_integer(C_word **ptr, C_word x, C_word y);
+static C_word rat_times_rat(C_word **ptr, C_word x, C_word y);
+static C_word cplx_times(C_word **ptr, C_word rx, C_word ix, C_word ry, C_word iy);
static C_word bignum_minus_unsigned(C_word **ptr, C_word x, C_word y);
static C_regparm void basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_r, C_word return_q) C_noret;
static C_regparm void integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) C_noret;
@@ -586,7 +589,6 @@ static C_PTABLE_ENTRY *create_initial_ptable();
static void dload_2(void *dummy) C_noret;
#endif
-
static void
C_dbg(C_char *prefix, C_char *fstr, ...)
{
@@ -843,7 +845,7 @@ static C_PTABLE_ENTRY *create_initial_ptable()
{
/* IMPORTANT: hardcoded table size -
this must match the number of C_pte calls + 1 (NULL terminator)! */
- C_PTABLE_ENTRY *pt = (C_PTABLE_ENTRY *)C_malloc(sizeof(C_PTABLE_ENTRY) * 72);
+ C_PTABLE_ENTRY *pt = (C_PTABLE_ENTRY *)C_malloc(sizeof(C_PTABLE_ENTRY) * 71);
int i = 0;
if(pt == NULL)
@@ -867,7 +869,6 @@ static C_PTABLE_ENTRY *create_initial_ptable()
C_pte(C_set_dlopen_flags);
C_pte(C_become);
C_pte(C_apply_values);
- /* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
C_pte(C_times);
C_pte(C_minus);
C_pte(C_plus);
@@ -912,7 +913,6 @@ static C_PTABLE_ENTRY *create_initial_ptable()
C_pte(C_flonum_to_string);
/* IMPORTANT: have you read the comments at the start and the end of this function? */
C_pte(C_signum);
- C_pte(C_2_basic_times);
C_pte(C_basic_quotient);
C_pte(C_basic_remainder);
C_pte(C_basic_divrem);
@@ -7382,60 +7382,342 @@ void C_ccall values_continuation(C_word c, C_word closure, C_word arg0, ...)
C_do_apply(n - 2, kont, k);
}
+static C_word rat_times_integer(C_word **ptr, C_word rat, C_word i)
+{
+ C_word ab[C_SIZEOF_FIX_BIGNUM*7], *a = ab,
+ num, denom, d_big, gcd, gcd_big, i_big, tmp_r, a_div_g, size;
-void C_ccall
-C_2_basic_times(C_word c, C_word self, C_word k, C_word x, C_word y)
+ switch (i) {
+ case C_fix(0): return C_fix(0);
+ case C_fix(1): return rat;
+ case C_fix(-1):
+ num = C_s_a_u_i_integer_negate(ptr, 1, C_block_item(rat, 1));
+ return C_ratnum(ptr, num , C_block_item(rat, 2));
+ /* default: CONTINUE BELOW */
+ }
+
+ num = C_block_item(rat, 1);
+ denom = C_block_item(rat, 2);
+
+ /* a/b * c/d = a*c / b*d [with b = 1] */
+ /* = ((a / g) * c) / (d / g) */
+ /* With g = gcd(a, d) and a = x [Knuth, 4.5.1] */
+ gcd = C_s_a_u_i_integer_gcd(&a, 2, i, denom);
+
+ /* Ensure bignums to simplify logic */
+ gcd_big = (gcd & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, gcd) : gcd;
+ i_big = (i & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, i) : i;
+ d_big = (denom & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, denom) : denom;
+
+ /* Calculate a/g (= i/gcd), which will later be multiplied by y */
+ switch(bignum_cmp_unsigned(i_big, gcd_big)) {
+ case 0:
+ a_div_g = C_bignum_negativep(i_big) ? C_fix(-1) : C_fix(1);
+ break;
+ case -1:
+ clear_buffer_object(ab, gcd);
+ return C_fix(0); /* Save some work */
+ break;
+ case 1:
+ default:
+ size = C_bignum_size(i_big) + 1 - C_bignum_size(gcd_big);
+ a_div_g = C_allocate_scratch_bignum(
+ &a, C_fix(size), C_mk_bool(C_bignum_negativep(i_big)),
+ C_SCHEME_FALSE);
+ size = C_bignum_size(i_big) + 1;
+ tmp_r = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(i_big, gcd_big, a_div_g,
+ tmp_r, C_SCHEME_FALSE);
+ free_tmp_bignum(tmp_r);
+ a_div_g = C_bignum_simplify(a_div_g);
+ }
+
+ /* Final numerator = a/g * c (= a_div_g * num) */
+ num = C_s_a_u_i_integer_times(&a, 2, a_div_g, num);
+ num = move_buffer_object(ptr, ab, num);
+
+ /* Calculate d/g (= denom/gcd). We already know |d| >= |g| */
+ if(bignum_cmp_unsigned(d_big, gcd_big) == 0) {
+ denom = C_bignum_negativep(d_big) ? C_fix(-1) : C_fix(1);
+ } else {
+ size = C_bignum_size(d_big) + 1 - C_bignum_size(gcd_big);
+ denom = C_allocate_scratch_bignum(
+ &a, C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ size = C_bignum_size(d_big) + 1;
+ tmp_r = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(d_big, gcd_big, denom,
+ tmp_r, C_SCHEME_FALSE);
+ free_tmp_bignum(tmp_r);
+ denom = move_buffer_object(ptr, ab, C_bignum_simplify(denom));
+ }
+
+ clear_buffer_object(ab, gcd);
+ clear_buffer_object(ab, a_div_g);
+
+ if (denom == C_fix(1)) return num;
+ else return C_ratnum(ptr, num, denom);
+}
+
+/* This is truly wretched */
+static C_word rat_times_rat(C_word **ptr, C_word x, C_word y)
+{
+ C_word ab[C_SIZEOF_FIX_BIGNUM*14], *a = ab,
+ xnum, xdenom, ynum, ydenom, xnum_big, ynum_big, xdenom_big, ydenom_big,
+ g1, g1_big, g2, g2_big, tmp_r, num, denom,
+ a_div_g1, b_div_g2, c_div_g2, d_div_g1, size;
+
+ xnum = C_block_item(x, 1);
+ xdenom = C_block_item(x, 2);
+ ynum = C_block_item(y, 1);
+ ydenom = C_block_item(y, 2);
+
+ /* a/b * c/d = a*c / b*d [generic] */
+ /* = ((a / g1) * (c / g2)) / ((b / g2) * (d / g1)) */
+ /* With g1 = gcd(a, d) and g2 = gcd(b, c) [Knuth, 4.5.1] */
+ g1 = C_s_a_u_i_integer_gcd(&a, 2, xnum, ydenom);
+ g2 = C_s_a_u_i_integer_gcd(&a, 2, ynum, xdenom);
+
+ /* Ensure bignums to simplify logic */
+ g1_big = (g1 & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, g1) : g1;
+ g2_big = (g2 & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, g2) : g2;
+ xnum_big = (xnum & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, xnum) : xnum;
+ xdenom_big = (xdenom & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, xdenom) : xdenom;
+ ynum_big = (ynum & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, ynum) : ynum;
+ ydenom_big = (ydenom & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, ydenom) : ydenom;
+
+ /* Calculate a/g1 (= xnum/g1), which will later be multiplied by c/g2 */
+ if (bignum_cmp_unsigned(xnum_big, g1_big) == 0) {
+ a_div_g1 = C_bignum_negativep(xnum_big) ? C_fix(-1) : C_fix(1);
+ } else { /* We know |xnum| >= |g1| */
+ size = C_bignum_size(xnum_big) + 1 - C_bignum_size(g1_big);
+ a_div_g1 = C_allocate_scratch_bignum(
+ &a, C_fix(size), C_mk_bool(C_bignum_negativep(xnum_big)),
+ C_SCHEME_FALSE);
+ size = C_bignum_size(xnum_big) + 1;
+ tmp_r = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(xnum_big, g1_big, a_div_g1,
+ tmp_r, C_SCHEME_FALSE);
+ free_tmp_bignum(tmp_r);
+ a_div_g1 = C_bignum_simplify(a_div_g1);
+ }
+
+ /* Calculate c/g2 (= ynum/g2), which will later be multiplied by a/g1 */
+ if (bignum_cmp_unsigned(ynum_big, g2_big) == 0) {
+ c_div_g2 = C_bignum_negativep(ynum_big) ? C_fix(-1) : C_fix(1);
+ } else { /* We know |ynum| >= |g2| */
+ size = C_bignum_size(ynum_big) + 1 - C_bignum_size(g2_big);
+ c_div_g2 = C_allocate_scratch_bignum(
+ &a, C_fix(size), C_mk_bool(C_bignum_negativep(ynum_big)),
+ C_SCHEME_FALSE);
+ size = C_bignum_size(ynum_big) + 1;
+ tmp_r = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(ynum_big, g2_big, c_div_g2,
+ tmp_r, C_SCHEME_FALSE);
+ free_tmp_bignum(tmp_r);
+ c_div_g2 = C_bignum_simplify(c_div_g2);
+ }
+
+ /* Final numerator = a/g1 * c/g2 */
+ num = C_s_a_u_i_integer_times(&a, 2, a_div_g1, c_div_g2);
+ num = move_buffer_object(ptr, ab, num);
+
+ /* Now, do the same for the denominator.... */
+
+ /* Calculate b/g2 (= xdenom/g2), which will later be multiplied by d/g1 */
+ if (bignum_cmp_unsigned(xdenom_big, g2_big) == 0) {
+ b_div_g2 = C_bignum_negativep(xdenom_big) ? C_fix(-1) : C_fix(1);
+ } else { /* We know |xdenom| >= |g2| */
+ size = C_bignum_size(xdenom_big) + 1 - C_bignum_size(g2_big);
+ b_div_g2 = C_allocate_scratch_bignum(
+ &a, C_fix(size), C_mk_bool(C_bignum_negativep(xdenom_big)),
+ C_SCHEME_FALSE);
+ size = C_bignum_size(xdenom_big) + 1;
+ tmp_r = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(xdenom_big, g2_big, b_div_g2,
+ tmp_r, C_SCHEME_FALSE);
+ free_tmp_bignum(tmp_r);
+ b_div_g2 = C_bignum_simplify(b_div_g2);
+ }
+
+ /* Calculate d/g1 (= ydenom/g1), which will later be multiplied by b/g2 */
+ if (bignum_cmp_unsigned(ydenom_big, g1_big) == 0) {
+ d_div_g1 = C_bignum_negativep(ydenom_big) ? C_fix(-1) : C_fix(1);
+ } else { /* We know |ydenom| >= |g1| */
+ size = C_bignum_size(ydenom_big) + 1 - C_bignum_size(g1_big);
+ d_div_g1 = C_allocate_scratch_bignum(
+ &a, C_fix(size), C_mk_bool(C_bignum_negativep(ydenom_big)),
+ C_SCHEME_FALSE);
+ size = C_bignum_size(ydenom_big) + 1;
+ tmp_r = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
+ bignum_destructive_divide_full(ydenom_big, g1_big, d_div_g1,
+ tmp_r, C_SCHEME_FALSE);
+ free_tmp_bignum(tmp_r);
+ d_div_g1 = C_bignum_simplify(d_div_g1);
+ }
+
+ /* Final denominator = b/g2 * d/g1 */
+ denom = C_s_a_u_i_integer_times(&a, 2, b_div_g2, d_div_g1);
+ denom = move_buffer_object(ptr, ab, denom);
+
+ clear_buffer_object(ab, g1);
+ clear_buffer_object(ab, g2);
+ clear_buffer_object(ab, a_div_g1);
+ clear_buffer_object(ab, b_div_g2);
+ clear_buffer_object(ab, c_div_g2);
+ clear_buffer_object(ab, d_div_g1);
+
+ if (denom == C_fix(1)) return num;
+ else return C_ratnum(ptr, num, denom);
+}
+
+static C_word
+cplx_times(C_word **ptr, C_word rx, C_word ix, C_word ry, C_word iy)
+{
+ /* Allocation here is kind of tricky: Each intermediate result can
+ * be at most a ratnum consisting of two bignums (2 digits), so
+ * C_SIZEOF_STRUCTURE(3) + C_SIZEOF_BIGNUM(2) = 11 words
+ */
+ C_word ab[(C_SIZEOF_STRUCTURE(3) + C_SIZEOF_BIGNUM(2))*6], *a = ab,
+ r1, r2, i1, i2, r, i;
+
+ /* a+bi * c+di = (a*c - b*d) + (a*d + b*c)i */
+ /* We call these: r1 = a*c, r2 = b*d, i1 = a*d, i2 = b*c */
+ r1 = C_s_a_i_times(&a, 2, rx, ry);
+ r2 = C_s_a_i_times(&a, 2, ix, iy);
+ i1 = C_s_a_i_times(&a, 2, rx, iy);
+ i2 = C_s_a_i_times(&a, 2, ix, ry);
+
+ r = C_s_a_i_minus(ptr, 2, r1, r2);
+ i = C_s_a_i_plus(ptr, 2, i1, i2);
+
+ r = move_buffer_object(ptr, ab, r);
+ i = move_buffer_object(ptr, ab, i);
+
+ clear_buffer_object(ab, r1);
+ clear_buffer_object(ab, r2);
+ clear_buffer_object(ab, i1);
+ clear_buffer_object(ab, i2);
+
+ if (C_truep(C_u_i_zerop(i))) return r;
+ else return C_cplxnum(ptr, r, i);
+}
+
+/* The maximum size this needs is that required to store a complex
+ * number result, where both real and imag parts consist of ratnums.
+ * The maximum size of those ratnums is if they consist of two bignums
+ * from a fixnum multiplication (2 digits each), so we're looking at
+ * C_SIZEOF_STRUCTURE(3) * 3 + C_SIZEOF_BIGNUM(2) * 4 = 40 words!
+ */
+C_regparm C_word C_fcall
+C_s_a_i_times(C_word **ptr, C_word n, C_word x, C_word y)
{
if (x & C_FIXNUM_BIT) {
if (y & C_FIXNUM_BIT) {
- C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
- C_kontinue(k, C_a_i_fixnum_times(&a, 2, x, y));
+ return C_a_i_fixnum_times(ptr, 2, x, y);
} else if (C_immediatep(y)) {
barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
} else if (C_block_header(y) == C_FLONUM_TAG) {
- C_word *a = C_alloc(C_SIZEOF_FLONUM);
- C_kontinue(k, C_flonum(&a, (double)C_unfix(x) * C_flonum_magnitude(y)));
+ return C_flonum(ptr, (double)C_unfix(x) * C_flonum_magnitude(y));
} else if (C_truep(C_bignump(y))) {
- C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
- C_kontinue(k, C_s_a_u_i_integer_times(&a, 2, x, y));
+ return C_s_a_u_i_integer_times(ptr, 2, x, y);
+ } else if (C_block_header(y) == C_STRUCTURE3_TAG) {
+ if (C_block_item(y, 0) == C_ratnum_type_tag) {
+ return rat_times_integer(ptr, y, x);
+ } else if (C_block_item(y, 0) == C_cplxnum_type_tag) {
+ return cplx_times(ptr, x, C_fix(0),
+ C_block_item(y, 1), C_block_item(y, 2));
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ }
} else {
- try_extended_number("\003sysextended-times", 3, k, x, y);
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
}
} else if (C_immediatep(x)) {
barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
} else if (C_block_header(x) == C_FLONUM_TAG) {
- C_word *a = C_alloc(C_SIZEOF_FLONUM);
if (y & C_FIXNUM_BIT) {
- C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) * (double)C_unfix(y)));
+ return C_flonum(ptr, C_flonum_magnitude(x) * (double)C_unfix(y));
} else if (C_immediatep(y)) {
barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
} else if (C_block_header(y) == C_FLONUM_TAG) {
- C_kontinue(k, C_a_i_flonum_times(&a, 2, x, y));
+ return C_a_i_flonum_times(ptr, 2, x, y);
} else if (C_truep(C_bignump(y))) {
- C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)*C_bignum_to_double(y)));
+ return C_flonum(ptr, C_flonum_magnitude(x) * C_bignum_to_double(y));
+ } else if (C_block_header(y) == C_STRUCTURE3_TAG) {
+ if (C_block_item(y, 0) == C_ratnum_type_tag) {
+ return C_s_a_i_times(ptr, 2, x, C_a_i_exact_to_inexact(ptr, 1, y));
+ } else if (C_block_item(y, 0) == C_cplxnum_type_tag) {
+ C_word ab[C_SIZEOF_FLONUM], *a = ab;
+ return cplx_times(ptr, x, C_flonum(&a, 0.0),
+ C_block_item(y, 1), C_block_item(y, 2));
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ }
} else {
- try_extended_number("\003sysextended-times", 3, k, x, y);
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
}
} else if (C_truep(C_bignump(x))) {
if (y & C_FIXNUM_BIT) {
- C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
- C_kontinue(k, C_s_a_u_i_integer_times(&a, 2, x, y));
+ return C_s_a_u_i_integer_times(ptr, 2, x, y);
} else if (C_immediatep(y)) {
barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
} else if (C_block_header(y) == C_FLONUM_TAG) {
- C_word *a = C_alloc(C_SIZEOF_FLONUM);
- C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)*C_flonum_magnitude(y)));
+ return C_flonum(ptr, C_bignum_to_double(x) * C_flonum_magnitude(y));
} else if (C_truep(C_bignump(y))) {
- C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
- C_kontinue(k, C_s_a_u_i_integer_times(&a, 2, x, y));
+ return C_s_a_u_i_integer_times(ptr, 2, x, y);
+ } else if (C_block_header(y) == C_STRUCTURE3_TAG) {
+ if (C_block_item(y, 0) == C_ratnum_type_tag) {
+ return rat_times_integer(ptr, y, x);
+ } else if (C_block_item(y, 0) == C_cplxnum_type_tag) {
+ return cplx_times(ptr, x, C_fix(0),
+ C_block_item(y, 1), C_block_item(y, 2));
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ }
} else {
- try_extended_number("\003sysextended-times", 3, k, x, y);
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ }
+ } else if (C_block_header(x) == C_STRUCTURE3_TAG) {
+ if (C_block_item(x, 0) == C_ratnum_type_tag) {
+ if (y & C_FIXNUM_BIT) {
+ return rat_times_integer(ptr, x, y);
+ } else if (C_immediatep(y)) {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ } else if (C_block_header(y) == C_FLONUM_TAG) {
+ return C_s_a_i_times(ptr, 2, C_a_i_exact_to_inexact(ptr, 1, x), y);
+ } else if (C_truep(C_bignump(y))) {
+ return rat_times_integer(ptr, x, y);
+ } else if (C_block_header(y) == C_STRUCTURE3_TAG) {
+ if (C_block_item(y, 0) == C_ratnum_type_tag) {
+ return rat_times_rat(ptr, x, y);
+ } else if (C_block_item(y, 0) == C_cplxnum_type_tag) {
+ return cplx_times(ptr, x, C_fix(0),
+ C_block_item(y, 1),C_block_item(y, 2));
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ }
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
+ }
+ } else if (C_block_item(x, 0) == C_cplxnum_type_tag) {
+ if (!C_immediatep(y) && C_block_header(y) == C_STRUCTURE3_TAG &&
+ C_block_item(y, 0) == C_cplxnum_type_tag) {
+ return cplx_times(ptr, C_block_item(x, 1), C_block_item(x, 2),
+ C_block_item(y, 1), C_block_item(y, 2));
+ } else {
+ C_word ab[C_SIZEOF_FLONUM], *a = ab, yi;
+ yi = C_truep(C_i_flonump(y)) ? C_flonum(&a,0) : C_fix(0);
+ return cplx_times(ptr, C_block_item(x, 1), C_block_item(x, 2), y, yi);
+ }
+ } else {
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
}
} else {
- try_extended_number("\003sysextended-times", 3, k, x, y);
+ barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
}
}
+
C_regparm C_word C_fcall
C_s_a_u_i_integer_times(C_word **ptr, C_word n, C_word x, C_word y)
{
@@ -7567,7 +7849,6 @@ bignum_times_bignum_karatsuba(C_word **ptr, C_word x, C_word y, C_word negp)
}
-/* XXX TODO OBSOLETE: This can be removed after recompiling c-platform.scm */
void C_ccall C_times(C_word c, C_word closure, C_word k, ...)
{
va_list v;
diff --git a/types.db b/types.db
index 10c1bc99..22ddaac8 100644
--- a/types.db
+++ b/types.db
@@ -374,7 +374,7 @@
((integer integer) (integer)
(##core#inline_allocate ("C_s_a_u_i_integer_times" 7) #(1) #(2)))
((* *) (number)
- (##sys#*-2 #(1) #(2))))
+ (##core#inline_allocate ("C_s_a_i_times" 40) #(1) #(2))))
(/ (#(procedure #:clean #:enforce #:foldable) / (number #!rest number) number)
((float fixnum) (float)
Trap