irrat.lisp 55.2 KB
Newer Older
wlott's avatar
wlott committed
1
2
3
;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
;;;
;;; **********************************************************************
4
5
6
7
;;; This code was written as part of the CMU Common Lisp project at
;;; Carnegie Mellon University, and has been placed in the public domain.
;;;
(ext:file-comment
8
  "$Header: src/code/irrat.lisp $")
9
;;;
wlott's avatar
wlott committed
10
11
12
13
14
15
16
17
18
;;; **********************************************************************
;;;
;;; This file contains all the irrational functions.  Actually, most of the
;;; work is done by calling out to C...
;;;
;;; Author: William Lott.
;;; 

(in-package "KERNEL")
19
(intl:textdomain "cmucl")
wlott's avatar
wlott committed
20
21
22
23
24


;;;; Random constants, utility functions, and macros.

(defconstant pi 3.14159265358979323846264338327950288419716939937511L0)
25
;(defconstant e 2.71828182845904523536028747135266249775724709369996L0)
wlott's avatar
wlott committed
26

27
28
29
;;; Make these INLINE, since the call to C is at least as compact as a Lisp
;;; call, and saves number consing to boot.
;;;
wlott's avatar
wlott committed
30
31
32
33
(defmacro def-math-rtn (name num-args)
  (let ((function (intern (concatenate 'simple-string
				       "%"
				       (string-upcase name)))))
34
    `(progn
pw's avatar
pw committed
35
       (declaim (inline ,function))
36
       (export ',function)
wlott's avatar
wlott committed
37
       (alien:def-alien-routine (,name ,function) double-float
38
39
40
41
42
	 ,@(let ((results nil))
	     (dotimes (i num-args (nreverse results))
	       (push (list (intern (format nil "ARG-~D" i))
			   'double-float)
		     results)))))))
wlott's avatar
wlott committed
43
44
45
46
47
48
49

(eval-when (compile load eval)

(defun handle-reals (function var)
  `((((foreach fixnum single-float bignum ratio))
     (coerce (,function (coerce ,var 'double-float)) 'single-float))
    ((double-float)
50
51
52
53
     (,function ,var))
    #+double-double
    ((double-double-float)
     (,(symbolicate "DD-" function) ,var))))
wlott's avatar
wlott committed
54
55
56
57
58
59
60
61
62

); eval-when (compile load eval)


;;;; Stubs for the Unix math library.

;;; Please refer to the Unix man pages for details about these routines.

;;; Trigonometric.
63
64
65
66
67
68
69
70
71
72
#-(and x86 (not sse2))
(progn
  ;; For x86 (without sse2), we can use x87 instructions to implement
  ;; these.  With sse2, we don't currently support that, so these
  ;; should be disabled.
  (def-math-rtn "sin" 1)
  (def-math-rtn "cos" 1)
  (def-math-rtn "tan" 1)
  (def-math-rtn "atan" 1)
  (def-math-rtn "atan2" 2))
wlott's avatar
wlott committed
73
74
75
76
77
(def-math-rtn "asin" 1)
(def-math-rtn "acos" 1)
(def-math-rtn "sinh" 1)
(def-math-rtn "cosh" 1)
(def-math-rtn "tanh" 1)
pw's avatar
pw committed
78
79
80
(def-math-rtn "asinh" 1)
(def-math-rtn "acosh" 1)
(def-math-rtn "atanh" 1)
wlott's avatar
wlott committed
81
82

;;; Exponential and Logarithmic.
83
84
85
86
87
88
#-(and x86 (not sse2))
(progn
  (def-math-rtn "exp" 1)
  (def-math-rtn "log" 1)
  (def-math-rtn "log10" 1))

wlott's avatar
wlott committed
89
(def-math-rtn "pow" 2)
90
91
#-(or x86 sparc-v7 sparc-v8 sparc-v9)
(def-math-rtn "sqrt" 1)
wlott's avatar
wlott committed
92
(def-math-rtn "hypot" 2)
ram's avatar
ram committed
93

94
95
96
97
98
99
100
101
;; Don't want log1p to use the x87 instruction.
#-(or hpux (and x86 (not sse2)))
(def-math-rtn "log1p" 1)

;; These are needed for use by byte-compiled files.  But don't use
;; these with sse2 since we don't support using the x87 instructions
;; here.
#+(and x86 (not sse2))
ram's avatar
ram committed
102
(progn
103
  #+nil
ram's avatar
ram committed
104
  (defun %sin (x)
105
    (declare (double-float x)
ram's avatar
ram committed
106
107
108
109
110
111
	     (values double-float))
    (%sin x))
  (defun %sin-quick (x)
    (declare (double-float x)
	     (values double-float))
    (%sin-quick x))
112
  #+nil
ram's avatar
ram committed
113
  (defun %cos (x)
114
    (declare (double-float x)
ram's avatar
ram committed
115
116
117
118
119
120
	     (values double-float))
    (%cos x))
  (defun %cos-quick (x)
    (declare (double-float x)
	     (values double-float))
    (%cos-quick x))
121
  #+nil
ram's avatar
ram committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
  (defun %tan (x)
    (declare (double-float x)
	     (values double-float))
    (%tan x))
  (defun %tan-quick (x)
    (declare (double-float x)
	     (values double-float))
    (%tan-quick x))
  (defun %atan (x)
    (declare (double-float x)
	     (values double-float))
    (%atan x))
  (defun %atan2 (x y)
    (declare (double-float x y)
	     (values double-float))
    (%atan2 x y))
  (defun %exp (x)
139
    (declare (double-float x)
ram's avatar
ram committed
140
141
142
	     (values double-float))
    (%exp x))
  (defun %log (x)
143
    (declare (double-float x)
ram's avatar
ram committed
144
145
146
	     (values double-float))
    (%log x))
  (defun %log10 (x)
147
    (declare (double-float x)
ram's avatar
ram committed
148
149
150
151
152
	     (values double-float))
    (%log10 x))
  #+nil ;; notyet
  (defun %pow (x y)
    (declare (type (double-float 0d0) x)
153
	     (double-float y)
ram's avatar
ram committed
154
155
156
	     (values (double-float 0d0)))
    (%pow x y))
  (defun %sqrt (x)
157
    (declare (double-float x)
ram's avatar
ram committed
158
159
160
	     (values double-float))
    (%sqrt x))
  (defun %scalbn (f ex)
161
    (declare (double-float f)
ram's avatar
ram committed
162
163
164
165
166
167
168
169
	     (type (signed-byte 32) ex)
	     (values double-float))
    (%scalbn f ex))
  (defun %scalb (f ex)
    (declare (double-float f ex)
	     (values double-float))
    (%scalb f ex))
  (defun %logb (x)
170
171
172
    (declare (double-float x)
	     (values double-float))
    (%logb x))
173
174
175
176
  (defun %log1p (x)
    (declare (double-float x)
	     (values double-float))
    (%log1p x))
ram's avatar
ram committed
177
  ) ; progn
wlott's avatar
wlott committed
178

179
180
181

;; As above for x86.  It also seems to be needed to handle
;; constant-folding in the compiler.
182
#+(or sparc (and x86 sse2))
183
184
185
186
187
188
189
(progn
  (defun %sqrt (x)
    (declare (double-float x)
	     (values double-float))
    (%sqrt x))
  )

190
191
192
193
194
195
196
197
198
199
200
;;; The standard libm routines for sin, cos, and tan on x86 (Linux,
;;; 32-bit.  64-bit is apparently ok) and ppc are not very accurate
;;; for large arguments when compared to sparc (and maxima).  This is
;;; basically caused by the fact that those libraries do not do an
;;; accurate argument reduction.  The following functions use some
;;; routines Sun's free fdlibm library to do accurate reduction.  Then
;;; we call the standard C functions (or vops for x86) on the reduced
;;; argument.  This produces much more accurate values.
;;;
;;; You can test this by computing (cos (scale-float 1d0 120)).  The
;;; true answer is -0.9258790228548379d0.
201

202
#+(or ppc x86)
203
(progn
204
(declaim (inline %%ieee754-rem-pi/2))
205
206
207
208
209
;; Basic argument reduction routine.  It returns two values: n and y
;; such that (n + 8*k)*pi/2+y = x where |y|<pi/4 and n indicates in
;; which octant the arg lies.  Y is actually computed in two parts,
;; y[0] and y[1] such that the sum is y, for accuracy.

210
(alien:def-alien-routine ("__ieee754_rem_pio2" %%ieee754-rem-pi/2) c-call:int
211
212
  (x double-float)
  (y (* double-float)))
213
214
215
216
217
218
219
220
221
222
223

;; Same as above, but instead of needing to pass an array in, the
;; output array is broken up into two output values instead.  This is
;; easier for the user, and we don't have to wrap calls with
;; without-gcing.
(declaim (inline %ieee754-rem-pi/2))
(alien:def-alien-routine ("ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
  (x double-float)
  (y0 double-float :out)
  (y1 double-float :out))

224
)
225

226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
;; If the C library is accurate, use %trig as the Lisp name.
#-(or ppc (and sse2 (not darwin)))
(progn
(declaim (inline %sin %cos %tan))
(macrolet ((frob (alien-name lisp-name)
	     `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
		(x double-float))))
  (frob "sin" %sin)
  (frob "cos" %cos)
  (frob "tan" %tan))
)

;; Make %%trig be the C library routines that don't do accurate
;; reduction.  This is for PPC and for any SSE2 build except on
;; Darwin. Darwin has accurate C library routines.
#+(or ppc (and sse2 (not darwin)))
242
243
244
245
246
247
248
249
(progn
(declaim (inline %%sin %%cos %%tan))
(macrolet ((frob (alien-name lisp-name)
	     `(alien:def-alien-routine (,alien-name ,lisp-name) double-float
		(x double-float))))
  (frob "sin" %%sin)
  (frob "cos" %%cos)
  (frob "tan" %%tan))
250
251
)

252
253
254
255
;; When the C library is not accurate, define %trig to do accurate
;; argument reduction and call the appropriate C function on the
;; reduced arg.  For x87, we can use the x87 FPU trig instructions.
#+(or ppc (and x86 (not darwin)))
256
257
(macrolet
    ((frob (sin cos tan)
258
       `(progn
259
	  ;; In all of the routines below, we just compute the sum of
rtoy's avatar
rtoy committed
260
261
262
263
264
265
266
267
268
	  ;; y0 and y1 and use that as the (reduced) argument for the
	  ;; trig functions.  This is slightly less accurate than what
	  ;; fdlibm does, which calls special functions using y0 and
	  ;; y1 separately, for greater accuracy.  This isn't
	  ;; implemented, and some spot checks indicate that what we
	  ;; have here is accurate.
	  ;;
	  ;; For x86 with an fsin/fcos/fptan instruction, the pi/4 is
	  ;; probably too restrictive.
269
270
271
272
273
	  (defun %sin (x)
	    (declare (double-float x))
	    (if (< (abs x) (/ pi 4))
		(,sin x)
		;; Argument reduction needed
274
275
276
277
278
279
280
281
		(multiple-value-bind (n y0 y1)
		    (%ieee754-rem-pi/2 x)
		  (let ((reduced (+ y0 y1)))
		    (case (logand n 3)
		      (0 (,sin reduced))
		      (1 (,cos reduced))
		      (2 (- (,sin reduced)))
		      (3 (- (,cos reduced))))))))
282
283
284
285
286
	  (defun %cos (x)
	    (declare (double-float x))
	    (if (< (abs x) (/ pi 4))
		(,cos x)
		;; Argument reduction needed
287
288
289
290
291
292
293
294
		(multiple-value-bind (n y0 y1)
		    (%ieee754-rem-pi/2 x)
		  (let ((reduced (+ y0 y1)))
		    (case (logand n 3)
		      (0 (,cos reduced))
		      (1 (- (,sin reduced)))
		      (2 (- (,cos reduced)))
		      (3 (,sin reduced)))))))
295
296
297
298
299
	  (defun %tan (x)
	    (declare (double-float x))
	    (if (< (abs x) (/ pi 4))
		(,tan x)
		;; Argument reduction needed
300
301
302
303
304
305
		(multiple-value-bind (n y0 y1)
		    (%ieee754-rem-pi/2 x)
		  (let ((reduced (+ y0 y1)))
		    (if (evenp n)
			(,tan reduced)
			(- (/ (,tan reduced)))))))))))
306
307
  ;; Don't want %sin-quick and friends with sse2.
  #+(and x86 (not sse2))
308
  (frob %sin-quick %cos-quick %tan-quick)
309
  #+(or ppc sse2)
310
311
  (frob %%sin %%cos %%tan))

312
313
314
315
316
;; Linux and sparc have a sincos function in the C library. Use it.
;; But on linux we need to do pi reduction ourselves because the C
;; library doesn't do accurate reduction.  Sparc does accurate pi
;; reduction, so we don't need to do it ourselves.
#+(or (and linux x86) sparc)
317
318
319
320
321
322
323
(progn
(declaim (inline %%sincos))
(export '%%sincos)
(alien:def-alien-routine ("sincos" %%sincos) c-call:void
  (x double-float)
  (sin double-float :out)
  (cos double-float :out))
324

325
#+(and linux x86)
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
(defun %sincos (theta)
  (declare (double-float theta))
  ;; Accurately reduce theta.
  (multiple-value-bind (n y0 y1)
      (%ieee754-rem-pi/2 theta)
    (multiple-value-bind (ignore s c)
	(%%sincos y0)
      (declare (ignore ignore))
      ;; Figure out which quadrant to use, and finish out the
      ;; computation using y1. This is done by using a 1st-order
      ;; Taylor expansion about y0.
      (flet ((sin2 (s c y)
	       ;; sin(x+y) = sin(x) + cos(x)*y
	       (+ s (* c y)))
	     (cos2 (s c y)
	       ;; cos(x+y) = cos(x) - sin(x)*y
	       (- c (* s y))))
	(case (logand n 3)
	  (0
345
	   (values (sin2 s c y1)
346
		   (cos2 s c y1)))
347
348
	  (1
	   (values (cos2 s c y1)
349
		   (- (sin2 s c y1))))
350
351
352
	  (2
	   (values (- (sin2 s c y1))
		   (- (cos2 s c y1))))
353
	  (3
354
355
	   (values (- (cos2 s c y1))
		   (sin2 s c y1))))))))
356
357
358
359
360
361
#+sparc
(declaim (inline %sinccos))
#+sparc
(defun %sincos (theta)
  (multiple-value-bind (ignore s c)
      (%%sincos theta)
362
363
    (declare (ignore ignore))
    (values s c)))
364
)
365
366


wlott's avatar
wlott committed
367
368
369
370

;;;; Power functions.

(defun exp (number)
371
  "Return e raised to the power NUMBER."
wlott's avatar
wlott committed
372
373
374
375
376
377
378
379
380
381
  (number-dispatch ((number number))
    (handle-reals %exp number)
    ((complex)
     (* (exp (realpart number))
	(cis (imagpart number))))))

;;; INTEXP -- Handle the rational base, integer power case.

(defparameter *intexp-maximum-exponent* 10000)

382
383
384
385
(define-condition intexp-limit-error (error)
  ((base :initarg :base :reader intexp-base)
   (power :initarg :power :reader intexp-power))
  (:report (lambda (condition stream)
386
	     (format stream (intl:gettext "The absolute value of ~S exceeds limit ~S.")
387
388
389
		     (intexp-power condition)
		     *intexp-maximum-exponent*))))

390
391
392
393
394
;;; This function precisely calculates base raised to an integral power.  It
;;; separates the cases by the sign of power, for efficiency reasons, as powers
;;; can be calculated more efficiently if power is a positive integer.  Values
;;; of power are calculated as positive integers, and inverted if negative.
;;;
wlott's avatar
wlott committed
395
(defun intexp (base power)
396
397
398
  ;; Handle the special case of 1^power and (-1)^power.  Maxima
  ;; sometimes does this, and there's no need to cause a continuable
  ;; error in this case.
399
400
  (when (eql base 1)
    (return-from intexp base))
401
402
  (when (eql base -1)
    (return-from intexp (if (oddp power) -1 1)))
403
  
wlott's avatar
wlott committed
404
  (when (> (abs power) *intexp-maximum-exponent*)
405
406
407
408
409
410
411
    ;; Allow user the option to continue with calculation, possibly
    ;; increasing the limit to the given power.
    (restart-case
	(error 'intexp-limit-error
	       :base base
	       :power power)
      (continue ()
412
	:report (lambda (stream)
413
		  (write-string (intl:gettext "Continue with calculation") stream)))
414
      (new-limit ()
415
	:report (lambda (stream)
416
		  (write-string (intl:gettext "Continue with calculation, update limit") stream))
417
	(setq *intexp-maximum-exponent* (abs power)))))
wlott's avatar
wlott committed
418
419
420
421
422
423
424
425
426
427
428
429
430
  (cond ((minusp power)
	 (/ (intexp base (- power))))
	((eql base 2)
	 (ash 1 power))
	(t
	 (do ((nextn (ash power -1) (ash power -1))
	      (total (if (oddp power) base 1)
		     (if (oddp power) (* base total) total)))
	     ((zerop nextn) total)
	   (setq base (* base base))
	   (setq power nextn)))))


431
432
433
434
435
436
437
438
;;; EXPT  --  Public
;;;
;;;    If an integer power of a rational, use INTEXP above.  Otherwise, do
;;; floating point stuff.  If both args are real, we try %POW right off,
;;; assuming it will return 0 if the result may be complex.  If so, we call
;;; COMPLEX-POW which directly computes the complex result.  We also separate
;;; the complex-real and real-complex cases from the general complex case.
;;;
wlott's avatar
wlott committed
439
(defun expt (base power)
440
  "Returns BASE raised to the POWER."
wlott's avatar
wlott committed
441
  (if (zerop power)
442
443
444
445
446
447
448
449
450
451
452
      ;; CLHS says that if the power is 0, the result is 1, subject to
      ;; numeric contagion.  But what happens if base is infinity or
      ;; NaN?  Do we silently return 1?  For now, I think we should
      ;; signal an error if the FP modes say so.
      (let ((result (1+ (* base power))))
	;; If we get an NaN here, that means base*power above didn't
	;; produce 0 and FP traps were disabled, so we handle that
	;; here.  Should this be a continuable restart?
	(if (and (floatp result) (float-nan-p result))
	    (float 1 result)
	    result))
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
    (labels (;; determine if the double float is an integer.
	     ;;  0 - not an integer
	     ;;  1 - an odd int
	     ;;  2 - an even int
	     (isint (ihi lo)
	       (declare (type (unsigned-byte 31) ihi)
			(type (unsigned-byte 32) lo)
			(optimize (speed 3) (safety 0)))
	       (let ((isint 0))
		 (declare (type fixnum isint))
		 (cond ((>= ihi #x43400000)	; exponent >= 53
			(setq isint 2))
		       ((>= ihi #x3ff00000)
			(let ((k (- (ash ihi -20) #x3ff)))	; exponent
			  (declare (type (mod 53) k))
			  (cond ((> k 20)
				 (let* ((shift (- 52 k))
					(j (logand (ash lo (- shift))))
					(j2 (ash j shift)))
				   (declare (type (mod 32) shift)
					    (type (unsigned-byte 32) j j2))
				   (when (= j2 lo)
				     (setq isint (- 2 (logand j 1))))))
				((= lo 0)
				 (let* ((shift (- 20 k))
					(j (ash ihi (- shift)))
					(j2 (ash j shift)))
				   (declare (type (mod 32) shift)
					    (type (unsigned-byte 31) j j2))
				   (when (= j2 ihi)
				     (setq isint (- 2 (logand j 1))))))))))
		 isint))
	     (real-expt (x y rtype)
	       (let ((x (coerce x 'double-float))
		     (y (coerce y 'double-float)))
		 (declare (double-float x y))
		 (let* ((x-hi (kernel:double-float-high-bits x))
			(x-lo (kernel:double-float-low-bits x))
			(x-ihi (logand x-hi #x7fffffff))
			(y-hi (kernel:double-float-high-bits y))
			(y-lo (kernel:double-float-low-bits y))
			(y-ihi (logand y-hi #x7fffffff)))
		   (declare (type (signed-byte 32) x-hi y-hi)
			    (type (unsigned-byte 31) x-ihi y-ihi)
			    (type (unsigned-byte 32) x-lo y-lo))
		   ;; y==zero: x**0 = 1
		   (when (zerop (logior y-ihi y-lo))
		     (return-from real-expt (coerce 1d0 rtype)))
		   ;; +-NaN return x+y
		   (when (or (> x-ihi #x7ff00000)
			     (and (= x-ihi #x7ff00000) (/= x-lo 0))
			     (> y-ihi #x7ff00000)
			     (and (= y-ihi #x7ff00000) (/= y-lo 0)))
		     (return-from real-expt (coerce (+ x y) rtype)))
		   (let ((yisint (if (< x-hi 0) (isint y-ihi y-lo) 0)))
		     (declare (type fixnum yisint))
		     ;; special value of y
		     (when (and (zerop y-lo) (= y-ihi #x7ff00000))
		       ;; y is +-inf
		       (return-from real-expt
			 (cond ((and (= x-ihi #x3ff00000) (zerop x-lo))
				;; +-1**inf is NaN
				(coerce (- y y) rtype))
			       ((>= x-ihi #x3ff00000)
				;; (|x|>1)**+-inf = inf,0
				(if (>= y-hi 0)
				    (coerce y rtype)
				    (coerce 0 rtype)))
			       (t
				;; (|x|<1)**-,+inf = inf,0
				(if (< y-hi 0)
				    (coerce (- y) rtype)
				    (coerce 0 rtype))))))

		     (let ((abs-x (abs x)))
		       (declare (double-float abs-x))
		       ;; special value of x
		       (when (and (zerop x-lo)
				  (or (= x-ihi #x7ff00000) (zerop x-ihi)
				      (= x-ihi #x3ff00000)))
			 ;; x is +-0,+-inf,+-1
			 (let ((z (if (< y-hi 0)
				      (/ 1 abs-x)	; z = (1/|x|)
				      abs-x)))
			   (declare (double-float z))
			   (when (< x-hi 0)
			     (cond ((and (= x-ihi #x3ff00000) (zerop yisint))
				    ;; (-1)**non-int
				    (let ((y*pi (* y pi)))
				      (declare (double-float y*pi))
				      (return-from real-expt
544
				        (complex
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
					 (coerce (%cos y*pi) rtype)
					 (coerce (%sin y*pi) rtype)))))
				   ((= yisint 1)
				    ;; (x<0)**odd = -(|x|**odd)
				    (setq z (- z)))))
			   (return-from real-expt (coerce z rtype))))
		       
		       (if (>= x-hi 0)
			   ;; x>0
			   (coerce (kernel::%pow x y) rtype)
			   ;; x<0
			   (let ((pow (kernel::%pow abs-x y)))
			     (declare (double-float pow))
			     (case yisint
			       (1 ; Odd
				(coerce (* -1d0 pow) rtype))
			       (2 ; Even
				(coerce pow rtype))
			       (t ; Non-integer
				(let ((y*pi (* y pi)))
				  (declare (double-float y*pi))
566
				  (complex
567
				   (coerce (* pow (%cos y*pi)) rtype)
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
				   (coerce (* pow (%sin y*pi)) rtype))))))))))))
	     (expt-xfrm (b p)
	       ;; Apply the same transformation as in the deftransform
	       ;; for expt in compiler/srctran.lisp.  Only call this
	       ;; if B is more contagious than P.  Otherwise, the type
	       ;; of the result will be wrong which will confuse the
	       ;; compiler!  Return NIL if the transform can't be
	       ;; applied.
	       (cond
		 ((= p 2) (* b b))
		 ((= p -2) (/ (* b b)))
		 ((= p 3) (* b b b))
		 ((= p -3) (/ (* b b b)))
		 ((= p 1/2) (sqrt b))
		 ((= p -1/2) (/ (sqrt b)))
		 (t nil))))
rtoy's avatar
rtoy committed
584
585
586
587
588
589
590
      ;; This is really messy and should be cleaned up.  The easiest
      ;; way to see if we're doing what we should is the macroexpand
      ;; the number-dispatch and check each branch.
      ;;
      ;; We try to apply the rule of float precision contagion (CLHS
      ;; 12.1.4.4): the result has the same precision has the most
      ;; precise argument.
591
      (number-dispatch ((base number) (power number))
rtoy's avatar
rtoy committed
592
593
        (((foreach fixnum (or bignum ratio) (complex rational))
	  integer)
594
	 (intexp base power))
rtoy's avatar
rtoy committed
595
596
	(((foreach single-float double-float)
	  rational)
Raymond Toy's avatar
Raymond Toy committed
597
598
	 (or (expt-xfrm base power)
	     (real-expt base power '(dispatch-type base))))
599
600
	(((foreach fixnum (or bignum ratio) single-float)
	  (foreach ratio single-float))
Raymond Toy's avatar
Raymond Toy committed
601
602
	 (or (expt-xfrm (coerce base 'single-float) power)
	     (real-expt base power 'single-float)))
603
604
	(((foreach fixnum (or bignum ratio) single-float double-float)
	  double-float)
Raymond Toy's avatar
Raymond Toy committed
605
606
	 (or (expt-xfrm (coerce base 'double-float) power)
	     (real-expt base power 'double-float)))
607
	((double-float single-float)
Raymond Toy's avatar
Raymond Toy committed
608
609
	 (or (expt-xfrm (coerce base 'double-float) power)
	     (real-expt base power 'double-float)))
610
	#+double-double
rtoy's avatar
rtoy committed
611
612
	(((foreach fixnum (or bignum ratio) single-float double-float
		   double-double-float)
613
	  double-double-float)
Raymond Toy's avatar
Raymond Toy committed
614
615
	 (or (expt-xfrm (coerce base 'double-double-float) power)
	     (dd-%pow (coerce base 'double-double-float) power)))
616
617
618
	#+double-double
	((double-double-float
	  (foreach fixnum (or bignum ratio) single-float double-float))
Raymond Toy's avatar
Raymond Toy committed
619
620
	 (or (expt-xfrm base power)
	     (dd-%pow base (coerce power 'double-double-float))))
rtoy's avatar
rtoy committed
621
622
623
	(((foreach (complex rational) (complex single-float) (complex double-float)
		   #+double-double (complex double-double-float))
	  rational)
624
625
626
	 (or (expt-xfrm base power)
	     (* (expt (abs base) power)
		(cis (* power (phase base))))))
rtoy's avatar
rtoy committed
627
628
	#+double-double
	((double-double-float
629
	  complex)
rtoy's avatar
rtoy committed
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
	     (exp (* power (* (log2 base 1w0) (log 2w0))))))
	(((foreach fixnum (or bignum ratio) single-float double-float)
	  (foreach (complex double-float)))
	 ;; Result should have double-float accuracy.  Use log2 in
	 ;; case the base won't fit in a double-float.
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
	     (exp (* power (* (log2 base) (log 2d0))))))
	((double-float
	  (foreach (complex rational) (complex single-float)))
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
	     (exp (* power (log base)))))
	#+double-double
	(((foreach fixnum (or bignum ratio) single-float double-float)
	  (foreach (complex double-double-float)))
	 ;; Result should have double-double-float accuracy.  Use log2
	 ;; in case the base won't fit in a double-float.
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
	     (exp (* power (* (log2 base 1w0) (log 2w0))))))
	(((foreach fixnum (or bignum ratio) single-float)
	  (foreach (complex single-float)))
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
	     (exp (* power (log base)))))
	(((foreach (complex rational) (complex single-float))
	  (foreach single-float (complex single-float)))
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
662
663
	     (or (expt-xfrm (coerce base '(complex single-float)) power)
		 (exp (* power (log base))))))
rtoy's avatar
rtoy committed
664
665
666
667
	(((foreach (complex rational) (complex single-float))
	  (foreach double-float (complex double-float)))
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
668
669
670
	     (or (expt-xfrm (coerce base '(complex double-float))
			    power)
		 (exp (* power (log (coerce base '(complex double-float))))))))
rtoy's avatar
rtoy committed
671
672
673
674
675
	#+double-double
	(((foreach (complex rational) (complex single-float))
	  (foreach double-double-float (complex double-double-float)))
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
676
677
678
	     (or (expt-xfrm (coerce base '(complex double-double-float))
			    power)
		 (exp (* power (log (coerce base '(complex double-double-float))))))))
rtoy's avatar
rtoy committed
679
	(((foreach (complex double-float))
680
681
	  (foreach single-float double-float
		   (complex single-float) (complex double-float)))
682
683
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
684
685
	     (or (expt-xfrm base power)
		 (exp (* power (log base))))))
rtoy's avatar
rtoy committed
686
687
688
689
690
	#+double-double
	(((foreach (complex double-float))
	  (foreach double-double-float (complex double-double-float)))
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
691
692
693
	     (or (expt-xfrm (coerce base '(complex double-double-float))
			    power)
		 (exp (* power (log (coerce base '(complex double-double-float))))))))
rtoy's avatar
rtoy committed
694
695
696
	#+double-double
	(((foreach (complex double-double-float))
	  (foreach float (complex float)))
697
698
	 (if (and (zerop base) (plusp (realpart power)))
	     (* base power)
699
700
	     (or (expt-xfrm base power)
		 (exp (* power (log base))))))))))
wlott's avatar
wlott committed
701

702
703
704
705
706
707
708
709
;; Log base 2 of a real number.  The result is a either a double-float
;; or double-double-float number (real or complex, as appropriate),
;; depending on the type of FLOAT-TYPE.
(defun log2 (x &optional (float-type 1d0))
  (labels ((log-of-2 (f)
	     ;; log(2), with the precision specified by the type of F
	     (number-dispatch ((f real))
	       ((double-float)
710
		#.(log 2d0))
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
	       #+double-double
	       ((double-double-float)
		#.(log 2w0))))
	   (log-2-pi (f)
	     ;; log(pi), with the precision specified by the type of F
	     (number-dispatch ((f real))
	       ((double-float)
		#.(/ pi (log 2d0)))
	       #+double-double
	       ((double-double-float)
		#.(/ dd-pi (log 2w0)))))
	   (log1p (x)
	     ;; log(1+x), with the precision specified by the type of
	     ;; X
	     (number-dispatch ((x real))
	       (((foreach single-float double-float))
		(%log1p (float x 1d0)))
	       #+double-double
	       ((double-double-float)
		(dd-%log1p x))))
	   (log2-bignum (bignum)
732
733
734
735
736
737
738
739
740
741
	     ;; Write x = 2^n*f where 1/2 < f <= 1.  Then log2(x) = n
	     ;; + log2(f).
	     ;;
	     ;; So we grab the top few bits of x and scale that
	     ;; appropriately, take the log of it and add it to n.
	     ;;
	     ;; Return n and log2(f) separately.
	     (if (minusp bignum)
		 (multiple-value-bind (n frac)
		     (log2-bignum (abs bignum))
742
743
744
745
746
747
748
749
750
		   (values n (complex frac (log-2-pi float-type))))
		 (let ((n (integer-length bignum))
		       (float-bits (float-digits float-type)))
		   (if (< n float-bits)
		       (values 0 (log (float bignum float-type)
				      (float 2 float-type)))
		       (let ((exp (min float-bits n))
			     (f (ldb (byte float-bits
					   (max 0 (- n float-bits)))
751
				     bignum)))
752
753
			 (values n (log (scale-float (float f float-type) (- exp))
					(float 2 float-type)))))))))
754
755
    (etypecase x
      (float
756
       (/ (log (float x float-type)) (log-of-2 float-type)))
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
      (ratio
       (let ((top (numerator x))
	     (bot (denominator x)))
	 ;; If the number of bits in the numerator and
	 ;; denominator are different, just use the fact
	 ;; log(x/y) = log(x) - log(y).  But to preserve
	 ;; accuracy, we actually do
	 ;; (log2(x)-log2(y))/log2(e)).
	 ;;
	 ;; However, if the numerator and denominator have the
	 ;; same number of bits, implying the quotient is near
	 ;; one, we use log1p(x) = log(1+x). Since the number is
	 ;; rational, we don't lose precision subtracting 1 from
	 ;; it, and converting it to double-float is accurate.
	 (if (= (integer-length top)
		(integer-length bot))
773
774
	     (/ (log1p (float (- x 1) float-type))
		(log-of-2 float-type))
775
776
777
778
779
780
781
782
783
784
	     (multiple-value-bind (top-n top-frac)
		 (log2-bignum top)
	       (multiple-value-bind (bot-n bot-frac)
		   (log2-bignum bot)
		 (+ (- top-n bot-n)
		    (- top-frac bot-frac)))))))
      (integer
       (multiple-value-bind (n frac)
	   (log2-bignum x)
	 (+ n frac))))))
785

wlott's avatar
wlott committed
786
(defun log (number &optional (base nil base-p))
787
  "Return the logarithm of NUMBER in the base BASE, which defaults to e."
wlott's avatar
wlott committed
788
  (if base-p
789
790
791
      (cond ((zerop base)
	     ;; ANSI spec
	     base)
792
793
794
795
796
797
798
799
800
801
802
803
804
805
	    ((and (realp number) (realp base))
	     ;; CLHS 12.1.4.1 says
	     ;;
	     ;;   When rationals and floats are combined by a
	     ;;   numerical function, the rational is first converted
	     ;;   to a float of the same format.
	     ;;
	     ;; So assume this applies to floats as well convert all
	     ;; numbers to the largest float format before computing
	     ;; the log.
	     ;;
	     ;; This makes (log 17 10.0) = (log 17.0 10) and so on.
	     (number-dispatch ((number real) (base real))
	       ((double-float
806
807
808
809
810
811
812
813
814
815
		 (foreach double-float single-float))
		(/ (log2 number) (log2 base)))
	       (((foreach fixnum bignum ratio)
		 (foreach fixnum bignum ratio single-float))
		(let* ((result (/ (log2 number) (log2 base))))
		  ;; Figure out the right result type
		  (if (realp result)
		      (coerce result 'single-float)
		      (coerce result '(complex single-float)))))
	       (((foreach fixnum bignum ratio)
816
		 double-float)
817
818
819
820
821
822
823
824
825
826
827
828
		(/ (log2 number) (log2 base)))
	       ((single-float
		 (foreach fixnum bignum ratio))
		(let* ((result (/ (log2 number) (log2 base))))
		  ;; Figure out the right result type
		  (if (realp result)
		      (coerce result 'single-float)
		      (coerce result '(complex single-float)))))
	       ((double-float
		 (foreach fixnum bignum ratio))
		(/ (log2 number) (log2 base)))
	       ((single-float double-float)
829
		(/ (log (coerce number 'double-float)) (log base)))
830
831
	       #+double-double
	       ((double-double-float
832
		 (foreach fixnum bignum ratio))
833
		(/ (log2 number 1w0) (log2 base 1w0)))
834
835
836
	       #+double-double
	       ((double-double-float
		 (foreach double-double-float double-float single-float))
837
838
		(/ (log number) (log (coerce base 'double-double-float))))
	       #+double-double
839
840
	       (((foreach fixnum bignum ratio)
		 double-double-float)
841
		(/ (log2 number 1w0) (log2 base 1w0)))
842
843
	       #+double-double
	       (((foreach double-float single-float)
844
845
		 double-double-float)
		(/ (log (coerce number 'double-double-float)) (log base)))
846
847
	       (((foreach single-float)
		 (foreach single-float))
848
849
850
851
852
853
854
855
856
857
858
859
860
		;; Converting everything to double-float helps the
		;; cases like (log 17 10) = (/ (log 17) (log 10)).
		;; This is usually handled above, but if we compute (/
		;; (log 17) (log 10)), we get a slightly different
		;; answer due to roundoff.  This makes it a bit more
		;; consistent.
		;;
		;; FIXME: This probably needs more work.
		(let ((result (/ (log (float number 1d0))
				 (log (float base 1d0)))))
		  (if (realp result)
		      (coerce result 'single-float)
		      (coerce result '(complex single-float)))))))
861
	    (t
862
	     ;; FIXME:  This probably needs some work as well.
863
	     (/ (log number) (log base))))
wlott's avatar
wlott committed
864
      (number-dispatch ((number number))
865
	(((foreach fixnum bignum))
866
	 (if (minusp number)
867
868
	     (complex (coerce (log (- number)) 'single-float)
		      (coerce pi 'single-float))
869
	     (coerce (/ (log2 number) #.(log (exp 1d0) 2d0)) 'single-float)))
870
871
872
873
	((ratio)
	 (if (minusp number)
	     (complex (coerce (log (- number)) 'single-float)
		      (coerce pi 'single-float))
874
875
876
877
878
879
	     ;; What happens when the ratio is close to 1?  We need to
	     ;; be careful to preserve accuracy.
	     (let ((top (numerator number))
		   (bot (denominator number)))
	       ;; If the number of bits in the numerator and
	       ;; denominator are different, just use the fact
880
881
882
883
884
885
886
887
888
	       ;; log(x/y) = log(x) - log(y).  But to preserve
	       ;; accuracy, we actually do
	       ;; (log2(x)-log2(y))/log2(e)).
	       ;;
	       ;; However, if the numerator and denominator have the
	       ;; same number of bits, implying the quotient is near
	       ;; one, we use log1p(x) = log(1+x). Since the number is
	       ;; rational, we don't lose precision subtracting 1 from
	       ;; it, and converting it to double-float is accurate.
889
	       (if (= (integer-length top)
890
		      (integer-length bot))
891
892
		   (coerce (%log1p (coerce (- number 1) 'double-float))
			   'single-float)
893
894
895
		   (coerce (/ (- (log2 top) (log2 bot))
			      #.(log (exp 1d0) 2d0))
			   'single-float)))))
ram's avatar
ram committed
896
897
898
899
900
901
902
903
904
	(((foreach single-float double-float))
	 ;; Is (log -0) -infinity (libm.a) or -infinity + i*pi (Kahan)?
	 ;; Since this doesn't seem to be an implementation issue
	 ;; I (pw) take the Kahan result.
	 (if (< (float-sign number)
		(coerce 0 '(dispatch-type number)))
	     (complex (log (- number)) (coerce pi '(dispatch-type number)))
	     (coerce (%log (coerce number 'double-float))
		     '(dispatch-type number))))
905
906
907
908
909
910
	#+double-double
	((double-double-float)
	 (let ((hi (kernel:double-double-hi number)))
	   (if (< (float-sign hi) 0d0)
	       (complex (dd-%log (- number)) dd-pi)
	       (dd-%log number))))
ram's avatar
ram committed
911
912
	((complex)
	 (complex-log number)))))
wlott's avatar
wlott committed
913
914

(defun sqrt (number)
915
  "Return the square root of NUMBER."
wlott's avatar
wlott committed
916
  (number-dispatch ((number number))
ram's avatar
ram committed
917
    (((foreach fixnum bignum ratio))
918
     (if (minusp number)
ram's avatar
ram committed
919
	 (complex-sqrt number)
920
	 (coerce (%sqrt (coerce number 'double-float)) 'single-float)))
ram's avatar
ram committed
921
    (((foreach single-float double-float))
dtc's avatar
dtc committed
922
     (if (minusp number)
ram's avatar
ram committed
923
924
925
	 (complex-sqrt number)
	 (coerce (%sqrt (coerce number 'double-float))
		 '(dispatch-type number))))
926
927
928
929
930
931
    #+double-double
    ((double-double-float)
     (if (minusp number)
	 (dd-complex-sqrt number)
	 (multiple-value-bind (hi lo)
	     (c::sqrt-dd (kernel:double-double-hi number) (kernel:double-double-lo number))
932
	   (kernel:%make-double-double-float hi lo))))
933
934
    ((complex)
     (complex-sqrt number))))
wlott's avatar
wlott committed
935
936
937
938


;;;; Trigonometic and Related Functions

939
(defun abs (number)
940
  "Returns the absolute value of the number."
941
  (number-dispatch ((number number))
942
943
    (((foreach single-float double-float fixnum rational
	       #+double-double double-double-float))
944
945
946
947
948
949
950
951
952
953
954
955
     (abs number))
    ((complex)
     (let ((rx (realpart number))
	   (ix (imagpart number)))
       (etypecase rx
	 (rational
	  (sqrt (+ (* rx rx) (* ix ix))))
	 (single-float
	  (coerce (%hypot (coerce rx 'double-float)
			  (coerce ix 'double-float))
		  'single-float))
	 (double-float
956
957
958
	  (%hypot rx ix))
	 #+double-double
	 (double-double-float
959
960
961
	  (multiple-value-bind (abs^2 scale)
	      (dd-cssqs number)
	    (scale-float (sqrt abs^2) scale))))))))
wlott's avatar
wlott committed
962
963

(defun phase (number)
964
  "Returns the angle part of the polar representation of a complex number.
965
  For complex numbers, this is (atan (imagpart number) (realpart number)).
wlott's avatar
wlott committed
966
967
  For non-complex positive numbers, this is 0.  For non-complex negative
  numbers this is PI."
968
  (etypecase number
ram's avatar
ram committed
969
    (rational
970
971
972
     (if (minusp number)
	 (coerce pi 'single-float)
	 0.0f0))
ram's avatar
ram committed
973
974
975
976
    (single-float
     (if (minusp (float-sign number))
	 (coerce pi 'single-float)
	 0.0f0))
977
    (double-float
ram's avatar
ram committed
978
     (if (minusp (float-sign number))
979
980
	 (coerce pi 'double-float)
	 0.0d0))
981
982
983
984
985
    #+double-double
    (double-double-float
     (if (minusp (float-sign number))
	 dd-pi
	 0w0))
986
987
988
    (complex
     (atan (imagpart number) (realpart number)))))

wlott's avatar
wlott committed
989
990

(defun sin (number)  
991
  "Return the sine of NUMBER."
wlott's avatar
wlott committed
992
993
994
995
996
  (number-dispatch ((number number))
    (handle-reals %sin number)
    ((complex)
     (let ((x (realpart number))
	   (y (imagpart number)))
ram's avatar
ram committed
997
998
       (complex (* (sin x) (cosh y))
		(* (cos x) (sinh y)))))))
wlott's avatar
wlott committed
999
1000

(defun cos (number)
1001
  "Return the cosine of NUMBER."
wlott's avatar
wlott committed
1002
1003
1004
1005
1006
  (number-dispatch ((number number))
    (handle-reals %cos number)
    ((complex)
     (let ((x (realpart number))
	   (y (imagpart number)))
ram's avatar
ram committed
1007
1008
       (complex (* (cos x) (cosh y))
		(- (* (sin x) (sinh y))))))))
wlott's avatar
wlott committed
1009
1010

(defun tan (number)
1011
  "Return the tangent of NUMBER."
wlott's avatar
wlott committed
1012
1013
1014
  (number-dispatch ((number number))
    (handle-reals %tan number)
    ((complex)
ram's avatar
ram committed
1015
     (complex-tan number))))
wlott's avatar
wlott committed
1016

1017
(defun cis (theta)
1018
  "Return cos(Theta) + i sin(Theta), AKA exp(i Theta)."
1019
  (if (complexp theta)
1020
      (error (intl:gettext "Argument to CIS is complex: ~S") theta)
1021
      #-(or (and linux x86) sparc)
1022
      (complex (cos theta) (sin theta))
1023
      #+(or (and linux x86) sparc)
1024
1025
1026
1027
1028
      (number-dispatch ((theta real))
	((rational)
	 (let ((arg (coerce theta 'double-float)))
	   (multiple-value-bind (s c)
	       (%sincos arg)
1029
1030
	     (complex (coerce c 'single-float)
		      (coerce s 'single-float)))))
1031
1032
1033
	(((foreach single-float double-float))
	 (multiple-value-bind (s c)
	     (%sincos (coerce theta 'double-float))
1034
1035
	   (complex (coerce c '(dispatch-type theta))
		    (coerce s '(dispatch-type theta)))))
1036
1037
1038
	#+double-double
	((double-double-float)
	 (complex (cos theta) (sin theta))))))
wlott's avatar
wlott committed
1039
1040

(defun asin (number)
1041
  "Return the arc sine of NUMBER."
wlott's avatar
wlott committed
1042
  (number-dispatch ((number number))
1043
1044
1045
1046
1047
    ((rational)
     (if (or (> number 1) (< number -1))
	 (complex-asin number)
	 (coerce (%asin (coerce number 'double-float)) 'single-float)))
    (((foreach single-float double-float))
1048
1049
1050
     (if (or (float-nan-p number)
	     (and (<= number (coerce 1 '(dispatch-type number)))
		  (>= number (coerce -1 '(dispatch-type number)))))
1051
	 (coerce (%asin (coerce number 'double-float))
1052
		 '(dispatch-type number))
1053
1054
1055
1056
1057
1058
1059
1060
	 (complex-asin number)))
    #+double-double
    ((double-double-float)
     (if (or (float-nan-p number)
	     (and (<= number 1w0)
		  (>= number -1w0)))
	 (dd-%asin number)
	 (dd-complex-asin number)))
wlott's avatar
wlott committed
1061
    ((complex)
1062
     (complex-asin number))))
wlott's avatar
wlott committed
1063
1064

(defun acos (number)
1065
  "Return the arc cosine of NUMBER."
wlott's avatar
wlott committed
1066
  (number-dispatch ((number number))
1067
1068
1069
1070
1071
    ((rational)
     (if (or (> number 1) (< number -1))
	 (complex-acos number)
	 (coerce (%acos (coerce number 'double-float)) 'single-float)))
    (((foreach single-float double-float))
1072
1073
1074
     (if (or (float-nan-p number)
	     (and (<= number (coerce 1 '(dispatch-type number)))
		  (>= number (coerce -1 '(dispatch-type number)))))
1075
	 (coerce (%acos (coerce number 'double-float))
1076
		 '(dispatch-type number))
1077
1078
1079
1080
1081
1082
1083
1084
	 (complex-acos number)))
    #+double-double
    ((double-double-float)
     (if (or (float-nan-p number)
	     (and (<= number 1w0)
		  (>= number -1w0)))
	 (dd-%acos number)
	 (complex-acos number)))
wlott's avatar
wlott committed
1085
    ((complex)
1086
     (complex-acos number))))
wlott's avatar
wlott committed
1087
1088
1089


(defun atan (y &optional (x nil xp))
1090
  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
wlott's avatar
wlott committed
1091
  (if xp
1092
      (flet ((atan2 (y x)
1093
1094
1095
1096
1097
1098
1099
1100
1101
	       (declare (type double-float y x)
			(values double-float))
	       (if (zerop x)
		   (if (zerop y)
		       (if (plusp (float-sign x))
			   y
			   (float-sign y pi))
		       (float-sign y (/ pi 2)))
		   (%atan2 y x))))
1102
1103
	;; If X is given, both X and Y must be real numbers.
	(number-dispatch ((y real) (x real))
1104
1105
1106
1107
1108
1109
1110
1111
1112
	  ((double-float
	    (foreach double-float single-float fixnum bignum ratio))
	   (atan2 y (coerce x 'double-float)))
	  (((foreach single-float fixnum bignum ratio)
	    double-float)
	   (atan2 (coerce y 'double-float) x))
	  (((foreach single-float fixnum bignum ratio)
	    (foreach single-float fixnum bignum ratio))
	   (coerce (atan2 (coerce y 'double-float) (coerce x 'double-float))
1113
1114
1115
1116
1117
1118
1119
1120
1121
		   'single-float))
	  #+double-double
	  ((double-double-float
	    (foreach double-double-float double-float single-float fixnum bignum ratio))
	   (dd-%atan2 y (coerce x 'double-double-float)))
	  #+double-double
	  (((foreach double-float single-float fixnum bignum ratio)
	    double-double-float)
	   (dd-%atan2 (coerce y 'double-double-float) x))))
wlott's avatar
wlott committed
1122
1123
1124
      (number-dispatch ((y number))
	(handle-reals %atan y)
	((complex)
ram's avatar
ram committed
1125
	 (complex-atan y)))))
1126

ram's avatar
ram committed
1127
(defun sinh (number)
1128
  "Return the hyperbolic sine of NUMBER."
ram's avatar
ram committed
1129
1130
1131
1132
1133
1134
1135
1136
1137
  (number-dispatch ((number number))
    (handle-reals %sinh number)
    ((complex)
     (let ((x (realpart number))
	   (y (imagpart number)))
       (complex (* (sinh x) (cos y))
		(* (cosh x) (sin y)))))))

(defun cosh (number)
1138
  "Return the hyperbolic cosine of NUMBER."
ram's avatar
ram committed
1139
1140
1141
1142
1143
1144
1145
1146
  (number-dispatch ((number number))
    (handle-reals %cosh number)
    ((complex)
     (let ((x (realpart number))
	   (y (imagpart number)))
       (complex (* (cosh x) (cos y))
		(* (sinh x) (sin y)))))))

1147
(defun tanh (number)
1148
  "Return the hyperbolic tangent of NUMBER."
ram's avatar
ram committed
1149
1150
1151
1152
  (number-dispatch ((number number))
    (handle-reals %tanh number)
    ((complex)
     (complex-tanh number))))
1153
1154

(defun asinh (number)
1155
  "Return the hyperbolic arc sine of NUMBER."
ram's avatar
ram committed
1156
1157
1158
1159
  (number-dispatch ((number number))
    (handle-reals %asinh number)
    ((complex)
     (complex-asinh number))))
1160
1161

(defun acosh (number)
1162
  "Return the hyperbolic arc cosine of NUMBER."
ram's avatar
ram committed
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
  (number-dispatch ((number number))
    ((rational)
     ;; acosh is complex if number < 1
     (if (< number 1)
	 (complex-acosh number)
	 (coerce (%acosh (coerce number 'double-float)) 'single-float)))
    (((foreach single-float double-float))
     (if (< number (coerce 1 '(dispatch-type number)))
	 (complex-acosh number)
	 (coerce (%acosh (coerce number 'double-float))
		 '(dispatch-type number))))
1174
1175
1176
1177
1178
    #+double-double
    ((double-double-float)
     (if (< number 1w0)
	 (complex-acosh number)
	 (dd-%acosh number)))
ram's avatar
ram committed
1179
1180
    ((complex)
     (complex-acosh number))))
1181
1182

(defun atanh (number)
1183
  "Return the hyperbolic arc tangent of NUMBER."
ram's avatar
ram committed
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
  (number-dispatch ((number number))
    ((rational)
     ;; atanh is complex if |number| > 1
     (if (or (> number 1) (< number -1))
	 (complex-atanh number)
	 (coerce (%atanh (coerce number 'double-float)) 'single-float)))
    (((foreach single-float double-float))
     (if (or (> number (coerce 1 '(dispatch-type number)))
	     (< number (coerce -1 '(dispatch-type number))))
	 (complex-atanh number)
	 (coerce (%atanh (coerce number 'double-float))
		 '(dispatch-type number))))
1196
1197
1198
1199
1200
1201
    #+double-double
    ((double-double-float)
     (if (or (> number 1w0)
	     (< number -1w0))
	 (complex-atanh number)
	 (dd-%atanh (coerce number 'double-double-float))))
ram's avatar
ram committed
1202
1203
    ((complex)
     (complex-atanh number))))
1204

1205
1206
1207
;;; HP-UX does not supply a C version of log1p, so use the definition.
;;; We really need to fix this.  The definition really loses big-time
;;; in roundoff as x gets small.
1208
1209

#+hpux
pw's avatar
pw committed
1210
(declaim (inline %log1p))
1211
#+hpux
pw's avatar
pw committed
1212
1213
1214
1215
(defun %log1p (number)
  (declare (double-float number)
	   (optimize (speed 3) (safety 0)))
  (the double-float (log (the (double-float 0d0) (+ number 1d0)))))
ram's avatar
ram committed
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267


;;;;
;;;; This is a set of routines that implement many elementary
;;;; transcendental functions as specified by ANSI Common Lisp.  The
;;;; implementation is based on Kahan's paper.
;;;;
;;;; I believe I have accurately implemented the routines and are
;;;; correct, but you may want to check for your self.
;;;;
;;;; These functions are written for CMU Lisp and take advantage of
;;;; some of the features available there.  It may be possible,
;;;; however, to port this to other Lisps.
;;;;
;;;; Some functions are significantly more accurate than the original
;;;; definitions in CMU Lisp.  In fact, some functions in CMU Lisp
;;;; give the wrong answer like (acos #c(-2.0 0.0)), where the true
;;;; answer is pi + i*log(2-sqrt(3)).
;;;;
;;;; All of the implemented functions will take any number for an
;;;; input, but the result will always be a either a complex
;;;; single-float or a complex double-float.
;;;;
;;;; General functions
;;;;   complex-sqrt
;;;;   complex-log
;;;;   complex-atanh
;;;;   complex-tanh
;;;;   complex-acos
;;;;   complex-acosh
;;;;   complex-asin
;;;;   complex-asinh
;;;;   complex-atan
;;;;   complex-tan
;;;;
;;;; Utility functions:
;;;;   scalb logb
;;;;
;;;; Internal functions:
;;;;    square coerce-to-complex-type cssqs complex-log-scaled
;;;;
;;;;
;;;; Please send any bug reports, comments, or improvements to Raymond
;;;; Toy at toy@rtp.ericsson.se.
;;;;
;;;; References
;;;;
;;;; Kahan, W. "Branch Cuts for Complex Elementary Functions, or Much
;;;; Ado About Nothing's Sign Bit" in Iserles and Powell (eds.) "The
;;;; State of the Art in Numerical Analysis", pp. 165-211, Clarendon
;;;; Press, 1987
;;;;
1268

ram's avatar
ram committed
1269
1270
(declaim (inline square))
(defun square (x)
1271
  (declare (float x))
ram's avatar
ram committed
1272
1273
1274
1275
1276
1277
  (* x x))

;; If you have these functions in libm, perhaps they should be used
;; instead of these Lisp versions.  These versions are probably good
;; enough, especially since they are portable.

1278
(declaim (inline scalb))
ram's avatar
ram committed
1279
(defun scalb (x n)
1280
  "Compute 2^N * X without compute 2^N first (use properties of the
ram's avatar
ram committed
1281
underlying floating-point format"
1282
  (declare (type float x)
1283
	   (type double-float-exponent n))
ram's avatar
ram committed
1284
1285
  (scale-float x n))

1286
1287
(declaim (inline logb-finite))