Skip to content

Commit

Permalink
math, dragonbox: introduce smallnum optimizations
Browse files Browse the repository at this point in the history 
Mostly ported from the reference implementation of Dragonbox.

They are activated by default only on SBCL, as on other implementations
they tend to be slower than bignum operations.

* code/math/implementation.lisp (+/64-64, +/128-64, */32-64/hi64)
(*/32-64/lo64, */64-64, */64-64/hi64, */64-64/lo64, */64-128/hi128)
(*/64-128/lo128): New functions, mostly ported from reference Dragonbox.
(hi/hi64/128): New smallnum optimization and compiler macro.
(round-to-odd/32-64, round-to-odd/64-128): New smallnum optimizations.
(floor-multiply/32-64q64, floor-multiply/64-128q128)
(floor-multiply/evenp/32-64q64, floor-multiply/evenp/64-128q128): New
smallnum optimizations from reference Dragonbox.
(*expt10/values/64*): Add smallnum version.
* code/math/packages.lisp: Push QUAVIVER/MATH/SMALLNUM feature on SBCL.
* code/dragonbox/implementation.lisp (floor-by-expt10): New smallnum
optimization from reference Dragonbox.
  • Loading branch information
@paulapatience
paulapatience committed Jul 4, 2024
1 parent 3a38565 commit 3b06535
Show file tree
Hide file tree
Showing 3 changed files with 191 additions and 2 deletions.
6 changes: 6 additions & 0 deletions code/dragonbox/implementation.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,10 +114,16 @@
(cond ((and (= size 32) (= power 1) (<= max-number 1073741828))
`(ldb (byte 32 32) (* ,number 429496730)))
((and (= size 64) (= power 1) (<= max-number 4611686018427387908))
#+quaviver/math/smallnum
`(quaviver/math::*/64-64/hi64 ,number 1844674407370955162)
#-quaviver/math/smallnum
`(ldb (byte 64 64) (* ,number 1844674407370955162)))
((and (= size 32) (= power 2))
`(ldb (byte 27 37) (* ,number 1374389535)))
((and (= size 64) (= power 3) (<= max-number 15534100272597517998))
#+quaviver/math/smallnum
`(ash (quaviver/math::*/64-64/hi64 ,number 4722366482869645214) -8)
#-quaviver/math/smallnum
`(ldb (byte 56 72) (* ,number 4722366482869645214)))
(t
`(floor ,number ,(expt 10 power))))))
Expand Down
184 changes: 182 additions & 2 deletions code/math/implementation.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,128 @@
#-quaviver/math/smallnum
`(unsigned-byte ,(* arithmetic-size count)))

;;; Low-level, smallnum-optimized arithmetic operations
;;;
;;; These operations represent bignums as individual ub64 values.
;;; They are not exported from quaviver/math.
;;;
;;; Some of these operations are ported from Dragonbox.

(declaim (ftype (function ((unsigned-byte 64) (unsigned-byte 64))
(values (unsigned-byte 64) &optional))
+/64-64)
(ftype (function ((unsigned-byte 64) (unsigned-byte 64) (unsigned-byte 64))
(values (unsigned-byte 64) (unsigned-byte 64) &optional))
+/128-64)
(ftype (function ((unsigned-byte 32) (unsigned-byte 64))
(values (unsigned-byte 64) &optional))
*/32-64/hi64
*/32-64/lo64)
(ftype (function ((unsigned-byte 64) (unsigned-byte 64))
(values (unsigned-byte 64) (unsigned-byte 64) &optional))
*/64-64)
(ftype (function ((unsigned-byte 64) (unsigned-byte 64))
(values (unsigned-byte 64) &optional))
*/64-64/hi64
*/64-64/lo64)
(ftype (function ((unsigned-byte 64) (unsigned-byte 64) (unsigned-byte 64))
(values (unsigned-byte 64) (unsigned-byte 64) &optional))
*/64-128/hi128
*/64-128/lo128)
(inline +/64-64
+/128-64
*/32-64/hi64
*/32-64/lo64
*/64-64
*/64-64/hi64
*/64-64/lo64
*/64-128/hi128
*/64-128/lo128))

(defun +/64-64 (x y)
(ldb (byte 64 0) (+ x y)))

;;; Based on https://github.com/jk-jeon/dragonbox/blob/04bc662afe22576fd0aa740c75dca63609297f19/include/dragonbox/dragonbox.h#L646-L650
(defun +/128-64 (xh xl y)
(let ((rl (+/64-64 xl y)))
(values (+/64-64 xh (if (< rl xl) 1 0))
rl)))

;;; Based on https://github.com/jk-jeon/dragonbox/blob/04bc662afe22576fd0aa740c75dca63609297f19/include/dragonbox/dragonbox.h#L835-L848
(defun */32-64/hi64 (x y)
(+/64-64 (* x (ldb (byte 32 32) y))
(ash (* x (ldb (byte 32 0) y))
-32)))

(defun */32-64/lo64 (x y)
;; SBCL elides bignums.
(ldb (byte 64 0) (* x y)))

;;; Based on https://github.com/jk-jeon/dragonbox/blob/04bc662afe22576fd0aa740c75dca63609297f19/include/dragonbox/dragonbox.h#L737-L753
(defun */64-64 (x y)
(let* ((a (ldb (byte 32 32) x))
(b (ldb (byte 32 0) x))
(c (ldb (byte 32 32) y))
(d (ldb (byte 32 0) y))
(ac (* a c))
(bc (* b c))
(ad (* a d))
(bd (* b d))
(u (+ (ldb (byte 32 32) bd)
(ldb (byte 32 0) ad)
(ldb (byte 32 0) bc))))
(declare ((unsigned-byte 32) a b c d)
((unsigned-byte 64) ac bc ad bd u))
;; (+ (ash ac 64) (ash ad 32) (ash bc 32) bd)
(values (ldb (byte 64 0)
(+ ac
(ldb (byte 32 32) u)
(ldb (byte 32 32) ad)
(ldb (byte 32 32) bc)))
(+/64-64 (ash (ldb (byte 32 0) u)
32)
(ldb (byte 32 0) bd)))))

;;; Based on https://github.com/jk-jeon/dragonbox/blob/04bc662afe22576fd0aa740c75dca63609297f19/include/dragonbox/dragonbox.h#L783-L798
(defun */64-64/hi64 (x y)
(let* ((a (ldb (byte 32 32) x))
(b (ldb (byte 32 0) x))
(c (ldb (byte 32 32) y))
(d (ldb (byte 32 0) y))
(ac (* a c))
(bc (* b c))
(ad (* a d))
(bd (* b d))
(u (+ (ldb (byte 32 32) bd)
(ldb (byte 32 0) ad)
(ldb (byte 32 0) bc))))
(declare ((unsigned-byte 32) a b c d)
((unsigned-byte 64) ac bc ad bd u))
;; (+ (ash ac 64) (ash ad 32) (ash bc 32) bd)
(ldb (byte 64 0)
(+ ac
(ldb (byte 32 32) u)
(ldb (byte 32 32) ad)
(ldb (byte 32 32) bc)))))

(defun */64-64/lo64 (x y)
;; SBCL elides bignums.
(ldb (byte 64 0) (* x y)))

;;; Based on https://github.com/jk-jeon/dragonbox/blob/04bc662afe22576fd0aa740c75dca63609297f19/include/dragonbox/dragonbox.h#L826-L831
(defun */64-128/hi128 (x yh yl)
(multiple-value-bind (rh rl) (*/64-64 x yh)
(+/128-64 rh rl (*/64-64/hi64 x yl))))

;;; Based on https://github.com/jk-jeon/dragonbox/blob/04bc662afe22576fd0aa740c75dca63609297f19/include/dragonbox/dragonbox.h#L852-L857
(defun */64-128/lo128 (x yh yl)
(multiple-value-bind (rh rl) (*/64-64 x yl)
(values (+/64-64 (*/64-64/lo64 x yh)
rh)
rl)))

;;; High-level arithmetic operations

(declaim (ftype (function ((arithmetic-word 64) (integer 0 64))
(values (arithmetic-word 64) &optional))
hi/64)
Expand Down Expand Up @@ -89,8 +211,20 @@
(otherwise whole)))

