/[cmucl]/src/code/irrat.lisp
ViewVC logotype

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.45.2.1 - (hide annotations)
Fri Jun 9 16:04:57 2006 UTC (7 years, 10 months ago) by rtoy
Branch: double-double-branch
CVS Tags: double-double-init-sparc-2, double-double-init-sparc, double-double-init-ppc, double-double-init-%make-sparc, double-double-init-checkpoint-1, double-double-reader-base, double-double-init-x86
Branch point for: double-double-reader-branch
Changes since 1.45: +49 -5 lines
Add basic support for kernel:double-double-float type.  The primitive
type is there, and basic arithmetic operations work as well as PRINT.
But many things do not work: No reader, formatted output, many mixed
type arithmetic operations, special functions are just double-float
values, coerced to double-double-float.

compiler/generic/interr.lisp:
o Add new error

compiler/generic/new-genesis.lisp:
o Dump double-double-float objects (barely tested)

compiler/generic/primtype.lisp:
o Tell compiler about the new primitive type double-double-float.

compiler/generic/vm-fndb.lisp:
o Make double-double-float-p a known function.

compiler/generic/vm-type.lisp:
o Update FLOAT-FORMAT-NAME to include double-double-float

compiler/generic/vm-typetran.lisp:
o Tell compiler about double-double-float type predicate.

compiler/sparc/float.lisp:
o Add necessary vops to move double-double-float args, store and load
  double-double-floats to/from the double-double-stack,
  double-double-reg moves, box and unbox double-double-floats, move
  double-double-floats to and from args
o Add necessary vops to create a double-double-float and to extract
  the high and low parts out of a double-double-float.

compiler/sparc/parms.lisp:
o Define double-double-float-digits

compiler/sparc/type-vops.lisp:
o Define double-double-float type vop
o Adjust number hierarchy to include double-double-float

compiler/sparc/vm.lisp:
o Define the necessary storage class and storage base for the
  double-double-reg and double-double-stack.

lisp/gencgc.c:
o Tell GC about double-double-float objects.

lisp/purify.c:
o Tell purify about double-double-float objects.

code/class.lisp:
o Add the new double-double-float class.

code/exports.lisp:
o Add the necessary symbols to the various packages.  (This is
  important to get right otherwise there's confusion on what symbol
  really represents double-double-float stuff.)

code/float.lisp:
o Implement some of the necessary functions to support
  double-double-float.

code/hash-new.lisp:
o Hash double-double-floats by xor'ing the hashes of each double-float
  part.  (Is that good enough?)

code/irrat.lisp:
o Implement the special functions by calling the double-float versions
  and coercing the result to a double-double-float.  This is needed to
  get type-derivation working, but the precise value isn't that
  important right now.  We'll have to implement them later.

code/kernel.lisp:
o Make make-double-double-float, double-double-hi, and
  double-double-lo known functions.

code/lispinit.lisp:
o Register the :double-double float feature.

code/load.lisp:
o Add FOP for reading double-double-float values from fasls.  (Barely
  tested, if at all.)

code/numbers.lisp:
o Implement basic arithmetic operations for double-double-floats.
  This needs quite a bit of work to clean up, but most things work.

code/pred.lisp:
o Tell the type system about double-double-float type.

code/print.lisp:
o Add very rudimentary printing for double-double-float.  Basically
  copied from code written by Richard Fateman, with permission.

code/seq.lisp:
o Tell coerce how to coerce things to a double-double-float.

code/type.lisp:
o Tell type system about the new float format double-double-float and
  how numeric contagion works with double-double-float.

code/dump.lisp:
o Tell dumper how to dump double-double-float values to a fasl.

compiler/float-tran.lisp:
o Add appropriate deftransforms to handle conversion of things to
  double-double-float and from from double-double-float to other float
  types.
o The basic implmentation of double-double-float arithmetic is also
  here.
o Add deftransforms to tell the compiler how to do basic arithmetic
  and comparisions on double-double-float numbers.

compiler/srctran.lisp:
o Fix a bug in interval-range-info when :low is 0dd0 and :high is 0.
  We weren't returning the right value, which confuses interval-div,
  which caused an infinite loop because we couldn't split the
  numerator into positive and negative parts.
