Ticket #869: ccl-math-reference.lisp

File ccl-math-reference.lisp, 53.2 KB (added by dfindlay, 3 years ago)

"Reference" implementation of math functions

Line 
1;;; Reference math function implementations
2;;; =======================================
3;;; Hopefully, all cases where the ANSI function definitions are "terrible methods
4;;; for floating point computation" have been addressed :-)
5;;; Care has been taken for rational (and complex rational) arguments to maintain
6;;; precision by not converting to (short) float too early.
7;;;
8;;; In most cases (exp and exp-like functions excepted) the common lisp functions
9;;; are only called for known floating point arguments. It is assumed that they
10;;; return accurate values in these cases!
11;;;
12;;; This lot should be applicable in any Lisp that (just) uses IEEE single and
13;;; double floating point formats.
14;;; (This code does assume that short-float = single-float)
15;;;
16;;; It should be "obvious" how this lot can be (back) ported to CCL...
17;;; (Replace calls to common lisp functions by #| ... |# commented-out blocks.)
18
19(cl:defpackage "MATH-REFERENCE"
20  (:use "COMMON-LISP")
21  (:nicknames "MATH")
22  (:shadow "SIGNUM" "ABS" "PHASE" "CIS" "SIN" "COS" "TAN" "ASIN" "ACOS" "ATAN"
23           "LOG" "EXP" "EXPT" "SQRT" "SINH" "COSH" "TANH" "ASINH" "ACOSH" "ATANH")
24  (:export "SIGNUM" "ABS" "PHASE" "CIS" "SIN" "COS" "TAN" "ASIN" "ACOS" "ATAN"
25           "LOG" "EXP" "EXPT" "SQRT" "SINH" "COSH" "TANH" "ASINH" "ACOSH" "ATANH"))
26
27(in-package :math)
28
29;;; === Stubs for low level CCL functions
30
31(defmacro number-case (keyform &body clauses)
32  (let ((value (gensym)))
33    `(let ((,value ,keyform))
34       (typecase ,value
35         ,@clauses
36         (otherwise
37          (error "value ~S is not of the expected type NUMBER." ,value))))))
38
39(defmacro %realpart (x)
40  `(realpart ,x))
41
42(defmacro %imagpart (x)
43  `(imagpart ,x))
44
45(defmacro %numerator (x)
46  `(numerator ,x))
47
48(defmacro %denominator (x)
49  `(denominator ,x))
50
51(defmacro %short-float (x)
52  `(float ,x 1s0))
53
54(defmacro %double-float (x)
55  `(float ,x 1d0))
56
57(defmacro with-stack-double-floats (vars &body body)
58  `(let ,(mapcar (lambda (item)
59                   `(,(first item) (%double-float ,(second item))))
60                 vars)
61     ,@body))
62
63(defmacro report-bad-arg (arg type)
64  `(error "~S: not of type ~S" ,arg ,type))
65
66(defmacro ratiop (x)
67  `(typep ,x 'ratio))
68
69(defmacro bignump (x)
70  `(typep ,x 'bignum))
71
72(defmacro fsqrt (x)
73  `(cl:sqrt ,x))
74
75
76(defconstant IEEE-single-float-bias 126)
77(defconstant IEEE-double-float-bias 1022)
78
79
80(defun %double-float-sign (x)
81  (or (minusp x) (eql x -0.0d0)))
82
83(defun %short-float-sign (x)
84  (or (minusp x) (eql x -0.0s0)))
85
86;;; These next two aren't quite accurate, but I don't think that matters :-)
87(defun %short-float-exp (x)
88  (+ (nth-value 1 (decode-float x)) IEEE-single-float-bias))
89
90(defun %double-float-exp (x)
91  (+ (nth-value 1 (decode-float x)) IEEE-double-float-bias))
92
93
94(defun %integer-power (b e)
95  (declare (type unsigned-byte e))
96  (if (zerop e)
97    (+ 1 (* b 0))
98    (if (eql b 2)
99      (ash 1 e)
100      (do* ((next (ash e -1) (ash e -1))
101            (oddexp (oddp e) (oddp e))
102            (total (if oddexp b 1) (if oddexp (* b total) total)))
103           ((zerop next) total)
104        (declare (type unsigned-byte next))
105        (setq b (* b b) e next)))))
106
107(defun %hypot (x y)
108  (let ((a (abs x))
109        (b (abs y)))
110    (when (> a b)
111      (psetq a b b a))
112    (if (= b 0.0d0)
113      0.0d0
114      (let ((c (/ a b)))
115        (* b (sqrt (+ 1.0d0 (* c c))))))))
116
117
118;;; === Unmodified CCL function definitions
119;;; (Copies needed to pick up shadowed versions of called functions)
120
121(defun abs (number)
122  "Return the absolute value of the number."
123  (number-case number
124   (real (cl:abs number))
125   #|
126   (fixnum
127    (locally (declare (fixnum number))
128      (if (minusp number) (- number) number)))
129   (double-float
130    (%double-float-abs number))
131   (short-float
132    (%short-float-abs number))
133   (bignum
134    (if (bignum-minusp number)(negate-bignum number) number))
135   (ratio
136    (if (minusp number) (- number) number))
137   |#
138   (complex
139    (let ((rx (%realpart number))
140          (ix (%imagpart number)))
141      (number-case rx
142        (rational
143         (sqrt (+ (* rx rx) (* ix ix))))
144        (short-float
145         (%short-float (%hypot (%double-float rx)
146                               (%double-float ix))))
147        (double-float
148         (%hypot rx ix)))))))
149
150(defun phase (number)
151  "Return the angle part of the polar representation of a complex number.
152  For complex numbers, this is (atan (imagpart number) (realpart number)).
153  For non-complex positive numbers, this is 0. For non-complex negative
154  numbers this is PI."
155  (number-case number
156    (rational
157     (if (minusp number)
158       (%short-float pi)
159       0.0f0))
160    (double-float
161     (if (minusp number)
162       (%double-float pi)
163       0.0d0))
164    (complex
165     (atan (%imagpart number) (%realpart number)))
166    (short-float
167     (if (minusp number)
168       (%short-float pi)
169       0.0s0))))
170
171(defun exp (x)
172  "Return e raised to the power NUMBER."
173  (typecase x
174    (complex (* (exp (realpart x)) (cis (imagpart x))))
175    #|
176    (double-float (%double-float-exp! x (%make-dfloat)))
177    (t
178     #+32-bit-target
179     (target::with-stack-short-floats ((sx x))
180       (%single-float-exp! sx (%make-sfloat)))
181     #+64-bit-target
182     (%single-float-exp (%short-float x)))
183    |#
184    (t (cl:exp x))))
185
186;;; === New function definitions ===
187;;; 1. in l0-float.lisp
188
189(eval-when (:execute :compile-toplevel)
190  (defconstant two^23 (ash 1 23))
191)
192
193(defun sin (x)
194  "Return the sine of NUMBER."
195  (cond ((complexp x)
196         (let* ((r (realpart x))
197                (i (imagpart x)))
198           (complex (* (sin r) (cosh i))
199                    (* (cos r) (sinh i)))))
200        ((or (typep x 'ratio)
201             (> (abs x) two^23))
202         (if (typep x 'double-float)
203           (imagpart (%extended-cis x))
204           (%short-float (imagpart (%extended-cis x)))))
205        #|
206        ((typep x 'double-float)
207         (%double-float-sin! x (%make-dfloat)))
208        |#
209        ((typep x 'double-float)
210         (cl:sin x))
211        #|
212        (t
213         #+32-bit-target
214         (target::with-stack-short-floats ((sx x))
215           (%single-float-sin! sx (%make-sfloat)))
216         #+64-bit-target
217         (%single-float-sin (%short-float x)))
218        |#
219        (t (cl:sin (%short-float x)))))
220
221(defun cos (x)
222  "Return the cosine of NUMBER."
223  (cond ((complexp x)
224         (let* ((r (realpart x))
225                (i (imagpart x)))
226           (complex (* (cos r) (cosh i))
227                    (- (* (sin r) (sinh i))))))
228        ((or (typep x 'ratio)
229             (> (abs x) two^23))
230         (if (typep x 'double-float)
231           (realpart (%extended-cis x))
232           (%short-float (realpart (%extended-cis x)))))
233        #|
234        ((typep x 'double-float)
235         (%double-float-cos! x (%make-dfloat)))
236        |#
237        ((typep x 'double-float)
238         (cl:cos x))
239        #|
240        (t
241         #+32-bit-target
242         (target::with-stack-short-floats ((sx x))
243           (%single-float-cos! sx (%make-sfloat)))
244         #+64-bit-target
245         (%single-float-cos (%short-float x)))
246        |#
247        (t (cl:cos (%short-float x)))))
248
249(defun tan (x)
250  "Return the tangent of NUMBER."
251  (cond ((complexp x)
252         (let ((r (realpart x))
253               (i (imagpart x)))
254           (if (zerop i)
255             (complex (tan r) i)
256             (let* ((tx (tan r))
257                    (ty (tanh i))
258                    (tx2 (* tx tx))
259                    (d (1+ (* tx2 (* ty ty))))
260                    (n (if (> (abs i) 20)
261                         (* 4 (exp (* -2 (abs i))))
262                         (let ((c (cosh i)))
263                           (/ (* c c))))))
264               (complex (/ (* n tx) d)
265                        (/ (* ty (1+ tx2)) d))))))
266        ((or (typep x 'ratio)
267             (> (abs x) two^23))
268         (let ((c (%extended-cis x)))
269           (if (typep x 'double-float)
270             (/ (imagpart c) (realpart c))
271             (%short-float (/ (imagpart c) (realpart c))))))
272        #|
273        ((typep x 'double-float)
274         (%double-float-tan! x (%make-dfloat)))
275        |#
276        ((typep x 'double-float)
277         (cl:tan x))
278        #|
279        (t
280         #+32-bit-target
281         (target::with-stack-short-floats ((sx x))
282           (%single-float-tan! sx (%make-sfloat)))
283         #+64-bit-target
284         (%single-float-tan (%short-float x))
285         )
286        |#
287        (t (cl:tan (%short-float x)))))
288
289;;; Helper function for sin/cos/tan for ratio or large arguments
290;;; (Will become inaccurate for ridiculously large arguments though)
291;;; Does not assume that float library returns accurate values for large arguments
292;;; (many don't)
293
294;;; hexdecimal representations of pi at various precisions
295(defconstant pi-vector
296  #(#x3243F6A8885A308D313198A2E0
297    #x3243F6A8885A308D313198A2E03707344A4093822299F31D008
298    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D
299    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC
300    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B5470
301    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310B
302    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB
303    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045
304    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70
305    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871
306    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D9
307    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD588215
308    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D
309    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085
310    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB
311    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3E
312    #x3243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3ED71577C1BD314B2778AF2FDA5
313    ))
314
315(defun %extended-cis (x)
316  (let (size pi-size)
317    (typecase x
318      (integer (setq size (1- (integer-length (abs x)))))
319      (ratio (setq size (- (integer-length (abs (numerator x)))
320                           (integer-length (denominator x)))))
321      (float (multiple-value-bind (mantissa exponent sign)
322                                  (integer-decode-float x)
323               (setq size (+ (1- (integer-length mantissa)) exponent))
324               (setq x (* sign mantissa (expt 2 exponent))))))
325    (setq pi-size (ceiling (+ size 64) 100))
326    (cond ((< pi-size 1) (setq pi-size 1))
327          ((> pi-size 17) (setq pi-size 17)))
328    (let* ((2pi-approx (/ (aref pi-vector (1- pi-size))
329                          (ash 1 (1- (* 100 pi-size)))))
330           (reduced-x (rem x 2pi-approx))
331           (x0 (float reduced-x 1.0d0))
332           (x1 (- reduced-x (rational x0))))
333      (* (cis x0) (cis (float x1 1.0d0))))))
334
335;;; Multiply by i (with additional optional scale factor)
336;;; Does the "right thing" with minus zeroes (see CLTL2)
337(defun i* (number &optional (scale 1))
338  (complex (* (- scale) (imagpart number))
339           (* scale (realpart number))))
340
341
342
343;;; >>> move these values of constants up from below <<<
344(eval-when (:execute :compile-toplevel)
345  (defconstant single-float-pi (coerce pi 'single-float))
346  (defconstant double-float-half-pi (/ pi 2))
347  (defconstant single-float-half-pi (coerce (/ pi 2) 'single-float))
348)
349
350(defun atan (y &optional (x nil x-p))
351  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
352  (cond (x-p
353         (cond ((or (typep x 'double-float)
354                    (typep y 'double-float))
355                (with-stack-double-floats ((dy y)
356                                           (dx x))
357                  #|
358                  (%df-atan2 dy dx)
359                  |#
360                  (cl:atan (%double-float dy) (%double-float dx))))
361               (t
362                (when (and (rationalp x) (rationalp y))
363                  ;; rescale arguments so that the maximum absolute value is 1
364                  (let ((x1 (abs x)) (y1 (abs y)))
365                    (cond ((> y1 x1)
366                           (setf x (/ x y1))
367                           (setf y (signum y)))
368                          ((not (zerop x))
369                           (setf y (/ y x1))
370                           (setf x (signum x))))))
371                #|
372                #+32-bit-target
373                (target::with-stack-short-floats ((sy y)
374                                                  (sx x))
375                  (%sf-atan2! sy sx))
376                #+64-bit-target
377                (%sf-atan2 (%short-float y) (%short-float x))
378                |#
379                (cl:atan (%short-float y) (%short-float x)))))
380        ((typep y 'float)
381         (cl:atan y))
382        #|
383        ((typep y 'double-float)
384         (%double-float-atan! y (%make-dfloat)))
385        ((typep y 'single-float)
386         #+32-bit-target
387         (%single-float-atan! y (%make-sfloat))
388         #+64-bit-target
389         (%single-float-atan y))
390        |#
391        ((typep y 'rational)
392         (cond ((<= (abs y) most-positive-short-float)
393                #|
394                #+32-bit-target
395                (target::with-stack-short-floats ((sy y))
396                  (%single-float-atan! sy (%make-sfloat)))
397                #+64-bit-target
398                (%single-float-atan (%short-float y))
399                |#
400                (cl:atan (%short-float y)))
401               ((minusp y)
402                #.(- single-float-half-pi))
403               (t
404                single-float-half-pi)))
405        (t
406         (let ((r (realpart y))
407               (i (imagpart y)))
408           (if (zerop i)
409             (complex (atan r) i)
410             (i* (%complex-atanh (complex (- i) r)) -1))))
411        ))
412
413
414;;; useful constants
415(eval-when (:execute :compile-toplevel)
416  (defconstant single-float-log2 0.6931472)                ; (log 2)
417  (defconstant double-float-log2 0.6931471805599453d0)     ; (log 2.0d0)
418  (defconstant double-float-log2^23 15.942385152878742d0)  ; (log (expt 2 23))
419)
420
421(defun log (x &optional (b nil b-p))
422  "Return the logarithm of NUMBER in the base BASE, which defaults to e."
423  (if b-p
424    (if (zerop b)
425      (if (zerop x)
426        (report-bad-arg x '(not (satisfies zerop) ))
427        ;(if (floatp x) (float 0.0d0 x) 0)
428        ; ** CORRECT THE CONTAGION for complex args **
429        (+ 0 (* x b)))
430      ; do the float/rational contagion before the division
431      ; but take care with negative zeroes
432      (let ((x1 (realpart x))
433            (b1 (realpart b)))
434        (if (and (typep x1 'float)
435                 (typep b1 'float))
436          (/ (log-e (* x (float 1.0 b1)))
437             (log-e (* b (float 1.0 x1))))
438          (let ((r (/ (cond ((typep x 'rational)
439                             (%rational-log x 1.0d0))
440                            ((typep x1 'rational)
441                             (%rational-complex-log x 1.0d0))
442                            (t
443                             (log-e (* x 1.0d0))))
444                      (cond ((typep b 'rational)
445                             (%rational-log b 1.0d0))
446                            ((typep b1 'rational)
447                             (%rational-complex-log b 1.0d0))
448                            (t
449                             (log-e (* b 1.0d0)))))))
450            (cond ((or (typep x1 'double-float)
451                       (typep b1 'double-float))
452                   r)
453                  ((complexp r)
454                   (complex (%short-float (realpart r))
455                            (%short-float (imagpart r))))
456                  (t
457                   (%short-float r)))))))
458    (log-e x)))
459
460;;; Cases rearranged to put more esoteric cases later
461(defun log-e (x)
462  (cond
463    ((typep x 'double-float)
464     (if (%double-float-sign x)
465       (complex (cl:log (- x)) pi)
466       (cl:log x)))
467    ((typep x 'short-float)
468     (if (%short-float-sign x)
469       (complex (cl:log (- x)) single-float-pi)
470       (cl:log x)))
471    #|
472    ((typep x 'double-float)
473     (if (%double-float-sign x)
474       (with-stack-double-floats ((dx x))
475         (complex (%double-float-log! (%%double-float-abs! dx dx) (%make-dfloat)) pi))
476       (%double-float-log! x (%make-dfloat))))
477    ((typep x 'short-float)
478     #+32-bit-target
479     (if (%short-float-sign x)
480       (target::with-stack-short-floats ((sx x))
481         (complex (%single-float-log! (%%short-float-abs! sx sx) (%make-sfloat))
482                  single-float-pi))
483       (%single-float-log! x (%make-sfloat)))
484     #+64-bit-target
485     (if (%short-float-sign x)
486       (complex (%single-float-log (%short-float-abs (%short-float x)))
487                single-float-pi)
488       (%single-float-log (%short-float x))))
489    |#
490    ((typep x 'complex)
491     (if (typep (realpart x) 'rational)
492       (%rational-complex-log x 1.0s0)
493       ;; take care that intermediate results do not overflow/underflow:
494       ;; pre-scale argument by an appropriate power of two;
495       ;; we only need to scale for very large and very small values -
496       ;;  hence the various 'magic' numbers (values may not be too
497       ;;  critical but do depend on the sqrt update to fix abs's operation)
498       (let ((m (max (abs (realpart x))
499                     (abs (imagpart x))))
500             (log-s 0)
501             (s 1))
502         (if (typep m 'short-float)
503           (let ((expon (- (%short-float-exp m) IEEE-single-float-bias)))
504             (cond ((> expon 126)
505                    (setq log-s double-float-log2^23)
506                    (setq s #.(ash 1 23)))
507                   ((< expon -124)
508                    (setq log-s #.(- double-float-log2^23))
509                    (setq s #.(/ 1.0s0 (ash 1 23))))))
510           (let ((expon (- (%double-float-exp m) IEEE-double-float-bias)))
511             (cond ((> expon 1022)
512                    (setq log-s double-float-log2^23)
513                    (setq s #.(ash 1 23)))
514                   ((< expon -1020)
515                    (setq log-s #.(- double-float-log2^23))
516                    (setq s #.(/ 1.0d0 (ash 1 23)))))))
517         (if (eql s 1)
518           (complex (log-abs x) (phase x))
519           (let ((temp (log-abs (/ x s))))
520             (complex (float (+ log-s temp) temp) (phase x)))))))
521    (t
522     (%rational-log x 1.0s0))))
523
524;;; (log1+ x) = (log (1+ x))
525;;; but is much more precise near x = 0
526(defun log1+ (x)
527  ;(cond ((typep x 'complex)
528  ;       (let ((r (realpart x))
529  ;             (i (imagpart x)))
530  ;         (if (and (< (abs r) 0.5)
531  ;                  (< (abs i) 3))
532  ;           (let* ((n (+ (* r (+ 2 r)) (* i i)))
533  ;                  (d (1+ (sqrt (1+ n)))))
534  ;             (complex (log1+ (/ n d)) (atan i (1+ r))))
535  ;          (log (1+ x)))))
536  ;      (t
537  (if (and (typep x 'ratio)
538           (< -0.5 x 2))
539    (setq x (%short-float x)))
540  (let ((y (1+ x)))
541    (if (eql y x)
542      (log-e y)
543      (let ((e (1- y)))
544        (if (zerop e)
545          (* x 1.0)
546          (- (log-e y) (/ (- e x) y)))))))
547  ;))
548
549;;; helper for complex log
550;;; uses more careful approach when (abs x) is near 1
551(defun log-abs (x)
552  (let ((a (abs x)))
553    (if (< 0.5 a 3)
554      (let* ((r (realpart x))
555             (i (imagpart x))
556             (n (if (> (abs r) (abs i))
557                  (+ (* (1+ r) (1- r)) (* i i))
558                  (+ (* r r) (* (1+ i) (1- i))))))
559        (log1+ (/ n (1+ a))))
560      (log-e a))))
561
562;;; helpers for rational log
563(defun %rational-log (x prototype)
564  (cond ((minusp x)
565         (complex (%rational-log (- x) prototype)
566                  (if (typep prototype 'short-float)
567                    single-float-pi
568                    pi)))
569        ((bignump x)
570         ;(let* ((base1 3)
571         ;       (guess (floor (1- (integer-length x))
572         ;                     (log base1 2)))
573         ;       (guess1 (* guess (log-e base1))))
574         ;  (+ guess1 (log-e (/ x (expt base1 guess)))))
575         ; Using base1 = 2 is *much* faster. Was there a reason for 3?
576         (let* ((guess (1- (integer-length x)))
577                (guess1 (* guess double-float-log2)))
578           (float (+ guess1 (log-e (float (/ x (ash 1 guess)) 1.0d0))) prototype)))
579        ((and (ratiop x)
580              ;; Rational arguments near +1 can be specified with great precision: don't lose it
581              (cond ((< 0.5 x 3)
582                     (log1+ (float (- x 1) prototype)))
583                    (
584                     ;; Escape out small values as well as large
585                     (or (> x most-positive-short-float)
586                         (< x least-positive-normalized-short-float))
587                     ;; avoid loss of precision due to subtracting logs of numerator and denominator
588                     (let* ((n (%numerator x))
589                            (d (%denominator x))
590                            (sn (1- (integer-length n)))
591                            (sd (1- (integer-length d))))
592                       (float (+ (* (- sn sd) double-float-log2)
593                                 (- (log1+ (float (1- (/ n (ash 1 sn))) 1.0d0))
594                                    (log1+ (float (1- (/ d (ash 1 sd))) 1.0d0))))
595                              prototype))))))
596        ((typep prototype 'short-float)
597         #|
598         #+32-bit-target
599         (target::with-stack-short-floats ((sx x))
600           (%single-float-log! sx (%make-sfloat)))
601         #+64-bit-target
602         (%single-float-log (%short-float x))
603         |#
604         (cl:log (%short-float x)))
605        (t
606         #|
607         (with-stack-double-floats ((dx x))
608           (%double-float-log! dx (%make-dfloat)))
609         |#
610         (cl:log (%double-float x)))))
611
612(defun %rational-complex-log (x prototype &aux ra ia)
613  (let* ((rx (realpart x))
614         (ix (imagpart x))
615         (x (abs rx))
616         (y (abs ix)))
617    (if (> y x)
618      (let ((r (float (/ rx y) 1.0d0)))
619        (setq ra (+ (%rational-log y 1.0d0)
620                    (/ (log1+ (* r r)) 2)))
621        (setq ia (atan (if (minusp ix) -1.0d0 1.0d0) r)))
622      (let ((r (float (/ ix x) 1.0d0)))
623        (setq ra (+ (%rational-log x 1.0d0)
624                    (/ (log1+ (* r r)) 2)))
625        (setq ia (atan r (if (minusp rx) -1.0d0 1.0d0)))))
626    (if (typep prototype 'short-float)
627      (complex (%short-float ra) (%short-float ia))
628      (complex ra ia))))
629
630
631
632;;; (It may be possible to do something with rational exponents, e.g. so that
633;;;       (expt x 1/2) == (sqrt x)
634;;;       (expt x 3/2) == (expt (sqrt x) 3)      ** NOT (sqrt (expt x 3)) !! **
635;;;                      or (* x (sqrt x))
636;;;       (expt x 1/8) == (sqrt (sqrt (sqrt x)))
637;;;    even, in the case of rational x, returning a rational result if possible.)
638;;;
639(defun expt (b e)
640  "Return BASE raised to the POWER."
641  (cond ((zerop e) (1+ (* b e)))
642        ((integerp e)
643         (if (minusp e) (/ 1 (%integer-power b (- e))) (%integer-power b e)))
644        ((zerop b)
645         (if (plusp (realpart e)) (* b e) (report-bad-arg e '(number (0) *))))
646        ((and (realp b) (plusp b) (realp e)
647              ; escape out very small or very large rationals
648              ; - avoid problems converting to float
649              (typecase b
650                (bignum (<= b most-positive-short-float))
651                (ratio (cond ((< b 0.5)
652                              (>= b least-positive-normalized-short-float))
653                             ((> b 3)
654                              (<= b most-positive-short-float))))
655                (t t)))
656         ;; assumes that the library routines are accurate
657         ;; (if not, just excise this whole clause)
658         #|
659         (if (or (typep b 'double-float)
660                 (typep e 'double-float))
661           (with-stack-double-floats ((b1 b)
662                                      (e1 e))
663             (%double-float-expt! b1 e1 (%make-dfloat)))
664           #+32-bit-target
665           (target::with-stack-short-floats ((b1 b)
666                                             (e1 e))
667             (%single-float-expt! b1 e1 (%make-sfloat)))
668           #+64-bit-target
669           (%single-float-expt (%short-float b) (%short-float e))
670           )
671         |#
672         (if (or (typep b 'double-float)
673                 (typep e 'double-float))
674           (cl:expt (%double-float b) (%double-float e))
675           (cl:expt (%short-float b) (%short-float e))))
676        ((typep b 'rational)
677         (let ((r (exp (* e (%rational-log b 1.0d0)))))
678           (cond ((typep (realpart e) 'double-float)
679                  r)
680                 ((typep r 'complex)
681                  (complex (%short-float (realpart r)) (%short-float (imagpart r))))
682                 (t
683                  (%short-float r)))))
684        ((typep (realpart b) 'rational)
685         (let ((r (exp (* e (%rational-complex-log b 1.0d0)))))
686           (if (typep (realpart e) 'double-float)
687             r
688             (complex (%short-float (realpart r)) (%short-float (imagpart r))))))
689        (t
690         ;; type upgrade b without losing -0.0 ...
691         (let ((r (exp (* e (log-e (* b 1.0d0))))))
692           (cond ((or (typep (realpart b) 'double-float)
693                      (typep (realpart e) 'double-float))
694                  r)
695                 ((typep r 'complex)
696                  (complex (%short-float (realpart r)) (%short-float (imagpart r))))
697                 (t
698                  (%short-float r)))))))
699
700
701
702;;; For complex square root,
703;;;  if (a + i*b)^2 = rx + i * ix
704;;;  then, the basic solution is
705;;;   a = sqrt((rx + sqrt(rx^2 + ix^2))/2)
706;;;   b = ix / (2 * a)
707;;;  but, in practice, there are various cases that must be handled specially
708;;; Note: for rational complex arguments, we can search for (and return) exact results
709(defun sqrt (x &aux a b)
710  "Return the square root of NUMBER."
711  (cond ((zerop x) x)
712        ((complexp x)
713         (let ((rx (realpart x))
714               (ix (imagpart x)))
715           (cond ((rationalp rx)
716                  (if (zerop rx)
717                    (let ((s (sqrt (/ (abs ix) 2))))
718                      (complex s (if (minusp ix) (- s) s)))
719                    (let* ((s (+ (* rx rx) (* ix ix)))
720                           (d (if (ratiop s)
721                                (/ (isqrt (%numerator s))
722                                   (isqrt (%denominator s)))
723                                (isqrt s))))
724                      (unless (eql s (* d d))
725                        (setf d (%hypot (%double-float rx)
726                                        (%double-float ix))))
727                      (cond ((minusp rx)
728                             (setq b (sqrt (/ (- d rx) 2)))
729                             (when (minusp ix)
730                               (setq b (- b)))
731                             (setq a (/ ix (+ b b))))
732                            (t
733                             (setq a (sqrt (/ (+ rx d) 2)))
734                             (setq b (/ ix (+ a a)))))
735                      (if (rationalp a)
736                        (complex a b)
737                        (complex (%short-float a) (%short-float b))))))
738                 ((minusp rx)
739                  (if (zerop ix)
740                    (complex 0 (float-sign ix (sqrt (- rx))))
741                    (let ((shift (cond ((< rx -1) -3)
742                                       ((and (> rx -5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
743                                       (t -1))))
744                      (setq rx (scale-float rx shift))
745                      (let ((s (fsqrt (- (abs (complex rx (scale-float ix shift))) rx))))
746                        (setq b (scale-float s (ash (- -1 shift) -1)))
747                        (when (minusp ix)
748                          (setq b (- b)))
749                        (setq a (/ ix (scale-float b 1)))
750                        (complex a b)))))
751                 (t
752                  (if (zerop ix)
753                    (complex (sqrt rx) ix)
754                    (let ((shift (cond ((> rx 1) -3)
755                                       ((and (< rx 5.9604645E-8) (< (abs ix) 5.9604645E-8)) 25)
756                                       (t -1))))
757                      (setq rx (scale-float rx shift))
758                      (let ((s (fsqrt (+ rx (abs (complex rx (scale-float ix shift)))))))
759                        (setq a (scale-float s (ash (- -1 shift) -1)))
760                        (setq b (/ ix (scale-float a 1)))
761                        (complex a b))))))))
762        ((minusp x) (complex 0 (sqrt (- x))))
763        ((floatp x)
764         (fsqrt x))
765        ((and (integerp x) (eql x (* (setq a (isqrt x)) a))) a)
766        ((and (ratiop x)
767              (let ((n (numerator x))
768                    d)
769                (and (eql n (* (setq a (isqrt n)) a))
770                     (eql (setq d (denominator x))
771                          (* (setq b (isqrt d)) b)))))
772         (/ a b))         
773        (t
774         ;; use double float (sometimes x can't be represented in single)
775         (with-stack-double-floats ((f1 x))
776           (%short-float (fsqrt f1))))))
777
778
779
780
781(defun asin (x)
782  "Return the arc sine of NUMBER."
783  (cond ((and (typep x 'double-float)
784              (locally (declare (type double-float x))
785                (and (<= -1.0d0 x)
786                     (<= x 1.0d0))))
787         #|
788         (%double-float-asin! x (%make-dfloat))
789         |#
790         (cl:asin x))
791        ((and (typep x 'single-float)
792              (locally (declare (type single-float x))
793                (and (<= -1.0s0 x)
794                     (<= x 1.0s0))))
795         #|
796         #+32-bit-target
797         (%single-float-asin! x (%make-sfloat))
798         #+64-bit-target
799         (%single-float-asin x)
800         |#
801         (cl:asin x))
802        ((and (typep x 'rational)
803              (<= (abs x) 1))
804         (if (integerp x)
805           (case x
806             (0 0.0s0)                          ; or simply 0 ??
807             (1 single-float-half-pi)
808             (-1 #.(- single-float-half-pi)))
809           (atan x (sqrt (- 1 (* x x))))))
810        (t
811         (%complex-asin/acos x nil))
812        ))
813
814
815(defun acos (x)
816  "Return the arc cosine of NUMBER."
817  (cond ((and (typep x 'double-float)
818              (locally (declare (type double-float x))
819                (and (<= -1.0d0 x)
820                     (<= x 1.0d0))))
821         #|
822         (%double-float-acos! x (%make-dfloat))
823         |#
824         (cl:acos x))
825        ((and (typep x 'short-float)
826              (locally (declare (type short-float x))
827                (and (<= -1.0s0 x)
828                     (<= x 1.0s0))))
829         #|
830         #+32-bit-target
831         (%single-float-acos! x (%make-sfloat))
832         #+64-bit-target
833         (%single-float-acos x)
834         |#
835         (cl:acos x))
836        ((and (typep x 'rational)
837              (<= (abs x) 1))
838         (if (integerp x)
839           (case x
840             (0 single-float-half-pi)
841             (1 0.0s0)                          ; or simply 0 ??
842             (-1 single-float-pi))
843           (atan (sqrt (- 1 (* x x))) x)))
844        (t
845         (%complex-asin/acos x t))
846        ))
847
848;;; combined complex asin/acos routine
849;;; argument acos is true for acos(z); false for asin(z)
850;;;
851;;; based on Hull, Fairgrieve & Tang, ACM TMS 23, 3, 299-335 (Sept. 1997)
852(defun %complex-asin/acos (z acos)
853  (let* ((rx (realpart z))
854         (ix (imagpart z))
855         (x (abs rx))
856         (y (abs ix))
857         (m (max x y)))
858    (if (> m 1.8014399E+16)
859      ;; Large argument escape
860      (let ((log-s 0))
861        (if (typep m 'double-float)
862          (if (> m #.(/ most-positive-double-float 2))
863            (setq log-s double-float-log2)
864            (setq z (* 2 z)))
865          (if (> m #.(/ most-positive-short-float 2))
866            (setq log-s single-float-log2)
867            (setq z (* 2 z))))
868        (if acos
869          (i* (+ log-s (log-e z))
870              (if (minusp ix) +1 -1))
871          (if (minusp ix)
872            (i* (+ log-s (log-e (i* z 1))) -1)
873            (i* (+ log-s (log-e (i* z -1))) 1))))
874      (let ((qrx rx)
875            (qx x)
876            x-1 y2 s)
877        (cond ((rationalp rx)
878               (setq x-1 (float (abs (- x 1))))
879               (setq rx (float rx))
880               (setq x (abs rx))
881               (setq y (float y))
882               (setq y2 (* y y))
883               (setq s (cond ((zerop x-1)
884                              y)
885                             ((> y x-1)
886                              (let ((c (/ x-1 y)))
887                                (* y (sqrt (1+ (* c c))))))
888                             (t
889                              (let ((c (/ y x-1)))
890                                (* x-1 (sqrt (1+ (* c c)))))))))
891              (t
892               (setq x-1 (abs (- x 1)))
893               (setq y2 (* y y))
894               (setq s (if (zerop x-1)
895                         y
896                         (sqrt (+ (* x-1 x-1) y2))))))
897        (let* ((x+1 (+ x 1))
898               (r (sqrt (+ (* x+1 x+1) y2)))
899               (a (/ (+ r s) 2))
900               (b (/ rx a))
901               (ra (if (<= (abs b) 0.6417)
902                     (if acos (acos b) (asin b))
903                     (let* ((r+x+1 (+ r x+1))
904                            (s+x-1 (+ s x-1))
905                            (a+x (+ a x))
906                            (ry (if (<= qx 1)
907                                  (let ((aa (+ (/ y2 r+x+1) s+x-1)))
908                                    (sqrt (/ (* a+x aa) 2)))
909                                  (let ((aa (+ (/ a+x r+x+1) (/ a+x s+x-1))))
910                                    (* y (sqrt (/ aa 2)))))))
911                       (if acos (atan ry rx) (atan rx ry)))))
912               (ia (if (<= a 1.5)
913                     (let* ((r+x+1 (+ r x+1))
914                            (s+x-1 (+ s x-1))
915                            (ll (if (< qx 1)
916                                  (let* ((aa (/ (+ (/ 1 r+x+1) (/ 1 s+x-1)) 2)))
917                                    (+ (* aa y2) (* y (sqrt (* aa (1+ a))))))
918                                  (let* ((a-1 (/ (+ (/ y2 r+x+1) s+x-1) 2)))
919                                    (+ a-1 (sqrt (* a-1 (1+ a))))))))
920                       (log1+ ll))
921                     (log (+ a (sqrt (1- (* a a))))))))
922          ;; final fixup of signs
923          (if acos
924            (if (complexp z)
925              (if (typep ix 'float)
926                (setq ia (float-sign (- ix) ia))
927                (if (plusp ix)
928                  (setq ia (- ia))))
929              (if (< qrx -1)
930                (setq ia (- ia))))
931            (if (complexp z)
932              (if (typep ix 'float)
933                (setq ia (float-sign ix ia))
934                (if (minusp ix)
935                  (setq ia (- ia))))
936              (if (> qrx 1)
937                (setq ia (- ia)))))
938          (complex ra ia))))))
939
940;;; complex atanh
941(defun %complex-atanh (z)
942  (let* ((rx (realpart z))
943         (ix (imagpart z))
944         (sign (typecase rx
945                 (double-float (%double-float-sign rx))
946                 (short-float (%short-float-sign rx))
947                 (t (minusp rx))))
948         (x rx)
949         (y ix)
950         (y1 (abs y))
951         ra ia)
952    ;; following code requires non-negative x
953    (when sign
954      (setf x (- x))
955      (setf y (- y)))
956    (cond ((> (max x y1) 1.8014399e+16)
957           ;; large value escape
958           (setq ra (if (> x y1)
959                      (let ((r (/ y x)))
960                        (/ (/ x) (1+ (* r r))))
961                      (let ((r (/ x y)))
962                        (/ (/ r y) (1+ (* r r))))))
963           (setq ia (typecase y
964                      (double-float (float-sign y double-float-half-pi))
965                      (single-float (float-sign y single-float-half-pi))
966                      (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi)))))
967          ((= x 1)
968           (setq ra (if (< y1 1e-9)
969                      (/ (log-e (/ 2 y1)) 2)
970                      (/ (log1+ (/ 4 (* y y))) 4)))
971           (setq ia (/ (atan (/ 2 y) -1) 2)))
972          (t
973           (let ((r2 (+ (* x x) (* y y))))
974             (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4))
975             (setq ia (/ (atan (* 2 y) (- 1 r2)) 2)))))
976    ;; fixup signs, with special case for real arguments
977    (cond (sign
978           (setq ra (- ra))
979           (when (typep z 'complex)
980             (setq ia (- ia))))
981          (t
982           (unless (typep z 'complex)
983             (setq ia (- ia)))))
984    (complex ra ia)))
985
986
987;;;
988;;; 2. In l0-numbers.lisp
989
990(defun cis (theta)
991  "Return cos(Theta) + i sin(Theta), i.e. exp(i Theta)."
992  (cond ((complexp theta)
993         (error "Argument to CIS is complex: ~S" theta))
994        ((or (typep theta 'ratio)
995             (> (abs theta) #.(ash 1 23)))
996         (if (typep theta 'double-float)
997           (%extended-cis theta)
998           (coerce (%extended-cis theta) '(complex single-float))))
999        (t
1000         (complex (cos theta) (sin theta)))))
1001
1002;;;
1003;;; 3. In numbers.lisp
1004
1005(defun signum (x)
1006  "If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."
1007  (cond ((complexp x) (if (zerop x)
1008                        x
1009                        (let ((m (max (abs (realpart x))
1010                                      (abs (imagpart x)))))
1011                          (cond ((rationalp m)
1012                                 ;; rescale to avoid intermediate under/overflow
1013                                 (setq x (/ x m)))
1014                                ((> m #.(ash 1 23))
1015                                 ;; ensure no overflow for abs
1016                                 (setq x (/ x 2))))
1017                          (/ x (abs x)))))
1018        ((rationalp x) (if (plusp x) 1 (if (zerop x) 0 -1)))
1019        ((zerop x) (float 0.0 x))
1020        (t (float-sign x))))
1021
1022(defun sinh (x)
1023  "Return the hyperbolic sine of NUMBER."
1024  (if (complexp x)
1025    (let ((r (realpart x))
1026          (i (imagpart x)))
1027      (complex (* (sinh r) (cos i))
1028               (* (cosh r) (sin i))))
1029    #|
1030    (if (typep x 'double-float)
1031      (%double-float-sinh! x (%make-dfloat))
1032      #+32-bit-target
1033      (target::with-stack-short-floats ((sx x))
1034        (%single-float-sinh! sx (%make-sfloat)))
1035      #+64-bit-target
1036        (%single-float-sinh (%short-float x)))
1037    |#
1038    (cl:sinh x)))
1039
1040(defun cosh (x)
1041  "Return the hyperbolic cosine of NUMBER."
1042  (if (complexp x)
1043    (let ((r (realpart x))
1044          (i (imagpart x)))
1045      (complex (* (cosh r) (cos i))
1046               (* (sinh r) (sin i))))
1047    #|
1048    (if (typep x 'double-float)
1049      (%double-float-cosh! x (%make-dfloat))
1050      #+32-bit-target
1051      (target::with-stack-short-floats ((sx x))
1052        (%single-float-cosh! sx (%make-sfloat)))
1053      #+64-bit-target
1054      (%single-float-cosh (%short-float x)))
1055    |#
1056    (cl:cosh x)))
1057
1058(defun tanh (x)
1059  "Return the hyperbolic tangent of NUMBER."
1060  (cond ((complexp x)
1061         (let ((r (realpart x))
1062               (i (imagpart x)))
1063           (if (zerop r)
1064             (complex r (tan i))
1065             (let* ((tx (tanh r))
1066                    (ty (tan i))
1067                    (ty2 (* ty ty))
1068                    (d (1+ (* (* tx tx) ty2)))
1069                    (n (if (> (abs r) 20)
1070                         (* 4 (exp (- (* 2 (abs r)))))
1071                         (let ((c (cosh r)))
1072                           (/ (* c c))))))
1073               (complex (/ (* tx (1+ ty2)) d)
1074                        (/ (* n ty) d))))))
1075        ((typep x 'double-float)
1076         (cl:tanh x))
1077        #|
1078        ((typep x 'double-float)
1079         (%double-float-tanh! x (%make-dfloat)))
1080        |#
1081        ((and (typep x 'rational)
1082              (> (abs x) 12))
1083         (if (plusp x) 1.0s0 -1.0s0))
1084        #|
1085        (t
1086         #+32-bit-target
1087         (target::with-stack-short-floats ((sx x))
1088           (%single-float-tanh! sx (%make-sfloat)))
1089         #+64-bit-target
1090         (%single-float-tanh (%short-float x)))
1091        |#
1092        (t (cl:tanh (%short-float x)))))
1093
1094(defun asinh (x)
1095  "Return the hyperbolic arc sine of NUMBER."
1096  (cond ((typep x 'double-float)
1097         #|
1098         (%double-float-asinh! x (%make-dfloat))
1099         |#
1100         (cl:asinh x))
1101        ((typep x 'short-float)
1102         #|
1103         #+32-bit-target
1104         (%single-float-asinh! x (%make-sfloat))
1105         #+64-bit-target
1106         (%single-float-asinh (%short-float x))
1107         |#
1108         (cl:asinh x))
1109        ((typep x 'rational)
1110         (if (< (abs x) most-positive-short-float)
1111           #|
1112           #+32-bit-target
1113           (target::with-stack-short-floats ((sx x))
1114             (%single-float-asinh! sx (%make-sfloat)))
1115           #+64-bit-target
1116           (%single-float-asinh (%short-float x))
1117           |#
1118           (cl:asinh (%short-float x))
1119           (* (signum x) (log-e (* 2 (abs x))))))
1120        (t
1121         (i* (%complex-asin/acos (i* x) nil) -1))))
1122
1123;;; for complex case, use acos and post-fix the branch cut
1124(defun acosh (x)
1125  "Return the hyperbolic arc cosine of NUMBER."
1126  (cond ((and (typep x 'double-float)
1127              (locally (declare (type double-float x))
1128                (<= 1.0d0 x)))
1129         #|
1130         (%double-float-acosh! x (%make-dfloat))
1131         |#
1132         (cl:acosh x))
1133        ((and (typep x 'short-float)
1134              (locally (declare (type short-float x))
1135                (<= 1.0s0 x)))
1136         #|
1137         #+32-bit-target
1138         (%single-float-acosh! x (%make-sfloat))
1139         #+64-bit-target
1140         (%single-float-acosh (%short-float x))
1141         |#
1142         (cl:acosh x))
1143        ((and (typep x 'rational)
1144              (<= 1 x))
1145         (cond ((< x 2)
1146                (log1+ (+ (- x 1) (sqrt (- (* x x) 1)))))
1147               ((<= x most-positive-short-float)
1148                #|
1149                #+32-bit-target
1150                (target::with-stack-short-floats ((x1 x))
1151                  (%single-float-acosh! x1 (%make-sfloat)))
1152                #+64-bit-target
1153                (%single-float-acosh (%short-float x))
1154                |#
1155                (cl:acosh (%short-float x)))
1156               (t
1157                (log-e (* 2 x)))))
1158        (t
1159         (let ((sign (and (typep x 'complex)
1160                          (let ((ix (imagpart x)))
1161                            (typecase ix
1162                              (double-float (%double-float-sign ix))
1163                              (single-float (%short-float-sign ix))
1164                              (t (minusp ix)))))))
1165           (i* (%complex-asin/acos x t) (if sign -1 1))))))
1166
1167(defun atanh (x)
1168  "Return the hyperbolic arc tangent of NUMBER."
1169  (cond ((and (typep x 'double-float)
1170              (locally (declare (type double-float x))
1171                (and (<= -1.0d0 x)
1172                     (<= x 1.0d0))))
1173         #|
1174         (%double-float-atanh! x (%make-dfloat))
1175         |#
1176         (cl:atanh x))
1177        ((and (typep x 'short-float)
1178              (locally (declare (type short-float x))
1179                (and (<= -1.0s0 x)
1180                     (<= x 1.0s0))))
1181         #|
1182         #+32-bit-target
1183         (%single-float-atanh! x (%make-sfloat))
1184         #+64-bit-target
1185         (%single-float-atanh x)
1186         |#
1187         (cl:atanh x))
1188        ((and (typep x 'rational)
1189              (<= (abs x) 1))
1190         (let ((n (numerator x))
1191               (d (denominator x)))
1192           (/ (log-e (/ (+ d n) (- d n))) 2)))
1193        (t
1194         (let ((r (realpart x)))
1195           (if (zerop r)
1196             (complex r (atan (imagpart x)))
1197             (%complex-atanh x))))))
1198
1199
1200(defun approx= (x y)
1201  (labels ((ulp (x)
1202             (multiple-value-bind (m e s) (integer-decode-float x)
1203               (if (eql m 0)
1204                 (values x x)
1205                 (let ((x2 (scale-float (float (1+ m) x) e))
1206                       (x1 (scale-float (float (1- m) x) e)))
1207                   (if (eql s 1)
1208                     (values x1 x2)
1209                     (values (- x2) (- x1)))))))
1210           (match (x y)
1211             (handler-case
1212                 (or (and (zerop x) (zerop y))
1213                     (multiple-value-bind (x1 x2) (ulp x)
1214                       (multiple-value-bind (y1 y2) (ulp y)
1215                         (and (> y2 x1) (> x2 y1)))))
1216               (error () nil))))
1217    (and (numberp x)
1218         (numberp y)
1219         (cond ((eql x y) t)
1220               ((or (and (typep x 'single-float) (typep y 'single-float))
1221                    (and (typep x 'double-float) (typep y 'double-float)))
1222                (and (match x y) :approx))
1223               ((or (and (typep x '(complex single-float)) (typep y '(complex single-float)))
1224                    (and (typep x '(complex double-float)) (typep y '(complex double-float))))
1225                (and (match (realpart x) (realpart y))
1226                     (match (imagpart x) (imagpart y))
1227                     :approx))))))
1228           
1229(defun test (list-of-forms)
1230  (flet ((eval-form (fn form)
1231           (handler-case
1232               (eval (cons fn (rest form)))
1233             (error () :Error))))
1234    (format t "~&; Expression~54TReference Value~108TCL Value~%")
1235    (format t   "; ----------~54T---------------~108T--------~%")
1236    (dolist (form list-of-forms)
1237      (let* ((label (princ-to-string form))
1238             (fn-name (symbol-name (car form)))
1239             (ref-fn (find-symbol fn-name :math))
1240             (cl-fn (find-symbol fn-name :cl))
1241             (ref-r (eval-form ref-fn form))
1242             (cl-r (eval-form cl-fn form))
1243             (mark (case (approx= ref-r cl-r)
1244                     ((nil) #\#)
1245                     (:approx #\~)
1246                     (t #\Space))))
1247        (format t
1248                "~&;~A~A~A~54T~A~108T~A~159T~A~%"
1249                mark
1250                (if (> (length label) 50)
1251                  (subseq label 0 50)
1252                  label)
1253                (if (> (length label) 50) #\\ #\Space)
1254                ref-r
1255                cl-r
1256                mark))))
1257  (values))
1258
1259(test '(
1260        (math:expt 1.0 0)
1261        (math:expt (complex 2.0 3.0) 0)
1262        (math:expt 0 1.2)
1263        (math:expt 0 (complex 2.0d0 3.0d0))
1264
1265        (math:log (expt 10 -66))
1266        (math:log (expt 10 66))
1267        (math:log (- (expt 10 -66)))
1268        (math:log (- (expt 10 66)))
1269        (math:log (complex (expt 10 66) (expt 10 68)))
1270        (math:log (complex (expt 10 -66) (expt 10 -68)))
1271        (math:log 8.0d0 2)
1272        (math:log #c(2 3) 0)
1273        (math:log #c(2.5 3.6) 0)
1274        (math:log #c(2.5d0 3.6d0) 0)
1275        (math:log (complex most-positive-single-float most-positive-single-float))
1276        (math:log (complex least-positive-single-float most-positive-single-float))
1277        (math:log (complex least-positive-single-float least-positive-single-float))
1278        (math:log (complex most-positive-double-float most-positive-double-float))
1279        (math:log (complex least-positive-double-float least-positive-double-float))
1280        (math:log (rational pi))
1281        (math:log (/ (expt 400 66) (expt 27 66)))
1282        (math:log (/ (expt 27 66) (expt 400 66)))
1283
1284        (math:sqrt (expt 10 47))
1285        (math:sqrt (/ (expt 10 47) 3))
1286        (math:sqrt (complex (expt 10 46) (expt 10 47)))
1287        (math:sqrt (complex (expt 10 -46) (expt 10 -47)))
1288        (math:sqrt (complex (expt 10 65) (expt 10 67)))
1289        (math:sqrt (complex (expt 10 -65) (expt 10 -67)))
1290        (math:sqrt (complex (/ (expt 10 65) 3) (expt 10 67)))
1291        (math:sqrt (complex most-positive-single-float most-positive-single-float))
1292        (math:sqrt (complex least-positive-single-float most-positive-single-float))
1293        (math:sqrt (complex least-positive-single-float least-positive-single-float))
1294        (math:sqrt (complex most-positive-double-float most-positive-double-float))
1295        (math:sqrt (complex least-positive-double-float least-positive-double-float))
1296
1297        (math:expt (expt 10 47) 1/2)
1298        (math:expt (/ (expt 10 47) 3) 1/2)
1299        (math:expt (complex (expt 10 46) (expt 10 47)) 1/2)
1300        (math:expt (complex (expt 10 -46) (expt 10 -47)) 1/2)
1301        (math:expt (complex (expt 10 65) (expt 10 67)) 1/2)
1302        (math:expt (complex (expt 10 -65) (expt 10 -67)) 1/2)
1303        (math:expt (complex (/ (expt 10 65) 3) (expt 10 67)) 1/2)
1304        (math:expt (expt 10 600) 0.5d0)
1305
1306        (math:expt (expt 10 47) 0.5)
1307        (math:expt (/ (expt 10 47) 3) 0.5)
1308        (math:expt (complex (expt 10 46) (expt 10 47)) 0.5)
1309        (math:expt (complex (expt 10 -46) (expt 10 -47)) 0.5)
1310        (math:expt (complex (expt 10 65) (expt 10 67)) 0.5)
1311        (math:expt (complex (expt 10 -65) (expt 10 -67)) 0.5)
1312        (math:expt (complex (/ (expt 10 65) 3) (expt 10 67)) 0.5)
1313
1314        (math:expt (complex -2.0 -0.0) (complex 0.5d0 0.0d0))
1315
1316        (math:atan (expt 10 46) (expt 10 47))
1317        (math:atan (expt 10 -46) (expt 10 -47))
1318        (math:atan most-positive-single-float most-positive-single-float)
1319        (math:atan least-positive-single-float least-positive-single-float)
1320
1321        (math:sin (expt 3 20))
1322        (math:cos (expt 3 20))
1323        (math:tan (expt 3 20))
1324        (math:sin (expt 3 40))
1325        (math:cos (expt 3 40))
1326        (math:tan (expt 3 40))
1327        (math:sin (expt 10 47))
1328        (math:cos (expt 10 47))
1329        (math:tan (expt 10 47))
1330        (math:tan (complex 3 4))
1331        (math:tan (complex 1.2 40))
1332
1333        (math:sinh (complex 1e-5 1e-6))
1334        (math:cosh (complex 1e-5 1e-6))
1335        (math:tanh (complex 1e6 10.0))
1336
1337        (math:tanh (complex 3 4))
1338        (math:tanh (complex 40 1.2))
1339
1340        (math:asin (complex 1.0 0.0))
1341        (math:asin (complex 0.1 1e-23))
1342        (math:asin (complex 0.99 1e-23))
1343        (math:asin (complex 1.01 1e-23))
1344        (math:asin (complex 10.0 1e-23))
1345        (math:asin (complex 1e-23 1.0))
1346        (math:asin (complex 1e-24 1e16))
1347        (math:asin (+ 1 (expt 10 -10)))
1348        (math:asin (- 1 (expt 10 -10)))
1349        (math:asin (complex 1 (expt 10 -10)))
1350        (math:asin (complex 1 (- (expt 10 -10))))
1351        (math:asin -10000)
1352        (math:asin (complex most-positive-single-float most-positive-single-float))
1353
1354        (math:acos (complex -1.0 0.0))
1355        (math:acos (complex 0.1 1e-23))
1356        (math:acos (complex 0.99 1e-23))
1357        (math:acos (complex 1.01 1e-23))
1358        (math:acos (complex 10.0 1e-23))
1359        (math:acos (complex 1e-23 1.0))
1360        (math:acos (complex 1e-24 1e16))
1361        (math:acos (+ 1 (expt 10 -10)))
1362        (math:acos (- 1 (expt 10 -10)))
1363        (math:acos (complex 1 (expt 10 -10)))
1364        (math:acos (complex 1 (- (expt 10 -10))))
1365        (math:acos -10000)
1366        (math:acos (complex most-positive-single-float most-positive-single-float))
1367
1368        (math:acosh (+ 1 (expt 10 -10)))
1369        (math:acosh (- 1 (expt 10 -10)))
1370        (math:log (+ 1 (expt 10 -10)))
1371        (math:log (- 1 (expt 10 -10)))
1372        (math:log 2 (+ 1 (expt 10 -10)))
1373        (math:log 2.0 (+ 1 (expt 10 -10)))
1374        (math:log 2.0d0 (+ 1 (expt 10 -10)))
1375
1376        (math:atan (complex 0 (+ 1 (expt 10 -15))))
1377        (math:atan (complex 0 (- 1 (expt 10 -15))))
1378        (math:atan (complex (expt 10 -10) 1))
1379        (math:atan (complex (- (expt 10 -10)) 1))
1380        (math:atanh (+ 1 (expt 10 -10)))
1381        (math:atanh (- 1 (expt 10 -10)))
1382
1383        (math:asinh (expt 10 46))
1384        (math:acosh (expt 10 46))
1385        (math:atanh (expt 10 46))
1386        (math:asinh (complex 1e-5 1e-6))
1387        (math:acosh (complex 1e-5 1e-6))
1388        (math:atanh (complex 1e-5 1e-6))
1389        (math:asinh (complex 1e-15 1e-6))
1390        (math:acosh (complex 1e-15 1e-6))
1391        (math:atanh (complex 1e-15 1e-6))
1392        (math:asinh (complex 1e-3 -1e6))
1393        (math:acosh (complex 1e-3 -1e6))
1394        (math:atanh (complex 1e-3 -1e6))
1395        (math:asinh (complex 1d-30 1d-12))
1396        (math:acosh (complex 1d-30 1d-12))
1397        (math:atanh (complex 1d-30 1d-12))
1398        (math:asinh (complex 1d-6 -1d12))
1399        (math:acosh (complex 1d-6 -1d12))
1400        (math:atanh (complex 1d-6 -1d12))
1401        (math:asinh (complex 1e-5 1e-16))
1402        (math:acosh (complex 1e-5 1e-16))
1403        (math:atanh (complex 1e-5 1e-16))
1404        (math:atanh (complex 1.0 1e-32))
1405))
1406