(defun hi/hi64/128 (x count)
#+quaviver/math/smallnum
(hi/64 (aref x 0) count)
#-quaviver/math/smallnum
(ash x (- count 128)))

(define-compiler-macro hi/hi64/128 (&whole whole x count)
#-quaviver/math/smallnum
(declare (ignore x))
(case count
(0 0)
#+quaviver/math/smallnum
(64 `(aref (the (arithmetic-word 64 2) ,x) 0))
(otherwise whole)))

(defmacro %round-to-odd-1 (cp g size)
`(let ((p (* ,cp ,g)))
(logior (ldb (byte ,size ,(ash size 1)) p)
Expand All @@ -103,10 +237,24 @@
(ash p ,(- size)))))

(defun round-to-odd/32-64 (cp g)
#-(or ecl cmucl) (%round-to-odd-1 cp g 32)
#+(or ecl cmucl) (%round-to-odd-2 cp g 32))
#+quaviver/math/smallnum
(let ((p (*/32-64/hi64 cp g)))
(if (ldb-test (byte 31 1) p)
(logior (ash p -32) 1)
(ash p -32)))
#+(and (not quaviver/math/smallnum) (not (or ecl cmucl)))
(%round-to-odd-1 cp g 32)
#+(and (not quaviver/math/smallnum) (or ecl cmucl))
(%round-to-odd-2 cp g 32))

(defun round-to-odd/64-128 (cp g)
#+quaviver/math/smallnum
(multiple-value-bind (ph pl)
(*/64-128/hi128 cp (aref g 0) (aref g 1))
(if (ldb-test (byte 63 1) pl)
(logior ph 1)
ph))
#-quaviver/math/smallnum
(%round-to-odd-2 cp g 64))

(defun round-to-odd/128-256 (cp g)
Expand All @@ -132,9 +280,19 @@
(zerop (ldb (byte ,size ,size) r))))) ; integer-p