o Tell other parts about double-double-float.
1 wlott 1.1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2     ;;;
3     ;;; **********************************************************************
4 ram 1.7 ;;; This code was written as part of the CMU Common Lisp project at
5     ;;; Carnegie Mellon University, and has been placed in the public domain.
6     ;;;
7     (ext:file-comment
8 rtoy 1.45.2.1 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.45.2.1 2006/06/09 16:04:57 rtoy Exp $")
9 ram 1.7 ;;;
10 wlott 1.1 ;;; **********************************************************************
11     ;;;
12     ;;; This file contains all the irrational functions. Actually, most of the
13     ;;; work is done by calling out to C...
14     ;;;
15     ;;; Author: William Lott.
16     ;;;
17    
18     (in-package "KERNEL")
19    
20    
21     ;;;; Random constants, utility functions, and macros.
22    
23     (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
24 wlott 1.2 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
25 wlott 1.1
26 ram 1.5 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
27     ;;; call, and saves number consing to boot.
28     ;;;
29 wlott 1.1 (defmacro def-math-rtn (name num-args)
30     (let ((function (intern (concatenate 'simple-string
31     "%"
32     (string-upcase name)))))
33 ram 1.4 `(progn
34 pw 1.31 (declaim (inline ,function))
35 ram 1.4 (export ',function)
36 wlott 1.10 (alien:def-alien-routine (,name ,function) double-float
37 ram 1.4 ,@(let ((results nil))
38     (dotimes (i num-args (nreverse results))
39     (push (list (intern (format nil "ARG-~D" i))
40     'double-float)
41     results)))))))
42 wlott 1.1
43     (eval-when (compile load eval)
44    
45     (defun handle-reals (function var)
46     `((((foreach fixnum single-float bignum ratio))
47     (coerce (,function (coerce ,var 'double-float)) 'single-float))
48     ((double-float)
49 rtoy 1.45.2.1 (,function ,var))
50     #+double-double
51     ((double-double-float)
52     ;; A hack until we write double-double-float versions of these
53     ;; special functions.
54     (kernel:make-double-double-float (,function (kernel:double-double-hi ,var))
55     0d0))))
56 wlott 1.1
57     ); eval-when (compile load eval)
58    
59    
60     ;;;; Stubs for the Unix math library.
61    
62     ;;; Please refer to the Unix man pages for details about these routines.
63    
64     ;;; Trigonometric.
65 ram 1.17 #-x86 (def-math-rtn "sin" 1)
66     #-x86 (def-math-rtn "cos" 1)
67     #-x86 (def-math-rtn "tan" 1)
68 wlott 1.1 (def-math-rtn "asin" 1)
69     (def-math-rtn "acos" 1)
70 ram 1.17 #-x86 (def-math-rtn "atan" 1)
71     #-x86 (def-math-rtn "atan2" 2)
72 wlott 1.1 (def-math-rtn "sinh" 1)
73     (def-math-rtn "cosh" 1)
74     (def-math-rtn "tanh" 1)
75 pw 1.18 (def-math-rtn "asinh" 1)
76     (def-math-rtn "acosh" 1)
77     (def-math-rtn "atanh" 1)
78 wlott 1.1
79     ;;; Exponential and Logarithmic.
80 dtc 1.19 #-x86 (def-math-rtn "exp" 1)
81     #-x86 (def-math-rtn "log" 1)
82     #-x86 (def-math-rtn "log10" 1)
83 wlott 1.1 (def-math-rtn "pow" 2)
84 dtc 1.30 #-(or x86 sparc-v7 sparc-v8 sparc-v9) (def-math-rtn "sqrt" 1)
85 wlott 1.1 (def-math-rtn "hypot" 2)
86 dtc 1.27 #-(or hpux x86) (def-math-rtn "log1p" 1)
87 ram 1.17
88     #+x86 ;; These are needed for use by byte-compiled files.
89     (progn
90     (defun %sin (x)
91 dtc 1.19 (declare (double-float x)
92 ram 1.17 (values double-float))
93     (%sin x))
94     (defun %sin-quick (x)
95     (declare (double-float x)
96     (values double-float))
97     (%sin-quick x))
98     (defun %cos (x)
99 dtc 1.19 (declare (double-float x)
100 ram 1.17 (values double-float))
101     (%cos x))
102     (defun %cos-quick (x)
103     (declare (double-float x)
104     (values double-float))
105     (%cos-quick x))
106     (defun %tan (x)
107     (declare (double-float x)
108     (values double-float))
109     (%tan x))
110     (defun %tan-quick (x)
111     (declare (double-float x)
112     (values double-float))
113     (%tan-quick x))
114     (defun %atan (x)
115     (declare (double-float x)
116     (values double-float))
117     (%atan x))
118     (defun %atan2 (x y)
119     (declare (double-float x y)
120     (values double-float))
121     (%atan2 x y))
122     (defun %exp (x)
123 dtc 1.19 (declare (double-float x)
124 ram 1.17 (values double-float))
125     (%exp x))
126     (defun %log (x)
127 dtc 1.19 (declare (double-float x)
128 ram 1.17 (values double-float))
129     (%log x))
130     (defun %log10 (x)
131 dtc 1.19 (declare (double-float x)
132 ram 1.17 (values double-float))
133     (%log10 x))
134     #+nil ;; notyet
135     (defun %pow (x y)
136     (declare (type (double-float 0d0) x)
137 dtc 1.19 (double-float y)
138 ram 1.17 (values (double-float 0d0)))
139     (%pow x y))
140     (defun %sqrt (x)
141 dtc 1.19 (declare (double-float x)
142 ram 1.17 (values double-float))
143     (%sqrt x))
144     (defun %scalbn (f ex)
145 dtc 1.19 (declare (double-float f)
146 ram 1.17 (type (signed-byte 32) ex)
147     (values double-float))
148     (%scalbn f ex))
149     (defun %scalb (f ex)
150     (declare (double-float f ex)
151     (values double-float))
152     (%scalb f ex))
153     (defun %logb (x)
154 dtc 1.19 (declare (double-float x)
155     (values double-float))
156     (%logb x))
157 dtc 1.27 (defun %log1p (x)
158     (declare (double-float x)
159     (values double-float))
160     (%log1p x))
161 ram 1.17 ) ; progn
162 wlott 1.1
163 rtoy 1.44
164     ;; As above for x86. It also seems to be needed to handle
165     ;; constant-folding in the compiler.
166     #+sparc
167     (progn
168     (defun %sqrt (x)
169     (declare (double-float x)
170     (values double-float))
171     (%sqrt x))
172     )
173    
174 wlott 1.1
175     ;;;; Power functions.
176    
177     (defun exp (number)
178     "Return e raised to the power NUMBER."
179     (number-dispatch ((number number))
180     (handle-reals %exp number)
181     ((complex)
182     (* (exp (realpart number))
183     (cis (imagpart number))))))
184    
185     ;;; INTEXP -- Handle the rational base, integer power case.
186    
187     (defparameter *intexp-maximum-exponent* 10000)
188    
189 ram 1.6 ;;; This function precisely calculates base raised to an integral power. It
190     ;;; separates the cases by the sign of power, for efficiency reasons, as powers
191     ;;; can be calculated more efficiently if power is a positive integer. Values
192     ;;; of power are calculated as positive integers, and inverted if negative.
193     ;;;
194 wlott 1.1 (defun intexp (base power)
195 rtoy 1.45 ;; Handle the special case of 1^power. Maxima sometimes does this,
196     ;; and there's no need to cause a continuable error in this case.
197     ;; Should we also handle (-1)^power?
198     (when (eql base 1)
199     (return-from intexp base))
200    
201 wlott 1.1 (when (> (abs power) *intexp-maximum-exponent*)
202     (cerror "Continue with calculation."
203     "The absolute value of ~S exceeds ~S."
204     power '*intexp-maximum-exponent* base power))
205     (cond ((minusp power)
206     (/ (intexp base (- power))))
207     ((eql base 2)
208     (ash 1 power))
209     (t
210     (do ((nextn (ash power -1) (ash power -1))
211     (total (if (oddp power) base 1)
212     (if (oddp power) (* base total) total)))
213     ((zerop nextn) total)
214     (setq base (* base base))
215     (setq power nextn)))))
216    
217    
218 ram 1.6 ;;; EXPT -- Public
219     ;;;
220     ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
221     ;;; floating point stuff. If both args are real, we try %POW right off,
222     ;;; assuming it will return 0 if the result may be complex. If so, we call
223     ;;; COMPLEX-POW which directly computes the complex result. We also separate
224     ;;; the complex-real and real-complex cases from the general complex case.
225     ;;;
226 wlott 1.1 (defun expt (base power)
227 wlott 1.3 "Returns BASE raised to the POWER."
228 wlott 1.1 (if (zerop power)
229 rtoy 1.40 ;; CLHS says that if the power is 0, the result is 1, subject to
230     ;; numeric contagion. But what happens if base is infinity or
231     ;; NaN? Do we silently return 1? For now, I think we should
232     ;; signal an error if the FP modes say so.
233     (let ((result (1+ (* base power))))
234     ;; If we get an NaN here, that means base*power above didn't
235     ;; produce 0 and FP traps were disabled, so we handle that
236     ;; here. Should this be a continuable restart?
237     (if (and (floatp result) (float-nan-p result))
238     (float 1 result)
239     result))
240 dtc 1.22 (labels (;; determine if the double float is an integer.
241     ;; 0 - not an integer
242     ;; 1 - an odd int
243     ;; 2 - an even int
244     (isint (ihi lo)
245     (declare (type (unsigned-byte 31) ihi)
246     (type (unsigned-byte 32) lo)
247     (optimize (speed 3) (safety 0)))
248     (let ((isint 0))
249     (declare (type fixnum isint))
250     (cond ((>= ihi #x43400000) ; exponent >= 53
251     (setq isint 2))
252     ((>= ihi #x3ff00000)
253     (let ((k (- (ash ihi -20) #x3ff))) ; exponent
254     (declare (type (mod 53) k))
255     (cond ((> k 20)
256     (let* ((shift (- 52 k))
257     (j (logand (ash lo (- shift))))
258     (j2 (ash j shift)))
259     (declare (type (mod 32) shift)
260     (type (unsigned-byte 32) j j2))
261     (when (= j2 lo)
262     (setq isint (- 2 (logand j 1))))))
263     ((= lo 0)
264     (let* ((shift (- 20 k))
265     (j (ash ihi (- shift)))
266     (j2 (ash j shift)))
267     (declare (type (mod 32) shift)
268     (type (unsigned-byte 31) j j2))
269     (when (= j2 ihi)
270     (setq isint (- 2 (logand j 1))))))))))
271     isint))
272     (real-expt (x y rtype)
273     (let ((x (coerce x 'double-float))
274     (y (coerce y 'double-float)))
275     (declare (double-float x y))
276     (let* ((x-hi (kernel:double-float-high-bits x))
277     (x-lo (kernel:double-float-low-bits x))
278     (x-ihi (logand x-hi #x7fffffff))
279     (y-hi (kernel:double-float-high-bits y))
280     (y-lo (kernel:double-float-low-bits y))
281     (y-ihi (logand y-hi #x7fffffff)))
282     (declare (type (signed-byte 32) x-hi y-hi)
283     (type (unsigned-byte 31) x-ihi y-ihi)
284     (type (unsigned-byte 32) x-lo y-lo))
285     ;; y==zero: x**0 = 1
286     (when (zerop (logior y-ihi y-lo))
287     (return-from real-expt (coerce 1d0 rtype)))
288     ;; +-NaN return x+y
289     (when (or (> x-ihi #x7ff00000)
290     (and (= x-ihi #x7ff00000) (/= x-lo 0))
291     (> y-ihi #x7ff00000)
292     (and (= y-ihi #x7ff00000) (/= y-lo 0)))
293     (return-from real-expt (coerce (+ x y) rtype)))
294     (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
295     (declare (type fixnum yisint))
296     ;; special value of y
297     (when (and (zerop y-lo) (= y-ihi #x7ff00000))
298     ;; y is +-inf
299     (return-from real-expt
300     (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
301     ;; +-1**inf is NaN
302     (coerce (- y y) rtype))
303     ((>= x-ihi #x3ff00000)
304     ;; (|x|>1)**+-inf = inf,0
305     (if (>= y-hi 0)
306     (coerce y rtype)
307     (coerce 0 rtype)))
308     (t
309     ;; (|x|<1)**-,+inf = inf,0
310     (if (< y-hi 0)
311     (coerce (- y) rtype)
312     (coerce 0 rtype))))))
313    
314     (let ((abs-x (abs x)))
315     (declare (double-float abs-x))
316     ;; special value of x
317     (when (and (zerop x-lo)
318     (or (= x-ihi #x7ff00000) (zerop x-ihi)
319     (= x-ihi #x3ff00000)))
320     ;; x is +-0,+-inf,+-1
321     (let ((z (if (< y-hi 0)
322     (/ 1 abs-x) ; z = (1/|x|)
323     abs-x)))
324     (declare (double-float z))
325     (when (< x-hi 0)
326     (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
327     ;; (-1)**non-int
328     (let ((y*pi (* y pi)))
329     (declare (double-float y*pi))
330     (return-from real-expt
331 dtc 1.24 (complex
332 dtc 1.22 (coerce (%cos y*pi) rtype)
333     (coerce (%sin y*pi) rtype)))))
334     ((= yisint 1)
335     ;; (x<0)**odd = -(|x|**odd)
336     (setq z (- z)))))
337     (return-from real-expt (coerce z rtype))))
338    
339     (if (>= x-hi 0)
340     ;; x>0
341     (coerce (kernel::%pow x y) rtype)
342     ;; x<0
343     (let ((pow (kernel::%pow abs-x y)))
344     (declare (double-float pow))
345     (case yisint
346     (1 ; Odd
347     (coerce (* -1d0 pow) rtype))
348     (2 ; Even
349     (coerce pow rtype))
350     (t ; Non-integer
351     (let ((y*pi (* y pi)))
352     (declare (double-float y*pi))
353 dtc 1.24 (complex
354 dtc 1.22 (coerce (* pow (%cos y*pi)) rtype)
355     (coerce (* pow (%sin y*pi)) rtype)))))))))))))
356     (declare (inline real-expt))
357     (number-dispatch ((base number) (power number))
358     (((foreach fixnum (or bignum ratio) (complex rational)) integer)
359     (intexp base power))
360     (((foreach single-float double-float) rational)
361     (real-expt base power '(dispatch-type base)))
362     (((foreach fixnum (or bignum ratio) single-float)
363     (foreach ratio single-float))
364     (real-expt base power 'single-float))
365     (((foreach fixnum (or bignum ratio) single-float double-float)
366     double-float)
367     (real-expt base power 'double-float))
368     ((double-float single-float)
369     (real-expt base power 'double-float))
370     (((foreach (complex rational) (complex float)) rational)
371     (* (expt (abs base) power)
372     (cis (* power (phase base)))))
373     (((foreach fixnum (or bignum ratio) single-float double-float)
374     complex)
375     (if (and (zerop base) (plusp (realpart power)))
376     (* base power)
377     (exp (* power (log base)))))
378     (((foreach (complex float) (complex rational))
379     (foreach complex double-float single-float))
380     (if (and (zerop base) (plusp (realpart power)))
381     (* base power)
382     (exp (* power (log base)))))))))
383 wlott 1.1
384 toy 1.34 ;; Compute the base 2 log of an integer
385     (defun log2 (x)
386     ;; Write x = 2^n*f where 1/2 < f <= 1. Then log2(x) = n + log2(f).
387     ;;
388     ;; So we grab the top few bits of x and scale that appropriately,
389     ;; take the log of it and add it to n.
390     (let ((n (integer-length x)))
391     (if (< n vm:double-float-digits)
392     (log (coerce x 'double-float) 2d0)
393     (let ((exp (min vm:double-float-digits n))
394     (f (ldb (byte vm:double-float-digits
395     (max 0 (- n vm:double-float-digits)))
396     x)))
397     (+ n (log (scale-float (float f 1d0) (- exp))
398     2d0))))))
399 rtoy 1.43
400 wlott 1.1 (defun log (number &optional (base nil base-p))
401 wlott 1.3 "Return the logarithm of NUMBER in the base BASE, which defaults to e."
402 wlott 1.1 (if base-p
403 toy 1.34 (cond ((zerop base)
404     ;; ANSI spec
405     base)
406     ((and (integerp number) (integerp base)
407     (plusp number) (plusp base))
408     ;; Let's try to do something nice when both the number
409     ;; and the base are positive integers. Use the rule that
410     ;; log_b(x) = log_2(x)/log_2(b)
411     (coerce (/ (log2 number) (log2 base)) 'single-float))
412 rtoy 1.43 ((and (realp number) (realp base))
413     ;; CLHS 12.1.4.1 says
414     ;;
415     ;; When rationals and floats are combined by a
416     ;; numerical function, the rational is first converted
417     ;; to a float of the same format.
418     ;;
419     ;; So assume this applies to floats as well convert all
420     ;; numbers to the largest float format before computing
421     ;; the log.
422     ;;
423     ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
424     (number-dispatch ((number real) (base real))
425     ((double-float
426     (foreach double-float single-float fixnum bignum ratio))
427     (/ (log number) (log (coerce base 'double-float))))
428     (((foreach single-float fixnum bignum ratio)
429     double-float)
430     (/ (log (coerce number 'double-float)) (log base)))
431     (((foreach single-float fixnum bignum ratio)
432     (foreach single-float fixnum bignum ratio))
433     ;; Converting everything to double-float helps the
434     ;; cases like (log 17 10) = (/ (log 17) (log 10)).
435     ;; This is usually handled above, but if we compute (/
436     ;; (log 17) (log 10)), we get a slightly different
437     ;; answer due to roundoff. This makes it a bit more
438     ;; consistent.
439     ;;
440     ;; FIXME: This probably needs more work.
441     (let ((result (/ (log (float number 1d0))
442     (log (float base 1d0)))))
443     (if (realp result)
444     (coerce result 'single-float)
445     (coerce result '(complex single-float)))))))
446 toy 1.34 (t
447 rtoy 1.43 ;; FIXME: This probably needs some work as well.
448 toy 1.34 (/ (log number) (log base))))
449 wlott 1.1 (number-dispatch ((number number))
450 toy 1.34 (((foreach fixnum bignum))
451 wlott 1.3 (if (minusp number)
452 toy 1.34 (complex (coerce (log (- number)) 'single-float)
453     (coerce pi 'single-float))
454 toy 1.37 (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
455 toy 1.34 ((ratio)
456     (if (minusp number)
457     (complex (coerce (log (- number)) 'single-float)
458     (coerce pi 'single-float))
459 toy 1.35 ;; What happens when the ratio is close to 1? We need to
460     ;; be careful to preserve accuracy.
461     (let ((top (numerator number))
462     (bot (denominator number)))
463     ;; If the number of bits in the numerator and
464     ;; denominator are different, just use the fact
465 toy 1.37 ;; log(x/y) = log(x) - log(y). But to preserve
466     ;; accuracy, we actually do
467     ;; (log2(x)-log2(y))/log2(e)).
468     ;;
469     ;; However, if the numerator and denominator have the
470     ;; same number of bits, implying the quotient is near
471     ;; one, we use log1p(x) = log(1+x). Since the number is
472     ;; rational, we don't lose precision subtracting 1 from
473     ;; it, and converting it to double-float is accurate.
474 toy 1.35 (if (= (integer-length top)
475 toy 1.36 (integer-length bot))
476 toy 1.35 (coerce (%log1p (coerce (- number 1) 'double-float))
477     'single-float)
478 toy 1.37 (coerce (/ (- (log2 top) (log2 bot))
479     #.(log (exp 1d0) 2d0))
480     'single-float)))))
481 ram 1.17 (((foreach single-float double-float))
482     ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
483     ;; Since this doesn't seem to be an implementation issue
484     ;; I (pw) take the Kahan result.
485     (if (< (float-sign number)
486     (coerce 0 '(dispatch-type number)))
487     (complex (log (- number)) (coerce pi '(dispatch-type number)))
488     (coerce (%log (coerce number 'double-float))
489     '(dispatch-type number))))
490 rtoy 1.45.2.1 #+double-double
491     ((double-double-float)
492     ;; Hack!
493     (let ((hi (kernel:double-double-hi number)))
494     (if (< (float-sign hi)
495     (coerce 0 '(dispatch-type number)))
496     (complex (coerce (log (- hi)) 'kernel:double-double-float)
497     (coerce pi '(dispatch-type number)))
498     (coerce (%log hi) '(dispatch-type number)))))
499 ram 1.17 ((complex)
500     (complex-log number)))))
501 wlott 1.1
502     (defun sqrt (number)
503     "Return the square root of NUMBER."
504     (number-dispatch ((number number))
505 ram 1.17 (((foreach fixnum bignum ratio))
506 wlott 1.3 (if (minusp number)
507 ram 1.17 (complex-sqrt number)
508 wlott 1.3 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
509 ram 1.17 (((foreach single-float double-float))
510 dtc 1.28 (if (minusp number)
511 ram 1.17 (complex-sqrt number)
512     (coerce (%sqrt (coerce number 'double-float))
513     '(dispatch-type number))))
514 rtoy 1.45.2.1 #+double-double
515     ((double-double-float)
516     (multiple-value-bind (hi lo)
517     (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
518     (kernel:make-double-double-float hi lo)))
519     ((complex)
520     (complex-sqrt number))))
521 wlott 1.1
522    
523     ;;;; Trigonometic and Related Functions
524    
525 wlott 1.2 (defun abs (number)
526     "Returns the absolute value of the number."
527     (number-dispatch ((number number))
528     (((foreach single-float double-float fixnum rational))
529     (abs number))
530 rtoy 1.45.2.1 #+(and nil double-double)
531     ((double-double-float)
532     ;; This is a hack until abs deftransform is working
533     (multiple-value-bind (hi lo)
534     (c::abs-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
535     (kernel:make-double-double-float hi lo)))
536     #+double-double
537     ((double-double-float)
538     ;; This is a hack until abs deftransform is working
539     (let ((hi (kernel:double-double-hi number))
540     (lo (kernel:double-double-lo number)))
541     (declare (double-float hi lo))
542     (when (minusp hi)
543     (setf hi (- hi))
544     (setf lo (- lo)))
545     (kernel:make-double-double-float hi lo)))
546 wlott 1.2 ((complex)
547     (let ((rx (realpart number))
548     (ix (imagpart number)))
549     (etypecase rx
550     (rational
551     (sqrt (+ (* rx rx) (* ix ix))))
552     (single-float
553     (coerce (%hypot (coerce rx 'double-float)
554     (coerce ix 'double-float))
555     'single-float))
556     (double-float
557 rtoy 1.45.2.1 (%hypot rx ix))
558     #+double-double
559     (double-double-float
560     (error "abs complex double-double-float not implemented!")))))))
561 wlott 1.1
562     (defun phase (number)
563     "Returns the angle part of the polar representation of a complex number.
564 wlott 1.3 For complex numbers, this is (atan (imagpart number) (realpart number)).
565 wlott 1.1 For non-complex positive numbers, this is 0. For non-complex negative
566     numbers this is PI."
567 wlott 1.3 (etypecase number
568 ram 1.17 (rational
569 wlott 1.3 (if (minusp number)
570     (coerce pi 'single-float)
571     0.0f0))
572 ram 1.17 (single-float
573     (if (minusp (float-sign number))
574     (coerce pi 'single-float)
575     0.0f0))
576 wlott 1.3 (double-float
577 ram 1.17 (if (minusp (float-sign number))
578 wlott 1.3 (coerce pi 'double-float)
579     0.0d0))
580 rtoy 1.45.2.1 #+double-double
581     (double-double-float
582     (if (minusp (float-sign number))
583     (coerce pi 'double-double-float)
584     (kernel:make-double-double-float 0d0 0d0)))
585 wlott 1.3 (complex
586     (atan (imagpart number) (realpart number)))))
587 wlott 1.1
588    
589     (defun sin (number)
590 wlott 1.3 "Return the sine of NUMBER."
591 wlott 1.1 (number-dispatch ((number number))
592     (handle-reals %sin number)
593     ((complex)
594     (let ((x (realpart number))
595     (y (imagpart number)))
596 ram 1.17 (complex (* (sin x) (cosh y))
597     (* (cos x) (sinh y)))))))
598 wlott 1.1
599     (defun cos (number)
600 wlott 1.3 "Return the cosine of NUMBER."
601 wlott 1.1 (number-dispatch ((number number))
602     (handle-reals %cos number)
603     ((complex)
604     (let ((x (realpart number))
605     (y (imagpart number)))
606 ram 1.17 (complex (* (cos x) (cosh y))
607     (- (* (sin x) (sinh y))))))))
608 wlott 1.1
609     (defun tan (number)
610 wlott 1.3 "Return the tangent of NUMBER."
611 wlott 1.1 (number-dispatch ((number number))
612     (handle-reals %tan number)
613     ((complex)
614 ram 1.17 (complex-tan number))))
615 wlott 1.1
616 wlott 1.2 (defun cis (theta)
617 wlott 1.3 "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
618 wlott 1.2 (if (complexp theta)
619     (error "Argument to CIS is complex: ~S" theta)
620     (complex (cos theta) (sin theta))))
621 wlott 1.1
622     (defun asin (number)
623 wlott 1.3 "Return the arc sine of NUMBER."
624 wlott 1.1 (number-dispatch ((number number))
625 wlott 1.3 ((rational)
626     (if (or (> number 1) (< number -1))
627     (complex-asin number)
628     (coerce (%asin (coerce number 'double-float)) 'single-float)))
629     (((foreach single-float double-float))
630 rtoy 1.42 (if (or (float-nan-p number)
631     (and (<= number (coerce 1 '(dispatch-type number)))
632     (>= number (coerce -1 '(dispatch-type number)))))
633 wlott 1.3 (coerce (%asin (coerce number 'double-float))
634 rtoy 1.42 '(dispatch-type number))
635     (complex-asin number)))
636 wlott 1.1 ((complex)
637 wlott 1.3 (complex-asin number))))
638 wlott 1.1
639     (defun acos (number)
640 wlott 1.3 "Return the arc cosine of NUMBER."
641 wlott 1.1 (number-dispatch ((number number))
642 wlott 1.3 ((rational)
643     (if (or (> number 1) (< number -1))
644     (complex-acos number)
645     (coerce (%acos (coerce number 'double-float)) 'single-float)))
646     (((foreach single-float double-float))
647 rtoy 1.42 (if (or (float-nan-p number)
648     (and (<= number (coerce 1 '(dispatch-type number)))
649     (>= number (coerce -1 '(dispatch-type number)))))
650 wlott 1.3 (coerce (%acos (coerce number 'double-float))
651 rtoy 1.42 '(dispatch-type number))
652     (complex-acos number)))
653 wlott 1.1 ((complex)
654 wlott 1.3 (complex-acos number))))
655 wlott 1.1
656    
657     (defun atan (y &optional (x nil xp))
658 wlott 1.3 "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
659 wlott 1.1 (if xp
660 wlott 1.12 (flet ((atan2 (y x)
661 wlott 1.13 (declare (type double-float y x)
662     (values double-float))
663     (if (zerop x)
664     (if (zerop y)
665     (if (plusp (float-sign x))
666     y
667     (float-sign y pi))
668     (float-sign y (/ pi 2)))
669     (%atan2 y x))))
670 toy 1.33 ;; If X is given, both X and Y must be real numbers.
671     (number-dispatch ((y real) (x real))
672 wlott 1.13 ((double-float
673     (foreach double-float single-float fixnum bignum ratio))
674     (atan2 y (coerce x 'double-float)))
675     (((foreach single-float fixnum bignum ratio)
676     double-float)
677     (atan2 (coerce y 'double-float) x))
678     (((foreach single-float fixnum bignum ratio)
679     (foreach single-float fixnum bignum ratio))
680     (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
681     'single-float))))
682 wlott 1.1 (number-dispatch ((y number))
683     (handle-reals %atan y)
684     ((complex)
685 ram 1.17 (complex-atan y)))))
686 wlott 1.3
687 ram 1.17 (defun sinh (number)
688     "Return the hyperbolic sine of NUMBER."
689     (number-dispatch ((number number))
690     (handle-reals %sinh number)
691     ((complex)
692     (let ((x (realpart number))
693     (y (imagpart number)))
694     (complex (* (sinh x) (cos y))
695     (* (cosh x) (sin y)))))))
696    
697     (defun cosh (number)
698     "Return the hyperbolic cosine of NUMBER."
699     (number-dispatch ((number number))
700     (handle-reals %cosh number)
701     ((complex)
702     (let ((x (realpart number))
703     (y (imagpart number)))
704     (complex (* (cosh x) (cos y))
705     (* (sinh x) (sin y)))))))
706    
707 wlott 1.3 (defun tanh (number)
708     "Return the hyperbolic tangent of NUMBER."
709 ram 1.17 (number-dispatch ((number number))
710     (handle-reals %tanh number)
711     ((complex)
712     (complex-tanh number))))
713 wlott 1.3
714     (defun asinh (number)
715     "Return the hyperbolic arc sine of NUMBER."
716 ram 1.17 (number-dispatch ((number number))
717     (handle-reals %asinh number)
718     ((complex)
719     (complex-asinh number))))
720 wlott 1.3
721     (defun acosh (number)
722     "Return the hyperbolic arc cosine of NUMBER."
723 ram 1.17 (number-dispatch ((number number))
724     ((rational)
725     ;; acosh is complex if number < 1
726     (if (< number 1)
727     (complex-acosh number)
728     (coerce (%acosh (coerce number 'double-float)) 'single-float)))
729     (((foreach single-float double-float))
730     (if (< number (coerce 1 '(dispatch-type number)))
731     (complex-acosh number)
732     (coerce (%acosh (coerce number 'double-float))
733     '(dispatch-type number))))
734     ((complex)
735     (complex-acosh number))))
736 wlott 1.3
737     (defun atanh (number)
738     "Return the hyperbolic arc tangent of NUMBER."
739 ram 1.17 (number-dispatch ((number number))
740     ((rational)
741     ;; atanh is complex if |number| > 1
742     (if (or (> number 1) (< number -1))
743     (complex-atanh number)
744     (coerce (%atanh (coerce number 'double-float)) 'single-float)))
745     (((foreach single-float double-float))
746     (if (or (> number (coerce 1 '(dispatch-type number)))
747     (< number (coerce -1 '(dispatch-type number))))
748     (complex-atanh number)
749     (coerce (%atanh (coerce number 'double-float))
750     '(dispatch-type number))))
751     ((complex)
752     (complex-atanh number))))
753 wlott 1.14
754 toy 1.32 ;;; HP-UX does not supply a C version of log1p, so use the definition.
755     ;;; We really need to fix this. The definition really loses big-time
756     ;;; in roundoff as x gets small.
757 wlott 1.14
758     #+hpux
759 pw 1.18 (declaim (inline %log1p))
760 wlott 1.14 #+hpux
761 pw 1.18 (defun %log1p (number)
762     (declare (double-float number)
763     (optimize (speed 3) (safety 0)))
764     (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
765 ram 1.17
766    
767     ;;;;
768     ;;;; This is a set of routines that implement many elementary
769     ;;;; transcendental functions as specified by ANSI Common Lisp. The
770     ;;;; implementation is based on Kahan's paper.
771     ;;;;
772     ;;;; I believe I have accurately implemented the routines and are
773     ;;;; correct, but you may want to check for your self.
774     ;;;;
775     ;;;; These functions are written for CMU Lisp and take advantage of
776     ;;;; some of the features available there. It may be possible,
777     ;;;; however, to port this to other Lisps.
778     ;;;;
779     ;;;; Some functions are significantly more accurate than the original
780     ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
781     ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
782     ;;;; answer is pi + i*log(2-sqrt(3)).
783     ;;;;
784     ;;;; All of the implemented functions will take any number for an
785     ;;;; input, but the result will always be a either a complex
786     ;;;; single-float or a complex double-float.
787     ;;;;
788     ;;;; General functions
789     ;;;; complex-sqrt
790     ;;;; complex-log
791     ;;;; complex-atanh
792     ;;;; complex-tanh
793     ;;;; complex-acos
794     ;;;; complex-acosh
795     ;;;; complex-asin
796     ;;;; complex-asinh
797     ;;;; complex-atan
798     ;;;; complex-tan
799     ;;;;
800     ;;;; Utility functions:
801     ;;;; scalb logb
802     ;;;;
803     ;;;; Internal functions:
804     ;;;; square coerce-to-complex-type cssqs complex-log-scaled
805     ;;;;
806     ;;;;
807     ;;;; Please send any bug reports, comments, or improvements to Raymond
808     ;;;; Toy at toy@rtp.ericsson.se.
809     ;;;;
810     ;;;; References
811     ;;;;
812     ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
813     ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
814     ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
815     ;;;; Press, 1987
816     ;;;;
817 toy 1.32
818 ram 1.17 (declaim (inline square))
819     (defun square (x)
820 toy 1.32 (declare (double-float x))
821 ram 1.17 (* x x))
822    
823     ;; If you have these functions in libm, perhaps they should be used
824     ;; instead of these Lisp versions. These versions are probably good
825     ;; enough, especially since they are portable.
826    
827 dtc 1.19 (declaim (inline scalb))
828 ram 1.17 (defun scalb (x n)
829     "Compute 2^N * X without compute 2^N first (use properties of the
830     underlying floating-point format"
831     (declare (type double-float x)
832 dtc 1.19 (type double-float-exponent n))
833 ram 1.17 (scale-float x n))
834    
835 toy 1.32 (declaim (inline logb-finite))
836     (defun logb-finite (x)
837     "Same as logb but X is not infinity and non-zero and not a NaN, so
838     that we can always return an integer"
839     (declare (type double-float x))
840     (multiple-value-bind (signif expon sign)
841     (decode-float x)
842     (declare (ignore signif sign))
843     ;; decode-float is almost right, except that the exponent
844     ;; is off by one
845     (1- expon)))
846    
847 ram 1.17 (defun logb (x)
848 toy 1.32 "Compute an integer N such that 1 <= |2^(-N) * x| < 2.
849 ram 1.17 For the special cases, the following values are used:
850    
851     x logb
852     NaN NaN
853     +/- infinity +infinity
854     0 -infinity
855     "
856     (declare (type double-float x))
857 dtc 1.26 (cond ((float-nan-p x)
858 ram 1.17 x)
859     ((float-infinity-p x)
860     #.ext:double-float-positive-infinity)
861     ((zerop x)
862     ;; The answer is negative infinity, but we are supposed to
863 toy 1.32 ;; signal divide-by-zero, so do the actual division
864 ram 1.17 (/ -1.0d0 x)
865     )
866     (t
867 toy 1.32 (logb-finite x))))
868    
869    
870 ram 1.17
871     ;; This function is used to create a complex number of the appropriate
872     ;; type.
873    
874     (declaim (inline coerce-to-complex-type))
875     (defun coerce-to-complex-type (x y z)
876     "Create complex number with real part X and imaginary part Y such that
877     it has the same type as Z. If Z has type (complex rational), the X
878     and Y are coerced to single-float."
879     (declare (double-float x y)
880 toy 1.32 (number z)
881     (optimize (extensions:inhibit-warnings 3)))
882     (if (typep (realpart z) 'double-float)
883 ram 1.17 (complex x y)
884     ;; Convert anything that's not a double-float to a single-float.
885 toy 1.32 (complex (float x 1f0)
886     (float y 1f0))))
887 ram 1.17
888     (defun cssqs (z)
889 dtc 1.21 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
890     ;; result is r + i*k, where k is an integer.
891 ram 1.17
892     ;; Save all FP flags
893 dtc 1.21 (let ((x (float (realpart z) 1d0))
894 toy 1.32 (y (float (imagpart z) 1d0)))
895 dtc 1.21 ;; Would this be better handled using an exception handler to
896     ;; catch the overflow or underflow signal? For now, we turn all
897     ;; traps off and look at the accrued exceptions to see if any
898     ;; signal would have been raised.
899     (with-float-traps-masked (:underflow :overflow)
900 toy 1.32 (let ((rho (+ (square x) (square y))))
901     (declare (optimize (speed 3) (space 0)))
902     (cond ((and (or (float-nan-p rho)
903     (float-infinity-p rho))
904     (or (float-infinity-p (abs x))
905     (float-infinity-p (abs y))))
906     (values ext:double-float-positive-infinity 0))
907     ((let ((threshold #.(/ least-positive-double-float
908     double-float-epsilon))
909     (traps (ldb vm::float-sticky-bits
910     (vm:floating-point-modes))))
911     ;; Overflow raised or (underflow raised and rho <
912     ;; lambda/eps)
913     (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
914     (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
915     (< rho threshold))))
916     ;; If we're here, neither x nor y are infinity and at
917     ;; least one is non-zero.. Thus logb returns a nice
918     ;; integer.
919     (let ((k (- (logb-finite (max (abs x) (abs y))))))
920     (values (+ (square (scalb x k))
921     (square (scalb y k)))
922     (- k))))
923     (t
924     (values rho 0)))))))
925 ram 1.17
926     (defun complex-sqrt (z)
927     "Principle square root of Z
928    
929     Z may be any number, but the result is always a complex."
930     (declare (number z))
931     (multiple-value-bind (rho k)
932     (cssqs z)
933 toy 1.32 (declare (type (or (member 0d0) (double-float 0d0)) rho)
934     (type fixnum k))
935 ram 1.17 (let ((x (float (realpart z) 1.0d0))
936     (y (float (imagpart z) 1.0d0))
937     (eta 0d0)
938     (nu 0d0))
939     (declare (double-float x y eta nu))
940    
941 toy 1.32 (locally
942     ;; space 0 to get maybe-inline functions inlined.
943     (declare (optimize (speed 3) (space 0)))
944    
945     (if (not (float-nan-p x))
946     (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
947    
948     (cond ((oddp k)
949     (setf k (ash k -1)))
950     (t
951     (setf k (1- (ash k -1)))
952     (setf rho (+ rho rho))))
953    
954     (setf rho (scalb (sqrt rho) k))
955    
956     (setf eta rho)
957     (setf nu y)
958    
959     (when (/= rho 0d0)
960     (when (not (float-infinity-p (abs nu)))
961     (setf nu (/ (/ nu rho) 2d0)))
962     (when (< x 0d0)
963     (setf eta (abs nu))
964     (setf nu (float-sign y rho))))
965     (coerce-to-complex-type eta nu z)))))
966 ram 1.17
967     (defun complex-log-scaled (z j)
968     "Compute log(2^j*z).
969    
970     This is for use with J /= 0 only when |z| is huge."
971     (declare (number z)
972     (fixnum j))
973     ;; The constants t0, t1, t2 should be evaluated to machine
974     ;; precision. In addition, Kahan says the accuracy of log1p
975     ;; influences the choices of these constants but doesn't say how to
976     ;; choose them. We'll just assume his choices matches our
977     ;; implementation of log1p.
978     (let ((t0 #.(/ 1 (sqrt 2.0d0)))
979     (t1 1.2d0)
980     (t2 3d0)
981     (ln2 #.(log 2d0))
982     (x (float (realpart z) 1.0d0))
983     (y (float (imagpart z) 1.0d0)))
984     (multiple-value-bind (rho k)
985     (cssqs z)
986 toy 1.32 (declare (optimize (speed 3)))
987 ram 1.17 (let ((beta (max (abs x) (abs y)))
988     (theta (min (abs x) (abs y))))
989 toy 1.32 (coerce-to-complex-type (if (and (zerop k)
990     (< t0 beta)
991     (or (<= beta t1)
992     (< rho t2)))
993     (/ (%log1p (+ (* (- beta 1.0d0)
994     (+ beta 1.0d0))
995     (* theta theta)))
996     2d0)
997     (+ (/ (log rho) 2d0)
998     (* (+ k j) ln2)))
999     (atan y x)
1000     z)))))
1001 ram 1.17
1002     (defun complex-log (z)
1003     "Log of Z = log |Z| + i * arg Z
1004    
1005     Z may be any number, but the result is always a complex."
1006     (declare (number z))
1007     (complex-log-scaled z 0))
1008    
1009     ;; Let us note the following "strange" behavior. atanh 1.0d0 is
1010     ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1011     ;; reason for the imaginary part is caused by the fact that arg i*y is
1012     ;; never 0 since we have positive and negative zeroes.
1013    
1014     (defun complex-atanh (z)
1015     "Compute atanh z = (log(1+z) - log(1-z))/2"
1016     (declare (number z))
1017 rtoy 1.41 (if (and (realp z) (< z -1))
1018     ;; atanh is continuous in quadrant III in this case.
1019     (complex-atanh (complex z -0f0))
1020     (let* ( ;; Constants
1021     (theta (/ (sqrt most-positive-double-float) 4.0d0))
1022     (rho (/ 4.0d0 (sqrt most-positive-double-float)))
1023     (half-pi (/ pi 2.0d0))
1024     (rp (float (realpart z) 1.0d0))
1025     (beta (float-sign rp 1.0d0))
1026     (x (* beta rp))
1027     (y (* beta (- (float (imagpart z) 1.0d0))))
1028     (eta 0.0d0)
1029     (nu 0.0d0))
1030     ;; Shouldn't need this declare.
1031     (declare (double-float x y))
1032     (locally
1033     (declare (optimize (speed 3)))
1034     (cond ((or (> x theta)
1035     (> (abs y) theta))
1036     ;; To avoid overflow...
1037     (setf nu (float-sign y half-pi))
1038     ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
1039     ;; which can cause overflow. Arrange this computation so
1040     ;; that it won't overflow.
1041     (setf eta (let* ((x-bigger (> x (abs y)))
1042     (r (if x-bigger (/ y x) (/ x y)))
1043     (d (+ 1.0d0 (* r r))))
1044     (if x-bigger
1045     (/ (/ x) d)
1046     (/ (/ r y) d)))))
1047     ((= x 1.0d0)
1048     ;; Should this be changed so that if y is zero, eta is set
1049     ;; to +infinity instead of approx 176? In any case
1050     ;; tanh(176) is 1.0d0 within working precision.
1051     (let ((t1 (+ 4d0 (square y)))
1052     (t2 (+ (abs y) rho)))
1053     (setf eta (log (/ (sqrt (sqrt t1))
1054     (sqrt t2))))
1055     (setf nu (* 0.5d0
1056     (float-sign y
1057     (+ half-pi (atan (* 0.5d0 t2))))))))
1058     (t
1059     (let ((t1 (+ (abs y) rho)))
1060     ;; Normal case using log1p(x) = log(1 + x)
1061     (setf eta (* 0.25d0
1062     (%log1p (/ (* 4.0d0 x)
1063     (+ (square (- 1.0d0 x))
1064     (square t1))))))
1065     (setf nu (* 0.5d0
1066     (atan (* 2.0d0 y)
1067     (- (* (- 1.0d0 x)
1068     (+ 1.0d0 x))
1069     (square t1))))))))
1070     (coerce-to-complex-type (* beta eta)
1071     (- (* beta nu))
1072     z)))))
1073 ram 1.17
1074     (defun complex-tanh (z)
1075     "Compute tanh z = sinh z / cosh z"
1076     (declare (number z))
1077     (let ((x (float (realpart z) 1.0d0))
1078     (y (float (imagpart z) 1.0d0)))
1079 toy 1.32 (locally
1080     ;; space 0 to get maybe-inline functions inlined
1081     (declare (optimize (speed 3) (space 0)))
1082     (cond ((> (abs x)
1083     #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1084     ;; This is more accurate under linux.
1085     #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1086     (%log most-positive-double-float)) 4d0))
1087     (coerce-to-complex-type (float-sign x)
1088     (float-sign y) z))
1089     (t
1090     (let* ((tv (%tan y))
1091     (beta (+ 1.0d0 (* tv tv)))
1092     (s (sinh x))
1093     (rho (sqrt (+ 1.0d0 (* s s)))))
1094     (if (float-infinity-p (abs tv))
1095     (coerce-to-complex-type (/ rho s)
1096     (/ tv)
1097     z)
1098     (let ((den (+ 1.0d0 (* beta s s))))
1099     (coerce-to-complex-type (/ (* beta rho s)
1100     den)
1101     (/ tv den)
1102     z)))))))))
1103 ram 1.17
1104     ;; Kahan says we should only compute the parts needed. Thus, the
1105     ;; realpart's below should only compute the real part, not the whole
1106     ;; complex expression. Doing this can be important because we may get
1107     ;; spurious signals that occur in the part that we are not using.
1108     ;;
1109     ;; However, we take a pragmatic approach and just use the whole
1110     ;; expression.
1111    
1112     ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1113     ;; it's the conjugate of the square root or the square root of the
1114     ;; conjugate. This needs to be checked.
1115    
1116     ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1117     ;; same as (sqrt (conjugate z)) for all z. This follows because
1118     ;;
1119     ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1120     ;;
1121     ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1122     ;;
1123 toy 1.35 ;; and these two expressions are equal if and only if arg conj z =
1124     ;; -arg z, which is clearly true for all z.
1125 ram 1.17
1126 rtoy 1.38 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1127     ;; complex, the real is converted to a complex before performing the
1128     ;; operation. However, Kahan says in this paper (pg 176):
1129     ;;
1130     ;; (iii) Careless handling can turn infinity or the sign of zero into
1131     ;; misinformation that subsequently disappears leaving behind
1132     ;; only a plausible but incorrect result. That is why compilers
1133     ;; must not transform z-1 into z-(1+i*0), as we have seen above,
1134     ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1135     ;; subsequent logarithm or square root produce a non-zero
1136     ;; imaginary part whose sign is opposite to what was intended.
1137     ;;
1138     ;; The interesting examples are too long and complicated to reproduce
1139     ;; here. We refer the reader to his paper.
1140     ;;
1141     ;; The functions below are intended to handle the cases where a real
1142     ;; is mixed with a complex and we don't want CL complex contagion to
1143     ;; occur..
1144    
1145     (declaim (inline 1+z 1-z z-1 z+1))
1146     (defun 1+z (z)
1147     (complex (+ 1 (realpart z)) (imagpart z)))
1148     (defun 1-z (z)
1149     (complex (- 1 (realpart z)) (- (imagpart z))))
1150     (defun z-1 (z)
1151     (complex (- (realpart z) 1) (imagpart z)))
1152     (defun z+1 (z)
1153     (complex (+ (realpart z) 1) (imagpart z)))
1154    
1155 ram 1.17 (defun complex-acos (z)
1156     "Compute acos z = pi/2 - asin z
1157    
1158     Z may be any number, but the result is always a complex."
1159     (declare (number z))
1160 rtoy 1.41 (if (and (realp z) (> z 1))
1161     ;; acos is continuous in quadrant IV in this case.
1162     (complex-acos (complex z -0f0))
1163     (let ((sqrt-1+z (complex-sqrt (1+z z)))
1164     (sqrt-1-z (complex-sqrt (1-z z))))
1165     (with-float-traps-masked (:divide-by-zero)
1166     (complex (* 2 (atan (/ (realpart sqrt-1-z)
1167     (realpart sqrt-1+z))))
1168     (asinh (imagpart (* (conjugate sqrt-1+z)
1169     sqrt-1-z))))))))
1170 ram 1.17
1171     (defun complex-acosh (z)
1172     "Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1173    
1174     Z may be any number, but the result is always a complex."
1175     (declare (number z))
1176 rtoy 1.38 (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1177     (sqrt-z+1 (complex-sqrt (z+1 z))))
1178 dtc 1.21 (with-float-traps-masked (:divide-by-zero)
1179     (complex (asinh (realpart (* (conjugate sqrt-z-1)
1180     sqrt-z+1)))
1181     (* 2 (atan (/ (imagpart sqrt-z-1)
1182     (realpart sqrt-z+1))))))))
1183 ram 1.17
1184    
1185     (defun complex-asin (z)
1186     "Compute asin z = asinh(i*z)/i
1187    
1188     Z may be any number, but the result is always a complex."
1189     (declare (number z))
1190 rtoy 1.41 (if (and (realp z) (> z 1))
1191     ;; asin is continuous in quadrant IV in this case.
1192     (complex-asin (complex z -0f0))
1193     (let ((sqrt-1-z (complex-sqrt (1-z z)))
1194     (sqrt-1+z (complex-sqrt (1+z z))))
1195     (with-float-traps-masked (:divide-by-zero)
1196     (complex (atan (/ (realpart z)
1197     (realpart (* sqrt-1-z sqrt-1+z))))
1198     (asinh (imagpart (* (conjugate sqrt-1-z)
1199     sqrt-1+z))))))))
1200 ram 1.17
1201     (defun complex-asinh (z)
1202     "Compute asinh z = log(z + sqrt(1 + z*z))
1203    
1204     Z may be any number, but the result is always a complex."
1205     (declare (number z))
1206     ;; asinh z = -i * asin (i*z)
1207     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1208     (result (complex-asin iz)))
1209     (complex (imagpart result)
1210     (- (realpart result)))))
1211    
1212     (defun complex-atan (z)
1213     "Compute atan z = atanh (i*z) / i
1214    
1215     Z may be any number, but the result is always a complex."
1216     (declare (number z))
1217     ;; atan z = -i * atanh (i*z)
1218     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1219     (result (complex-atanh iz)))
1220     (complex (imagpart result)
1221     (- (realpart result)))))
1222    
1223     (defun complex-tan (z)
1224     "Compute tan z = -i * tanh(i * z)
1225    
1226     Z may be any number, but the result is always a complex."
1227     (declare (number z))
1228     ;; tan z = -i * tanh(i*z)
1229     (let* ((iz (complex (- (imagpart z)) (realpart z)))
1230     (result (complex-tanh iz)))
1231     (complex (imagpart result)
1232     (- (realpart result)))))
1233 toy 1.32

  ViewVC Help
Powered by ViewVC 1.1.5