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 (defun bit-accuracy (est true
)
30 (let* ((diff (abs (- est true
)))
31 (err (float (if (zerop true
)
39 (defun print-result (est true
)
40 (format t
"est: ~A~%" est
)
41 (format t
"tru: ~A~%" true
)
42 (format t
"err: ~A~%" (float (- est true
) 1d0
))
43 (format t
"bits: ~,1f~%" (bit-accuracy est true
)))
46 (make-instance 'qd-real
:value qdi
::+qd-e
+))
49 (make-instance 'qd-real
:value qdi
::+qd-log2
+))
52 ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
54 ;; Arctan is computed using the Taylor series
56 ;; arctan(x) = x - x^3/3 + x^5/5 - x^7/7
57 (flet ((atan-series (x)
59 (eps (float (scale-float 1d0 -
212) #q1
))
65 (loop while
(> tmp eps
) do
82 (format t
"~2&pi via Machin's atan formula~%")
87 (declare (optimize (speed 3)))
88 ;; Salamin-Brent Quadratic formula for pi
96 (declare (double-float m
))
99 (let* ((a-new (* (+ a b
) .5d0
))
100 (b-new (sqrt (* a b
)))
102 (* (- (* a-new a-new
)
108 (setf p
(/ (* (* a a
) 2d0
)
110 (format t
"~2&Salamin-Brent Quadratic formula for pi~%")
111 (print-result p
+pi
+)
115 (declare (optimize (speed 3)))
116 ;; Borwein Quartic formula for pi
124 (declare (double-float m
))
127 (let ((r (expt (- 1 (expt y
4))
138 (format t
"~2&Borwein's Quartic formula for pi~%")
139 (print-result p
+pi
+)
143 ;; Taylor series for e
149 (loop while
(> tmp
1d-100
) do
154 (format t
"~2&e via Taylor series~%")
159 ;; Taylor series for log 2
161 ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
163 ;; with x = 1/2 to get log(1/2) = -log(2)
168 (loop while
(> tt
1d-100
) do
171 (setf tt
(* tt
.5d0
))
174 (format t
"~2&log(2) via Taylor series~%")
175 (print-result s
+log2
+)
179 (let* ((arg (/ (sqrt #q3
)))
180 (y (/ (atan arg
) +pi
+))
182 (format t
"~2&atan for special args~%")
183 (format t
"atan(1/sqrt(3))/pi = 1/6~%")
184 (print-result y true
))
185 ;; atan(sqrt(3)) = %pi/3
186 (let* ((arg (sqrt #q3
))
187 (y (/ (atan arg
) +pi
+))
189 (format t
"atan(sqrt(3))/pi = 1/3~%")
190 (print-result y true
))
193 (y (/ (atan arg
) +pi
+))
195 (format t
"atan(1)/pi = 1/4~%")
196 (print-result y true
))
198 (y (/ (atan arg
) +pi
+))
200 (format t
"atan(1q100)/pi = 1/2~%")
201 (print-result y true
))
202 (let* ((arg #q-1q100
)
203 (y (/ (atan arg
) +pi
+))
205 (format t
"atan(-1q100)/pi = -1/2~%")
206 (print-result y true
)))
209 (format t
"~2&sin for special args~%")
210 (let* ((arg (/ +pi
+ 6))
213 (format t
"sin(pi/6) = 1/2~%")
214 (print-result y true
))
215 (let* ((arg (/ +pi
+ 4))
218 (format t
"sin(pi/4) = 1/sqrt(2)~%")
219 (print-result y true
))
220 (let* ((arg (/ +pi
+ 3))
222 (true (/ (sqrt #q3
) 2)))
223 (format t
"sin(pi/3) = sqrt(3)/2~%")
224 (print-result y true
)))
227 (format t
"~2&tan for special args~%")
228 (let* ((arg (/ +pi
+ 6))
230 (true (/ (sqrt #q3
))))
231 (format t
"tan(pi/6) = 1/sqrt(3)~%")
232 (print-result y true
))
233 (let* ((arg (/ +pi
+ 4))
236 (format t
"tan(pi/4) = 1~%")
237 (print-result y true
))
238 (let* ((arg (/ +pi
+ 3))
241 (format t
"tan(pi/3) = sqrt(3)~%")
242 (print-result y true
)))
245 (format t
"~2&asin for special args~%")
249 (format t
"asin(1/2) = pi/6~%")
250 (print-result y true
))
251 (let* ((arg (sqrt #q
.5))
254 (format t
"asin(1/sqrt(2) = pi/4~%")
255 (print-result y true
))
256 (let* ((arg (/ (sqrt #q3
) 2))
259 (format t
"asin(sqrt(3)/2) = pi/3~%")
260 (print-result y true
)))
263 (format t
"~2&Log for special args~%")
267 (format t
"log(2)~%")
268 (print-result y true
))
271 (true (make-instance 'qd-real
:value qdi
::+qd-log10
+)))
272 (format t
"log(10)~%")
273 (print-result y true
))
274 (let* ((arg (+ 1 (scale-float #q1 -
80)))
276 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25
))
277 (format t
"log(1+2^-80)~%")
278 (print-result y true
)))
281 (format t
"~2&Sqrt for special args~%")
282 (dolist (f '((#q1q200
#q1q100
)
286 (#q1q-250
#q1q-125
)))
287 (destructuring-bind (arg true
)
289 (let ((y (sqrt arg
)))
290 (format t
"sqrt(~/qdi::qd-format/)~%" (qd-value arg
))
291 (print-result y true
)))))