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

Contents of /src/code/irrat.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.60.4.1 - (show annotations)
Thu Feb 25 20:34:50 2010 UTC (4 years, 1 month ago) by rtoy
Branch: intl-2-branch
Changes since 1.60: +44 -39 lines
Restart internalization work.  This new branch starts with code from
the intl-branch on date 2010-02-12 18:00:00+0500.  This version works
and

LANG=en@piglatin bin/lisp

works (once the piglatin translation is added).
1 ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
2 ;;;
3 ;;; **********************************************************************
4 ;;; 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 "$Header: /tiger/var/lib/cvsroots/cmucl/src/code/irrat.lisp,v 1.60.4.1 2010/02/25 20:34:50 rtoy Exp $")
9 ;;;
10 ;;; **********************************************************************
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 (intl:textdomain "cmucl")
20
21
22 ;;;; Random constants, utility functions, and macros.
23
24 (defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
25 ;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
26
27 ;;; Make these INLINE, since the call to C is at least as compact as a Lisp
28 ;;; call, and saves number consing to boot.
29 ;;;
30 (defmacro def-math-rtn (name num-args)
31 (let ((function (intern (concatenate 'simple-string
32 "%"
33 (string-upcase name)))))
34 `(progn
35 (declaim (inline ,function))
36 (export ',function)
37 (alien:def-alien-routine (,name ,function) double-float
38 ,@(let ((results nil))
39 (dotimes (i num-args (nreverse results))
40 (push (list (intern (format nil "ARG-~D" i))
41 'double-float)
42 results)))))))
43
44 (eval-when (compile load eval)
45
46 (defun handle-reals (function var)
47 `((((foreach fixnum single-float bignum ratio))
48 (coerce (,function (coerce ,var 'double-float)) 'single-float))
49 ((double-float)
50 (,function ,var))
51 #+double-double
52 ((double-double-float)
53 (,(symbolicate "DD-" function) ,var))))
54
55 ); eval-when (compile load eval)
56
57
58 ;;;; Stubs for the Unix math library.
59
60 ;;; Please refer to the Unix man pages for details about these routines.
61
62 ;;; Trigonometric.
63 #-(and x86 (not sse2))
64 (progn
65 ;; For x86 (without sse2), we can use x87 instructions to implement
66 ;; these. With sse2, we don't currently support that, so these
67 ;; should be disabled.
68 (def-math-rtn "sin" 1)
69 (def-math-rtn "cos" 1)
70 (def-math-rtn "tan" 1)
71 (def-math-rtn "atan" 1)
72 (def-math-rtn "atan2" 2))
73 (def-math-rtn "asin" 1)
74 (def-math-rtn "acos" 1)
75 (def-math-rtn "sinh" 1)
76 (def-math-rtn "cosh" 1)
77 (def-math-rtn "tanh" 1)
78 (def-math-rtn "asinh" 1)
79 (def-math-rtn "acosh" 1)
80 (def-math-rtn "atanh" 1)
81
82 ;;; Exponential and Logarithmic.
83 #-(and x86 (not sse2))
84 (progn
85 (def-math-rtn "exp" 1)
86 (def-math-rtn "log" 1)
87 (def-math-rtn "log10" 1))
88
89 (def-math-rtn "pow" 2)
90 #-(or x86 sparc-v7 sparc-v8 sparc-v9)
91 (def-math-rtn "sqrt" 1)
92 (def-math-rtn "hypot" 2)
93
94 ;; Don't want log1p to use the x87 instruction.
95 #-(or hpux (and x86 (not sse2)))
96 (def-math-rtn "log1p" 1)
97
98 ;; These are needed for use by byte-compiled files. But don't use
99 ;; these with sse2 since we don't support using the x87 instructions
100 ;; here.
101 #+(and x86 (not sse2))
102 (progn
103 #+nil
104 (defun %sin (x)
105 (declare (double-float x)
106 (values double-float))
107 (%sin x))
108 (defun %sin-quick (x)
109 (declare (double-float x)
110 (values double-float))
111 (%sin-quick x))
112 #+nil
113 (defun %cos (x)
114 (declare (double-float x)
115 (values double-float))
116 (%cos x))
117 (defun %cos-quick (x)
118 (declare (double-float x)
119 (values double-float))
120 (%cos-quick x))
121 #+nil
122 (defun %tan (x)
123 (declare (double-float x)
124 (values double-float))
125 (%tan x))
126 (defun %tan-quick (x)
127 (declare (double-float x)
128 (values double-float))
129 (%tan-quick x))
130 (defun %atan (x)
131 (declare (double-float x)
132 (values double-float))
133 (%atan x))
134 (defun %atan2 (x y)
135 (declare (double-float x y)
136 (values double-float))
137 (%atan2 x y))
138 (defun %exp (x)
139 (declare (double-float x)
140 (values double-float))
141 (%exp x))
142 (defun %log (x)
143 (declare (double-float x)
144 (values double-float))
145 (%log x))
146 (defun %log10 (x)
147 (declare (double-float x)
148 (values double-float))
149 (%log10 x))
150 #+nil ;; notyet
151 (defun %pow (x y)
152 (declare (type (double-float 0d0) x)
153 (double-float y)
154 (values (double-float 0d0)))
155 (%pow x y))
156 (defun %sqrt (x)
157 (declare (double-float x)
158 (values double-float))
159 (%sqrt x))
160 (defun %scalbn (f ex)
161 (declare (double-float f)
162 (type (signed-byte 32) ex)
163 (values double-float))
164 (%scalbn f ex))
165 (defun %scalb (f ex)
166 (declare (double-float f ex)
167 (values double-float))
168 (%scalb f ex))
169 (defun %logb (x)
170 (declare (double-float x)
171 (values double-float))
172 (%logb x))
173 (defun %log1p (x)
174 (declare (double-float x)
175 (values double-float))
176 (%log1p x))
177 ) ; progn
178
179
180 ;; As above for x86. It also seems to be needed to handle
181 ;; constant-folding in the compiler.
182 #+sparc
183 (progn
184 (defun %sqrt (x)
185 (declare (double-float x)
186 (values double-float))
187 (%sqrt x))
188 )
189
190 ;;; The standard libm routines for sin, cos, and tan on x86 (Linux)
191 ;;; and ppc are not very accurate for large arguments when compared to
192 ;;; sparc (and maxima). This is basically caused by the fact that
193 ;;; those libraries do not do an accurate argument reduction. The
194 ;;; following functions use some routines Sun's free fdlibm library to
195 ;;; do accurate reduction. Then we call the standard C functions (or
196 ;;; vops for x86) on the reduced argument. This produces much more
197 ;;; accurate values.
198
199 #+(or ppc x86)
200 (progn
201 (declaim (inline %%ieee754-rem-pi/2))
202 ;; Basic argument reduction routine. It returns two values: n and y
203 ;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
204 ;; which octant the arg lies. Y is actually computed in two parts,
205 ;; y[0] and y[1] such that the sum is y, for accuracy.
206
207 (alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
208 (x double-float)
209 (y (* double-float)))
210
211 ;; Same as above, but instead of needing to pass an array in, the
212 ;; output array is broken up into two output values instead. This is
213 ;; easier for the user, and we don't have to wrap calls with
214 ;; without-gcing.
215 (declaim (inline %ieee754-rem-pi/2))
216 (alien:def-alien-routine ("ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
217 (x double-float)
218 (y0 double-float :out)
219 (y1 double-float :out))
220
221 )
222
223 #+(or ppc sse2)
224 (progn
225 (declaim (inline %%sin %%cos %%tan))
226 (macrolet ((frob (alien-name lisp-name)
227 `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
228 (x double-float))))
229 (frob "sin" %%sin)
230 (frob "cos" %%cos)
231 (frob "tan" %%tan))
232 )
233
234 #+(or ppc x86)
235 (macrolet
236 ((frob (sin cos tan)
237 `(progn
238 ;; In all of the routines below, we just compute the sum of
239 ;; y0 and y1 and use that as the (reduced) argument for the
240 ;; trig functions. This is slightly less accurate than what
241 ;; fdlibm does, which calls special functions using y0 and
242 ;; y1 separately, for greater accuracy. This isn't
243 ;; implemented, and some spot checks indicate that what we
244 ;; have here is accurate.
245 ;;
246 ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
247 ;; probably too restrictive.
248 (defun %sin (x)
249 (declare (double-float x))
250 (if (< (abs x) (/ pi 4))
251 (,sin x)
252 ;; Argument reduction needed
253 (multiple-value-bind (n y0 y1)
254 (%ieee754-rem-pi/2 x)
255 (let ((reduced (+ y0 y1)))
256 (case (logand n 3)
257 (0 (,sin reduced))
258 (1 (,cos reduced))
259 (2 (- (,sin reduced)))
260 (3 (- (,cos reduced))))))))
261 (defun %cos (x)
262 (declare (double-float x))
263 (if (< (abs x) (/ pi 4))
264 (,cos x)
265 ;; Argument reduction needed
266 (multiple-value-bind (n y0 y1)
267 (%ieee754-rem-pi/2 x)
268 (let ((reduced (+ y0 y1)))
269 (case (logand n 3)
270 (0 (,cos reduced))
271 (1 (- (,sin reduced)))
272 (2 (- (,cos reduced)))
273 (3 (,sin reduced)))))))
274 (defun %tan (x)
275 (declare (double-float x))
276 (if (< (abs x) (/ pi 4))
277 (,tan x)
278 ;; Argument reduction needed
279 (multiple-value-bind (n y0 y1)
280 (%ieee754-rem-pi/2 x)
281 (let ((reduced (+ y0 y1)))
282 (if (evenp n)
283 (,tan reduced)
284 (- (/ (,tan reduced)))))))))))
285 ;; Don't want %sin-quick and friends with sse2.
286 #+(and x86 (not sse2))
287 (frob %sin-quick %cos-quick %tan-quick)
288 #+(or ppc sse2)
289 (frob %%sin %%cos %%tan))
290
291
292
293 ;;;; Power functions.
294
295 (defun exp (number)
296 _N"Return e raised to the power NUMBER."
297 (number-dispatch ((number number))
298 (handle-reals %exp number)
299 ((complex)
300 (* (exp (realpart number))
301 (cis (imagpart number))))))
302
303 ;;; INTEXP -- Handle the rational base, integer power case.
304
305 (defparameter *intexp-maximum-exponent* 10000)
306
307 (define-condition intexp-limit-error (error)
308 ((base :initarg :base :reader intexp-base)
309 (power :initarg :power :reader intexp-power))
310 (:report (lambda (condition stream)
311 (format stream _"The absolute value of ~S exceeds limit ~S."
312 (intexp-power condition)
313 *intexp-maximum-exponent*))))
314
315 ;;; This function precisely calculates base raised to an integral power. It
316 ;;; separates the cases by the sign of power, for efficiency reasons, as powers
317 ;;; can be calculated more efficiently if power is a positive integer. Values
318 ;;; of power are calculated as positive integers, and inverted if negative.
319 ;;;
320 (defun intexp (base power)
321 ;; Handle the special case of 1^power. Maxima sometimes does this,
322 ;; and there's no need to cause a continuable error in this case.
323 ;; Should we also handle (-1)^power?
324 (when (eql base 1)
325 (return-from intexp base))
326
327 (when (> (abs power) *intexp-maximum-exponent*)
328 ;; Allow user the option to continue with calculation, possibly
329 ;; increasing the limit to the given power.
330 (restart-case
331 (error 'intexp-limit-error
332 :base base
333 :power power)
334 (continue ()
335 :report (lambda (condition stream)
336 (declare (ignore condition))
337 (write-string _"Continue with calculation" stream)))
338 (new-limit ()
339 :report (lambda (condition stream)
340 (declare (ignore condition))
341 (write-string _"Continue with calculation, update limit" stream))
342 (setq *intexp-maximum-exponent* (abs power)))))
343 (cond ((minusp power)
344 (/ (intexp base (- power))))
345 ((eql base 2)
346 (ash 1 power))
347 (t
348 (do ((nextn (ash power -1) (ash power -1))
349 (total (if (oddp power) base 1)
350 (if (oddp power) (* base total) total)))
351 ((zerop nextn) total)
352 (setq base (* base base))
353 (setq power nextn)))))
354
355
356 ;;; EXPT -- Public
357 ;;;
358 ;;; If an integer power of a rational, use INTEXP above. Otherwise, do
359 ;;; floating point stuff. If both args are real, we try %POW right off,
360 ;;; assuming it will return 0 if the result may be complex. If so, we call
361 ;;; COMPLEX-POW which directly computes the complex result. We also separate
362 ;;; the complex-real and real-complex cases from the general complex case.
363 ;;;
364 (defun expt (base power)
365 _N"Returns BASE raised to the POWER."
366 (if (zerop power)
367 ;; CLHS says that if the power is 0, the result is 1, subject to
368 ;; numeric contagion. But what happens if base is infinity or
369 ;; NaN? Do we silently return 1? For now, I think we should
370 ;; signal an error if the FP modes say so.
371 (let ((result (1+ (* base power))))
372 ;; If we get an NaN here, that means base*power above didn't
373 ;; produce 0 and FP traps were disabled, so we handle that
374 ;; here. Should this be a continuable restart?
375 (if (and (floatp result) (float-nan-p result))
376 (float 1 result)
377 result))
378 (labels (;; determine if the double float is an integer.
379 ;; 0 - not an integer
380 ;; 1 - an odd int
381 ;; 2 - an even int
382 (isint (ihi lo)
383 (declare (type (unsigned-byte 31) ihi)
384 (type (unsigned-byte 32) lo)
385 (optimize (speed 3) (safety 0)))
386 (let ((isint 0))
387 (declare (type fixnum isint))
388 (cond ((>= ihi #x43400000) ; exponent >= 53
389 (setq isint 2))
390 ((>= ihi #x3ff00000)
391 (let ((k (- (ash ihi -20) #x3ff))) ; exponent
392 (declare (type (mod 53) k))
393 (cond ((> k 20)
394 (let* ((shift (- 52 k))
395 (j (logand (ash lo (- shift))))
396 (j2 (ash j shift)))
397 (declare (type (mod 32) shift)
398 (type (unsigned-byte 32) j j2))
399 (when (= j2 lo)
400 (setq isint (- 2 (logand j 1))))))
401 ((= lo 0)
402 (let* ((shift (- 20 k))
403 (j (ash ihi (- shift)))
404 (j2 (ash j shift)))
405 (declare (type (mod 32) shift)
406 (type (unsigned-byte 31) j j2))
407 (when (= j2 ihi)
408 (setq isint (- 2 (logand j 1))))))))))
409 isint))
410 (real-expt (x y rtype)
411 (let ((x (coerce x 'double-float))
412 (y (coerce y 'double-float)))
413 (declare (double-float x y))
414 (let* ((x-hi (kernel:double-float-high-bits x))
415 (x-lo (kernel:double-float-low-bits x))
416 (x-ihi (logand x-hi #x7fffffff))
417 (y-hi (kernel:double-float-high-bits y))
418 (y-lo (kernel:double-float-low-bits y))
419 (y-ihi (logand y-hi #x7fffffff)))
420 (declare (type (signed-byte 32) x-hi y-hi)
421 (type (unsigned-byte 31) x-ihi y-ihi)
422 (type (unsigned-byte 32) x-lo y-lo))
423 ;; y==zero: x**0 = 1
424 (when (zerop (logior y-ihi y-lo))
425 (return-from real-expt (coerce 1d0 rtype)))
426 ;; +-NaN return x+y
427 (when (or (> x-ihi #x7ff00000)
428 (and (= x-ihi #x7ff00000) (/= x-lo 0))
429 (> y-ihi #x7ff00000)
430 (and (= y-ihi #x7ff00000) (/= y-lo 0)))
431 (return-from real-expt (coerce (+ x y) rtype)))
432 (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
433 (declare (type fixnum yisint))
434 ;; special value of y
435 (when (and (zerop y-lo) (= y-ihi #x7ff00000))
436 ;; y is +-inf
437 (return-from real-expt
438 (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
439 ;; +-1**inf is NaN
440 (coerce (- y y) rtype))
441 ((>= x-ihi #x3ff00000)
442 ;; (|x|>1)**+-inf = inf,0
443 (if (>= y-hi 0)
444 (coerce y rtype)
445 (coerce 0 rtype)))
446 (t
447 ;; (|x|<1)**-,+inf = inf,0
448 (if (< y-hi 0)
449 (coerce (- y) rtype)
450 (coerce 0 rtype))))))
451
452 (let ((abs-x (abs x)))
453 (declare (double-float abs-x))
454 ;; special value of x
455 (when (and (zerop x-lo)
456 (or (= x-ihi #x7ff00000) (zerop x-ihi)
457 (= x-ihi #x3ff00000)))
458 ;; x is +-0,+-inf,+-1
459 (let ((z (if (< y-hi 0)
460 (/ 1 abs-x) ; z = (1/|x|)
461 abs-x)))
462 (declare (double-float z))
463 (when (< x-hi 0)
464 (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
465 ;; (-1)**non-int
466 (let ((y*pi (* y pi)))
467 (declare (double-float y*pi))
468 (return-from real-expt
469 (complex
470 (coerce (%cos y*pi) rtype)
471 (coerce (%sin y*pi) rtype)))))
472 ((= yisint 1)
473 ;; (x<0)**odd = -(|x|**odd)
474 (setq z (- z)))))
475 (return-from real-expt (coerce z rtype))))
476
477 (if (>= x-hi 0)
478 ;; x>0
479 (coerce (kernel::%pow x y) rtype)
480 ;; x<0
481 (let ((pow (kernel::%pow abs-x y)))
482 (declare (double-float pow))
483 (case yisint
484 (1 ; Odd
485 (coerce (* -1d0 pow) rtype))
486 (2 ; Even
487 (coerce pow rtype))
488 (t ; Non-integer
489 (let ((y*pi (* y pi)))
490 (declare (double-float y*pi))
491 (complex
492 (coerce (* pow (%cos y*pi)) rtype)
493 (coerce (* pow (%sin y*pi)) rtype)))))))))))))
494 ;; This is really messy and should be cleaned up. The easiest
495 ;; way to see if we're doing what we should is the macroexpand
496 ;; the number-dispatch and check each branch.
497 ;;
498 ;; We try to apply the rule of float precision contagion (CLHS
499 ;; 12.1.4.4): the result has the same precision has the most
500 ;; precise argument.
501 (number-dispatch ((base number) (power number))
502 (((foreach fixnum (or bignum ratio) (complex rational))
503 integer)
504 (intexp base power))
505 (((foreach single-float double-float)
506 rational)
507 (real-expt base power '(dispatch-type base)))
508 (((foreach fixnum (or bignum ratio) single-float)
509 (foreach ratio single-float))
510 (real-expt base power 'single-float))
511 (((foreach fixnum (or bignum ratio) single-float double-float)
512 double-float)
513 (real-expt base power 'double-float))
514 ((double-float single-float)
515 (real-expt base power 'double-float))
516 #+double-double
517 (((foreach fixnum (or bignum ratio) single-float double-float
518 double-double-float)
519 double-double-float)
520 (dd-%pow (coerce base 'double-double-float) power))
521 #+double-double
522 ((double-double-float
523 (foreach fixnum (or bignum ratio) single-float double-float))
524 (dd-%pow base (coerce power 'double-double-float)))
525 (((foreach (complex rational) (complex single-float) (complex double-float)
526 #+double-double (complex double-double-float))
527 rational)
528 (* (expt (abs base) power)
529 (cis (* power (phase base)))))
530 #+double-double
531 ((double-double-float
532 complex)
533 (if (and (zerop base) (plusp (realpart power)))
534 (* base power)
535 (exp (* power (* (log2 base 1w0) (log 2w0))))))
536 (((foreach fixnum (or bignum ratio) single-float double-float)
537 (foreach (complex double-float)))
538 ;; Result should have double-float accuracy. Use log2 in
539 ;; case the base won't fit in a double-float.
540 (if (and (zerop base) (plusp (realpart power)))
541 (* base power)
542 (exp (* power (* (log2 base) (log 2d0))))))
543 ((double-float
544 (foreach (complex rational) (complex single-float)))
545 (if (and (zerop base) (plusp (realpart power)))
546 (* base power)
547 (exp (* power (log base)))))
548 #+double-double
549 (((foreach fixnum (or bignum ratio) single-float double-float)
550 (foreach (complex double-double-float)))
551 ;; Result should have double-double-float accuracy. Use log2
552 ;; in case the base won't fit in a double-float.
553 (if (and (zerop base) (plusp (realpart power)))
554 (* base power)
555 (exp (* power (* (log2 base 1w0) (log 2w0))))))
556 (((foreach fixnum (or bignum ratio) single-float)
557 (foreach (complex single-float)))
558 (if (and (zerop base) (plusp (realpart power)))
559 (* base power)
560 (exp (* power (log base)))))
561 (((foreach (complex rational) (complex single-float))
562 (foreach single-float (complex single-float)))
563 (if (and (zerop base) (plusp (realpart power)))
564 (* base power)
565 (exp (* power (log base)))))
566 (((foreach (complex rational) (complex single-float))
567 (foreach double-float (complex double-float)))
568 (if (and (zerop base) (plusp (realpart power)))
569 (* base power)
570 (exp (* power (log (coerce base '(complex double-float)))))))
571 #+double-double
572 (((foreach (complex rational) (complex single-float))
573 (foreach double-double-float (complex double-double-float)))
574 (if (and (zerop base) (plusp (realpart power)))
575 (* base power)
576 (exp (* power (log (coerce base '(complex double-double-float)))))))
577 (((foreach (complex double-float))
578 (foreach single-float double-float (complex single-float)
579 (complex double-float)))
580 (if (and (zerop base) (plusp (realpart power)))
581 (* base power)
582 (exp (* power (log base)))))
583 #+double-double
584 (((foreach (complex double-float))
585 (foreach double-double-float (complex double-double-float)))
586 (if (and (zerop base) (plusp (realpart power)))
587 (* base power)
588 (exp (* power (log (coerce base '(complex double-double-float)))))))
589 #+double-double
590 (((foreach (complex double-double-float))
591 (foreach float (complex float)))
592 (if (and (zerop base) (plusp (realpart power)))
593 (* base power)
594 (exp (* power (log base)))))))))
595
596 ;; Log base 2 of a real number. The result is a either a double-float
597 ;; or double-double-float number (real or complex, as appropriate),
598 ;; depending on the type of FLOAT-TYPE.
599 (defun log2 (x &optional (float-type 1d0))
600 (labels ((log-of-2 (f)
601 ;; log(2), with the precision specified by the type of F
602 (number-dispatch ((f real))
603 ((double-float)
604 #.(log 2d0))
605 #+double-double
606 ((double-double-float)
607 #.(log 2w0))))
608 (log-2-pi (f)
609 ;; log(pi), with the precision specified by the type of F
610 (number-dispatch ((f real))
611 ((double-float)
612 #.(/ pi (log 2d0)))
613 #+double-double
614 ((double-double-float)
615 #.(/ dd-pi (log 2w0)))))
616 (log1p (x)
617 ;; log(1+x), with the precision specified by the type of
618 ;; X
619 (number-dispatch ((x real))
620 (((foreach single-float double-float))
621 (%log1p (float x 1d0)))
622 #+double-double
623 ((double-double-float)
624 (dd-%log1p x))))
625 (log2-bignum (bignum)
626 ;; Write x = 2^n*f where 1/2 < f <= 1. Then log2(x) = n
627 ;; + log2(f).
628 ;;
629 ;; So we grab the top few bits of x and scale that
630 ;; appropriately, take the log of it and add it to n.
631 ;;
632 ;; Return n and log2(f) separately.
633 (if (minusp bignum)
634 (multiple-value-bind (n frac)
635 (log2-bignum (abs bignum))
636 (values n (complex frac (log-2-pi float-type))))
637 (let ((n (integer-length bignum))
638 (float-bits (float-digits float-type)))
639 (if (< n float-bits)
640 (values 0 (log (float bignum float-type)
641 (float 2 float-type)))
642 (let ((exp (min float-bits n))
643 (f (ldb (byte float-bits
644 (max 0 (- n float-bits)))
645 bignum)))
646 (values n (log (scale-float (float f float-type) (- exp))
647 (float 2 float-type)))))))))
648 (etypecase x
649 (float
650 (/ (log (float x float-type)) (log-of-2 float-type)))
651 (ratio
652 (let ((top (numerator x))
653 (bot (denominator x)))
654 ;; If the number of bits in the numerator and
655 ;; denominator are different, just use the fact
656 ;; log(x/y) = log(x) - log(y). But to preserve
657 ;; accuracy, we actually do
658 ;; (log2(x)-log2(y))/log2(e)).
659 ;;
660 ;; However, if the numerator and denominator have the
661 ;; same number of bits, implying the quotient is near
662 ;; one, we use log1p(x) = log(1+x). Since the number is
663 ;; rational, we don't lose precision subtracting 1 from
664 ;; it, and converting it to double-float is accurate.
665 (if (= (integer-length top)
666 (integer-length bot))
667 (/ (log1p (float (- x 1) float-type))
668 (log-of-2 float-type))
669 (multiple-value-bind (top-n top-frac)
670 (log2-bignum top)
671 (multiple-value-bind (bot-n bot-frac)
672 (log2-bignum bot)
673 (+ (- top-n bot-n)
674 (- top-frac bot-frac)))))))
675 (integer
676 (multiple-value-bind (n frac)
677 (log2-bignum x)
678 (+ n frac))))))
679
680 (defun log (number &optional (base nil base-p))
681 _N"Return the logarithm of NUMBER in the base BASE, which defaults to e."
682 (if base-p
683 (cond ((zerop base)
684 ;; ANSI spec
685 base)
686 ((and (realp number) (realp base))
687 ;; CLHS 12.1.4.1 says
688 ;;
689 ;; When rationals and floats are combined by a
690 ;; numerical function, the rational is first converted
691 ;; to a float of the same format.
692 ;;
693 ;; So assume this applies to floats as well convert all
694 ;; numbers to the largest float format before computing
695 ;; the log.
696 ;;
697 ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
698 (number-dispatch ((number real) (base real))
699 ((double-float
700 (foreach double-float single-float))
701 (/ (log2 number) (log2 base)))
702 (((foreach fixnum bignum ratio)
703 (foreach fixnum bignum ratio single-float))
704 (let* ((result (/ (log2 number) (log2 base))))
705 ;; Figure out the right result type
706 (if (realp result)
707 (coerce result 'single-float)
708 (coerce result '(complex single-float)))))
709 (((foreach fixnum bignum ratio)
710 double-float)
711 (/ (log2 number) (log2 base)))
712 ((single-float
713 (foreach fixnum bignum ratio))
714 (let* ((result (/ (log2 number) (log2 base))))
715 ;; Figure out the right result type
716 (if (realp result)
717 (coerce result 'single-float)
718 (coerce result '(complex single-float)))))
719 ((double-float
720 (foreach fixnum bignum ratio))
721 (/ (log2 number) (log2 base)))
722 ((single-float double-float)
723 (/ (log (coerce number 'double-float)) (log base)))
724 #+double-double
725 ((double-double-float
726 (foreach fixnum bignum ratio))
727 (/ (log2 number 1w0) (log2 base 1w0)))
728 #+double-double
729 ((double-double-float
730 (foreach double-double-float double-float single-float))
731 (/ (log number) (log (coerce base 'double-double-float))))
732 #+double-double
733 (((foreach fixnum bignum ratio)
734 double-double-float)
735 (/ (log2 number 1w0) (log2 base 1w0)))
736 #+double-double
737 (((foreach double-float single-float)
738 double-double-float)
739 (/ (log (coerce number 'double-double-float)) (log base)))
740 (((foreach single-float)
741 (foreach single-float))
742 ;; Converting everything to double-float helps the
743 ;; cases like (log 17 10) = (/ (log 17) (log 10)).
744 ;; This is usually handled above, but if we compute (/
745 ;; (log 17) (log 10)), we get a slightly different
746 ;; answer due to roundoff. This makes it a bit more
747 ;; consistent.
748 ;;
749 ;; FIXME: This probably needs more work.
750 (let ((result (/ (log (float number 1d0))
751 (log (float base 1d0)))))
752 (if (realp result)
753 (coerce result 'single-float)
754 (coerce result '(complex single-float)))))))
755 (t
756 ;; FIXME: This probably needs some work as well.
757 (/ (log number) (log base))))
758 (number-dispatch ((number number))
759 (((foreach fixnum bignum))
760 (if (minusp number)
761 (complex (coerce (log (- number)) 'single-float)
762 (coerce pi 'single-float))
763 (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
764 ((ratio)
765 (if (minusp number)
766 (complex (coerce (log (- number)) 'single-float)
767 (coerce pi 'single-float))
768 ;; What happens when the ratio is close to 1? We need to
769 ;; be careful to preserve accuracy.
770 (let ((top (numerator number))
771 (bot (denominator number)))
772 ;; If the number of bits in the numerator and
773 ;; denominator are different, just use the fact
774 ;; log(x/y) = log(x) - log(y). But to preserve
775 ;; accuracy, we actually do
776 ;; (log2(x)-log2(y))/log2(e)).
777 ;;
778 ;; However, if the numerator and denominator have the
779 ;; same number of bits, implying the quotient is near
780 ;; one, we use log1p(x) = log(1+x). Since the number is
781 ;; rational, we don't lose precision subtracting 1 from
782 ;; it, and converting it to double-float is accurate.
783 (if (= (integer-length top)
784 (integer-length bot))
785 (coerce (%log1p (coerce (- number 1) 'double-float))
786 'single-float)
787 (coerce (/ (- (log2 top) (log2 bot))
788 #.(log (exp 1d0) 2d0))
789 'single-float)))))
790 (((foreach single-float double-float))
791 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
792 ;; Since this doesn't seem to be an implementation issue
793 ;; I (pw) take the Kahan result.
794 (if (< (float-sign number)
795 (coerce 0 '(dispatch-type number)))
796 (complex (log (- number)) (coerce pi '(dispatch-type number)))
797 (coerce (%log (coerce number 'double-float))
798 '(dispatch-type number))))
799 #+double-double
800 ((double-double-float)
801 (let ((hi (kernel:double-double-hi number)))
802 (if (< (float-sign hi) 0d0)
803 (complex (dd-%log (- number)) dd-pi)
804 (dd-%log number))))
805 ((complex)
806 (complex-log number)))))
807
808 (defun sqrt (number)
809 _N"Return the square root of NUMBER."
810 (number-dispatch ((number number))
811 (((foreach fixnum bignum ratio))
812 (if (minusp number)
813 (complex-sqrt number)
814 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
815 (((foreach single-float double-float))
816 (if (minusp number)
817 (complex-sqrt number)
818 (coerce (%sqrt (coerce number 'double-float))
819 '(dispatch-type number))))
820 #+double-double
821 ((double-double-float)
822 (if (minusp number)
823 (dd-complex-sqrt number)
824 (multiple-value-bind (hi lo)
825 (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
826 (kernel:%make-double-double-float hi lo))))
827 ((complex)
828 (complex-sqrt number))))
829
830
831 ;;;; Trigonometic and Related Functions
832
833 (defun abs (number)
834 _N"Returns the absolute value of the number."
835 (number-dispatch ((number number))
836 (((foreach single-float double-float fixnum rational
837 #+double-double double-double-float))
838 (abs number))
839 ((complex)
840 (let ((rx (realpart number))
841 (ix (imagpart number)))
842 (etypecase rx
843 (rational
844 (sqrt (+ (* rx rx) (* ix ix))))
845 (single-float
846 (coerce (%hypot (coerce rx 'double-float)
847 (coerce ix 'double-float))
848 'single-float))
849 (double-float
850 (%hypot rx ix))
851 #+double-double
852 (double-double-float
853 (multiple-value-bind (abs^2 scale)
854 (dd-cssqs number)
855 (scale-float (sqrt abs^2) scale))))))))
856
857 (defun phase (number)
858 _N"Returns the angle part of the polar representation of a complex number.
859 For complex numbers, this is (atan (imagpart number) (realpart number)).
860 For non-complex positive numbers, this is 0. For non-complex negative
861 numbers this is PI."
862 (etypecase number
863 (rational
864 (if (minusp number)
865 (coerce pi 'single-float)
866 0.0f0))
867 (single-float
868 (if (minusp (float-sign number))
869 (coerce pi 'single-float)
870 0.0f0))
871 (double-float
872 (if (minusp (float-sign number))
873 (coerce pi 'double-float)
874 0.0d0))
875 #+double-double
876 (double-double-float
877 (if (minusp (float-sign number))
878 dd-pi
879 0w0))
880 (complex
881 (atan (imagpart number) (realpart number)))))
882
883
884 (defun sin (number)
885 _N"Return the sine of NUMBER."
886 (number-dispatch ((number number))
887 (handle-reals %sin number)
888 ((complex)
889 (let ((x (realpart number))
890 (y (imagpart number)))
891 (complex (* (sin x) (cosh y))
892 (* (cos x) (sinh y)))))))
893
894 (defun cos (number)
895 _N"Return the cosine of NUMBER."
896 (number-dispatch ((number number))
897 (handle-reals %cos number)
898 ((complex)
899 (let ((x (realpart number))
900 (y (imagpart number)))
901 (complex (* (cos x) (cosh y))
902 (- (* (sin x) (sinh y))))))))
903
904 (defun tan (number)
905 _N"Return the tangent of NUMBER."
906 (number-dispatch ((number number))
907 (handle-reals %tan number)
908 ((complex)
909 (complex-tan number))))
910
911 (defun cis (theta)
912 _N"Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
913 (if (complexp theta)
914 (error _"Argument to CIS is complex: ~S" theta)
915 (complex (cos theta) (sin theta))))
916
917 (defun asin (number)
918 _N"Return the arc sine of NUMBER."
919 (number-dispatch ((number number))
920 ((rational)
921 (if (or (> number 1) (< number -1))
922 (complex-asin number)
923 (coerce (%asin (coerce number 'double-float)) 'single-float)))
924 (((foreach single-float double-float))
925 (if (or (float-nan-p number)
926 (and (<= number (coerce 1 '(dispatch-type number)))
927 (>= number (coerce -1 '(dispatch-type number)))))
928 (coerce (%asin (coerce number 'double-float))
929 '(dispatch-type number))
930 (complex-asin number)))
931 #+double-double
932 ((double-double-float)
933 (if (or (float-nan-p number)
934 (and (<= number 1w0)
935 (>= number -1w0)))
936 (dd-%asin number)
937 (dd-complex-asin number)))
938 ((complex)
939 (complex-asin number))))
940
941 (defun acos (number)
942 _N"Return the arc cosine of NUMBER."
943 (number-dispatch ((number number))
944 ((rational)
945 (if (or (> number 1) (< number -1))
946 (complex-acos number)
947 (coerce (%acos (coerce number 'double-float)) 'single-float)))
948 (((foreach single-float double-float))
949 (if (or (float-nan-p number)
950 (and (<= number (coerce 1 '(dispatch-type number)))
951 (>= number (coerce -1 '(dispatch-type number)))))
952 (coerce (%acos (coerce number 'double-float))
953 '(dispatch-type number))
954 (complex-acos number)))
955 #+double-double
956 ((double-double-float)
957 (if (or (float-nan-p number)
958 (and (<= number 1w0)
959 (>= number -1w0)))
960 (dd-%acos number)
961 (complex-acos number)))
962 ((complex)
963 (complex-acos number))))
964
965
966 (defun atan (y &optional (x nil xp))
967 _N"Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
968 (if xp
969 (flet ((atan2 (y x)
970 (declare (type double-float y x)
971 (values double-float))
972 (if (zerop x)
973 (if (zerop y)
974 (if (plusp (float-sign x))
975 y
976 (float-sign y pi))
977 (float-sign y (/ pi 2)))
978 (%atan2 y x))))
979 ;; If X is given, both X and Y must be real numbers.
980 (number-dispatch ((y real) (x real))
981 ((double-float
982 (foreach double-float single-float fixnum bignum ratio))
983 (atan2 y (coerce x 'double-float)))
984 (((foreach single-float fixnum bignum ratio)
985 double-float)
986 (atan2 (coerce y 'double-float) x))
987 (((foreach single-float fixnum bignum ratio)
988 (foreach single-float fixnum bignum ratio))
989 (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
990 'single-float))
991 #+double-double
992 ((double-double-float
993 (foreach double-double-float double-float single-float fixnum bignum ratio))
994 (dd-%atan2 y (coerce x 'double-double-float)))
995 #+double-double
996 (((foreach double-float single-float fixnum bignum ratio)
997 double-double-float)
998 (dd-%atan2 (coerce y 'double-double-float) x))))
999 (number-dispatch ((y number))
1000 (handle-reals %atan y)
1001 ((complex)
1002 (complex-atan y)))))
1003
1004 (defun sinh (number)
1005 _N"Return the hyperbolic sine of NUMBER."
1006 (number-dispatch ((number number))
1007 (handle-reals %sinh number)
1008 ((complex)
1009 (let ((x (realpart number))
1010 (y (imagpart number)))
1011 (complex (* (sinh x) (cos y))
1012 (* (cosh x) (sin y)))))))
1013
1014 (defun cosh (number)
1015 _N"Return the hyperbolic cosine of NUMBER."
1016 (number-dispatch ((number number))
1017 (handle-reals %cosh number)
1018 ((complex)
1019 (let ((x (realpart number))
1020 (y (imagpart number)))
1021 (complex (* (cosh x) (cos y))
1022 (* (sinh x) (sin y)))))))
1023
1024 (defun tanh (number)
1025 _N"Return the hyperbolic tangent of NUMBER."
1026 (number-dispatch ((number number))
1027 (handle-reals %tanh number)
1028 ((complex)
1029 (complex-tanh number))))
1030
1031 (defun asinh (number)
1032 _N"Return the hyperbolic arc sine of NUMBER."
1033 (number-dispatch ((number number))
1034 (handle-reals %asinh number)
1035 ((complex)
1036 (complex-asinh number))))
1037
1038 (defun acosh (number)
1039 _N"Return the hyperbolic arc cosine of NUMBER."
1040 (number-dispatch ((number number))
1041 ((rational)
1042 ;; acosh is complex if number < 1
1043 (if (< number 1)
1044 (complex-acosh number)
1045 (coerce (%acosh (coerce number 'double-float)) 'single-float)))
1046 (((foreach single-float double-float))
1047 (if (< number (coerce 1 '(dispatch-type number)))
1048 (complex-acosh number)
1049 (coerce (%acosh (coerce number 'double-float))
1050 '(dispatch-type number))))
1051 #+double-double
1052 ((double-double-float)
1053 (if (< number 1w0)
1054 (complex-acosh number)
1055 (dd-%acosh number)))
1056 ((complex)
1057 (complex-acosh number))))
1058
1059 (defun atanh (number)
1060 _N"Return the hyperbolic arc tangent of NUMBER."
1061 (number-dispatch ((number number))
1062 ((rational)
1063 ;; atanh is complex if |number| > 1
1064 (if (or (> number 1) (< number -1))
1065 (complex-atanh number)
1066 (coerce (%atanh (coerce number 'double-float)) 'single-float)))
1067 (((foreach single-float double-float))
1068 (if (or (> number (coerce 1 '(dispatch-type number)))
1069 (< number (coerce -1 '(dispatch-type number))))
1070 (complex-atanh number)
1071 (coerce (%atanh (coerce number 'double-float))
1072 '(dispatch-type number))))
1073 #+double-double
1074 ((double-double-float)
1075 (if (or (> number 1w0)
1076 (< number -1w0))
1077 (complex-atanh number)
1078 (dd-%atanh (coerce number 'double-double-float))))
1079 ((complex)
1080 (complex-atanh number))))
1081
1082 ;;; HP-UX does not supply a C version of log1p, so use the definition.
1083 ;;; We really need to fix this. The definition really loses big-time
1084 ;;; in roundoff as x gets small.
1085
1086 #+hpux
1087 (declaim (inline %log1p))
1088 #+hpux
1089 (defun %log1p (number)
1090 (declare (double-float number)
1091 (optimize (speed 3) (safety 0)))
1092 (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
1093
1094
1095 ;;;;
1096 ;;;; This is a set of routines that implement many elementary
1097 ;;;; transcendental functions as specified by ANSI Common Lisp. The
1098 ;;;; implementation is based on Kahan's paper.
1099 ;;;;
1100 ;;;; I believe I have accurately implemented the routines and are
1101 ;;;; correct, but you may want to check for your self.
1102 ;;;;
1103 ;;;; These functions are written for CMU Lisp and take advantage of
1104 ;;;; some of the features available there. It may be possible,
1105 ;;;; however, to port this to other Lisps.
1106 ;;;;
1107 ;;;; Some functions are significantly more accurate than the original
1108 ;;;; definitions in CMU Lisp. In fact, some functions in CMU Lisp
1109 ;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
1110 ;;;; answer is pi + i*log(2-sqrt(3)).
1111 ;;;;
1112 ;;;; All of the implemented functions will take any number for an
1113 ;;;; input, but the result will always be a either a complex
1114 ;;;; single-float or a complex double-float.
1115 ;;;;
1116 ;;;; General functions
1117 ;;;; complex-sqrt
1118 ;;;; complex-log
1119 ;;;; complex-atanh
1120 ;;;; complex-tanh
1121 ;;;; complex-acos
1122 ;;;; complex-acosh
1123 ;;;; complex-asin
1124 ;;;; complex-asinh
1125 ;;;; complex-atan
1126 ;;;; complex-tan
1127 ;;;;
1128 ;;;; Utility functions:
1129 ;;;; scalb logb
1130 ;;;;
1131 ;;;; Internal functions:
1132 ;;;; square coerce-to-complex-type cssqs complex-log-scaled
1133 ;;;;
1134 ;;;;
1135 ;;;; Please send any bug reports, comments, or improvements to Raymond
1136 ;;;; Toy at toy@rtp.ericsson.se.
1137 ;;;;
1138 ;;;; References
1139 ;;;;
1140 ;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
1141 ;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
1142 ;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
1143 ;;;; Press, 1987
1144 ;;;;
1145
1146 (declaim (inline square))
1147 (defun square (x)
1148 (declare (float x))
1149 (* x x))
1150
1151 ;; If you have these functions in libm, perhaps they should be used
1152 ;; instead of these Lisp versions. These versions are probably good
1153 ;; enough, especially since they are portable.
1154
1155 (declaim (inline scalb))
1156 (defun scalb (x n)
1157 _N"Compute 2^N * X without compute 2^N first (use properties of the
1158 underlying floating-point format"
1159 (declare (type float x)
1160 (type double-float-exponent n))
1161 (scale-float x n))
1162
1163 (declaim (inline logb-finite))
1164 (defun logb-finite (x)
1165 _N"Same as logb but X is not infinity and non-zero and not a NaN, so
1166 that we can always return an integer"
1167 (declare (type float x))
1168 (multiple-value-bind (signif expon sign)
1169 (decode-float x)
1170 (declare (ignore signif sign))
1171 ;; decode-float is almost right, except that the exponent
1172 ;; is off by one
1173 (1- expon)))
1174
1175 (defun logb (x)
1176 _N"Compute an integer N such that 1 <= |2^(-N) * x| < 2.
1177 For the special cases, the following values are used:
1178
1179 x logb
1180 NaN NaN
1181 +/- infinity +infinity
1182 0 -infinity
1183 "
1184 (declare (type float x))
1185 (cond ((float-nan-p x)
1186 x)
1187 ((float-infinity-p x)
1188 #.ext:double-float-positive-infinity)
1189 ((zerop x)
1190 ;; The answer is negative infinity, but we are supposed to
1191 ;; signal divide-by-zero, so do the actual division
1192 (/ -1 x)
1193 )
1194 (t
1195 (logb-finite x))))
1196
1197
1198
1199 ;; This function is used to create a complex number of the appropriate
1200 ;; type.
1201
1202 (declaim (inline coerce-to-complex-type))
1203 (defun coerce-to-complex-type (x y z)
1204 _N"Create complex number with real part X and imaginary part Y such that
1205 it has the same type as Z. If Z has type (complex rational), the X
1206 and Y are coerced to single-float."
1207 (declare (double-float x y)
1208 (number z)
1209 (optimize (extensions:inhibit-warnings 3)))
1210 (if (typep (realpart z) 'double-float)
1211 (complex x y)
1212 ;; Convert anything that's not a double-float to a single-float.
1213 (complex (float x 1f0)
1214 (float y 1f0))))
1215
1216 (defun cssqs (z)
1217 ;; Compute |(x+i*y)/2^k|^2 scaled to avoid over/underflow. The
1218 ;; result is r + i*k, where k is an integer.
1219
1220 ;; Save all FP flags
1221 (let ((x (float (realpart z) 1d0))
1222 (y (float (imagpart z) 1d0)))
1223 ;; Would this be better handled using an exception handler to
1224 ;; catch the overflow or underflow signal? For now, we turn all
1225 ;; traps off and look at the accrued exceptions to see if any
1226 ;; signal would have been raised.
1227 (with-float-traps-masked (:underflow :overflow)
1228 (let ((rho (+ (square x) (square y))))
1229 (declare (optimize (speed 3) (space 0)))
1230 (cond ((and (or (float-nan-p rho)
1231 (float-infinity-p rho))
1232 (or (float-infinity-p (abs x))
1233 (float-infinity-p (abs y))))
1234 (values ext:double-float-positive-infinity 0))
1235 ((let ((threshold #.(/ least-positive-double-float
1236 double-float-epsilon))
1237 (traps (ldb vm::float-sticky-bits
1238 (vm:floating-point-modes))))
1239 ;; Overflow raised or (underflow raised and rho <
1240 ;; lambda/eps)
1241 (or (not (zerop (logand vm:float-overflow-trap-bit traps)))
1242 (and (not (zerop (logand vm:float-underflow-trap-bit traps)))
1243 (< rho threshold))))
1244 ;; If we're here, neither x nor y are infinity and at
1245 ;; least one is non-zero.. Thus logb returns a nice
1246 ;; integer.
1247 (let ((k (- (logb-finite (max (abs x) (abs y))))))
1248 (values (+ (square (scalb x k))
1249 (square (scalb y k)))
1250 (- k))))
1251 (t
1252 (values rho 0)))))))
1253
1254 (defun complex-sqrt (z)
1255 _N"Principle square root of Z
1256
1257 Z may be any number, but the result is always a complex."
1258 (declare (number z))
1259 #+double-double
1260 (when (typep z '(or double-double-float (complex double-double-float)))
1261 (return-from complex-sqrt (dd-complex-sqrt z)))
1262 (multiple-value-bind (rho k)
1263 (cssqs z)
1264 (declare (type (or (member 0d0) (double-float 0d0)) rho)
1265 (type fixnum k))
1266 (let ((x (float (realpart z) 1.0d0))
1267 (y (float (imagpart z) 1.0d0))
1268 (eta 0d0)
1269 (nu 0d0))
1270 (declare (double-float x y eta nu))
1271
1272 (locally
1273 ;; space 0 to get maybe-inline functions inlined.
1274 (declare (optimize (speed 3) (space 0)))
1275
1276 (if (not (locally (declare (optimize (inhibit-warnings 3)))
1277 (float-nan-p x)))
1278 (setf rho (+ (scalb (abs x) (- k)) (sqrt rho))))
1279
1280 (cond ((oddp k)
1281 (setf k (ash k -1)))
1282 (t
1283 (setf k (1- (ash k -1)))
1284 (setf rho (+ rho rho))))
1285
1286 (setf rho (scalb (sqrt rho) k))
1287
1288 (setf eta rho)
1289 (setf nu y)
1290
1291 (when (/= rho 0d0)
1292 (when (not (float-infinity-p (abs nu)))
1293 (setf nu (/ (/ nu rho) 2d0)))
1294 (when (< x 0d0)
1295 (setf eta (abs nu))
1296 (setf nu (float-sign y rho))))
1297 (coerce-to-complex-type eta nu z)))))
1298
1299 (defun complex-log-scaled (z j)
1300 _N"Compute log(2^j*z).
1301
1302 This is for use with J /= 0 only when |z| is huge."
1303 (declare (number z)
1304 (fixnum j))
1305 ;; The constants t0, t1, t2 should be evaluated to machine
1306 ;; precision. In addition, Kahan says the accuracy of log1p
1307 ;; influences the choices of these constants but doesn't say how to
1308 ;; choose them. We'll just assume his choices matches our
1309 ;; implementation of log1p.
1310 (let ((t0 #.(/ 1 (sqrt 2.0d0)))
1311 (t1 1.2d0)
1312 (t2 3d0)
1313 (ln2 #.(log 2d0))
1314 (x (float (realpart z) 1.0d0))
1315 (y (float (imagpart z) 1.0d0)))
1316 (multiple-value-bind (rho k)
1317 (cssqs z)
1318 (declare (optimize (speed 3)))
1319 (let ((beta (max (abs x) (abs y)))
1320 (theta (min (abs x) (abs y))))
1321 (coerce-to-complex-type (if (and (zerop k)
1322 (< t0 beta)
1323 (or (<= beta t1)
1324 (< rho t2)))
1325 (/ (%log1p (+ (* (- beta 1.0d0)
1326 (+ beta 1.0d0))
1327 (* theta theta)))
1328 2d0)
1329 (+ (/ (log rho) 2d0)
1330 (* (+ k j) ln2)))
1331 (atan y x)
1332 z)))))
1333
1334 (defun complex-log (z)
1335 _N"Log of Z = log |Z| + i * arg Z
1336
1337 Z may be any number, but the result is always a complex."
1338 (declare (number z))
1339 #+double-double
1340 (when (typep z '(or double-double-float (complex double-double-float)))
1341 (return-from complex-log (dd-complex-log-scaled z 0)))
1342 (complex-log-scaled z 0))
1343
1344 ;; Let us note the following "strange" behavior. atanh 1.0d0 is
1345 ;; +infinity, but the following code returns approx 176 + i*pi/4. The
1346 ;; reason for the imaginary part is caused by the fact that arg i*y is
1347 ;; never 0 since we have positive and negative zeroes.
1348
1349 (defun complex-atanh (z)
1350 _N"Compute atanh z = (log(1+z) - log(1-z))/2"
1351 (declare (number z))
1352 #+double-double
1353 (when (typep z '(or double-double-float (complex double-double-float)))
1354 (return-from complex-atanh (dd-complex-atanh z)))
1355
1356 (if (and (realp z) (< z -1))
1357 ;; atanh is continuous in quadrant III in this case.
1358 (complex-atanh (complex z -0f0))
1359 (let* ( ;; Constants
1360 (theta (/ (sqrt most-positive-double-float) 4.0d0))
1361 (rho (/ 4.0d0 (sqrt most-positive-double-float)))
1362 (half-pi (/ pi 2.0d0))
1363 (rp (float (realpart z) 1.0d0))
1364 (beta (float-sign rp 1.0d0))
1365 (x (* beta rp))
1366 (y (* beta (- (float (imagpart z) 1.0d0))))
1367 (eta 0.0d0)
1368 (nu 0.0d0))
1369 ;; Shouldn't need this declare.
1370 (declare (double-float x y))
1371 (locally
1372 (declare (optimize (speed 3)))
1373 (cond ((or (> x theta)
1374 (> (abs y) theta))
1375 ;; To avoid overflow...
1376 (setf nu (float-sign y half-pi))
1377 ;; eta is real part of 1/(x + iy). This is x/(x^2+y^2),
1378 ;; which can cause overflow. Arrange this computation so
1379 ;; that it won't overflow.
1380 (setf eta (let* ((x-bigger (> x (abs y)))
1381 (r (if x-bigger (/ y x) (/ x y)))
1382 (d (+ 1.0d0 (* r r))))
1383 (if x-bigger
1384 (/ (/ x) d)
1385 (/ (/ r y) d)))))
1386 ((= x 1.0d0)
1387 ;; Should this be changed so that if y is zero, eta is set
1388 ;; to +infinity instead of approx 176? In any case
1389 ;; tanh(176) is 1.0d0 within working precision.
1390 (let ((t1 (+ 4d0 (square y)))
1391 (t2 (+ (abs y) rho)))
1392 (setf eta (log (/ (sqrt (sqrt t1))
1393 (sqrt t2))))
1394 (setf nu (* 0.5d0
1395 (float-sign y
1396 (+ half-pi (atan (* 0.5d0 t2))))))))
1397 (t
1398 (let ((t1 (+ (abs y) rho)))
1399 ;; Normal case using log1p(x) = log(1 + x)
1400 (setf eta (* 0.25d0
1401 (%log1p (/ (* 4.0d0 x)
1402 (+ (square (- 1.0d0 x))
1403 (square t1))))))
1404 (setf nu (* 0.5d0
1405 (atan (* 2.0d0 y)
1406 (- (* (- 1.0d0 x)
1407 (+ 1.0d0 x))
1408 (square t1))))))))
1409 (coerce-to-complex-type (* beta eta)
1410 (- (* beta nu))
1411 z)))))
1412
1413 (defun complex-tanh (z)
1414 _N"Compute tanh z = sinh z / cosh z"
1415 (declare (number z))
1416 #+double-double
1417 (when (typep z '(or double-double-float (complex double-double-float)))
1418 (return-from complex-tanh (dd-complex-tanh z)))
1419
1420 (let ((x (float (realpart z) 1.0d0))
1421 (y (float (imagpart z) 1.0d0)))
1422 (locally
1423 ;; space 0 to get maybe-inline functions inlined
1424 (declare (optimize (speed 3) (space 0)))
1425 (cond ((> (abs x)
1426 #-(or linux hpux) #.(/ (%asinh most-positive-double-float) 4d0)
1427 ;; This is more accurate under linux.
1428 #+(or linux hpux) #.(/ (+ (%log 2.0d0)
1429 (%log most-positive-double-float)) 4d0))
1430 (coerce-to-complex-type (float-sign x)
1431 (float-sign y) z))
1432 (t
1433 (let* ((tv (%tan y))
1434 (beta (+ 1.0d0 (* tv tv)))
1435 (s (sinh x))
1436 (rho (sqrt (+ 1.0d0 (* s s)))))
1437 (if (float-infinity-p (abs tv))
1438 (coerce-to-complex-type (/ rho s)
1439 (/ tv)
1440 z)
1441 (let ((den (+ 1.0d0 (* beta s s))))
1442 (coerce-to-complex-type (/ (* beta rho s)
1443 den)
1444 (/ tv den)
1445 z)))))))))
1446
1447 ;; Kahan says we should only compute the parts needed. Thus, the
1448 ;; realpart's below should only compute the real part, not the whole
1449 ;; complex expression. Doing this can be important because we may get
1450 ;; spurious signals that occur in the part that we are not using.
1451 ;;
1452 ;; However, we take a pragmatic approach and just use the whole
1453 ;; expression.
1454
1455 ;; NOTE: The formula given by Kahan is somewhat ambiguous in whether
1456 ;; it's the conjugate of the square root or the square root of the
1457 ;; conjugate. This needs to be checked.
1458
1459 ;; I checked. It doesn't matter because (conjugate (sqrt z)) is the
1460 ;; same as (sqrt (conjugate z)) for all z. This follows because
1461 ;;
1462 ;; (conjugate (sqrt z)) = exp(0.5*log |z|)*exp(-0.5*j*arg z).
1463 ;;
1464 ;; (sqrt (conjugate z)) = exp(0.5*log|z|)*exp(0.5*j*arg conj z)
1465 ;;
1466 ;; and these two expressions are equal if and only if arg conj z =
1467 ;; -arg z, which is clearly true for all z.
1468
1469 ;; NOTE: The rules of Common Lisp says that if you mix a real with a
1470 ;; complex, the real is converted to a complex before performing the
1471 ;; operation. However, Kahan says in this paper (pg 176):
1472 ;;
1473 ;; (iii) Careless handling can turn infinity or the sign of zero into
1474 ;; misinformation that subsequently disappears leaving behind
1475 ;; only a plausible but incorrect result. That is why compilers
1476 ;; must not transform z-1 into z-(1+i*0), as we have seen above,
1477 ;; nor -(-x-x^2) into (x+x^2), as we shall see below, lest a
1478 ;; subsequent logarithm or square root produce a non-zero
1479 ;; imaginary part whose sign is opposite to what was intended.
1480 ;;
1481 ;; The interesting examples are too long and complicated to reproduce
1482 ;; here. We refer the reader to his paper.
1483 ;;
1484 ;; The functions below are intended to handle the cases where a real
1485 ;; is mixed with a complex and we don't want CL complex contagion to
1486 ;; occur..
1487
1488 (declaim (inline 1+z 1-z z-1 z+1))
1489 (defun 1+z (z)
1490 (complex (+ 1 (realpart z)) (imagpart z)))
1491 (defun 1-z (z)
1492 (complex (- 1 (realpart z)) (- (imagpart z))))
1493 (defun z-1 (z)
1494 (complex (- (realpart z) 1) (imagpart z)))
1495 (defun z+1 (z)
1496 (complex (+ (realpart z) 1) (imagpart z)))
1497
1498 (defun complex-acos (z)
1499 _N"Compute acos z = pi/2 - asin z
1500
1501 Z may be any number, but the result is always a complex."
1502 (declare (number z))
1503 #+double-double
1504 (when (typep z '(or double-double-float (complex double-double-float)))
1505 (return-from complex-acos (dd-complex-acos z)))
1506 (if (and (realp z) (> z 1))
1507 ;; acos is continuous in quadrant IV in this case.
1508 (complex-acos (complex z -0f0))
1509 (let ((sqrt-1+z (complex-sqrt (1+z z)))
1510 (sqrt-1-z (complex-sqrt (1-z z))))
1511 (with-float-traps-masked (:divide-by-zero)
1512 (complex (* 2 (atan (/ (realpart sqrt-1-z)
1513 (realpart sqrt-1+z))))
1514 (asinh (imagpart (* (conjugate sqrt-1+z)
1515 sqrt-1-z))))))))
1516
1517 (defun complex-acosh (z)
1518 _N"Compute acosh z = 2 * log(sqrt((z+1)/2) + sqrt((z-1)/2))
1519
1520 Z may be any number, but the result is always a complex."
1521 (declare (number z))
1522 (let ((sqrt-z-1 (complex-sqrt (z-1 z)))
1523 (sqrt-z+1 (complex-sqrt (z+1 z))))
1524 (with-float-traps-masked (:divide-by-zero)
1525 (complex (asinh (realpart (* (conjugate sqrt-z-1)
1526 sqrt-z+1)))
1527 (* 2 (atan (/ (imagpart sqrt-z-1)
1528 (realpart sqrt-z+1))))))))
1529
1530
1531 (defun complex-asin (z)
1532 _N"Compute asin z = asinh(i*z)/i
1533
1534 Z may be any number, but the result is always a complex."
1535 (declare (number z))
1536 #+double-double
1537 (when (typep z '(or double-double-float (complex double-double-float)))
1538 (return-from complex-asin (dd-complex-asin z)))
1539 (if (and (realp z) (> z 1))
1540 ;; asin is continuous in quadrant IV in this case.
1541 (complex-asin (complex z -0f0))
1542 (let ((sqrt-1-z (complex-sqrt (1-z z)))
1543 (sqrt-1+z (complex-sqrt (1+z z))))
1544 (with-float-traps-masked (:divide-by-zero)
1545 (complex (atan (/ (realpart z)
1546 (realpart (* sqrt-1-z sqrt-1+z))))
1547 (asinh (imagpart (* (conjugate sqrt-1-z)
1548 sqrt-1+z))))))))
1549
1550 (defun complex-asinh (z)
1551 _N"Compute asinh z = log(z + sqrt(1 + z*z))
1552
1553 Z may be any number, but the result is always a complex."
1554 (declare (number z))
1555 ;; asinh z = -i * asin (i*z)
1556 #+double-double
1557 (when (typep z '(or double-double-float (complex double-double-float)))
1558 (return-from complex-asinh (dd-complex-asinh z)))
1559 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1560 (result (complex-asin iz)))
1561 (complex (imagpart result)
1562 (- (realpart result)))))
1563
1564 (defun complex-atan (z)
1565 _N"Compute atan z = atanh (i*z) / i
1566
1567 Z may be any number, but the result is always a complex."
1568 (declare (number z))
1569 ;; atan z = -i * atanh (i*z)
1570 #+double-double
1571 (when (typep z '(or double-double-float (complex double-double-float)))
1572 (return-from complex-atan (dd-complex-atan z)))
1573 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1574 (result (complex-atanh iz)))
1575 (complex (imagpart result)
1576 (- (realpart result)))))
1577
1578 (defun complex-tan (z)
1579 _N"Compute tan z = -i * tanh(i * z)
1580
1581 Z may be any number, but the result is always a complex."
1582 (declare (number z))
1583 ;; tan z = -i * tanh(i*z)
1584 #+double-double
1585 (when (typep z '(or double-double-float (complex double-double-float)))
1586 (return-from complex-tan (dd-complex-tan z)))
1587 (let* ((iz (complex (- (imagpart z)) (realpart z)))
1588 (result (complex-tanh iz)))
1589 (complex (imagpart result)
1590 (- (realpart result)))))
1591

  ViewVC Help
Powered by ViewVC 1.1.5