qd-methods.lisp 32.2 KB
;;;; -*- Mode: lisp -*-
;;;;
;;;; 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.

(in-package #:oct)

(defmethod make-qd ((x cl:rational))
;; We should do something better than this.
(make-instance 'qd-real :value (rational-to-qd x)))

;; 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))))

(cl::1+ a))

(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)

(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))
(+ b a))

(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))
(* b a))

(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)
b)))

(defmethod two-arg-/ ((a cl:float) (b qd-real))
(make-instance 'qd-real :value (div-qd (make-qd-d (cl:float a 1d0))
(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))))
(cl-name (intern (symbol-name name) :cl))
(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))))))
(frob zerop number)
(frob plusp)
(frob minusp))

(defun bignum-to-qd (bignum)
(make-instance 'qd-real
:value (rational-to-qd bignum)))

(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))))

#+nil
(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)))

(defmethod qfloat ((x ratio) (num-type qd-real))
(make-instance 'qd-real :value (rational-to-qd x)))

#+cmu
(defmethod qfloat ((x ext:double-double-float) (num-type qd-real))
(make-instance 'qd-real :value (make-qd-dd x 0w0)))

(defmethod qfloat ((x qd-real) (num-type cl:float))
(multiple-value-bind (q0 q1 q2 q3)
(qd-parts (qd-value x))
(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)
(qd-parts (qd-value x))
(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 num-type)
(qfloat x num-type))

(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-qd 0d0))
(defmethod qimagpart ((x qd-complex))
(make-instance 'qd-real :value (qd-imag x)))

(defun imagpart (x)
(qimagpart 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))
'#:octi)))
(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)))))
(progn
(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))
(cl:sqrt x))

(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))

(declaim (inline qd-cssqs))
(defun qd-cssqs (z)
(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)))

#+nil
(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)))

(defmethod abs ((z qd-complex))
;; 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)
(if b
(cl:log a b)
(cl:log a)))

(defmethod log ((a qd-real) &optional b)
(if b
(/ (log a) (log 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))))))

(defmethod log1p ((a qd-real))
(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 real))
(exp (* y (log x))))

(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
:value (npow (qd-value x) y)))

(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)
(qd-= (qd-value a) (make-qd-d (cl:float b 1d0)))
nil))

(defmethod two-arg-= ((a number) (b qd-real))
(if (cl:realp a)
(qd-= (make-qd-d (cl:float a 1d0)) (qd-value b))
nil))

(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))
(qd-= (qd-value a) (qd-value b)))

(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)))
((atom nlist) t)
(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))
(nlist more-numbers (cdr nlist)))
((atom nlist) t)
(declare (list nlist))
(unless (do* ((nl nlist (cdr nl)))
((atom nl) t)
(declare (list 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)
:imag (make-qd-d y)))

(defun complex (x &optional (y 0))
(qcomplex x y))

(defmethod qinteger-decode-float ((f cl:float))
(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))

(defmethod qdecode-float ((f cl:float))
(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))))
(let ((f (ffloor-qd (qd-value x))))
(multiple-value-bind (int exp sign)
(integer-decode-qd f)
(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)
(values f
rem)
(values (+ f 1)
(- rem 1)))))

(defun fceiling (x &optional y)
(multiple-value-bind (f rem)
(ffloor x y)
(if (zerop rem)
(values f
rem)
(values (+ f 1)
(- rem 1)))))

(defun truncate (x &optional (y 1))
(if (minusp x)
(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))
(if (minusp x)
(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)
(oddp shifted)
(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))
(cl:asin x))

(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))
(cl:acos x))

(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))
(cl:acosh x))

(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))
(cl:atanh x))

