Newer
Older
;;;; Copyright (c) 2007, 2008, 2011 Raymond Toy
;;;;
;;;; Permission is hereby granted, free of charge, to any person
;;;; obtaining a copy of this software and associated documentation
;;;; files (the "Software"), to deal in the Software without
;;;; restriction, including without limitation the rights to use,
;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
;;;; copies of the Software, and to permit persons to whom the
;;;; Software is furnished to do so, subject to the following
;;;; conditions:
;;;;
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;;
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
(defmethod make-qd ((x cl:rational))
;; We should do something better than this.
(make-instance 'qd-real :value (rational-to-qd x)))
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
;; Determine which of x and y has the higher precision and return the
;; value of the higher precision number. If both x and y are
;; rationals, just return 1f0, for a single-float value.
(defun float-contagion-2 (x y)
(etypecase x
(cl:rational
(etypecase y
(cl:rational
1f0)
(cl:float
y)
(qd-real
y)))
(single-float
(etypecase y
((or cl:rational single-float)
x)
((or double-float qd-real)
y)))
(double-float
(etypecase y
((or cl:rational single-float double-float)
x)
(qd-real
y)))
(qd-real
x)))
;; Return a floating point (or complex) type of the highest precision
;; among all of the given arguments.
(defun float-contagion (&rest args)
;; It would be easy if we could just add the args together and let
;; normal contagion do the work, but we could easily introduce
;; overflows or other errors that way. So look at each argument and
;; determine the precision and choose the highest precision.
(etypecase (reduce #'float-contagion-2 (mapcar #'realpart (if (cdr args)
args
(list (car args) 0))))
(single-float 'single-float)
(double-float 'double-float)
(qd-real 'qd-real)))
(defun apply-contagion (number precision)
(etypecase number
((or cl:real qd-real)
(coerce number precision))
((or cl:complex qd-complex)
(complex (coerce (realpart number) precision)
(coerce (imagpart number) precision)))))
;; WITH-FLOATING-POINT-CONTAGION - macro
;;
;; Determines the highest precision of the variables in VARLIST and
;; converts each of the values to that precision.
(defmacro with-floating-point-contagion (varlist &body body)
(let ((precision (gensym "PRECISION-")))
`(let ((,precision (float-contagion ,@varlist)))
(let (,@(mapcar #'(lambda (v)
`(,v (apply-contagion ,v ,precision)))
varlist))
,@body))))
(defmethod add1 ((a number))
(cl::1+ a))
(defmethod add1 ((a qd-real))
(make-instance 'qd-real :value (add-qd-d (qd-value a) 1d0)))
(defmethod sub1 ((a number))
(cl::1- a))
(defmethod sub1 ((a qd-real))
(make-instance 'qd-real :value (sub-qd-d (qd-value a) 1d0)))
(declaim (inline 1+ 1-))
(defun 1+ (x)
(add1 x))
(defun 1- (x)
(sub1 x))
(defmethod two-arg-+ ((a qd-real) (b qd-real))
(make-instance 'qd-real :value (add-qd (qd-value a) (qd-value b))))
(defmethod two-arg-+ ((a qd-real) (b cl:float))
(make-instance 'qd-real :value (add-qd-d (qd-value a) (cl:float b 1d0))))
#+cmu
(defmethod two-arg-+ ((a qd-real) (b ext:double-double-float))
(make-instance 'qd-real :value (add-qd-dd (qd-value a) b)))
(defmethod two-arg-+ ((a real) (b qd-real))
(defmethod two-arg-+ ((a number) (b number))
(cl:+ a b))
(defun + (&rest args)
(if (null args)
0
(do ((args (cdr args) (cdr args))
(res (car args)
(two-arg-+ res (car args))))
((null args) res))))
(defmethod two-arg-- ((a qd-real) (b qd-real))
(make-instance 'qd-real :value (sub-qd (qd-value a) (qd-value b))))
(defmethod two-arg-- ((a qd-real) (b cl:float))
(make-instance 'qd-real :value (sub-qd-d (qd-value a) (cl:float b 1d0))))
#+cmu
(defmethod two-arg-- ((a qd-real) (b ext:double-double-float))
(make-instance 'qd-real :value (sub-qd-dd (qd-value a) b)))
(defmethod two-arg-- ((a cl:float) (b qd-real))
(make-instance 'qd-real :value (sub-d-qd (cl:float a 1d0) (qd-value b))))
(defmethod two-arg-- ((a number) (b number))
(cl:- a b))
(defmethod unary-minus ((a number))
(cl:- a))
(defmethod unary-minus ((a qd-real))
(make-instance 'qd-real :value (neg-qd (qd-value a))))
(defun - (number &rest more-numbers)
(if more-numbers
(do ((nlist more-numbers (cdr nlist))
(result number))
((atom nlist) result)
(declare (list nlist))
(setq result (two-arg-- result (car nlist))))
(unary-minus number)))
(defmethod two-arg-* ((a qd-real) (b qd-real))
(make-instance 'qd-real :value (mul-qd (qd-value a) (qd-value b))))
(defmethod two-arg-* ((a qd-real) (b cl:float))
(make-instance 'qd-real :value (mul-qd-d (qd-value a) (cl:float b 1d0))))
#+cmu
(defmethod two-arg-* ((a qd-real) (b ext:double-double-float))
;; We'd normally want to use mul-qd-dd, but mul-qd-dd is broken.
(make-instance 'qd-real :value (mul-qd (qd-value a)
(make-qd-dd b 0w0))))
(defmethod two-arg-* ((a real) (b qd-real))
(defmethod two-arg-* ((a number) (b number))
(cl:* a b))
(defun * (&rest args)
(if (null args)
1
(do ((args (cdr args) (cdr args))
(res (car args)
(two-arg-* res (car args))))
((null args) res))))
(defmethod two-arg-/ ((a qd-real) (b qd-real))
(make-instance 'qd-real :value (div-qd (qd-value a) (qd-value b))))
(defmethod two-arg-/ ((a qd-real) (b cl:float))
(make-instance 'qd-real :value (div-qd-d (qd-value a) (cl:float b 1d0))))
#+cmu
(defmethod two-arg-/ ((a qd-real) (b ext:double-double-float))
(make-instance 'qd-real :value (div-qd-dd (qd-value a)
(defmethod two-arg-/ ((a cl:float) (b qd-real))
(qd-value b))))
#+cmu
(defmethod two-arg-/ ((a ext:double-double-float) (b qd-real))
(make-instance 'qd-real :value (div-qd (make-qd-dd a 0w0)
(qd-value b))))
(defmethod two-arg-/ ((a number) (b number))
(cl:/ a b))
(defmethod unary-divide ((a number))
(cl:/ a))
(defmethod unary-divide ((a qd-real))
(make-instance 'qd-real :value (div-qd +qd-one+ (qd-value a))))
(defun / (number &rest more-numbers)
(if more-numbers
(do ((nlist more-numbers (cdr nlist))
(result number))
((atom nlist) result)
(declare (list nlist))
(setq result (two-arg-/ result (car nlist))))
(unary-divide number)))
(macrolet ((frob (name &optional (type 'real))
(let ((method-name (intern (concatenate 'string
(string '#:q)
(symbol-name name))))
(qd-name (intern (concatenate 'string
(symbol-name name)
(string '#:-qd)))))
`(progn
(defmethod ,method-name ((x ,type))
(,cl-name x))
(defmethod ,method-name ((x qd-real))
(,qd-name (qd-value x)))
(declaim (inline ,name))
(defun ,name (x)
(,method-name x))))))
(defun bignum-to-qd (bignum)
(make-instance 'qd-real
(defun floatp (x)
(typep x '(or short-float single-float double-float long-float qd-real)))
(defmethod qfloat ((x real) (num-type cl:float))
(cl:float x num-type))
(defmethod qfloat ((x cl:float) (num-type qd-real))
(make-instance 'qd-real :value (make-qd-d (cl:float x 1d0))))
(defmethod qfloat ((x integer) (num-type qd-real))
(cond ((typep x 'fixnum)
(make-instance 'qd-real :value (make-qd-d (cl:float x 1d0))))
(t
;; A bignum
(bignum-to-qd x))))
(defmethod qfloat ((x ratio) (num-type qd-real))
;; This probably has some issues with roundoff
(two-arg-/ (qfloat (numerator x) num-type)
(qfloat (denominator x) num-type)))
#+cmu
(defmethod qfloat ((x ext:double-double-float) (num-type qd-real))
(make-instance 'qd-real :value (make-qd-dd x 0w0)))
(multiple-value-bind (q0 q1 q2 q3)
(cl:float (cl:+ q0 q1 q2 q3) num-type)))
#+cmu
(defmethod qfloat ((x qd-real) (num-type ext:double-double-float))
(multiple-value-bind (q0 q1 q2 q3)
(cl:+ (cl:float q0 1w0)
(cl:float q1 1w0)
(cl:float q2 1w0)
(cl:float q3 1w0))))
(defmethod qfloat ((x qd-real) (num-type qd-real))
x)
(declaim (inline float))
(defun float (x &optional (other nil otherp))
(if otherp
(qfloat x other)
(if (floatp x)
x
(qfloat x 1.0))))
(defmethod qrealpart ((x number))
(cl:realpart x))
(defmethod qrealpart ((x qd-real))
x)
(defmethod qrealpart ((x qd-complex))
(make-instance 'qd-real :value (qd-real x)))
(defun realpart (x)
(qrealpart x))
(defmethod qimagpart ((x number))
(cl:imagpart x))
(defmethod qimagpart ((x qd-real))
(make-instance 'qd-real :value (qd-imag x)))
(defmethod qconjugate ((a number))
(cl:conjugate a))
(defmethod qconjugate ((a qd-real))
(make-instance 'qd-real :value (qd-value a)))
(defmethod qconjugate ((a qd-complex))
(make-instance 'qd-complex
:real (qd-real a)
:imag (neg-qd (qd-imag a))))
(defun conjugate (z)
(qconjugate z))
(defmethod qscale-float ((f cl:float) (n integer))
(cl:scale-float f n))
(defmethod qscale-float ((f qd-real) (n integer))
(make-instance 'qd-real :value (scale-float-qd (qd-value f) n)))
(declaim (inline scale-float))
(defun scale-float (f n)
(qscale-float f n))
(macrolet
((frob (op)
(let ((method (intern (concatenate 'string
(string '#:two-arg-)
(symbol-name op))))
(cl-fun (find-symbol (symbol-name op) :cl))
(qd-fun (intern (concatenate 'string (string '#:qd-) (symbol-name op))
`(progn
(defmethod ,method ((a real) (b real))
(,cl-fun a b))
(defmethod ,method ((a qd-real) (b real))
(,qd-fun (qd-value a) (make-qd-d (cl:float b 1d0))))
(defmethod ,method ((a real) (b qd-real))
;; This is not really right if A is a rational. We're
;; supposed to compare them as rationals.
(,qd-fun (make-qd-d (cl:float a 1d0)) (qd-value b)))
(defmethod ,method ((a qd-real) (b qd-real))
(,qd-fun (qd-value a) (qd-value b)))
(defun ,op (number &rest more-numbers)
"Returns T if its arguments are in strictly increasing order, NIL otherwise."
(declare (optimize (safety 2))
(dynamic-extent more-numbers))
(do* ((n number (car nlist))
(nlist more-numbers (cdr nlist)))
((atom nlist) t)
(declare (list nlist))
(if (not (,method n (car nlist))) (return nil))))))))
(frob <)
(frob >)
(frob <=)
(frob >=))
;; Handle the special functions for a real argument. Complex args are
;; handled elsewhere.
(macrolet
((frob (name)
(let ((method-name
(intern (concatenate 'string (string '#:q)
(symbol-name name))))
(cl-name (intern (symbol-name name) :cl))
(qd-name (intern (concatenate 'string (symbol-name name)
(string '#:-qd)))))
(defmethod ,name ((x number))
(,cl-name x))
(defmethod ,name ((x qd-real))
(make-instance 'qd-real :value (,qd-name (qd-value x))))))))
(frob abs)
(frob exp)
(frob sin)
(frob cos)
(frob tan)
;;(frob asin)
;;(frob acos)
(frob sinh)
(frob cosh)
(frob tanh)
(frob asinh)
;;(frob acosh)
;;(frob atanh)
)
(defmethod sqrt ((x number))
(defmethod sqrt ((x qd-real))
(if (minusp x)
(make-instance 'qd-complex
:real +qd-zero+
:imag (sqrt-qd (neg-qd (qd-value x))))
(make-instance 'qd-real :value (sqrt-qd (qd-value x)))))
(defun scalb (x n)
"Compute 2^N * X without compute 2^N first (use properties of the
underlying floating-point format"
(declare (type qd-real x))
(scale-float x n))
(multiple-value-bind (rho k)
(octi::hypot-aux-qd (qd-value (realpart z))
(qd-value (imagpart z)))
(values (make-instance 'qd-real :value rho)
k)))
(defmethod qabs ((z qd-complex))
;; sqrt(x^2+y^2)
;; If |x| > |y| then sqrt(x^2+y^2) = |x|*sqrt(1+(y/x)^2)
(multiple-value-bind (abs^2 rho)
(hypot-qd (qd-value (realpart z))
(qd-value (imagpart z)))
(scale-float (make-instance 'qd-real :value (sqrt abs^2))
rho)))
;; sqrt(x^2+y^2)
;; If |x| > |y| then sqrt(x^2+y^2) = |x|*sqrt(1+(y/x)^2)
(make-instance 'qd-real
:value (hypot-qd (qd-value (realpart z))
(qd-value (imagpart z)))))
(defmethod log ((a number) &optional b)
(defmethod log ((a qd-real) &optional b)
(if (minusp (float-sign a))
(make-instance 'qd-complex
:real (log-qd (abs-qd (qd-value a)))
:imag +qd-pi+)
(make-instance 'qd-real :value (log-qd (qd-value a))))))
(make-instance 'qd-real :value (log1p-qd (qd-value a))))
(defmethod atan ((y real) &optional x)
(cond (x
(cond ((typep x 'qd-real)
(make-instance 'qd-real
:value (atan2-qd (qd-value y) (qd-value x))))
(t
(cl:atan y x))))
(t
(cl:atan y))))
(defmethod atan ((y qd-real) &optional x)
(make-instance 'qd-real
:value
(if x
(atan2-qd (qd-value y) (qd-value x))
(atan-qd (qd-value y)))))
(defmethod qexpt ((x number) (y number))
(cl:expt x y))
(defmethod qexpt ((x number) (y qd-real))
(exp (* y (log (apply-contagion x 'qd-real)))))
(defmethod qexpt ((x number) (y qd-complex))
(exp (* y (log (apply-contagion x 'qd-real)))))
(defmethod qexpt ((x qd-real) (y cl:complex))
(exp (* (make-instance 'qd-complex
:real (qd-value (realpart y))
:imag (qd-value (imagpart y)))
(log x))))
(defmethod qexpt ((x qd-real) (y qd-real))
;; x^y = exp(y*log(x))
(exp (* y (log x))))
(defmethod qexpt ((x qd-real) (y integer))
(make-instance 'qd-real
(declaim (inline expt))
(defun expt (x y)
(qexpt x y))
(defmethod two-arg-= ((a number) (b number))
(cl:= a b))
(defmethod two-arg-= ((a qd-real) (b number))
(if (cl:realp b)
(defmethod two-arg-= ((a number) (b qd-real))
(if (cl:realp a)
(defmethod two-arg-= ((a qd-complex) b)
(and (two-arg-= (realpart a) (realpart b))
(two-arg-= (imagpart a) (imagpart b))))
(defmethod two-arg-= (a (b qd-complex))
(and (two-arg-= (realpart a) (realpart b))
(two-arg-= (imagpart a) (imagpart b))))
(defmethod two-arg-= ((a qd-real) (b qd-real))
(defun = (number &rest more-numbers)
"Returns T if all of its arguments are numerically equal, NIL otherwise."
(declare (optimize (safety 2))
(dynamic-extent more-numbers))
(do ((nlist more-numbers (cdr nlist)))
(declare (list nlist))
(if (not (two-arg-= (car nlist) number))
(return nil))))
(defun /= (number &rest more-numbers)
"Returns T if no two of its arguments are numerically equal, NIL otherwise."
(declare (optimize (safety 2))
(dynamic-extent more-numbers))
(do* ((head number (car nlist))
(nlist more-numbers (cdr nlist)))
((atom nlist) t)
(declare (list nlist))
(unless (do* ((nl nlist (cdr nl)))
(declare (list nl))
(if (two-arg-= head (car nl))
(return nil)))
(return nil))))
(defmethod qcomplex ((x cl:real) (y cl:real))
(cl:complex x y))
(defmethod qcomplex ((x cl:real) (y qd-real))
(qcomplex (make-qd x) y))
(defmethod qcomplex ((x qd-real) (y qd-real))
(make-instance 'qd-complex
:real (qd-value x)
:imag (qd-value y)))
(defmethod qcomplex ((x qd-real) (y cl:real))
(make-instance 'qd-complex
:real (qd-value x)
(defun complex (x &optional (y 0))
(qcomplex x y))
(cl:integer-decode-float f))
(defmethod qinteger-decode-float ((f qd-real))
(integer-decode-qd (qd-value f)))
(declaim (inline integer-decode-float))
(defun integer-decode-float (f)
(qinteger-decode-float f))
(cl:decode-float f))
(defmethod qdecode-float ((f qd-real))
(multiple-value-bind (frac exp s)
(decode-float-qd (qd-value f))
(values (make-instance 'qd-real :value frac)
exp
(make-instance 'qd-real :value s))))
(declaim (inline decode-float))
(defun decode-float (f)
(qdecode-float f))
(defmethod qfloor ((x real) &optional y)
(if y
(cl:floor x y)
(cl:floor x)))
(defmethod qfloor ((x qd-real) &optional y)
(if (and y (/= y 1))
(let ((f (qfloor (/ x y))))
(values f
(- x (* f y))))
(multiple-value-bind (int exp sign)
(integer-decode-qd f)
(values (ash (* sign int) exp)
(make-instance 'qd-real
:value (qd-value
(- x (make-instance 'qd-real
:value f)))))))))
(defun floor (x &optional y)
(qfloor x y))
(defmethod qffloor ((x real) &optional y)
(if y
(cl:ffloor x y)
(cl:ffloor x)))
(defmethod qffloor ((x qd-real) &optional y)
(if (and y (/= y 1))
(let ((f (qffloor (/ x y))))
(values f
(- x (* f y))))
(let ((f (make-instance 'qd-real :value (ffloor-qd (qd-value x)))))
(values f
(- x f)))))
(defun ffloor (x &optional y)
(qffloor x y))
(defun ceiling (x &optional y)
(multiple-value-bind (f rem)
(floor x y)
(if (zerop rem)
rem)
(values (+ f 1)
(- rem 1)))))
(defun fceiling (x &optional y)
(multiple-value-bind (f rem)
(ffloor x y)
(if (zerop rem)
rem)
(values (+ f 1)
(- rem 1)))))
(defun truncate (x &optional (y 1))
(ceiling x y)
(floor x y)))
(defun rem (x y)
(nth-value 1 (truncate x y)))
(defun mod (x y)
(nth-value 1 (floor x y)))
(defun ftruncate (x &optional (y 1))
(fceiling x y)
(ffloor x y)))
(defmethod %unary-round ((x real))
(cl::round x))
(defmethod %unary-round ((number qd-real))
(multiple-value-bind (bits exp)
(integer-decode-float number)
(let* ((shifted (ash bits exp))
(rounded (if (and (minusp exp)
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
(not (zerop (logand bits
(ash 1 (- -1 exp))))))
(1+ shifted)
shifted)))
(if (minusp number)
(- rounded)
rounded))))
(defun round (number &optional (divisor 1))
(if (eql divisor 1)
(let ((r (%unary-round number)))
(values r
(- number r)))
(multiple-value-bind (tru rem)
(truncate number divisor)
(if (zerop rem)
(values tru rem)
(let ((thresh (/ (abs divisor) 2)))
(cond ((or (> rem thresh)
(and (= rem thresh) (oddp tru)))
(if (minusp divisor)
(values (- tru 1) (+ rem divisor))
(values (+ tru 1) (- rem divisor))))
((let ((-thresh (- thresh)))
(or (< rem -thresh)
(and (= rem -thresh) (oddp tru))))
(if (minusp divisor)
(values (+ tru 1) (- rem divisor))
(values (- tru 1) (+ rem divisor))))
(t (values tru rem))))))))
(defun fround (number &optional (divisor 1))
"Same as ROUND, but returns first value as a float."
(multiple-value-bind (res rem)
(round number divisor)
(values (float res (if (floatp rem) rem 1.0)) rem)))
(defmethod qfloat-sign ((a real) &optional (f (float 1 a)))
(cl:float-sign a f))
(defmethod qfloat-sign ((a qd-real) &optional f)
(if f
(make-instance 'qd-real
:value (mul-qd-d (abs-qd (qd-value f))
(cl:float-sign (qd-0 (qd-value a)))))
(make-instance 'qd-real :value (make-qd-d (cl:float-sign (qd-0 (qd-value a)))))))
(declaim (inline float-sign))
(defun float-sign (n &optional (float2 nil float2p))
(if float2p
(qfloat-sign n float2)
(qfloat-sign n)))
(defun max (number &rest more-numbers)
"Returns the greatest of its arguments."
(declare (optimize (safety 2)) (type (or real qd-real) number)
(dynamic-extent more-numbers))
(dolist (real more-numbers)
(when (> real number)
(setq number real)))
number)
(defun min (number &rest more-numbers)
"Returns the least of its arguments."
(declare (optimize (safety 2)) (type (or real qd-real) number)
(dynamic-extent more-numbers))
(do ((nlist more-numbers (cdr nlist))
(result (the (or real qd-real) number)))
((null nlist) (return result))
(declare (list nlist))
(if (< (car nlist) result)
(setq result (car nlist)))))
(defmethod asin ((x number))
(defmethod asin ((x qd-real))
(if (<= -1 x 1)
(make-instance 'qd-real :value (asin-qd (qd-value x)))
(qd-complex-asin x)))
(defmethod acos ((x number))
(defmethod acos ((x qd-real))
(cond ((> (abs x) 1)
(qd-complex-acos x))
(t
(make-instance 'qd-real :value (acos-qd (qd-value x))))))
(defmethod acosh ((x number))
(defmethod acosh ((x qd-real))
(if (< x 1)
(qd-complex-acosh x)
(make-instance 'qd-real :value (acosh-qd (qd-value x)))))
(defmethod atanh ((x number))
(defmethod atanh ((x qd-real))
(if (> (abs x) 1)
(make-instance 'qd-real :value (atanh-qd (qd-value x)))))
(defmethod cis ((x real))
(defmethod cis ((x qd-real))
(multiple-value-bind (s c)
(sincos-qd (qd-value x))
(make-instance 'qd-complex
:real c
:imag s)))
(defmethod phase ((x number))
(defmethod phase ((x qd-real))
(defun signum (number)
"If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
(if (zerop number)
number
(if (rationalp number)
(if (plusp number) 1 -1)
(/ number (abs number)))))
(defmethod random ((x cl:real) &optional (state *random-state*))
(cl:random x state))
(defmethod random ((x qd-real) &optional (state *random-state*))
(* x (make-instance 'qd-real
:value (octi:random-qd state))))
(defmethod float-digits ((x cl:real))
(cl:float-digits x))
(defmethod float-digits ((x qd-real))
(* 4 (float-digits 1d0)))
(defmethod rational ((x real))
(cl:rational x))
(defmethod rational ((x qd-real))
(with-qd-parts (x0 x1 x2 x3)
(qd-value x)
(+ (cl:rational x0)
(cl:rational x1)
(cl:rational x2)
(cl:rational x3))))
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
(defmethod rationalize ((x cl:real))
(cl:rationalize x))
;;; The algorithm here is the method described in CLISP. Bruno Haible has
;;; graciously given permission to use this algorithm. He says, "You can use
;;; it, if you present the following explanation of the algorithm."
;;;
;;; Algorithm (recursively presented):
;;; If x is a rational number, return x.
;;; If x = 0.0, return 0.
;;; If x < 0.0, return (- (rationalize (- x))).
;;; If x > 0.0:
;;; Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
;;; exponent, sign).
;;; If m = 0 or e >= 0: return x = m*2^e.
;;; Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e
;;; with smallest possible numerator and denominator.
;;; Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
;;; But in this case the result will be x itself anyway, regardless of
;;; the choice of a. Therefore we can simply ignore this case.
;;; Note 2: At first, we need to consider the closed interval [a,b].
;;; but since a and b have the denominator 2^(|e|+1) whereas x itself
;;; has a denominator <= 2^|e|, we can restrict the seach to the open
;;; interval (a,b).
;;; So, for given a and b (0 < a < b) we are searching a rational number
;;; y with a <= y <= b.
;;; Recursive algorithm fraction_between(a,b):
;;; c := (ceiling a)
;;; if c < b
;;; then return c ; because a <= c < b, c integer
;;; else
;;; ; a is not integer (otherwise we would have had c = a < b)
;;; k := c-1 ; k = floor(a), k < a < b <= k+1
;;; return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
;;; ; note 1 <= 1/(b-k) < 1/(a-k)
;;;
;;; You can see that we are actually computing a continued fraction expansion.
;;;
;;; Algorithm (iterative):
;;; If x is rational, return x.
;;; Call (integer-decode-float x). It returns a m,e,s (mantissa,
;;; exponent, sign).
;;; If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.)
;;; Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1)
;;; (positive and already in lowest terms because the denominator is a
;;; power of two and the numerator is odd).
;;; Start a continued fraction expansion
;;; p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
;;; Loop
;;; c := (ceiling a)
;;; if c >= b
;;; then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
;;; goto Loop
;;; finally partial_quotient(c).
;;; Here partial_quotient(c) denotes the iteration
;;; i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].
;;; At the end, return s * (p[i]/q[i]).
;;; This rational number is already in lowest terms because
;;; p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.
;;;
(defmethod rationalize ((x qd-real))
;; This is a fairly straigtforward implementation of the iterative
;; algorithm above.
(multiple-value-bind (frac expo sign)
(integer-decode-float x)
(cond ((or (zerop frac) (>= expo 0))
(if (minusp sign)
(- (ash frac expo))
(ash frac expo)))
(t
;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e),
;; so build the fraction up immediately, without having to do
;; a gcd.
(let ((a (/ (- (* 2 frac) 1) (ash 1 (- 1 expo))))
(b (/ (+ (* 2 frac) 1) (ash 1 (- 1 expo))))
(p0 0)
(q0 1)
(p1 1)
(q1 0))
(do ((c (ceiling a) (ceiling a)))
((< c b)
(let ((top (+ (* c p1) p0))
(bot (+ (* c q1) q0)))
(/ (if (minusp sign)
(- top)
top)
bot)))
(let* ((k (- c 1))
(p2 (+ (* k p1) p0))
(q2 (+ (* k q1) q0)))
(psetf a (/ (- b k))
b (/ (- a k)))
(setf p0 p1
q0 q1
p1 p2
q1 q2))))))))
(define-compiler-macro + (&whole form &rest args)
(declare (ignore form))
(if (null args)
0
(do ((args (cdr args) (cdr args))
(res (car args)
`(two-arg-+ ,res ,(car args))))
((null args) res))))
(define-compiler-macro - (&whole form number &rest more-numbers)
(declare (ignore form))
(if more-numbers
(do ((nlist more-numbers (cdr nlist))
(result number))
((atom nlist) result)
(declare (list nlist))
(setq result `(two-arg-- ,result ,(car nlist))))
`(unary-minus ,number)))
(define-compiler-macro * (&whole form &rest args)
(declare (ignore form))
(if (null args)
1
(do ((args (cdr args) (cdr args))
(res (car args)
`(two-arg-* ,res ,(car args))))
((null args) res))))
(define-compiler-macro / (number &rest more-numbers)
(if more-numbers
(do ((nlist more-numbers (cdr nlist))
(result number))
((atom nlist) result)
(declare (list nlist))
(setq result `(two-arg-/ ,result ,(car nlist))))
`(unary-divide ,number)))
;; Compiler macros to convert <, >, <=, and >= into multiple calls of
;; the corresponding two-arg-<foo> function.
(macrolet
((frob (op)
(let ((method (intern (concatenate 'string
(string '#:two-arg-)
(symbol-name op)))))
`(define-compiler-macro ,op (number &rest more-numbers)
(do* ((n number (car nlist))
(nlist more-numbers (cdr nlist))