Even faster (truncate ratio ratio)

Do simpler divison arithmetic, avoding multiple GCDs.
This commit is contained in:
Stas Boukarev 2025-04-11 00:01:02 +03:00
parent 3b55917c64
commit 8b547bb076
2 changed files with 75 additions and 145 deletions

View file

@ -48,7 +48,7 @@
;;; our guess for the preferred order in which to do type tests
;;; (cheaper and/or more probable first.)
(defconstant-eqx +type-test-ordering+
'(fixnum single-float double-float integer #+long-float long-float
'((eql 1) fixnum single-float double-float integer #+long-float long-float
sb-vm:signed-word word bignum
complex ratio)
#'equal)
@ -56,8 +56,8 @@
;;; Should TYPE1 be tested before TYPE2?
(defun type-test-order (type1 type2)
(let ((o1 (position type1 +type-test-ordering+))
(o2 (position type2 +type-test-ordering+)))
(let ((o1 (position type1 +type-test-ordering+ :test #'equal))
(o2 (position type2 +type-test-ordering+ :test #'equal)))
(cond ((not o1) nil)
((not o2) t)
(t

View file

@ -312,6 +312,14 @@
(declare (type word x y))
(%multiply-high x y))
(defmacro dispatch-ratio ((ratio numerator denominator) &body body)
`(let ((,numerator (numerator ,ratio))
(,denominator (denominator ,ratio)))
(if (and (fixnump ,numerator)
(fixnump ,denominator))
(progn ,@body)
(progn ,@body))))
(defmacro dispatch-two-ratios ((numerator1-var numerator1)
(denominator1-var denominator1)
(numerator2-var numerator2)
@ -377,47 +385,6 @@
(declaim (maybe-inline truncate floor ceiling round fround))
;;; Do operations on ratios without consing intermediate ratios.
(defmacro with-ratio-ops (&body body)
`(flet ((integer*ratio (x y)
(let* ((ny (numerator y))
(dy (denominator y))
(gcd (gcd x dy)))
(if (eql gcd 1)
(values (* x ny) dy)
(let ((nn (* (truncate x gcd) ny))
(nd (truncate dy gcd)))
(values nn nd)))))
(truncate-ratio (number numerator denominator)
(if (= numerator 1)
(* number denominator)
(truncate (* number denominator)
numerator)))
(ratio-ratio (nx dx ny dy)
(if (= dy 1)
(%make-ratio (- nx (* ny dx)) dx)
(let ((g1 (gcd dx dy)))
(if (eql g1 1)
(%make-ratio (- (* nx dy) (* dx ny)) (* dx dy))
(let* ((t2 (truncate dx g1))
(t1 (- (* nx (truncate dy g1)) (* t2 ny)))
(g2 (gcd t1 g1)))
(cond ((eql t1 0) 0)
((eql g2 1)
(%make-ratio t1 (* t2 dy)))
(t
(let* ((nn (truncate t1 g2))
(t3 (truncate dy g2))
(nd
(if (eql t2 1)
t3
(* t2 t3))))
(if (eql nd 1)
nn
(%make-ratio nn nd)))))))))))
(declare (inline integer*ratio truncate-ratio ratio-ratio))
,@body))
(defun truncate (number &optional (divisor 1))
"Return number (or number/divisor) as an integer, rounded toward 0.
The second returned value is the remainder."
@ -450,27 +417,21 @@
(if (single-digit-bignum-p divisor)
(bignum-truncate-single-digit (make-small-bignum number) divisor)
(bignum-truncate (make-small-bignum number) divisor)))
((ratio integer)
(let ((numerator (numerator number))
(denominator (denominator number)))
(if (eql divisor 1)
(let ((q (truncate numerator denominator)))
(values q (- number q)))
(multiple-value-bind (q rem) (truncate numerator
(* denominator divisor))
(values q (%make-ratio rem denominator))))))
((ratio (foreach (eql 1) fixnum bignum))
(dispatch-ratio (number numerator denominator)
(multiple-value-bind (q rem) (truncate numerator
(* denominator divisor))
(values q (%make-ratio rem denominator)))))
((ratio ratio)
(let ((numerator (numerator number))
(denominator (denominator number)))
;; (let ((q (truncate (numerator number)
;; (* (denominator number) divisor))))
;; (values q (- number (* q divisor))))
(with-ratio-ops
(multiple-value-bind (div-num div-den) (integer*ratio denominator divisor)
(let ((q (truncate-ratio numerator
div-num div-den)))
(multiple-value-bind (rem-num rem-den) (integer*ratio q divisor)
(values q (ratio-ratio numerator denominator rem-num rem-den))))))))
(dispatch-two-ratios
(n-num (numerator number))
(n-den (denominator number))
(d-num (numerator divisor))
(d-den (denominator divisor))
(let ((q-num (* n-num d-den))
(q-den (* n-den d-num)))
(multiple-value-bind (q rem) (truncate q-num q-den)
(values q (/ rem (* n-den d-den)))))))
((bignum fixnum)
(bignum-truncate-single-digit number divisor))
((bignum bignum)
@ -534,43 +495,32 @@
(if (single-digit-bignum-p divisor)
(bignum-truncate-single-digit (make-small-bignum number) divisor)
(bignum-truncate (make-small-bignum number) divisor))))
((ratio integer)
(let ((numerator (numerator number)))
(if (eql divisor 1)
(let ((q (truncate numerator (denominator number))))
(if (< numerator 0)
(let ((q-1 (- q 1)))
(values q-1 (- number q-1)))
(values q (- number q))))
(let ((q (truncate numerator (* (denominator number) divisor))))
(if (if (minusp divisor)
(> numerator 0)
(< numerator 0))
(let ((q-1 (- q 1)))
(values q-1 (- number (* q-1 divisor))))
(values q (- number (* q divisor))))))))
((ratio (foreach (eql 1) fixnum bignum))
(dispatch-ratio (number numerator denominator)
(let* ((q-num numerator)
(q-den (* denominator divisor))
(rem-den denominator))
(multiple-value-bind (q rem) (truncate q-num q-den)
(if (if (minusp divisor)
(> rem 0)
(< rem 0))
(values (1- q) (%make-ratio (+ rem q-den) rem-den))
(values q (%make-ratio rem rem-den)))))))
((ratio ratio)
(let ((numerator (numerator number))
(denominator (denominator number)))
;; (let* ((q (truncate (numerator number) (* (denominator number) divisor)))
;; (mult (* q divisor)))
;; (if (if (minusp divisor)
;; (< mult number)
;; (> mult number))
;; (let ((q-1 (- q 1)))
;; (values q-1 (- number (* q-1 divisor))))
;; (values q (- number (* q divisor)))))
(with-ratio-ops
(multiple-value-bind (div-num div-den) (integer*ratio denominator divisor)
(let ((q (truncate-ratio numerator div-num div-den)))
(multiple-value-bind (mult-num mult-den) (integer*ratio q divisor)
(if (if (minusp divisor)
(compare-two-ratios < mult-num mult-den numerator denominator)
(compare-two-ratios > mult-num mult-den numerator denominator))
(let ((q-1 (- q 1)))
(multiple-value-bind (mult-num mult-den) (integer*ratio q-1 divisor)
(values q-1 (ratio-ratio numerator denominator mult-num mult-den))))
(values q (ratio-ratio numerator denominator mult-num mult-den)))))))))
(dispatch-two-ratios
(n-num (numerator number))
(n-den (denominator number))
(d-num (numerator divisor))
(d-den (denominator divisor))
(let ((q-num (* n-num d-den))
(q-den (* n-den d-num))
(rem-den (* n-den d-den)))
(multiple-value-bind (q rem) (truncate q-num q-den)
(if (if (minusp d-num)
(> rem 0)
(< rem 0))
(values (1- q) (/ (+ rem q-den) rem-den))
(values q (/ rem rem-den)))))))
((bignum fixnum)
(fixup (bignum-truncate-single-digit number divisor)))
((bignum bignum)
@ -631,44 +581,32 @@
(if (single-digit-bignum-p divisor)
(bignum-truncate-single-digit (make-small-bignum number) divisor)
(bignum-truncate (make-small-bignum number) divisor))))
((ratio integer)
(let ((numerator (numerator number)))
(if (eql divisor 1)
(let ((q (truncate numerator (denominator number))))
(if (> numerator 0)
(let ((q+1 (+ q 1)))
(values q+1 (- number q+1)))
(values q (- number q))))
(let ((q (truncate numerator (* (denominator number) divisor))))
(if (if (minusp divisor)
(< numerator 0)
(> numerator 0))
(let ((q+1 (+ q 1)))
(values q+1 (- number (* q+1 divisor))))
(values q (- number (* q divisor))))))))
((ratio (foreach (eql 1) fixnum bignum))
(dispatch-ratio (number numerator denominator)
(let* ((q-num numerator)
(q-den (* denominator divisor))
(rem-den denominator))
(multiple-value-bind (q rem) (truncate q-num q-den)
(if (if (minusp divisor)
(< rem 0)
(> rem 0))
(values (1+ q) (%make-ratio (- rem q-den) rem-den))
(values q (%make-ratio rem rem-den)))))))
((ratio ratio)
;; (let* ((q (truncate (numerator number)
;; (* (denominator number) divisor)))
;; (mult (* q divisor)))
;; (if (if (minusp divisor)
;; (> mult number)
;; (< mult number))
;; (let ((q+1 (+ q 1)))
;; (values q+1 (- number (* q+1 divisor))))
;; (values q (- number (* q divisor)))))
(let ((numerator (numerator number))
(denominator (denominator number)))
(with-ratio-ops
(multiple-value-bind (div-num div-den) (integer*ratio denominator divisor)
(let ((q (truncate-ratio numerator div-num div-den)))
(multiple-value-bind (mult-num mult-den) (integer*ratio q divisor)
(if (if (minusp divisor)
(compare-two-ratios > mult-num mult-den numerator denominator)
(compare-two-ratios < mult-num mult-den numerator denominator))
(let ((q+1 (+ q 1)))
(multiple-value-bind (mult-num mult-den) (integer*ratio q+1 divisor)
(values q+1 (ratio-ratio numerator denominator mult-num mult-den))))
(values q (ratio-ratio numerator denominator mult-num mult-den)))))))))
(dispatch-two-ratios
(n-num (numerator number))
(n-den (denominator number))
(d-num (numerator divisor))
(d-den (denominator divisor))
(let ((q-num (* n-num d-den))
(q-den (* n-den d-num))
(rem-den (* n-den d-den)))
(multiple-value-bind (q rem) (truncate q-num q-den)
(if (if (minusp d-num)
(< rem 0)
(> rem 0))
(values (1+ q) (/ (- rem q-den) rem-den))
(values q (/ rem rem-den)))))))
((bignum fixnum)
(fixup (bignum-truncate-single-digit number divisor)))
((bignum bignum)
@ -1202,14 +1140,6 @@ the first."
(,op x (rational y))))))
) ; EVAL-WHEN
(defmacro dispatch-ratio ((ratio numerator denominator) &body body)
`(let ((,numerator (numerator ,ratio))
(,denominator (denominator ,ratio)))
(if (and (fixnump ,numerator)
(fixnump ,denominator))
(progn ,@body)
(progn ,@body))))
(macrolet ((def-two-arg-</> (name op ratio-arg1 ratio-arg2 &rest cases)
`(defun ,name (x y)
(declare (explicit-check))