qd-complex.lisp 23.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
;;;; -*- Mode: lisp -*-
;;;;
;;;; Copyright (c) 2007 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.

26
(in-package #:oct)
toy's avatar
toy committed
27

28 29 30
(eval-when (:compile-toplevel :load-toplevel :execute)
  (setf *readtable* *oct-readtable*))

31 32
(defmethod add1 ((a qd-complex))
  (make-instance 'qd-complex
Raymond Toy's avatar
Raymond Toy committed
33 34 35 36 37 38
		 :real (add-qd-d (qd-value (realpart a)) 1d0)
		 :imag (qd-value (imagpart a))))

(defmethod sub1 ((a qd-complex))
  (make-instance 'qd-complex
		 :real (sub-qd-d (qd-value (realpart a)) 1d0)
39 40
		 :imag (qd-value (imagpart a))))
		 
41
(defmethod two-arg-/ ((a qd-real) (b cl:rational))
toy's avatar
toy committed
42 43 44
  (make-instance 'qd-real :value (div-qd (qd-value a)
					 (qd-value (float b #q0)))))

45
(defmethod two-arg-/ ((a cl:rational) (b qd-real))
toy's avatar
toy committed
46 47 48
  (make-instance 'qd-real :value (div-qd (qd-value (float a #q0))
					 (qd-value b))))

49
(defmethod two-arg-* ((a qd-real) (b cl:rational))
toy's avatar
toy committed
50 51
  (make-instance 'qd-real :value (mul-qd (qd-value a) (qd-value (float b #q0)))))

52
(defmethod two-arg-+ ((a qd-real) (b cl:rational))
toy's avatar
toy committed
53 54
  (make-instance 'qd-real :value (add-qd (qd-value a) (qd-value (float b #q0)))))

55
(defmethod two-arg-+ ((a cl:rational) (b qd-real))
toy's avatar
toy committed
56 57
  (make-instance 'qd-real :value (add-qd (qd-value b) (qd-value (float a #q0)))))

58
(defmethod two-arg-- ((a qd-real) (b cl:rational))
toy's avatar
toy committed
59 60
  (make-instance 'qd-real :value (sub-qd (qd-value a) (qd-value (float b #q0)))))

61
(defmethod two-arg-- ((a cl:rational) (b qd-real))
toy's avatar
toy committed
62 63
  (make-instance 'qd-real :value (sub-qd (qd-value (float a #q0)) (qd-value b))))

toy's avatar
toy committed
64 65 66 67
(defmethod unary-minus ((z qd-complex))
  (complex (- (realpart z))
	   (- (imagpart z))))

68 69 70 71
(defmethod qzerop ((z qd-complex))
  (and (zerop (realpart z))
       (zerop (imagpart z))))

72 73 74 75 76 77 78 79
(defmethod two-arg-+ ((a qd-complex) (b qd-complex))
  (complex (+ (realpart a) (realpart b))
	   (+ (imagpart a) (imagpart b))))

(defmethod two-arg-+ ((a qd-complex) (b real))
  (complex (+ (realpart a) b)
	   (imagpart a)))

toy's avatar
toy committed
80 81 82 83 84 85 86 87
(defmethod two-arg-+ ((a qd-complex) (b qd-real))
  (complex (+ (realpart a) b)
	   (imagpart a)))

(defmethod two-arg-+ ((a qd-real) (b qd-complex))
  (complex (+ a (realpart b))
	   (imagpart b)))

88 89 90 91 92 93 94 95 96 97 98 99 100
(defmethod two-arg-+ ((a qd-complex) (b cl:complex))
  (complex (+ (realpart a) (imagpart b))
	   (+ (imagpart a) (imagpart b))))

(defmethod two-arg-+ ((a number) (b qd-complex))
  (two-arg-+ b a))
  
(defmethod two-arg-- ((a qd-complex) (b qd-complex))
  (complex (- (realpart a) (realpart b))
	   (- (imagpart a) (imagpart b))))

(defmethod two-arg-- ((a qd-complex) (b real))
  (complex (- (realpart a) b)
toy's avatar
toy committed
101
	   (imagpart a)))
102 103

(defmethod two-arg-- ((a qd-complex) (b cl:complex))
toy's avatar
toy committed
104
  (complex (- (realpart a) (realpart b))
105 106
	   (- (imagpart a) (imagpart b))))

toy's avatar
toy committed
107 108 109 110
(defmethod two-arg-- ((a qd-complex) (b qd-real))
  (complex (- (realpart a) b)
	   (imagpart a)))

111 112 113
(defmethod two-arg-- ((a number) (b qd-complex))
  (complex (- (realpart a) (realpart b))
	   (- (imagpart a) (imagpart b))))
toy's avatar
toy committed
114 115 116 117

(defmethod two-arg-- ((a qd-real) (b qd-complex))
  (complex (- a (realpart b))
	   (- (imagpart b))))
118
  
119 120 121 122 123 124 125
(defmethod two-arg-* ((a qd-complex) (b qd-complex))
  (let* ((rx (realpart a))
	 (ix (imagpart a))
	 (ry (realpart b))
	 (iy (imagpart b)))
    (complex (- (* rx ry) (* ix iy))
	     (+ (* rx iy) (* ix ry)))))
126 127 128 129 130 131 132

(defmethod two-arg-* ((a qd-complex) (b real))
  (let* ((rx (realpart a))
	 (ix (imagpart a)))
    (complex (* rx b)
	     (* ix b))))

133 134 135 136 137 138 139 140 141 142 143
(defmethod two-arg-* ((a qd-complex) (b qd-real))
  (let* ((rx (realpart a))
	 (ix (imagpart a)))
    (complex (* rx b)
	     (* ix b))))

(defmethod two-arg-* ((a qd-real) (b qd-complex))
  (two-arg-* b a))



144
#+cmu
145
(defmethod two-arg-* ((a qd-complex) (b cl:complex))
146
  ;; For now, convert B into a qd-complex and use that.
147 148 149
  (let ((re (coerce (realpart b) 'ext:double-double-float))
	(im (coerce (imagpart b) 'ext:double-double-float)))
    (two-arg-* a (make-instance 'qd-complex
150 151
				:real (make-qd-dd re 0w0)
				:imag (make-qd-dd im 0w0)))))
152

153 154 155 156 157 158
#-cmu
(defmethod two-arg-* ((a qd-complex) (b cl:complex))
  ;; For now, convert B into a qd-complex and use that.
  (let ((re (coerce (realpart b) 'double-float))
	(im (coerce (imagpart b) 'double-float)))
    (two-arg-* a (make-instance 'qd-complex
159 160
				:real (make-qd-d re)
				:imag (make-qd-d im)))))
161

162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
(defmethod two-arg-* ((a number) (b qd-complex))
  (two-arg-* b a))

(defmethod two-arg-/ ((x qd-complex) (y qd-complex))
  (let* ((rx (realpart x))
	 (ix (imagpart x))
	 (ry (realpart y))
	 (iy (imagpart y)))
    (if (> (abs ry) (abs iy))
	(let* ((r (/ iy ry))
	       (dn (+ ry (* r iy))))
	  (complex (/ (+ rx (* ix r)) dn)
		   (/ (- ix (* rx r)) dn)))
	(let* ((r (/ ry iy))
	       (dn (+ iy (* r ry))))
	  (complex (/ (+ (* rx r) ix) dn)
		   (/ (- (* ix r) rx) dn))))))

(defmethod two-arg-/ ((x qd-complex) (y qd-real))
  (complex (/ (realpart x) y)
	   (/ (imagpart x) y)))

184 185 186 187
(defmethod two-arg-/ ((x qd-complex) (y number))
  (complex (/ (realpart x) y)
	   (/ (imagpart x) y)))

188
(defmethod two-arg-/ ((x number) (y qd-complex))
toy's avatar
toy committed
189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
  (let* ((rx (realpart x))
	 (ix (imagpart x))
	 (ry (realpart y))
	 (iy (imagpart y)))
    (if (> (abs ry) (abs iy))
	(let* ((r (/ iy ry))
	       (dn (+ ry (* r iy))))
	  (complex (/ (+ rx (* ix r)) dn)
		   (/ (- ix (* rx r)) dn)))
	(let* ((r (/ ry iy))
	       (dn (+ iy (* r ry))))
	  (complex (/ (+ (* rx r) ix) dn)
		   (/ (- (* ix r) rx) dn))))))

(defmethod two-arg-/ ((x qd-real) (y qd-complex))
  ;; This can be simplified since X is real.
205 206 207 208 209 210 211 212 213 214 215 216 217
  (let* ((rx (realpart x))
	 (ix (imagpart x))
	 (ry (realpart y))
	 (iy (imagpart y)))
    (if (> (abs ry) (abs iy))
	(let* ((r (/ iy ry))
	       (dn (+ ry (* r iy))))
	  (complex (/ (+ rx (* ix r)) dn)
		   (/ (- ix (* rx r)) dn)))
	(let* ((r (/ ry iy))
	       (dn (+ iy (* r ry))))
	  (complex (/ (+ (* rx r) ix) dn)
		   (/ (- (* ix r) rx) dn))))))
218

219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
(defmethod two-arg-+ ((a qd-real) (b cl:complex))
  (complex (+ a (realpart b))
	   (imagpart b)))

(defmethod two-arg-+ ((a cl:complex) (b qd-real))
  (complex (+ (realpart a) b)
	   (imagpart a)))

(defmethod two-arg-- ((a qd-real) (b cl:complex))
  (complex (- a (realpart b))
	   (- (imagpart b))))

(defmethod two-arg-- ((a cl:complex) (b qd-real))
  (complex (- (realpart a) b)
	   (imagpart a)))

(defmethod two-arg-* ((a qd-real) (b cl:complex))
  (complex (* a (realpart b))
	   (* a (imagpart b))))

(defmethod two-arg-* ((a cl:complex) (b qd-real))
  (complex (* (realpart a) b)
	   (* (imagpart a) b)))


(defmethod two-arg-/ ((a qd-real) (b cl:complex))
  (two-arg-/ a (coerce b 'qd-complex)))

(defmethod two-arg-/ ((a cl:complex) (b qd-real))
  (complex (/ (realpart a) b)
	   (/ (imagpart a) b)))


252 253 254
(defmethod unary-divide ((a qd-complex))
  (two-arg-/ #q1 a))

Raymond Toy's avatar
Raymond Toy committed
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
(defmethod coerce ((obj t) (type t))
  (cl:coerce obj type))

(defmethod coerce ((number cl:real) (type (eql 'qd-real)))
  (float number #q0))

(defmethod coerce ((number qd-real) (type (eql 'qd-real)))
  number)

(defmethod coerce ((number cl:number) (type (eql 'qd-complex)))
  (complex (float (realpart number) #q0)
	   (float (imagpart number) #q0)))

(defmethod coerce ((number qd-complex) (type (eql 'qd-complex)))
  number)

Raymond Toy's avatar
Raymond Toy committed
271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297
;; These two macros are borrowed from CMUCL.
(defmacro incf (place &optional (delta 1) &environment env)
  "The first argument is some location holding a number. This number is
  incremented by the second argument, DELTA, which defaults to 1."
  (multiple-value-bind (dummies vals newval setter getter)
      (get-setf-expansion place env)
    (let ((d (gensym)))
      `(let* (,@(mapcar #'list dummies vals)
              (,d ,delta)
              (,(car newval) (+ ,getter ,d)))
         ,setter))))

(defmacro decf (place &optional (delta 1) &environment env)
  "The first argument is some location holding a number. This number is
  decremented by the second argument, DELTA, which defaults to 1."
  (multiple-value-bind (dummies vals newval setter getter)
      (get-setf-expansion place env)
    (let ((d (gensym)))
      `(let* (,@(mapcar #'list dummies vals)
              (,d ,delta)
              (,(car newval) (- ,getter ,d)))
         ,setter))))


;; Most of this code taken from CMUCL and slightly modified to support
;; QD-COMPLEX.

toy's avatar
toy committed
298 299 300
(declaim (inline square))
(defun square (x)
  (declare (type qd-real x))
301
  (make-instance 'qd-real :value (sqr-qd (qd-value x))))
toy's avatar
toy committed
302 303 304 305 306 307 308 309 310 311

(defun qd-complex-sqrt (z)
  "Principle square root of Z

Z may be any number, but the result is always a complex."
  (declare (type qd-complex z))
  (multiple-value-bind (rho k)
      (qd-cssqs z)
    (declare (type qd-real rho)
	     (type fixnum k))
312 313
    (let ((x (realpart z))
	  (y (imagpart z))
toy's avatar
toy committed
314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352
	  (eta #q0.0)
	  (nu #q0.0))
      (declare (type qd-real x y eta nu))

      (locally
	  ;; space 0 to get maybe-inline functions inlined.
	  (declare (optimize (speed 3) (space 0)))

	(setf rho (+ (scalb (abs x) (- k)) (sqrt rho)))

	(cond ((oddp k)
	       (setf k (ash k -1)))
	      (t
	       (setf k (1- (ash k -1)))
	       (setf rho (+ rho rho))))

	(setf rho (scalb (sqrt rho) k))

	(setf eta rho)
	(setf nu y)

	(when (not (zerop rho))
	  (setf nu (/ (/ nu rho) 2d0)))
	(when (minusp x)
	  (setf eta (abs nu))
	  (setf nu (float-sign y rho))))
      (complex eta nu))))

(defun qd-complex-log-scaled (z j)
  "Compute log(2^j*z).

This is for use with J /= 0 only when |z| is huge."
  (declare (type qd-complex z)
	   (fixnum j))
  ;; The constants t0, t1, t2 should be evaluated to machine
  ;; precision.  In addition, Kahan says the accuracy of log1p
  ;; influences the choices of these constants but doesn't say how to
  ;; choose them.  We'll just assume his choices matches our
  ;; implementation of log1p.
353
  (let ((t0 (/ 1 (sqrt #q2.0q0)))
toy's avatar
toy committed
354 355
	(t1 #q1.2q0)
	(t2 #q3q0)
Raymond Toy's avatar
Raymond Toy committed
356
	(ln2 +log2+)
357 358
	(x (realpart z))
	(y (imagpart z)))
toy's avatar
toy committed
359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388
    (multiple-value-bind (rho k)
	(qd-cssqs z)
      (declare (optimize (speed 3)))
      (let ((beta (max (abs x) (abs y)))
	    (theta (min (abs x) (abs y))))
	(complex (if (and (zerop k)
			  (< t0 beta)
			  (or (<= beta t1)
			      (< rho t2)))
		     (/ (log1p (+ (* (- beta 1.0d0)
				     (+ beta 1.0d0))
				  (* theta theta)))
			2d0)
		     (+ (/ (log rho) 2d0)
			(* (+ k j) ln2)))
		 (atan y x))))))

(defun qd-complex-log (z)
  "Log of Z = log |Z| + i * arg Z

Z may be any number, but the result is always a complex."
  (declare (type qd-complex z))
  (qd-complex-log-scaled z 0))
	       

;; Let us note the following "strange" behavior.  atanh 1.0d0 is
;; +infinity, but the following code returns approx 176 + i*pi/4. The
;; reason for the imaginary part is caused by the fact that arg i*y is
;; never 0 since we have positive and negative zeroes.

389 390 391 392 393 394
;;
;; atanh(z) = (log(1+z) - log(1-z))/2
;;
;; The branch cut is on the real axis for |x| >= 1.  For x =< -1,
;; atanh is continuous with quadrant III; for x >= 1, continuous with
;; quadrant I.
toy's avatar
toy committed
395 396 397
(defun qd-complex-atanh (z)
  "Compute atanh z = (log(1+z) - log(1-z))/2"
  (declare (type (or qd-real qd-complex) z))
398 399
  (cond ((and (typep z 'qd-real) (< z -1))
	 (qd-complex-atanh (complex z -0d0)))
400
	(t
401 402 403 404 405 406 407 408 409 410 411 412
	 (flet ((careful-mul (a b)
		  ;; Carefully multiply a and b, taking care to handle
		  ;; signed zeroes.  Only need to handle the case of b
		  ;; being zero.
		  (if (zerop b)
		      (if (minusp (* (float-sign a) (float-sign b)))
			  #q-0
			  #q0)
		      (* a b))))
	   (let* ( ;; Constants
		  (theta (/ (sqrt most-positive-double-float) 4.0d0))
		  (rho (/ 4.0d0 (sqrt most-positive-double-float)))
Raymond Toy's avatar
Raymond Toy committed
413
		  (half-pi +pi/2+)
414 415
		  (rp (realpart z))
		  (beta (float-sign rp))
416
		  (x (* beta rp))
417
		  (y (careful-mul beta (- (imagpart z))))
418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462
		  (eta #q0.0q0)
		  (nu #q0.0q0))
	     ;; Shouldn't need this declare.
	     (declare (type qd-real x y))
	     (locally
		 (declare (optimize (speed 3)))
	       (cond ((or (> x theta)
			  (> (abs y) theta))
		      ;; To avoid overflow...
		      (setf nu (float-sign y half-pi))
		      ;; eta is real part of 1/(x + iy).  This is x/(x^2+y^2),
		      ;; which can cause overflow.  Arrange this computation so
		      ;; that it won't overflow.
		      (setf eta (let* ((x-bigger (> x (abs y)))
				       (r (if x-bigger (/ y x) (/ x y)))
				       (d (+ 1.0d0 (* r r))))
				  (if x-bigger
				      (/ (/ x) d)
				      (/ (/ r y) d)))))
		     ((= x #q1.0q0)
		      ;; Should this be changed so that if y is zero, eta is set
		      ;; to +infinity instead of approx 176?  In any case
		      ;; tanh(176) is 1.0d0 within working precision.
		      (let ((t1 (+ 4d0 (square y)))
			    (t2 (+ (abs y) rho)))
			(setf eta (log (/ (sqrt (sqrt t1))
					  (sqrt t2))))
			(setf nu (* 0.5d0
				    (float-sign y
						(+ half-pi (atan (* 0.5d0 t2))))))
			))
		     (t
		      (let ((t1 (+ (abs y) rho)))
			;; Normal case using log1p(x) = log(1 + x)
			(setf eta (* 0.25d0
				     (log1p (/ (* 4.0d0 x)
					       (+ (square (- 1.0d0 x))
						  (square t1))))))
			(setf nu (* 0.5d0
				    (atan (careful-mul 2.0d0 y)
					  (- (* (- 1.0d0 x)
						(+ 1.0d0 x))
					     (square t1))))))))
	       (complex (* beta eta)
			(- (* beta nu)))))))))
toy's avatar
toy committed
463 464


465 466
;; tanh(z) = sinh(z)/cosh(z)
;;
467 468 469
(defun qd-complex-tanh (z)
  "Compute tanh z = sinh z / cosh z"
  (declare (type (or qd-real qd-complex) z))
470 471
  (let ((x (realpart z))
	(y (imagpart z)))
472 473 474
    (locally
	;; space 0 to get maybe-inline functions inlined
	(declare (optimize (speed 3) (space 0)))
475 476 477 478 479 480 481 482 483
      (cond ((> (abs x) #.(/ (+ (log most-positive-double-float)
				(log 2d0))
			     4d0))
	     ;; The threshold above is
	     ;; asinh(most-positive-double-float)/4, but many Lisps
	     ;; cannot actually compute that.  Hence use the
	     ;; (accurate) approximation
	     ;; asinh(most-positive-double-float) =
	     ;; log(most-positive-double-float) + log(2)
484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504
	     (complex (float-sign x)
		      (float-sign y)))
	    (t
	     ;; With quad-double's it happens that tan(pi/2) will
	     ;; actually produce a division by zero error.  We need to
	     ;; handle that case carefully.
	     (let* ((tv (ignore-errors (tan y)))
		    (s (sinh x))
		    (rho (sqrt (+ 1.0d0 (* s s)))))
	       (cond (tv
		      (let* ((beta (+ 1.0d0 (* tv tv)))
			     (den (+ 1.0d0 (* beta s s))))
			(complex (/ (* beta rho s)
				    den)
				 (/ tv den))))
		     (t
		      ;; This means tan(y) produced some error.  We'll
		      ;; assume it's an overflow error because y is
		      ;; pi/2 + 2*k*pi.  But we need a value for tv to
		      ;; compute (/ tv).  This would be a signed-zero.
		      ;; For now, just return +0.
505 506
		      (complex (/ rho s)
			       #q0)))))))))
507

toy's avatar
toy committed
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549
;; Kahan says we should only compute the parts needed.  Thus, the
;; realpart's below should only compute the real part, not the whole
;; complex expression.  Doing this can be important because we may get
;; spurious signals that occur in the part that we are not using.
;;
;; However, we take a pragmatic approach and just use the whole
;; expression.

;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
;; it's the conjugate of the square root or the square root of the
;; conjugate.  This needs to be checked.

;; I checked.  It doesn't matter because (conjugate (sqrt z)) is the
;; same as (sqrt (conjugate z)) for all z.  This follows because
;;
;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
;;
;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
;;
;; and these two expressions are equal if and only if arg conj z =
;; -arg z, which is clearly true for all z.

;; NOTE: The rules of Common Lisp says that if you mix a real with a
;; complex, the real is converted to a complex before performing the
;; operation.  However, Kahan says in this paper (pg 176):
;;
;; (iii) Careless handling can turn infinity or the sign of zero into
;;       misinformation that subsequently disappears leaving behind
;;       only a plausible but incorrect result.  That is why compilers
;;       must not transform z-1 into z-(1+i*0), as we have seen above,
;;       nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
;;       subsequent logarithm or square root produce a non-zero
;;       imaginary part whose sign is opposite to what was intended.
;;
;; The interesting examples are too long and complicated to reproduce
;; here.  We refer the reader to his paper.
;;
;; The functions below are intended to handle the cases where a real
;; is mixed with a complex and we don't want CL complex contagion to
;; occur..

(declaim (inline 1+z 1-z z-1 z+1))
550
(defun 1+z (z)
toy's avatar
toy committed
551
  (complex (+ 1 (realpart z)) (imagpart z)))
552
(defun 1-z (z)
toy's avatar
toy committed
553
  (complex (- 1 (realpart z)) (- (imagpart z))))
554
(defun z-1 (z)
toy's avatar
toy committed
555
  (complex (- (realpart z) 1) (imagpart z)))
556
(defun z+1 (z)
toy's avatar
toy committed
557 558 559 560 561 562
  (complex (+ (realpart z) 1) (imagpart z)))

(defun qd-complex-acos (z)
  "Compute acos z = pi/2 - asin z

Z may be any number, but the result is always a complex."
563 564
  (declare (type (or qd-real qd-complex) z))
  (if (and (typep z 'qd-real) (> z 1))
toy's avatar
toy committed
565
      ;; acos is continuous in quadrant IV in this case.
566 567 568
      (qd-complex-acos (complex z -0f0))
      (let ((sqrt-1+z (qd-complex-sqrt (1+z z)))
	    (sqrt-1-z (qd-complex-sqrt (1-z z))))
toy's avatar
toy committed
569 570 571 572 573 574
	(cond ((zerop (realpart sqrt-1+z))
	       ;; Same as below, but we compute atan ourselves (because we
	       ;; have atan +/- infinity).
	       (complex 
			(if (minusp (float-sign (* (realpart sqrt-1-z)
						   (realpart sqrt-1+z))))
575 576
			    (- +pi+)
			    +pi+)
toy's avatar
toy committed
577 578 579 580 581 582 583 584 585 586 587 588
			(asinh (imagpart (* (conjugate sqrt-1+z)
					    sqrt-1-z)))))
	      (t
	       (complex (* 2 (atan (/ (realpart sqrt-1-z)
				      (realpart sqrt-1+z))))
			(asinh (imagpart (* (conjugate sqrt-1+z)
					    sqrt-1-z)))))))))

(defun qd-complex-acosh (z)
  "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))

Z may be any number, but the result is always a complex."
589 590 591
  (declare (type (or qd-real qd-complex) z))
  (let* ((sqrt-z-1 (qd-complex-sqrt (z-1 z)))
	 (sqrt-z+1 (qd-complex-sqrt (z+1 z))))
toy's avatar
toy committed
592 593 594 595 596 597 598 599 600 601
    ;; We need to handle the case where real part of sqrt-z+1 is zero,
    ;; because division by zero with double-double-floats doesn't
    ;; produce infinity.
    (cond ((zerop (realpart sqrt-z+1))
	   ;; Same as below, but we compute atan ourselves (because we
	   ;; have atan +/- infinity).
	   (complex (asinh (realpart (* (conjugate sqrt-z-1)
					sqrt-z+1)))
		    (if (minusp (float-sign (* (imagpart sqrt-z-1)
					       (realpart sqrt-z+1))))
602 603
			(- +pi+)
			+pi+)))
toy's avatar
toy committed
604 605 606 607 608 609
	  (t
	   (complex (asinh (realpart (* (conjugate sqrt-z-1)
					sqrt-z+1)))
		    (* 2 (atan (/ (imagpart sqrt-z-1)
				  (realpart sqrt-z+1)))))))))

toy's avatar
toy committed
610 611
;; asin(z) = asinh(i*z)/i
;;         = -i log(i*z + sqrt(1-z^2))
612 613 614 615 616
(defun qd-complex-asin (z)
  "Compute asin z = asinh(i*z)/i

Z may be any number, but the result is always a complex."
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
617 618
  (cond ((and (typep z 'qd-real) (> z 1))
	 (qd-complex-asin (complex z -0d0)))
619 620 621 622 623 624
	(t
	 (let* ((sqrt-1-z (qd-complex-sqrt (1-z z)))
		(sqrt-1+z (qd-complex-sqrt (1+z z)))
		(den (realpart (* sqrt-1-z sqrt-1+z))))
	   (cond ((zerop den)
		  ;; Like below but we handle atan part ourselves.
toy's avatar
toy committed
625 626 627 628
		  ;; Must be sure to take into account the sign of
		  ;; (realpart z) and den!
		  (complex (if (minusp (* (float-sign (realpart z))
					  (float-sign den)))
629 630 631 632 633
			       (- (/ +pi+ 2))
			       (/ +pi+ 2))
			   (asinh (imagpart (* (conjugate sqrt-1-z)
					       sqrt-1+z)))))
		 (t
634 635 636 637 638
		  ;; We get a invalid operation here when z is real and |z| > 1.
		  (complex (atan (/ (realpart z)
				    (realpart (* sqrt-1-z sqrt-1+z))))
			   (asinh (imagpart (* (conjugate sqrt-1-z)
					       sqrt-1+z))))))))))
639

toy's avatar
toy committed
640 641 642 643
(defun qd-complex-asinh (z)
  "Compute asinh z = log(z + sqrt(1 + z*z))

Z may be any number, but the result is always a complex."
644
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
645 646
  ;; asinh z = -i * asin (i*z)
  (let* ((iz (complex (- (imagpart z)) (realpart z)))
647
	 (result (qd-complex-asin iz)))
toy's avatar
toy committed
648 649 650 651 652 653 654
    (complex (imagpart result)
	     (- (realpart result)))))
	 
(defun qd-complex-atan (z)
  "Compute atan z = atanh (i*z) / i

Z may be any number, but the result is always a complex."
655
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
656 657
  ;; atan z = -i * atanh (i*z)
  (let* ((iz (complex (- (imagpart z)) (realpart z)))
658
	 (result (qd-complex-atanh iz)))
toy's avatar
toy committed
659 660 661 662 663 664 665
    (complex (imagpart result)
	     (- (realpart result)))))

(defun qd-complex-tan (z)
  "Compute tan z = -i * tanh(i * z)

Z may be any number, but the result is always a complex."
666
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
667 668
  ;; tan z = -i * tanh(i*z)
  (let* ((iz (complex (- (imagpart z)) (realpart z)))
669
	 (result (qd-complex-tanh iz)))
toy's avatar
toy committed
670 671
    (complex (imagpart result)
	     (- (realpart result)))))
Raymond Toy's avatar
Raymond Toy committed
672
;; End of implementation of complex functions from CMUCL.
673

674
(defmethod asin ((x qd-complex))
675 676
  (qd-complex-asin x))

677
(defmethod acos ((x qd-complex))
678 679
  (qd-complex-acos x))

680
(defmethod acosh ((x qd-complex))
681 682
  (qd-complex-acosh x))

683
(defmethod atanh ((x qd-complex))
684
  (qd-complex-atanh x))
685

686
(defmethod sin ((z qd-complex))
687 688 689 690 691
  (let ((x (realpart z))
	(y (imagpart z)))
    (complex (* (sin x) (cosh y))
	     (* (cos x) (sinh y)))))

692
(defmethod cos ((z qd-complex))
693 694 695 696 697
  (let ((x (realpart z))
	(y (imagpart z)))
    (complex (* (cos x) (cosh y))
	     (- (* (sin x) (sinh y))))))

698
(defmethod tan ((z qd-complex))
699 700
  (qd-complex-tan z))

701
(defmethod sinh ((z qd-complex))
702 703 704 705 706
  (let ((x (realpart z))
	(y (imagpart z)))
    (complex (* (sinh x) (cos y))
	     (* (cosh x) (sin y)))))

707
(defmethod cosh ((z qd-complex))
toy's avatar
toy committed
708 709
  (let ((x (realpart z))
	(y (imagpart z)))
710 711 712
    (complex (* (cosh x) (cos y))
	     (* (sinh x) (sin y)))))

713
(defmethod tanh ((z qd-complex))
714
  (qd-complex-tanh z))
715

716
(defmethod sqrt ((z qd-complex))
717
  (qd-complex-sqrt z))
toy's avatar
toy committed
718

719
(defmethod atan ((y qd-complex) &optional x)
720 721 722 723
  (if x
      (error "First arg of 2-arg ATAN must be real")
      (qd-complex-atan y)))

724
(defmethod atan ((y cl:complex) &optional x)
725 726 727
  (if x
      (error "First arg of 2-arg ATAN must be real")
      (cl:atan y)))
728

729
(defmethod exp ((z qd-complex))
730 731 732 733 734 735
  (let* ((x (realpart z))
	 (y (imagpart z))
	 (ex (exp x)))
    (complex (* ex (cos y))
	     (* ex (sin y)))))

736
(defmethod log ((a qd-complex) &optional b)
737
  (if b
Raymond Toy's avatar
Raymond Toy committed
738
      (/ (log a) (log b))
739 740 741 742 743 744
      (complex (log (abs a))
	       (atan (imagpart a) (realpart a)))))

(defmethod qexpt ((x qd-complex) (y number))
  (exp (* (float y #q0) (log x))))

745 746 747
(defmethod qexpt ((x qd-complex) (y qd-real))
  (exp (* y (log x))))

748 749 750
(defmethod qexpt ((x number) (y qd-complex))
  (exp (* y (log (float x #q0)))))

toy's avatar
toy committed
751 752 753
(defmethod qexpt ((x qd-complex) (y qd-complex))
  (exp (* y (log x))))

754
(defmethod phase ((z qd-complex))
toy's avatar
toy committed
755
  (atan (imagpart z) (realpart z)))
Raymond Toy's avatar
Raymond Toy committed
756 757 758 759 760 761 762 763 764 765 766 767

(defun realp (x)
  (or (typep x 'qd-real)
      (cl:realp x)))

(defun complexp (x)
  (or (typep x 'qd-complex)
      (cl:complexp x)))

(defun numberp (x)
  (or (realp x)
      (complexp x)))