Ticket #869: fn-ticket.lisp

File fn-ticket.lisp, 9.7 KB (added by dfindlay, 3 years ago)

sqrt/log/atan changes [demo]

Line 
1(in-package :cl-user)
2
3;;; This code is written in cl-user, not ccl
4;;; The code for new-sqrt, new-log and new-atan follows that for CCL's sqrt, log
5;;; and atan (in l0-float.lisp) in terms of the case structure; but where the branch
6;;; is not changed the CCL function is called (rather than copying the code in CCL's
7;;; branch, which probably wouldn't work here). I hope it shouldn't be too difficult
8;;; to cut and paste the appropriate code bits into CCL/s code.
9
10;;; Use of new %hypot is assumed:
11;;; (but this version isn't as efficient as the *real* one)
12(defun ccl::%hypot (x y)
13  (let ((a (abs x))
14        (b (abs y)))
15    (when (> a b)
16      (psetq a b b a))
17    (if (= b 0.0d0)
18      b
19      (let ((c (/ a b)))
20        (* b (sqrt (+ 1.0d0 (* c c))))))))
21 
22
23(defmacro single-float-sqrt (x)
24  ;; usual expression for single float sqrt
25  ;; (referred to several times in new-sqrt)
26  `(sqrt ,x))
27
28
29;; fix sqrt to return values for rationals whose square root can be represented
30;; as single float, even though they can't be themselves.
31;; also take care to avoid intermediate overflow for some complex arguments
32(defun new-sqrt (x)
33  "Return the square root of NUMBER."
34  (cond ((zerop x) x)
35        ((complexp x)
36         ;; take care that intermediate results do not overflow/underflow:
37         ;; pre-scale argument by an appropriate power of two
38         ;; (there may be more efficient ways of doing this; and strictly
39         ;; if doesn't need to be done for 'reasonably sized' arguments)
40         (let ((m (max (abs (realpart x))
41                       (abs (imagpart x))))
42               (s 1)
43               (s2 1))
44           (typecase m
45             (integer (let ((expon (floor (integer-length m) 2)))
46                        (setq s (ash 1 expon))
47                        (setq s2 (ash s expon))))
48             (ratio (let ((expon (floor (- (integer-length (numerator m))
49                                           (integer-length (denominator m)))
50                                        2)))
51                      (setq s (expt 2 expon))
52                      (setq s2 (* s s))))
53             (float (let ((expon (truncate (nth-value 1 (decode-float m)) 2)))
54                      (if (minusp expon)
55                        (incf expon)
56                        (decf expon))
57                      (setq s (scale-float (float 1 m) expon))
58                      (setq s2 (scale-float s expon)))))
59           (let ((z (/ x s2)))
60             (* s (* (sqrt (abs z)) (cis (/ (phase z) 2)))))))
61        ((minusp x) (complex 0 (new-sqrt (- x))))
62        ((floatp x)
63         ;; usual floating point sqrt
64         (sqrt x))
65        ((integerp x)
66         (let* ((a (isqrt x))
67                (test (* a a)))
68           (cond ((eql x test)
69                  a)
70                 ((<= x most-positive-single-float)
71                  (single-float-sqrt x))
72                 (t
73                  ;; handle case where (sqrt x) can be represented as single float
74                  ;; but x can't
75                  (let ((error (/ x test)))
76                    (* (float a 1.0)
77                       (single-float-sqrt error)))))))
78        ((ratiop x)
79         ;; similar to integer case - be careful when (sqrt x) can be represented
80         ;;  as single float but x can't
81         (let* ((n (numerator x))
82                (a (isqrt n))
83                (test1 (* a a)))
84           (unless (eql n test1)
85             (when (and (<= x most-positive-single-float)
86                        (>= x least-positive-normalized-single-float))
87               (return-from new-sqrt (single-float-sqrt x))))
88           (let* ((d (denominator x))
89                  (b (isqrt d))
90                  (test2 (* b b)))
91             (cond ((and (eql d test2)
92                         (eql n test1))
93                    (/ a b))
94                   ((and (<= x most-positive-single-float)
95                         (>= x least-positive-normalized-single-float))
96                    (single-float-sqrt x))
97                   (t
98                    (let ((error (/ x (/ test1 test2))))
99                      (* (float (/ a b) 1.0)
100                         (single-float-sqrt error))))))))
101        (t
102         (single-float-sqrt x))))
103
104
105;; makes log work for small rationals (it already worked for big ones)
106;; correct wrong return type for negative bignums
107;; fix float/rational contagion in two argument case
108;; propagate rational behaviour to complex rationals as well
109(defun new-log (x &optional (b nil b-p))
110  "Return the logarithm of NUMBER in the base BASE, which defaults to e."
111  (if b-p
112    (if (zerop b)
113      (if (zerop x)
114        ;(report-bad-arg x '(not (satisfies zerop) ))
115        (error "NEW-LOG: zero arguments")
116        ;(if (floatp x) (float 0.0d0 x) 0))   ** CORRECT THE CONTAGION for complex args **
117        ; (is this the best way to get a zero of the right type?)
118        (+ 0 (* x b)))
119      ;; ** Don't we have to do the float/rational contagion before the division?
120      ;; (is this the best way to do it?)
121      (/ (new-log-e (+ x (* b 0)))
122         (new-log-e (+ b (* x 0)))))
123    (new-log-e x)))
124
125
126;;; constant required: will need to be 'magic' number somewhere
127(defconstant +log2+ #.(log 2.0d0)) ; = 0.6931471805599453D0
128
129(defun new-log-e (x)
130  (cond
131    ((bignump x)
132     (if (minusp x)
133       (complex (new-log-e (- x)) #.(coerce pi 'single-float))   ; ** FIX RETURN TYPE ERROR **
134       (let* ((base1 3)
135              (guess (floor (1- (integer-length x))
136                            (log base1 2)))
137              (guess1 (* guess (new-log-e base1))))
138         (+ guess1 (new-log-e (/ x (expt base1 guess)))))))
139    ((and (ratiop x)                                ; Escape out small values as well as large
140          (if (minusp x)
141            (or (> x least-negative-normalized-short-float)
142                (< x most-negative-short-float))
143            (or (> x most-positive-short-float)
144                (< x least-positive-normalized-short-float))))
145     (- (new-log-e (numerator x)) (new-log-e (denominator x))))
146    ((typep x 'complex)
147     ;; take care that intermediate results do not overflow/underflow:
148     ;; pre-scale argument by an appropriate power of two
149     ;; we only need to scale for very large and very small values -
150     ;;  hence the various 'magic' numbers (values may not be too
151     ;;  critical but do depend on the sqrt update to fix abs's operation)
152     (let ((m (max (abs (realpart x))
153                   (abs (imagpart x))))
154           (log-s 0)
155           (s 1))
156       (typecase m
157         (integer (let ((expon (integer-length m)))
158                    (when (> expon 126)
159                      (setq log-s (float (* expon +log2+) 1.0))
160                      (setq s (ash 1 expon)))))
161         (ratio (let ((expon (- (integer-length (numerator m))
162                                (integer-length (denominator m)))))
163                  (when (or (> expon 125)
164                            (< expon -122))
165                    (setq log-s (float (* expon +log2+) 1.0))
166                    (setq s (expt 2 expon)))))
167         (single-float (let ((expon (nth-value 1 (decode-float m))))
168                         (when (or (> expon 126)
169                                   (< expon -124))
170                           (if (minusp expon)
171                             (incf expon)
172                             (decf expon))
173                           (setq log-s (float (* expon +log2+) 1.0))     ; assumes (float-radix m) = 2
174                           (setq s (scale-float 1.0 expon)))))
175         (double-float (let ((expon (nth-value 1 (decode-float m))))
176                         (when (or (> expon 1022)
177                                   (< expon -1020))
178                           (if (minusp expon)
179                             (incf expon)
180                             (decf expon))
181                           (setq log-s (* expon +log2+))     ; assumes (float-radix m) = 2
182                           (setq s (scale-float 1.0d0 expon))))))
183       (if (eql s 1)
184         (complex (new-log-e (abs x)) (phase x))
185         (let ((z (/ x s)))
186           (complex (+ log-s (new-log-e (abs z))) (phase z))))))
187    ((typep x 'double-float)
188     ;; usual double float log-e
189     (log x))
190    (t
191     ;; usual single float log-e
192     (log x)
193     )))
194
195
196;; atan needs to be more careful about small and large rational arguments
197;; (this will make phase as careful too for complex rationals)
198(defun new-atan (y &optional (x nil x-p))
199  "Return the arc tangent of Y if X is omitted or Y/X if X is supplied."
200  (cond (x-p
201         (when (or (rationalp x) (rationalp y))
202           (unless (zerop y)           ; don't tread on special zero handling
203             ;; rescale arguments so that the maximum absolute value is 1
204             (cond ((> (abs y) (abs x))
205                    (setf x (/ x (abs y)))
206                    (setf y (signum y)))
207                   (t
208                    (setf y (/ y (abs x)))
209                    (setf x (signum x))))))
210         ;; usual expression
211         (atan y x))
212        (t
213         ;; usual expressions
214         ;; [but (sqrt -1) in the complex case could be just #C(0 1)]
215         (atan y))))
216
217
218
219
220;;; FIXES
221#|
222
223(sqrt (expt 10 47)) => 3.1622778E+23
224(sqrt (/ (expt 10 47) 3)) => 1.8257418E+23
225(sqrt (complex (expt 10 46) (expt 10 47))) => #C(2.3505187E+23 2.12719E+23)
226(sqrt (complex most-positive-short-float most-positive-short-float)) => #C(2.0267142E+19 8.394926E+18)
227
228(log (expt 10 -66)) => -151.97063
229(log (- (expt 10 66))) => #C(151.97063 3.1415927) [corrects type]
230(log (complex (expt 10 65) (expt 10 66))) => #C(151.9756 1.4711276)
231(log (complex (expt 10 -65) (expt 10 -66))) => #C(-149.66307 0.09966865)
232(log 8.0d0 2) => 3.0D0 [corrects accuracy]
233(log #C(0.0 1.0) 0) => #C(0.0 0.0) [corrects type]
234
235(atan (expt 10 46) (expt 10 47)) => 0.09966865
236(atan (expt 10 -46) (expt 10 -47)) => 1.4711276 [corrects loss of accuracy]
237
238|#
239