1 ;;;; -*- Mode: lisp -*-
3 ;;;; Copyright (c) 2007 Raymond Toy
5 ;;;; Permission is hereby granted, free of charge, to any person
6 ;;;; obtaining a copy of this software and associated documentation
7 ;;;; files (the "Software"), to deal in the Software without
8 ;;;; restriction, including without limitation the rights to use,
9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10 ;;;; copies of the Software, and to permit persons to whom the
11 ;;;; Software is furnished to do so, subject to the following
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
26 ;; Most of this code taken from CMUCL and slightly modified to support
31 (defmethod two-arg-/ ((a qd-real
) (b rational
))
32 (make-instance 'qd-real
:value
(div-qd (qd-value a
)
33 (qd-value (float b
#q0
)))))
35 (defmethod two-arg-/ ((a rational
) (b qd-real
))
36 (make-instance 'qd-real
:value
(div-qd (qd-value (float a
#q0
))
39 (defmethod two-arg-* ((a qd-real
) (b rational
))
40 (make-instance 'qd-real
:value
(mul-qd (qd-value a
) (qd-value (float b
#q0
)))))
42 (defmethod two-arg-+ ((a qd-real
) (b rational
))
43 (make-instance 'qd-real
:value
(add-qd (qd-value a
) (qd-value (float b
#q0
)))))
45 (defmethod two-arg-+ ((a rational
) (b qd-real
))
46 (make-instance 'qd-real
:value
(add-qd (qd-value b
) (qd-value (float a
#q0
)))))
48 (defmethod two-arg-- ((a qd-real
) (b rational
))
49 (make-instance 'qd-real
:value
(sub-qd (qd-value a
) (qd-value (float b
#q0
)))))
51 (defmethod two-arg-- ((a rational
) (b qd-real
))
52 (make-instance 'qd-real
:value
(sub-qd (qd-value (float a
#q0
)) (qd-value b
))))
54 (defmethod unary-minus ((z qd-complex
))
55 (complex (- (realpart z
))
58 (defmethod qzerop ((z qd-complex
))
59 (and (zerop (realpart z
))
60 (zerop (imagpart z
))))
62 (defmethod two-arg-+ ((a qd-complex
) (b qd-complex
))
63 (complex (+ (realpart a
) (realpart b
))
64 (+ (imagpart a
) (imagpart b
))))
66 (defmethod two-arg-+ ((a qd-complex
) (b real
))
67 (complex (+ (realpart a
) b
)
70 (defmethod two-arg-+ ((a qd-complex
) (b qd-real
))
71 (complex (+ (realpart a
) b
)
74 (defmethod two-arg-+ ((a qd-real
) (b qd-complex
))
75 (complex (+ a
(realpart b
))
78 (defmethod two-arg-+ ((a qd-complex
) (b cl
:complex
))
79 (complex (+ (realpart a
) (imagpart b
))
80 (+ (imagpart a
) (imagpart b
))))
82 (defmethod two-arg-+ ((a number
) (b qd-complex
))
85 (defmethod two-arg-- ((a qd-complex
) (b qd-complex
))
86 (complex (- (realpart a
) (realpart b
))
87 (- (imagpart a
) (imagpart b
))))
89 (defmethod two-arg-- ((a qd-complex
) (b real
))
90 (complex (- (realpart a
) b
)
93 (defmethod two-arg-- ((a qd-complex
) (b cl
:complex
))
94 (complex (- (realpart a
) (realpart b
))
95 (- (imagpart a
) (imagpart b
))))
97 (defmethod two-arg-- ((a qd-complex
) (b qd-real
))
98 (complex (- (realpart a
) b
)
101 (defmethod two-arg-- ((a number
) (b qd-complex
))
102 (complex (- (realpart a
) (realpart b
))
103 (- (imagpart a
) (imagpart b
))))
105 (defmethod two-arg-- ((a qd-real
) (b qd-complex
))
106 (complex (- a
(realpart b
))
109 (defmethod two-arg-* ((a qd-complex
) (b qd-complex
))
110 (let* ((rx (realpart a
))
114 (complex (- (* rx ry
) (* ix iy
))
115 (+ (* rx iy
) (* ix ry
)))))
117 (defmethod two-arg-* ((a qd-complex
) (b real
))
118 (let* ((rx (realpart a
))
123 (defmethod two-arg-* ((a qd-complex
) (b qd-real
))
124 (let* ((rx (realpart a
))
129 (defmethod two-arg-* ((a qd-real
) (b qd-complex
))
135 (defmethod two-arg-* ((a qd-complex
) (b cl
:complex
))
136 ;; For now, convert B into a qd-complex and use that.
137 (let ((re (coerce (realpart b
) 'ext
:double-double-float
))
138 (im (coerce (imagpart b
) 'ext
:double-double-float
)))
139 (two-arg-* a
(make-instance 'qd-complex
140 :real
(make-qd-dd re
0w0
)
141 :imag
(make-qd-dd im
0w0
)))))
144 (defmethod two-arg-* ((a qd-complex
) (b cl
:complex
))
145 ;; For now, convert B into a qd-complex and use that.
146 (let ((re (coerce (realpart b
) 'double-float
))
147 (im (coerce (imagpart b
) 'double-float
)))
148 (two-arg-* a
(make-instance 'qd-complex
150 :imag
(make-qd-d im
)))))
152 (defmethod two-arg-* ((a number
) (b qd-complex
))
155 (defmethod two-arg-/ ((x qd-complex
) (y qd-complex
))
156 (let* ((rx (realpart x
))
160 (if (> (abs ry
) (abs iy
))
162 (dn (+ ry
(* r iy
))))
163 (complex (/ (+ rx
(* ix r
)) dn
)
164 (/ (- ix
(* rx r
)) dn
)))
166 (dn (+ iy
(* r ry
))))
167 (complex (/ (+ (* rx r
) ix
) dn
)
168 (/ (- (* ix r
) rx
) dn
))))))
170 (defmethod two-arg-/ ((x qd-complex
) (y qd-real
))
171 (complex (/ (realpart x
) y
)
174 (defmethod two-arg-/ ((x qd-complex
) (y number
))
175 (complex (/ (realpart x
) y
)
178 (defmethod two-arg-/ ((x number
) (y qd-complex
))
179 (let* ((rx (realpart x
))
183 (if (> (abs ry
) (abs iy
))
185 (dn (+ ry
(* r iy
))))
186 (complex (/ (+ rx
(* ix r
)) dn
)
187 (/ (- ix
(* rx r
)) dn
)))
189 (dn (+ iy
(* r ry
))))
190 (complex (/ (+ (* rx r
) ix
) dn
)
191 (/ (- (* ix r
) rx
) dn
))))))
193 (defmethod two-arg-/ ((x qd-real
) (y qd-complex
))
194 ;; This can be simplified since X is real.
195 (let* ((rx (realpart x
))
199 (if (> (abs ry
) (abs iy
))
201 (dn (+ ry
(* r iy
))))
202 (complex (/ (+ rx
(* ix r
)) dn
)
203 (/ (- ix
(* rx r
)) dn
)))
205 (dn (+ iy
(* r ry
))))
206 (complex (/ (+ (* rx r
) ix
) dn
)
207 (/ (- (* ix r
) rx
) dn
))))))
209 (defmethod unary-divide ((a qd-complex
))
212 (defmethod coerce ((obj t
) (type t
))
213 (cl:coerce obj type
))
215 (defmethod coerce ((number cl
:real
) (type (eql 'qd-real
)))
218 (defmethod coerce ((number qd-real
) (type (eql 'qd-real
)))
221 (defmethod coerce ((number cl
:number
) (type (eql 'qd-complex
)))
222 (complex (float (realpart number
) #q0
)
223 (float (imagpart number
) #q0
)))
225 (defmethod coerce ((number qd-complex
) (type (eql 'qd-complex
)))
228 (declaim (inline square
))
230 (declare (type qd-real x
))
231 (make-instance 'qd-real
:value
(sqr-qd (qd-value x
))))
233 (defun qd-complex-sqrt (z)
234 "Principle square root of Z
236 Z may be any number, but the result is always a complex."
237 (declare (type qd-complex z
))
238 (multiple-value-bind (rho k
)
240 (declare (type qd-real rho
)
242 (let ((x (realpart z
))
246 (declare (type qd-real x y eta nu
))
249 ;; space 0 to get maybe-inline functions inlined.
250 (declare (optimize (speed 3) (space 0)))
252 (setf rho
(+ (scalb (abs x
) (- k
)) (sqrt rho
)))
257 (setf k
(1- (ash k -
1)))
258 (setf rho
(+ rho rho
))))
260 (setf rho
(scalb (sqrt rho
) k
))
265 (when (not (zerop rho
))
266 (setf nu
(/ (/ nu rho
) 2d0
)))
269 (setf nu
(float-sign y rho
))))
272 (defun qd-complex-log-scaled (z j
)
275 This is for use with J /= 0 only when |z| is huge."
276 (declare (type qd-complex z
)
278 ;; The constants t0, t1, t2 should be evaluated to machine
279 ;; precision. In addition, Kahan says the accuracy of log1p
280 ;; influences the choices of these constants but doesn't say how to
281 ;; choose them. We'll just assume his choices matches our
282 ;; implementation of log1p.
283 (let ((t0 (/ 1 (sqrt #q2.0q0
)))
289 (multiple-value-bind (rho k
)
291 (declare (optimize (speed 3)))
292 (let ((beta (max (abs x
) (abs y
)))
293 (theta (min (abs x
) (abs y
))))
294 (complex (if (and (zerop k
)
298 (/ (log1p (+ (* (- beta
1.0d0
)
306 (defun qd-complex-log (z)
307 "Log of Z = log |Z| + i * arg Z
309 Z may be any number, but the result is always a complex."
310 (declare (type qd-complex z
))
311 (qd-complex-log-scaled z
0))
314 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
315 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
316 ;; reason for the imaginary part is caused by the fact that arg i*y is
317 ;; never 0 since we have positive and negative zeroes.
320 ;; atanh(z) = (log(1+z) - log(1-z))/2
322 ;; The branch cut is on the real axis for |x| >= 1. For x =< -1,
323 ;; atanh is continuous with quadrant III; for x >= 1, continuous with
325 (defun qd-complex-atanh (z)
326 "Compute atanh z = (log(1+z) - log(1-z))/2"
327 (declare (type (or qd-real qd-complex
) z
))
328 (cond ((and (typep z
'qd-real
) (< z -
1))
329 (qd-complex-atanh (complex z -
0d0
)))
331 (flet ((careful-mul (a b
)
332 ;; Carefully multiply a and b, taking care to handle
333 ;; signed zeroes. Only need to handle the case of b
336 (if (minusp (* (float-sign a
) (float-sign b
)))
341 (theta (/ (sqrt most-positive-double-float
) 4.0d0
))
342 (rho (/ 4.0d0
(sqrt most-positive-double-float
)))
343 (half-pi #.
(/ +pi
+ 2d0
))
345 (beta (float-sign rp
))
347 (y (careful-mul beta
(- (imagpart z
))))
350 ;; Shouldn't need this declare.
351 (declare (type qd-real x y
))
353 (declare (optimize (speed 3)))
354 (cond ((or (> x theta
)
356 ;; To avoid overflow...
357 (setf nu
(float-sign y half-pi
))
358 ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
359 ;; which can cause overflow. Arrange this computation so
360 ;; that it won't overflow.
361 (setf eta
(let* ((x-bigger (> x
(abs y
)))
362 (r (if x-bigger
(/ y x
) (/ x y
)))
363 (d (+ 1.0d0
(* r r
))))
368 ;; Should this be changed so that if y is zero, eta is set
369 ;; to +infinity instead of approx 176? In any case
370 ;; tanh(176) is 1.0d0 within working precision.
371 (let ((t1 (+ 4d0
(square y
)))
372 (t2 (+ (abs y
) rho
)))
373 (setf eta
(log (/ (sqrt (sqrt t1
))
377 (+ half-pi
(atan (* 0.5d0 t2
))))))
380 (let ((t1 (+ (abs y
) rho
)))
381 ;; Normal case using log1p(x) = log(1 + x)
383 (log1p (/ (* 4.0d0 x
)
384 (+ (square (- 1.0d0 x
))
387 (atan (careful-mul 2.0d0 y
)
391 (complex (* beta eta
)
392 (- (* beta nu
)))))))))
395 ;; tanh(z) = sinh(z)/cosh(z)
397 (defun qd-complex-tanh (z)
398 "Compute tanh z = sinh z / cosh z"
399 (declare (type (or qd-real qd-complex
) z
))
400 (let ((x (realpart z
))
403 ;; space 0 to get maybe-inline functions inlined
404 (declare (optimize (speed 3) (space 0)))
405 (cond ((> (abs x
) #.
(/ (+ (log most-positive-double-float
)
408 ;; The threshold above is
409 ;; asinh(most-positive-double-float)/4, but many Lisps
410 ;; cannot actually compute that. Hence use the
411 ;; (accurate) approximation
412 ;; asinh(most-positive-double-float) =
413 ;; log(most-positive-double-float) + log(2)
414 (complex (float-sign x
)
417 ;; With quad-double's it happens that tan(pi/2) will
418 ;; actually produce a division by zero error. We need to
419 ;; handle that case carefully.
420 (let* ((tv (ignore-errors (tan y
)))
422 (rho (sqrt (+ 1.0d0
(* s s
)))))
424 (let* ((beta (+ 1.0d0
(* tv tv
)))
425 (den (+ 1.0d0
(* beta s s
))))
426 (complex (/ (* beta rho s
)
430 ;; This means tan(y) produced some error. We'll
431 ;; assume it's an overflow error because y is
432 ;; pi/2 + 2*k*pi. But we need a value for tv to
433 ;; compute (/ tv). This would be a signed-zero.
434 ;; For now, just return +0.
438 ;; Kahan says we should only compute the parts needed. Thus, the
439 ;; realpart's below should only compute the real part, not the whole
440 ;; complex expression. Doing this can be important because we may get
441 ;; spurious signals that occur in the part that we are not using.
443 ;; However, we take a pragmatic approach and just use the whole
446 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
447 ;; it's the conjugate of the square root or the square root of the
448 ;; conjugate. This needs to be checked.
450 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
451 ;; same as (sqrt (conjugate z)) for all z. This follows because
453 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
455 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
457 ;; and these two expressions are equal if and only if arg conj z =
458 ;; -arg z, which is clearly true for all z.
460 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
461 ;; complex, the real is converted to a complex before performing the
462 ;; operation. However, Kahan says in this paper (pg 176):
464 ;; (iii) Careless handling can turn infinity or the sign of zero into
465 ;; misinformation that subsequently disappears leaving behind
466 ;; only a plausible but incorrect result. That is why compilers
467 ;; must not transform z-1 into z-(1+i*0), as we have seen above,
468 ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
469 ;; subsequent logarithm or square root produce a non-zero
470 ;; imaginary part whose sign is opposite to what was intended.
472 ;; The interesting examples are too long and complicated to reproduce
473 ;; here. We refer the reader to his paper.
475 ;; The functions below are intended to handle the cases where a real
476 ;; is mixed with a complex and we don't want CL complex contagion to
479 (declaim (inline 1+z
1-z z-1 z
+1))
481 (complex (+ 1 (realpart z
)) (imagpart z
)))
483 (complex (- 1 (realpart z
)) (- (imagpart z
))))
485 (complex (- (realpart z
) 1) (imagpart z
)))
487 (complex (+ (realpart z
) 1) (imagpart z
)))
489 (defun qd-complex-acos (z)
490 "Compute acos z = pi/2 - asin z
492 Z may be any number, but the result is always a complex."
493 (declare (type (or qd-real qd-complex
) z
))
494 (if (and (typep z
'qd-real
) (> z
1))
495 ;; acos is continuous in quadrant IV in this case.
496 (qd-complex-acos (complex z -
0f0
))
497 (let ((sqrt-1+z
(qd-complex-sqrt (1+z z
)))
498 (sqrt-1-z (qd-complex-sqrt (1-z z
))))
499 (cond ((zerop (realpart sqrt-1
+z
))
500 ;; Same as below, but we compute atan ourselves (because we
501 ;; have atan +/- infinity).
503 (if (minusp (float-sign (* (realpart sqrt-1-z
)
504 (realpart sqrt-1
+z
))))
507 (asinh (imagpart (* (conjugate sqrt-1
+z
)
510 (complex (* 2 (atan (/ (realpart sqrt-1-z
)
511 (realpart sqrt-1
+z
))))
512 (asinh (imagpart (* (conjugate sqrt-1
+z
)
515 (defun qd-complex-acosh (z)
516 "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
518 Z may be any number, but the result is always a complex."
519 (declare (type (or qd-real qd-complex
) z
))
520 (let* ((sqrt-z-1 (qd-complex-sqrt (z-1 z
)))
521 (sqrt-z+1 (qd-complex-sqrt (z+1 z
))))
522 ;; We need to handle the case where real part of sqrt-z+1 is zero,
523 ;; because division by zero with double-double-floats doesn't
525 (cond ((zerop (realpart sqrt-z
+1))
526 ;; Same as below, but we compute atan ourselves (because we
527 ;; have atan +/- infinity).
528 (complex (asinh (realpart (* (conjugate sqrt-z-1
)
530 (if (minusp (float-sign (* (imagpart sqrt-z-1
)
531 (realpart sqrt-z
+1))))
535 (complex (asinh (realpart (* (conjugate sqrt-z-1
)
537 (* 2 (atan (/ (imagpart sqrt-z-1
)
538 (realpart sqrt-z
+1)))))))))
540 ;; asin(z) = asinh(i*z)/i
541 ;; = -i log(i*z + sqrt(1-z^2))
542 (defun qd-complex-asin (z)
543 "Compute asin z = asinh(i*z)/i
545 Z may be any number, but the result is always a complex."
546 (declare (type (or qd-real qd-complex
) z
))
547 (cond ((and (typep z
'qd-real
) (> z
1))
548 (qd-complex-asin (complex z -
0d0
)))
550 (let* ((sqrt-1-z (qd-complex-sqrt (1-z z
)))
551 (sqrt-1+z
(qd-complex-sqrt (1+z z
)))
552 (den (realpart (* sqrt-1-z sqrt-1
+z
))))
554 ;; Like below but we handle atan part ourselves.
555 ;; Must be sure to take into account the sign of
556 ;; (realpart z) and den!
557 (complex (if (minusp (* (float-sign (realpart z
))
561 (asinh (imagpart (* (conjugate sqrt-1-z
)
564 ;; We get a invalid operation here when z is real and |z| > 1.
565 (complex (atan (/ (realpart z
)
566 (realpart (* sqrt-1-z sqrt-1
+z
))))
567 (asinh (imagpart (* (conjugate sqrt-1-z
)
570 (defun qd-complex-asinh (z)
571 "Compute asinh z = log(z + sqrt(1 + z*z))
573 Z may be any number, but the result is always a complex."
574 (declare (type (or qd-real qd-complex
) z
))
575 ;; asinh z = -i * asin (i*z)
576 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
577 (result (qd-complex-asin iz
)))
578 (complex (imagpart result
)
579 (- (realpart result
)))))
581 (defun qd-complex-atan (z)
582 "Compute atan z = atanh (i*z) / i
584 Z may be any number, but the result is always a complex."
585 (declare (type (or qd-real qd-complex
) z
))
586 ;; atan z = -i * atanh (i*z)
587 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
588 (result (qd-complex-atanh iz
)))
589 (complex (imagpart result
)
590 (- (realpart result
)))))
592 (defun qd-complex-tan (z)
593 "Compute tan z = -i * tanh(i * z)
595 Z may be any number, but the result is always a complex."
596 (declare (type (or qd-real qd-complex
) z
))
597 ;; tan z = -i * tanh(i*z)
598 (let* ((iz (complex (- (imagpart z
)) (realpart z
)))
599 (result (qd-complex-tanh iz
)))
600 (complex (imagpart result
)
601 (- (realpart result
)))))
603 (defmethod qasin ((x qd-complex
))
606 (defmethod qacos ((x qd-complex
))
609 (defmethod qacosh ((x qd-complex
))
610 (qd-complex-acosh x
))
612 (defmethod qatanh ((x qd-complex
))
613 (qd-complex-atanh x
))
615 (defmethod qsin ((z qd-complex
))
616 (let ((x (realpart z
))
618 (complex (* (sin x
) (cosh y
))
619 (* (cos x
) (sinh y
)))))
621 (defmethod qcos ((z qd-complex
))
622 (let ((x (realpart z
))
624 (complex (* (cos x
) (cosh y
))
625 (- (* (sin x
) (sinh y
))))))
627 (defmethod qtan ((z qd-complex
))
630 (defmethod qsinh ((z qd-complex
))
631 (let ((x (realpart z
))
633 (complex (* (sinh x
) (cos y
))
634 (* (cosh x
) (sin y
)))))
636 (defmethod qcosh ((z qd-complex
))
637 (let ((x (realpart z
))
639 (complex (* (cosh x
) (cos y
))
640 (* (sinh x
) (sin y
)))))
642 (defmethod qtanh ((z qd-complex
))
645 (defmethod qsqrt ((z qd-complex
))
648 (defmethod qatan ((y qd-complex
) &optional x
)
650 (error "First arg of 2-arg ATAN must be real")
651 (qd-complex-atan y
)))
653 (defmethod qatan ((y cl
:complex
) &optional x
)
655 (error "First arg of 2-arg ATAN must be real")
658 (defmethod qexp ((z qd-complex
))
659 (let* ((x (realpart z
))
662 (complex (* ex
(cos y
))
665 (defmethod qlog ((a qd-complex
) &optional b
)
667 (/ (qlog a
) (qlog b
))
668 (complex (log (abs a
))
669 (atan (imagpart a
) (realpart a
)))))
671 (defmethod qexpt ((x qd-complex
) (y number
))
672 (exp (* (float y
#q0
) (log x
))))
674 (defmethod qexpt ((x number
) (y qd-complex
))
675 (exp (* y
(log (float x
#q0
)))))
677 (defmethod qexpt ((x qd-complex
) (y qd-complex
))
680 (defmethod qphase ((z qd-complex
))
681 (atan (imagpart z
) (realpart z
)))
684 (or (typep x
'qd-real
)
688 (or (typep x
'qd-complex
)