1 ;;;; -*- Mode: lisp -*-
3 ;;;; Copyright (c) 2007 Raymond Toy
5 ;;;; Permission is hereby granted, free of charge, to any person
6 ;;;; obtaining a copy of this software and associated documentation
7 ;;;; files (the "Software"), to deal in the Software without
8 ;;;; restriction, including without limitation the rights to use,
9 ;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
10 ;;;; copies of the Software, and to permit persons to whom the
11 ;;;; Software is furnished to do so, subject to the following
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
17 ;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
18 ;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
19 ;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
20 ;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
21 ;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
22 ;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 ;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
24 ;;;; OTHER DEALINGS IN THE SOFTWARE.
29 ;; Compute to how many bits EST and TRUE are equal. If they are
30 ;; identical, return T.
31 (defun bit-accuracy (est true
)
32 (let* ((diff (abs-qd (sub-qd est true
)))
33 (err (if (zerop-qd true
)
35 (cl:/ (qd-0 diff
) (abs (qd-0 true
))))))
36 (if (zerop (qd-0 diff
))
38 (cl:-
(log err
2d0
)))))
40 (defun print-result (est true
)
41 (format t
"est: ~/qdi::qd-format/~%" est
)
42 (format t
"tru: ~/qdi::qd-format/~%" true
)
43 (format t
"err: ~A~%" (qd-0 (sub-qd est true
)))
44 (format t
"bits: ~,1f~%" (bit-accuracy est true
)))
46 ;; Machin's formula for pi
48 (defun atan-series (x)
50 (eps (make-qd-d (scale-float 1d0 -
212)
54 (s1 (make-qd-dd 0w0
0w0
))
57 (loop while
(qd-> tmp eps
) do
61 (sub-qd s1
(div-qd tmp
(make-qd-d d
0d0
0d0
0d0
)))
62 (add-qd s1
(div-qd tmp
(make-qd-d d
0d0
0d0
0d0
)))))
64 (setf tmp
(mul-qd tmp r
))
65 (setf sign
(cl:- sign
)))
69 ;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0
71 ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
73 ;; Arctan is computed using the Taylor series
75 ;; arctan(x) = x - x^3/3 + x^5/5 - x^7/7
76 (flet ((atan-series (x)
78 (eps (make-qd-d (scale-float 1d0 -
212)))
84 (loop while
(qd-> tmp eps
) do
88 (sub-qd s1
(div-qd tmp
(make-qd-d d
)))
89 (add-qd s1
(div-qd tmp
(make-qd-d d
)))))
91 (setf tmp
(mul-qd tmp r
))
92 (setf sign
(cl:- sign
)))
94 (let* ((x1 (div-qd +qd-one
+
100 (p (mul-qd-d (sub-qd (mul-qd-d s1
4d0
)
103 (format t
"~2&pi via Machin's atan formula~%")
104 (print-result p
+qd-pi
+)
108 (declare (optimize (speed 3)))
109 ;; Salamin-Brent Quadratic formula for pi
111 (b (sqrt-qd (make-qd-d 0.5d0
)))
112 (s (make-qd-d 0.5d0
))
114 (p (div-qd (mul-qd-d (sqr-qd a
) 2d0
)
116 (declare (double-float m
))
119 (let* ((a-new (mul-qd-d (add-qd a b
) .5d0
))
120 (b-new (sqrt-qd (mul-qd a b
)))
122 (mul-qd-d (sub-qd (sqr-qd a-new
)
128 (setf p
(div-qd (mul-qd-d (sqr-qd a
) 2d0
)
130 (format t
"~2&Salamin-Brent Quadratic formula for pi~%")
131 (print-result p
+qd-pi
+)
135 (declare (optimize (speed 3)))
136 ;; Borwein Quartic formula for pi
137 (let* ((a (sub-qd (make-qd-d 6d0
)
138 (mul-qd-d (sqrt-qd (make-qd-d 2d0
))
140 (y (sub-qd-d (sqrt-qd (make-qd-d 2d0
))
145 (declare (double-float m
))
148 (let ((r (nroot-qd (sub-qd +qd-one
+
151 (setf y
(div-qd (sub-d-qd 1d0
155 (setf a
(sub-qd (mul-qd a
156 (sqr-qd (sqr-qd (add-qd-d y
1d0
))))
158 (add-qd-d (add-qd y
(sqr-qd y
))
161 (setf p
(div-qd +qd-one
+
163 (format t
"~2&Borwein's Quartic formula for pi~%")
164 (print-result p
+qd-pi
+)
168 ;; 2.718281828459045235360287471352662497757247093699959574966967627724L0
170 ;; Taylor series for e
171 (let ((s (make-qd-d 2d0
))
176 (loop while
(qd-> tmp
(make-qd-d 1d-100
)) do
179 (setf tmp
(div-qd tmp
180 (make-qd-d (float n
1d0
))))
181 (setf s
(add-qd s tmp
)))
182 (format t
"~2&e via Taylor series~%")
183 (print-result s
+qd-e
+)
187 ;; 0.6931471805599453094172321214581765680755001343602552541206800094934L0
189 ;; Taylor series for log 2
191 ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
193 ;; with x = 1/2 to get log(1/2) = -log(2)
194 (let ((s (make-qd-d .5d0
))
195 (tt (make-qd-d .5d0
))
198 (loop while
(qd-> tt
(make-qd-d 1d-100
)) do
201 (setf tt
(mul-qd-d tt
.5d0
))
203 (div-qd tt
(make-qd-d (float n
1d0
))))))
204 (format t
"~2&log(2) via Taylor series~%")
205 (print-result s
+qd-log2
+)
208 (defun test-atan (&optional
(fun #'atan-qd
))
209 ;; Compute atan for known values
211 (format t
"~2&atan via ~S~%" fun
)
212 ;; atan(1/sqrt(3)) = pi/6
213 (let* ((arg (div-qd +qd-one
+ (sqrt-qd (make-qd-d 3d0
))))
214 (y (div-qd (funcall fun arg
) +qd-pi
+))
215 (true (div-qd +qd-one
+ (make-qd-d 6d0
))))
216 (format t
"atan(1/sqrt(3))/pi = ~/qdi::qd-format/~%" y
)
217 (format t
"1/6 = ~/qdi::qd-format/~%" true
)
218 (format t
"bits = ~,1f~%"
219 (bit-accuracy y true
)))
220 ;; atan(sqrt(3)) = %pi/3
221 (let* ((arg (sqrt-qd (make-qd-d 3d0
)))
222 (y (div-qd (funcall fun arg
) +qd-pi
+))
223 (true (div-qd +qd-one
+ (make-qd-d 3d0
))))
224 (format t
"atan(sqrt(3))/pi = ~/qdi::qd-format/~%" y
)
225 (format t
"1/3 = ~/qdi::qd-format/~%" true
)
226 (format t
"bits = ~,1f~%"
227 (bit-accuracy y true
)))
229 (let* ((arg +qd-one
+)
230 (y (div-qd (funcall fun arg
) +qd-pi
+))
231 (true (div-qd +qd-one
+ (make-qd-d 4d0
))))
232 (format t
"atan(1)/pi = ~/qdi::qd-format/~%" y
)
233 (format t
"1/4 = ~/qdi::qd-format/~%" true
)
234 (format t
"bits = ~,1f~%"
235 (bit-accuracy y true
))))
238 (format t
"~2&sin~%")
239 (let* ((arg (div-qd +qd-pi
+ (make-qd-d 6d0
)))
241 (true (make-qd-d 0.5d0
)))
242 (format t
"sin(pi/6) = ~/qdi::qd-format/~%" y
)
243 (format t
"1/2 = ~/qdi::qd-format/~%" true
)
244 (format t
"bits = ~,1f~%"
245 (bit-accuracy y true
)))
246 (let* ((arg (div-qd +qd-pi
+ (make-qd-d 4d0
)))
248 (true (sqrt-qd (make-qd-d 0.5d0
))))
249 (format t
"sin(pi/4) = ~/qdi::qd-format/~%" y
)
250 (format t
"1/sqrt(2) = ~/qdi::qd-format/~%" true
)
251 (format t
"bits = ~,1f~%"
252 (bit-accuracy y true
)))
253 (let* ((arg (div-qd +qd-pi
+ (make-qd-d 3d0
)))
255 (true (div-qd (sqrt-qd (make-qd-d 3d0
)) (make-qd-d 2d0
))))
256 (format t
"sin(pi/3) = ~/qdi::qd-format/~%" y
)
257 (format t
"sqrt(3)/2 = ~/qdi::qd-format/~%" true
)
258 (format t
"bits = ~,1f~%"
259 (bit-accuracy y true
)))
262 (defun test-tan (&optional
(f #'tan-qd
))
263 (format t
"~2&tan via ~S~%" f
)
264 (let* ((arg (div-qd +qd-pi
+ (make-qd-d 6d0
)))
266 (true (div-qd +qd-one
+ (sqrt-qd (make-qd-d 3d0
)))))
267 (format t
"tan(pi/6) = ~/qdi::qd-format/~%" y
)
268 (format t
"1/sqrt(3) = ~/qdi::qd-format/~%" true
)
269 (format t
"bits = ~,1f~%"
270 (bit-accuracy y true
)))
271 (let* ((arg (div-qd +qd-pi
+ (make-qd-d 4d0
)))
274 (format t
"tan(pi/4) = ~/qdi::qd-format/~%" y
)
275 (format t
"1 = ~/qdi::qd-format/~%" true
)
276 (format t
"bits = ~,1f~%"
277 (bit-accuracy y true
)))
278 (let* ((arg (div-qd +qd-pi
+ (make-qd-d 3d0
)))
280 (true (sqrt-qd (make-qd-d 3d0
))))
281 (format t
"tan(pi/3) = ~/qdi::qd-format/~%" y
)
282 (format t
"sqrt(3) = ~/qdi::qd-format/~%" true
)
283 (format t
"bits = ~,1f~%"
284 (bit-accuracy y true
)))
288 (format t
"~2&asin~%")
289 (let* ((arg (make-qd-d 0.5d0
))
291 (true (div-qd +qd-pi
+ (make-qd-d 6d0
))))
292 (format t
"asin(1/2) = ~/qdi::qd-format/~%" y
)
293 (format t
"pi/6 = ~/qdi::qd-format/~%" true
)
294 (format t
"bits = ~,1f~%"
295 (bit-accuracy y true
)))
296 (let* ((arg (sqrt-qd (make-qd-d 0.5d0
)))
298 (true (div-qd +qd-pi
+ (make-qd-d 4d0
))))
299 (format t
"asin(1/sqrt(2))= ~/qdi::qd-format/~%" y
)
300 (format t
"pi/4 = ~/qdi::qd-format/~%" true
)
301 (format t
"bits = ~,1f~%"
302 (bit-accuracy y true
)))
303 (let* ((arg (div-qd (sqrt-qd (make-qd-d 3d0
)) (make-qd-d 2d0
)))
305 (true (div-qd +qd-pi
+ (make-qd-d 3d0
))))
306 (format t
"asin(sqrt(3)/2)= ~/qdi::qd-format/~%" y
)
307 (format t
"pi/3 = ~/qdi::qd-format/~%" true
)
308 (format t
"bits = ~,1f~%"
309 (bit-accuracy y true
)))
313 (format t
"~2&Log via ~A~%" f
)
314 (let* ((arg (make-qd-d 2d0
))
317 (format t
"log(2) = ~/qdi::qd-format/~%" y
)
318 (format t
"log(2) = ~/qdi::qd-format/~%" true
)
319 (format t
"bits = ~,1f~%"
320 (bit-accuracy y true
)))
321 (let* ((arg (make-qd-d 10d0
))
324 (format t
"log(10) = ~/qdi::qd-format/~%" y
)
325 (format t
"log(10) = ~/qdi::qd-format/~%" true
)
326 (format t
"bits = ~,1f~%"
327 (bit-accuracy y true
)))
328 (let* ((arg (add-d-qd 1d0
(scale-float-qd (make-qd-d 1d0
) -
80)))
330 (true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
331 (format t
"log(1+2^-80) = ~/qdi::qd-format/~%" y
)
332 (format t
"log(1+2^-80) = ~/qdi::qd-format/~%" true
)
333 (format t
"bits = ~,1f~%"
334 (bit-accuracy y true
)))
337 (defun test-log1p (f)
338 (format t
"~2&Log1p via ~A~%" f
)
339 (let* ((arg (make-qd-d 1d0
))
342 (format t
"log1p(1) = ~/qdi::qd-format/~%" y
)
343 (format t
"log(2) = ~/qdi::qd-format/~%" true
)
344 (format t
"bits = ~,1f~%"
345 (bit-accuracy y true
)))
346 (let* ((arg (make-qd-d 9d0
))
348 (true (qd-from-string "2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0")))
349 (format t
"log1p(9) = ~/qdi::qd-format/~%" y
)
350 (format t
"log(10) = ~/qdi::qd-format/~%" true
)
351 (format t
"bits = ~,1f~%"
352 (bit-accuracy y true
)))
353 (let* ((arg (scale-float-qd (make-qd-d 1d0
) -
80))
355 (true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
356 (format t
"log1p(2^-80) = ~/qdi::qd-format/~%" y
)
357 (format t
"log(1+2^-80) = ~/qdi::qd-format/~%" true
)
358 (format t
"bits = ~,1f~%"
359 (bit-accuracy y true
)))
363 (format t
"~2&Exp via ~A~%" f
)
364 (let* ((arg +qd-zero
+)
367 (format t
"exp(0)-1 = ~/qdi::qd-format/~%" y
)
368 (format t
"0 = ~/qdi::qd-format/~%" true
)
369 (format t
"bits = ~,1f~%"
370 (bit-accuracy y true
)))
371 (let* ((arg +qd-one
+)
373 (true (qd-from-string "1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0")))
374 (format t
"exp(1)-1 = ~/qdi::qd-format/~%" y
)
375 (format t
"e-1 = ~/qdi::qd-format/~%" true
)
376 (format t
"bits = ~,1f~%"
377 (bit-accuracy y true
)))
379 (let* ((arg (scale-float-qd +qd-one
+ -
100))
381 (true (qd-from-string "7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31")))
382 (format t
"exp(2^-80)-1 = ~/qdi::qd-format/~%" y
)
383 (format t
"exp(2^-80)-1 = ~/qdi::qd-format/~%" true
)
384 (format t
"bits = ~,1f~%"
385 (bit-accuracy y true
)))
396 (dolist (f '(atan-qd/newton atan-qd
/cordic atan-qd
/duplication
))
398 (dolist (f '(tan-qd/sincos tan-qd
/cordic
))
400 (dolist (f '(log-qd/newton log-qd
/agm log-qd
/agm2 log-qd
/agm3 log-qd
/halley
))
402 (dolist (f '(log1p-qd/duplication
))
404 (dolist (f (list #'(lambda (x)
405 (sub-qd-d (exp-qd/reduce x
) 1d0
))
407 #'expm1-qd
/duplication
))