math: Add log-expt functions

code/benchmark/integer-float.lisp

@@ -16,13 +16,14 @@
 
 (defun random-float (type)
  (list (random (ash 1 (quaviver:significand-size type)))
- (quaviver/math:floor-log10-expt2 (+ (random (- (quaviver:max-exponent type)
- (quaviver:min-exponent type)
- (quaviver:significand-size type)
- 4))
- (quaviver:min-exponent type)
- (quaviver:significand-size type)
- 2))
+ (quaviver/math:floor-log-expt 10 2
+ (+ (random (- (quaviver:max-exponent type)
+ (quaviver:min-exponent type)
+ (quaviver:significand-size type)
+ 4))
+ (quaviver:min-exponent type)
+ (quaviver:significand-size type)
+ 2))
  (if (zerop (random 2)) 1 -1)))
 
 (defun integer-float (&key (base 10)

code/liebler/implementation.lisp

@@ -11,7 +11,7 @@
  (zerop ,significand))
  (quaviver:integer-float ,client ',result-type 2
  ,significand ,exponent ,sign)
- (let ((k (quaviver/math:floor-log2-expt10 ,exponent))
+ (let ((k (quaviver/math:floor-log-expt 2 10 ,exponent))
  (shift (- ,arithmetic-size (integer-length ,significand))))
  (setf ,significand (,round-to-odd (,expt10 (- ,exponent))
  (ash ,significand shift))

code/math/implementation.lisp

@@ -6,20 +6,25 @@
  (ftype (function ((unsigned-byte 128) (unsigned-byte 64))
  (unsigned-byte 64))
  round-to-odd/64)
+ (ftype (function ((unsigned-byte 256) (unsigned-byte 128))
+ (unsigned-byte 128))
+ round-to-odd/128)
  (ftype (function (fixnum) (unsigned-byte 64))
  expt10/32)
  (ftype (function (fixnum) (unsigned-byte 128))
  expt10/64)
- (ftype (function (fixnum) fixnum)
- floor-log2-expt10)
- (ftype (function (fixnum &optional boolean) fixnum)
- floor-log10-expt2)
+ (ftype (function (fixnum) (unsigned-byte 256))
+ expt10/128)
+ (ftype (function (fixnum fixnum fixnum &optional boolean) fixnum)
+ floor-log-expt ceiling-log-expt)
  (inline round-to-odd/32
  round-to-odd/64
+ round-to-odd/128
  expt10/32
  expt10/64
- floor-log2-expt10
- floor-log10-expt2))
+ expt10/128
+ floor-log-expt
+ ceiling-log-expt))
 
 (defmacro %round-to-odd-1 (g cp size)
  `(let ((p (* ,g ,cp)))
@@ -47,7 +52,7 @@
 (defconstant +expt10/max-exponent/32+ 53)
 
 (defvar *expt10/values/32*
- (compute-expt10 +expt10/min-exponent/32+ +expt10/max-exponent/32+ 64))
+ (compute-expt +expt10/min-exponent/32+ +expt10/max-exponent/32+ 64))
 
 (defun expt10/32 (power)
  (svref *expt10/values/32*
@@ -58,7 +63,7 @@
 (defconstant +expt10/max-exponent/64+ 342)
 
 (defvar *expt10/values/64*
- (compute-expt10 +expt10/min-exponent/64+ +expt10/max-exponent/64+ 128))
+ (compute-expt +expt10/min-exponent/64+ +expt10/max-exponent/64+ 128))
 
 (defun expt10/64 (power)
  (svref *expt10/values/64*
@@ -73,13 +78,80 @@
 (defun expt10/128 (power)
  (svref (or *expt10/values/128*
  (setf *expt10/values/128*
- (compute-expt10 +expt10/min-exponent/128+ +expt10/max-exponent/128+ 256)))
+ (compute-expt +expt10/min-exponent/128+ +expt10/max-exponent/128+ 256)))
  (- (- +expt10/min-exponent/128+) power)))
 
-(defun floor-log2-expt10 (e)
- (ash (* e 1741647) -19))
-
-(defun floor-log10-expt2 (e &optional three-quarters-p)
- (ash (- (* e 1262611)
- (if three-quarters-p 524031 0))
- -22))
+(defconstant +min-base+ 2)
+
+(defconstant +max-base+ 36)
+
+(defconstant +log-expt-shift+ 22)
+
+(defvar *log-expt*
+ (compute-log-expt +min-base+ +max-base+ +log-expt-shift+))
+
+(defvar *log-3/4*
+ (compute-log-3/4 +min-base+ +max-base+ +log-expt-shift+))
+
+(defun floor-log-expt (log-base expt-base exp &optional three-quarters-p)
+ (declare (optimize speed))
+ (ash (+ (* exp (aref *log-expt*
+ (- log-base +min-base+)
+ (- expt-base +min-base+)))
+ (if three-quarters-p
+ (svref *log-3/4* (- log-base +min-base+))
+ 0))
+ (- +log-expt-shift+)))
+
+(define-compiler-macro floor-log-expt
+ (&whole whole log-base expt-base exp &optional three-quarters-p)
+ (if (or (not (constantp log-base))
+ (not (constantp expt-base)))
+ whole
+ (let ((multiplier (aref *log-expt*
+ (- log-base +min-base+)
+ (- expt-base +min-base+)))
+ (offset (- +log-expt-shift+))
+ (shift (- +log-expt-shift+)))
+ (cond ((null three-quarters-p)
+ `(ash (* ,exp ,multiplier) ,shift))
+ ((constantp three-quarters-p)
+ `(ash (+ (* ,exp ,multiplier) ,offset)
+ ,shift))
+ (t
+ `(ash (+ (* ,exp ,multiplier)
+ (if ,three-quarters-p
+ ,offset
+ 0))
+ ,shift))))))
+
+(defun ceiling-log-expt (log-base expt-base exp &optional three-quarters-p)
+ (values (ceiling (+ (* exp (aref *log-expt*
+ (- log-base +min-base+)
+ (- expt-base +min-base+)))
+ (if three-quarters-p
+ (svref *log-3/4* (- log-base +min-base+))
+ 0))
+ (ash 1 +log-expt-shift+))))
+
+(define-compiler-macro ceiling-log-expt
+ (&whole whole log-base expt-base exp &optional three-quarters-p)
+ (if (or (not (constantp log-base))
+ (not (constantp expt-base)))
+ whole
+ (let ((multiplier (aref *log-expt*
+ (- log-base +min-base+)
+ (- expt-base +min-base+)))
+ (offset (- +log-expt-shift+))
+ (divisor (ash 1 +log-expt-shift+)))
+ (cond ((null three-quarters-p)
+ `(values (ceiling (* ,exp ,multiplier) ,divisor)))
+ ((constantp three-quarters-p)
+ `(values (ceiling (+ (* ,exp ,multiplier) ,offset)
+ ,divisor)))
+ (t
+ `(values (ceiling (+ (* ,exp ,multiplier)
+ (if ,three-quarters-p
+ ,offset
+ 0))
+ ,divisor)))))))

code/math/packages.lisp

@@ -6,5 +6,6 @@
  #:round-to-odd/32
  #:round-to-odd/64
  #:round-to-odd/128
+ #:floor-log-expt
  #:floor-log2-expt10
  #:floor-log10-expt2))

code/math/utility.lisp

@@ -1,6 +1,6 @@
 (in-package #:quaviver/math)
 
-(defun compute-expt10 (k-min k-max width &optional (base 10))
+(defun compute-expt (k-min k-max width &optional (base 10))
  (make-array (- k-max k-min -1)
  :initial-contents
  (loop for k from k-min upto k-max
@@ -14,3 +14,19 @@
  (ceiling (/ l (ash 1 (abs shift))))
  (ceiling (* (ash 1 shift)
  l)))))))))
+
+(defun compute-log-expt (min-base max-base shift)
+ (make-array (list (- max-base min-base -1)
+ (- max-base min-base -1))
+ :initial-contents
+ (loop for log-base from min-base upto max-base
+ collect (loop for expt-base from min-base upto max-base
+ collect (floor (* (log (coerce expt-base 'double-float)
+ log-base)
+ (ash 1 shift)))))))
+
+(defun compute-log-3/4 (min-base max-base shift)
+ (make-array (- max-base min-base -1)
+ :initial-contents
+ (loop for log-base from min-base upto max-base
+ collect (- (floor (* (log (coerce 4/3 'double-float) log-base) (ash 1 shift)))))))

code/schubfach/implementation.lisp

@@ -19,8 +19,8 @@
  (values significand exponent sign)
  (let* ((lower-boundary-is-closer (= significand ,(ash 1 (1- significand-size))))
  (is-even (evenp significand))
- (k (quaviver/math:floor-log10-expt2 exponent lower-boundary-is-closer))
- (h (+ exponent 1 (quaviver/math:floor-log2-expt10 (- k))))
+ (k (quaviver/math:floor-log-expt 10 2 exponent lower-boundary-is-closer))
+ (h (+ exponent 1 (quaviver/math:floor-log-expt 2 10 (- k))))
  (expt10 (,expt10 k)))
  (declare (type (unsigned-byte ,(ash arithmetic-size 1))
  expt10)