(defmethod atanh ((x qd-real))
(if (> (abs x) 1)
(qd-complex-atanh x)
(make-instance 'qd-real :value (atanh-qd (qd-value x)))))

(defmethod cis ((x real))
(cl:cis x))

(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))
(cl:phase x))

(defmethod phase ((x qd-real))
(if (minusp x)
(- +pi+)
(make-instance 'qd-real :value (make-qd-d 0d0))))

(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))))

(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))
(res nil))
((atom nlist)
(and ,@(nreverse res)))
(push (,',method ,n ,(car nlist)) res))))))
(frob <)
(frob >)
(frob <=)
(frob >=))

(define-compiler-macro /= (&whole form number &rest more-numbers)
;; Convert (/= x y) to (not (two-arg-= x y)).  Should we try to
;; handle a few more cases?
(if (cdr more-numbers)
form
(not (two-arg-= ,number ,(car more-numbers)))))

;; Define compiler macro the convert two-arg-foo into the appropriate
;; CL function or QD-REAL function so we don't have to do CLOS
;; dispatch.
#+(or)
(macrolet
((frob (name cl-op qd-op)
(define-compiler-macro ,name (&whole form x y &environment env)
(flet ((arg-type (arg)
(multiple-value-bind (def-type localp decl)
(ext:variable-information arg env)
(declare (ignore localp))
(when def-type
(cdr (assoc 'type decl))))))
(let ((x-type (arg-type x))
(y-type (arg-type y)))
(cond ((and (subtypep x-type 'cl:number)
(subtypep y-type 'cl:number))
(,',cl-op ,x ,y))
((and (subtypep x-type 'qd-real)
(subtypep y-type 'qd-real))
(make-instance 'qd-real :value (,',qd-op (qd-value ,x)
(qd-value ,y))))
(t
;; Don't know how to handle this, so give up.
form)))))))
(frob two-arg-- cl:- sub-qd)
(frob two-arg-* cl:* mul-qd)
(frob two-arg-/ cl:/ div-qd))

#+(or)
(macrolet
((frob (name cl-op qd-op cl-qd-op qd-cl-op)
(define-compiler-macro ,name (&whole form x y &environment env)
(flet ((arg-type (arg)
(multiple-value-bind (def-type localp decl)
(ext:variable-information arg env)
(declare (ignore localp))
(when def-type
(cdr (assoc 'type decl))))))
(let ((x-type (arg-type x))
(y-type (arg-type y)))
(cond ((subtypep x-type 'cl:float)
(cond ((subtypep y-type 'cl:number)
(,',cl-op ,x ,y))
((subtypep y-type 'qd-real)
(if ,cl-qd-op
(make-instance 'qd-real :value (,',cl-qd-op (cl:float ,x 1d0)
(qd-value ,y)))
form))
(t form)))
((subtypep x-type 'qd-real)
(cond ((subtypep y-type 'cl:float)
(if ,qd-cl-op
(make-instance 'qd-real :value (,',qd-cl-op (qd-value ,x)
(float ,y 1d0)))
form))
((subtypep y-type 'qd-real)
(make-instance 'qd-real :value (,',qd-op (qd-value ,x)
(qd-value ,y))))
(t form)))
(t
;; Don't know how to handle this, so give up.
form)))))))
(frob two-arg-- cl:- sub-qd sub-d-qd sub-qd-d)
(frob two-arg-* cl:* mul-qd mul-d-qd mul-qd-d)
(frob two-arg-/ cl:/ div-qd nil nil))

(defgeneric epsilon (m)
(:documentation
"Return an epsilon value of the same precision as the argument.  It is
the smallest number x such that 1+x /= x.  The argument can be
complex"))

(defmethod epsilon ((m cl:float))
(etypecase m
(single-float single-float-epsilon)
(double-float double-float-epsilon)))

(defmethod epsilon ((m cl:complex))
(epsilon (realpart m)))

(defmethod epsilon ((m qd-real))
;; What is the epsilon value for a quad-double?  This is complicated
;; by the fact that things like (+ #q1 #q1q-100) is representable as
;; a quad-double.  For most purposes we want epsilon to be close to
;; the 212 bits of precision (4*53 bits) that we normally have with
(scale-float +qd-real-one+ -212))

(defmethod epsilon ((m qd-complex))
(epsilon (realpart m)))

(defgeneric float-pi (x)
(:documentation
"Return a floating-point value of the mathematical constant pi that is
the same precision as the argument.  The argument can be complex."))

(defmethod float-pi ((x cl:rational))
(float pi 1f0))

(defmethod float-pi ((x cl:float))
(float pi x))

(defmethod float-pi ((x qd-real))
+pi+)

(defmethod float-pi ((z cl:complex))
(float pi (realpart z)))

(defmethod float-pi ((z qd-complex))
+pi+)

(define-condition domain-error (simple-error)
((function-name :accessor condition-function-name
:initarg :function-name))
(:report (lambda (condition stream)
(format stream "Domain Error for function ~S:~&"
(condition-function-name condition))
(pprint-logical-block (stream nil :per-line-prefix "  ")
(apply #'format stream
(simple-condition-format-control condition)
(simple-condition-format-arguments condition))))))`