oct-complex.lisp 23.1 KB
Newer Older
Raymond Toy's avatar
Raymond Toy committed
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.

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

Raymond Toy's avatar
Raymond Toy committed
28
(named-readtables:in-readtable oct-readtable))
29

Raymond Toy's avatar
Raymond Toy committed
30
31
(defmethod add1 ((a qd-complex))
  (make-instance 'qd-complex
Raymond Toy's avatar
Raymond Toy committed
32
33
34
35
36
37
		 :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)
Raymond Toy's avatar
Raymond Toy committed
38
39
		 :imag (qd-value (imagpart a))))
		 
40
(defmethod two-arg-/ ((a qd-real) (b cl:rational))
toy's avatar
toy committed
41
42
43
  (make-instance 'qd-real :value (div-qd (qd-value a)
					 (qd-value (float b #q0)))))

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

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

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

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

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

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

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

toy's avatar
toy committed
67
68
69
70
(defmethod qzerop ((z qd-complex))
  (and (zerop (realpart z))
       (zerop (imagpart z))))

71
72
73
74
75
76
77
78
(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
79
80
81
82
83
84
85
86
(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)))

87
88
89
90
91
92
93
94
95
96
97
98
99
(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
100
	   (imagpart a)))
101
102

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

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

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

(defmethod two-arg-- ((a qd-real) (b qd-complex))
  (complex (- a (realpart b))
	   (- (imagpart b))))
117
  
118
119
120
121
122
123
124
(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)))))
125
126
127
128
129
130
131

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

132
133
134
135
136
137
138
139
140
141
142
(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))



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

152
153
154
155
156
157
#-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
toy's avatar
toy committed
158
159
				:real (make-qd-d re)
				:imag (make-qd-d im)))))
160

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
(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)))

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

187
(defmethod two-arg-/ ((x number) (y qd-complex))
toy's avatar
toy committed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  (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.
204
205
206
207
208
209
210
211
212
213
214
215
216
  (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))))))
217

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


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

Raymond Toy's avatar
Raymond Toy committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
(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
270
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
;; 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
297
298
299
(declaim (inline square))
(defun square (x)
  (declare (type qd-real x))
300
  (make-instance 'qd-real :value (sqr-qd (qd-value x))))
toy's avatar
toy committed
301
302
303
304
305
306
307
308
309
310

(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))
311
312
    (let ((x (realpart z))
	  (y (imagpart z))
toy's avatar
toy committed
313
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
	  (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.
352
  (let ((t0 (/ 1 (sqrt #q2.0q0)))
toy's avatar
toy committed
353
354
	(t1 #q1.2q0)
	(t2 #q3q0)
Raymond Toy's avatar
Raymond Toy committed
355
	(ln2 +log2+)
356
357
	(x (realpart z))
	(y (imagpart z)))
toy's avatar
toy committed
358
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
    (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.

toy's avatar
toy committed
388
389
390
391
392
393
;;
;; 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
394
395
396
(defun qd-complex-atanh (z)
  "Compute atanh z = (log(1+z) - log(1-z))/2"
  (declare (type (or qd-real qd-complex) z))
397
398
  (cond ((and (typep z 'qd-real) (< z -1))
	 (qd-complex-atanh (complex z -0d0)))
399
	(t
400
401
402
403
404
405
406
407
408
409
410
411
	 (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
412
		  (half-pi +pi/2+)
413
414
		  (rp (realpart z))
		  (beta (float-sign rp))
415
		  (x (* beta rp))
416
		  (y (careful-mul beta (- (imagpart z))))
417
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
		  (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
462
463


toy's avatar
toy committed
464
465
;; tanh(z) = sinh(z)/cosh(z)
;;
466
467
468
(defun qd-complex-tanh (z)
  "Compute tanh z = sinh z / cosh z"
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
469
470
  (let ((x (realpart z))
	(y (imagpart z)))
471
472
473
    (locally
	;; space 0 to get maybe-inline functions inlined
	(declare (optimize (speed 3) (space 0)))
474
475
476
477
478
479
480
481
482
      (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)
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
	     (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.
toy's avatar
toy committed
504
505
		      (complex (/ rho s)
			       #q0)))))))))
506

toy's avatar
toy committed
507
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
;; 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))
549
(defun 1+z (z)
toy's avatar
toy committed
550
  (complex (+ 1 (realpart z)) (imagpart z)))
551
(defun 1-z (z)
toy's avatar
toy committed
552
  (complex (- 1 (realpart z)) (- (imagpart z))))
553
(defun z-1 (z)
toy's avatar
toy committed
554
  (complex (- (realpart z) 1) (imagpart z)))
555
(defun z+1 (z)
toy's avatar
toy committed
556
557
558
559
560
561
  (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."
562
563
  (declare (type (or qd-real qd-complex) z))
  (if (and (typep z 'qd-real) (> z 1))
toy's avatar
toy committed
564
      ;; acos is continuous in quadrant IV in this case.
565
566
567
      (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
568
569
570
571
572
573
	(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))))
574
575
			    (- +pi+)
			    +pi+)
toy's avatar
toy committed
576
577
578
579
580
581
582
583
584
585
586
587
			(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."
588
589
590
  (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
591
592
593
594
595
596
597
598
599
600
    ;; 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))))
601
602
			(- +pi+)
			+pi+)))
toy's avatar
toy committed
603
604
605
606
607
608
	  (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
609
610
;; asin(z) = asinh(i*z)/i
;;         = -i log(i*z + sqrt(1-z^2))
611
612
613
614
615
(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
616
617
  (cond ((and (typep z 'qd-real) (> z 1))
	 (qd-complex-asin (complex z -0d0)))
618
619
620
621
622
623
	(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
624
625
626
627
		  ;; Must be sure to take into account the sign of
		  ;; (realpart z) and den!
		  (complex (if (minusp (* (float-sign (realpart z))
					  (float-sign den)))
628
629
630
631
632
			       (- (/ +pi+ 2))
			       (/ +pi+ 2))
			   (asinh (imagpart (* (conjugate sqrt-1-z)
					       sqrt-1+z)))))
		 (t
633
634
635
636
637
		  ;; 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))))))))))
638

toy's avatar
toy committed
639
640
641
642
(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."
643
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
644
645
  ;; asinh z = -i * asin (i*z)
  (let* ((iz (complex (- (imagpart z)) (realpart z)))
646
	 (result (qd-complex-asin iz)))
toy's avatar
toy committed
647
648
649
650
651
652
653
    (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."
654
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
655
656
  ;; atan z = -i * atanh (i*z)
  (let* ((iz (complex (- (imagpart z)) (realpart z)))
657
	 (result (qd-complex-atanh iz)))
toy's avatar
toy committed
658
659
660
661
662
663
664
    (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."
665
  (declare (type (or qd-real qd-complex) z))
toy's avatar
toy committed
666
667
  ;; tan z = -i * tanh(i*z)
  (let* ((iz (complex (- (imagpart z)) (realpart z)))
668
	 (result (qd-complex-tanh iz)))
toy's avatar
toy committed
669
670
    (complex (imagpart result)
	     (- (realpart result)))))
Raymond Toy's avatar
Raymond Toy committed
671
;; End of implementation of complex functions from CMUCL.
672

673
(defmethod asin ((x qd-complex))
674
675
  (qd-complex-asin x))

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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