Ticket #1059: math-patches.lisp

File math-patches.lisp, 6.4 KB (added by dfindlay, 18 months ago)

Updated math patches

Line 
1;;; Math Patches
2
3;;; phase (l0-numbers.lisp)
4;;; Returns wrong value for negative zero
5;;;  - use %double-float-sign, %short-float-sign instead of minusp
6
7(defun phase (number)
8  "Return the angle part of the polar representation of a complex number.
9  For complex numbers, this is (atan (imagpart number) (realpart number)).
10  For non-complex positive numbers, this is 0. For non-complex negative
11  numbers this is PI."
12  (number-case number
13    (rational
14     (if (minusp number)
15       (%short-float pi)
16       0.0f0))
17    (double-float
18
19;;; Changed code (1)
20
21     (if (%double-float-sign number)
22
23;;; End of change (1)
24
25       (%double-float pi)
26       0.0d0))
27    (complex
28     (atan (%imagpart number) (%realpart number)))
29    (short-float
30
31;;; Changed code (2)
32
33     (if (%short-float-sign number)
34
35;;; End of change (2)
36
37       (%short-float pi)
38       0.0s0))))
39
40;;; exp (l0-float.lisp)
41;;; Allow ridiculously large negative real rationals
42;;;  - these just need to return 0.0
43
44(defun exp (x)
45  "Return e raised to the power NUMBER."
46  (typecase x
47    (complex (* (exp (realpart x)) (cis (imagpart x))))
48    (double-float (%double-float-exp! x (%make-dfloat)))
49    (t
50
51;;; Changed code
52
53     (if (and (typep x 'rational)
54              (< x -104))
55       0.0s0
56     #+32-bit-target
57     (target::with-stack-short-floats ((sx x))
58       (%single-float-exp! sx (%make-sfloat)))
59     #+64-bit-target
60     (%single-float-exp (%short-float x))))))
61;;; End of change
62
63;;; tan (l0-float.lisp)
64;;; Complex numbers with very large imaginary parts were overflowing in (* 2 (abs i))
65;;;  - recode to avoid this
66
67(defun tan (x)
68  "Return the tangent of NUMBER."
69  (cond ((complexp x)
70         (let ((r (realpart x))
71               (i (imagpart x)))
72           (if (zerop i)
73             (complex (tan r) i)
74             (let* ((tx (tan r))
75                    (ty (tanh i))
76                    (tx2 (* tx tx))
77                    (d (1+ (* tx2 (* ty ty))))
78
79;;; Changed code
80
81                    (c (if (> (abs i) 20)
82                         (* 2 (exp (- (abs i))))
83                         (/ (cosh i)))))
84               (complex (/ (* (* c c) tx) d)
85                        (/ (* ty (1+ tx2)) d))))))
86
87;;; End of change
88
89        ((or (typep x 'ratio)
90             (> (abs x) two^23))
91         (let ((c (%extended-cis x)))
92           (if (typep x 'double-float)
93             (/ (imagpart c) (realpart c))
94             (%short-float (/ (imagpart c) (realpart c))))))
95        ((typep x 'double-float)
96         (%double-float-tan! x (%make-dfloat)))
97        (t
98         #+32-bit-target
99         (target::with-stack-short-floats ((sx x))
100           (%single-float-tan! sx (%make-sfloat)))
101         #+64-bit-target
102         (%single-float-tan (%short-float x))
103         )))
104
105;;; %complex-atanh (l0-float.lisp)
106;;; The expression (/ 2 y) was overflowing for very small y
107;;;  - recode log(2/y) as log(2) - log(y)
108;;;    replace atan(2/y, -1) with atan(2, -y) of atan(-2, y) depending on sign of y
109;;; Argument to log1+ becoming -1 when x is near 1 and y is small
110;;;  - add additional cond clause to handle this
111
112(defun %complex-atanh (z)
113  (let* ((rx (realpart z))
114         (ix (imagpart z))
115         (sign (typecase rx
116                 (double-float (%double-float-sign rx))
117                 (short-float (%short-float-sign rx))
118                 (t (minusp rx))))
119         (x rx)
120         (y ix)
121         (y1 (abs y))
122         ra ia)
123    ;; following code requires non-negative x
124    (when sign
125      (setf x (- x))
126      (setf y (- y)))
127    (cond ((> (max x y1) 1.8014399e+16)
128           ;; large value escape
129           (setq ra (if (> x y1)
130                      (let ((r (/ y x)))
131                        (/ (/ x) (1+ (* r r))))
132                      (let ((r (/ x y)))
133                        (/ (/ r y) (1+ (* r r))))))
134           (setq ia (typecase y
135                      (double-float (float-sign y double-float-half-pi))
136                      (single-float (float-sign y single-float-half-pi))
137                      (t (if (minusp y) #.(- single-float-half-pi) single-float-half-pi)))))
138          ((= x 1)
139
140;;; Changed code
141
142           (cond ((< y1 1e-9)
143                  (setq ra (/ (- (if (typep y 'double-float) double-float-log2 single-float-log2)
144                                 (log-e y1))
145                              2))
146                  (setq ia (/ (if (minusp y) (atan -2 y) (atan 2 (- y))) 2)))
147                 (t
148                  (setq ra (/ (log1+ (/ 4 (* y y))) 4))
149                  (setq ia (/ (atan (/ 2 y) -1) 2)))))
150          ((and (< y1 1)
151                (< 0.5 x 2))
152           (let ((x-1 (- x 1))
153                 (x+1 (+ x 1))
154                 (y2 (* y y)))
155             (setq ra (/ (log-e (/ (+ (* x-1 x-1) y2) (+ (* x+1 x+1) y2))) -4))
156             (setq ia (/ (atan (* 2 y) (- 1 (+ (* x x) y2))) 2))))
157
158;;; End of change
159
160          (t
161           (let ((r2 (+ (* x x) (* y y))))
162             (setq ra (/ (log1+ (/ (* -4 x) (1+ (+ (* 2 x) r2)))) -4))
163             (setq ia (/ (atan (* 2 y) (- 1 r2)) 2)))))
164    ;; fixup signs, with special case for real arguments
165    (cond (sign
166           (setq ra (- ra))
167           (when (typep z 'complex)
168             (setq ia (- ia))))
169          (t
170           (unless (typep z 'complex)
171             (setq ia (- ia)))))
172    (complex ra ia)))
173
174;;; tanh (numbers.lisp)
175;;; Complex numbers with very large real parts were overflowing in (* 2 (abs r))
176;;;  - recode to avoid this
177;;; (analogous fix to that in tan)
178
179(defun tanh (x)
180  "Return the hyperbolic tangent of NUMBER."
181  (cond ((complexp x)
182         (let ((r (realpart x))
183               (i (imagpart x)))
184           (if (zerop r)
185             (complex r (tan i))
186             (let* ((tx (tanh r))
187                    (ty (tan i))
188                    (ty2 (* ty ty))
189                    (d (1+ (* (* tx tx) ty2)))
190
191;;; Changed code
192
193                    (c (if (> (abs r) 20)
194                         (* 2 (exp (- (abs r))))
195                         (/ (cosh r)))))
196               (complex (/ (* tx (1+ ty2)) d)
197                        (/ (* (* c c) ty) d))))))
198
199;;; End of change
200
201        ((typep x 'double-float)
202         (%double-float-tanh! x (%make-dfloat)))
203        ((and (typep x 'rational)
204              (> (abs x) 12))
205         (if (plusp x) 1.0s0 -1.0s0))
206        (t
207         #+32-bit-target
208         (target::with-stack-short-floats ((sx x))
209           (%single-float-tanh! sx (%make-sfloat)))
210         #+64-bit-target
211         (%single-float-tanh (%short-float x)))))
212