(defun floor-multiply/32-64q64 (x y &optional (pre-shift 0))
#+quaviver/math/smallnum
(let ((r (*/32-64/hi64 (ash x pre-shift) y)))
(values (ldb (byte 32 32) r) ; integer part
(not (ldb-test (byte 32 0) r)))) ; integer-p
#-quaviver/math/smallnum
(%floor-multiply/n-2nq2n x y pre-shift 32))

(defun floor-multiply/64-128q128 (x y &optional (pre-shift 0))
#+quaviver/math/smallnum
(multiple-value-bind (rh rl)
(*/64-128/hi128 (ash x pre-shift) (aref y 0) (aref y 1))
(values rh (zerop rl))) ; integer part, integer-p
#-quaviver/math/smallnum
(%floor-multiply/n-2nq2n x y pre-shift 64))

;;; The FLOOR-MULTIPLY/EVENP operations compute the same thing as
Expand All @@ -152,9 +310,21 @@
(zerop (ldb (byte ,size (- ,size ,pre-shift)) r))))) ; integer-p

(defun floor-multiply/evenp/32-64q64 (x y &optional (pre-shift 0))
#+quaviver/math/smallnum
(let ((r (*/32-64/lo64 x y)))
(values (logbitp (- 64 pre-shift) r) ; even-p
(not (ldb-test (byte 32 (- 32 pre-shift)) r)))) ; integer-p
#-quaviver/math/smallnum
(%floor-multiply/evenp/n-2nq2n x y pre-shift 32))

(defun floor-multiply/evenp/64-128q128 (x y &optional (pre-shift 0))
#+quaviver/math/smallnum
(multiple-value-bind (rh rl)
(*/64-128/lo128 x (aref y 0) (aref y 1))
(values (logbitp (- 64 pre-shift) rh) ; even-p
(and (zerop (ldb (byte (- 64 pre-shift) pre-shift) rh)) ; integer-p
(zerop (ldb (byte 64 (- 64 pre-shift)) rl)))))
#-quaviver/math/smallnum
(%floor-multiply/evenp/n-2nq2n x y pre-shift 64))

(defconstant +expt10/min-exponent/32+ -53)
Expand All @@ -173,6 +343,16 @@
(defconstant +expt10/max-exponent/64+ 342)

(defvar *expt10/values/64*
#+quaviver/math/smallnum
(let ((table (compute-expt +expt10/min-exponent/64+ +expt10/max-exponent/64+ 128)))
(make-array (length table)
:initial-contents
(loop for x across table
collect (make-array 2 :element-type '(unsigned-byte 64)
:initial-contents
(list (ldb (byte 64 64) x)
(ldb (byte 64 0) x))))))
#-quaviver/math/smallnum
(compute-expt +expt10/min-exponent/64+ +expt10/max-exponent/64+ 128))

(defun expt10/64 (power)
Expand Down
3 changes: 3 additions & 0 deletions code/math/packages.lisp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,3 +18,6 @@
#:floor-multiply/evenp/64-128q128
#:floor-log-expt
#:ceiling-log-expt))

#+sbcl
(pushnew :quaviver/math/smallnum *features*)

0 comments on commit 3b06535

Please sign in to comment.