/[oct]/oct/qd-test.lisp
ViewVC logotype

Contents of /oct/qd-test.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.21 - (show annotations)
Mon Oct 15 18:53:44 2007 UTC (6 years, 6 months ago) by rtoy
Branch: MAIN
CVS Tags: GIT_CONVERSION, three-arg-branch-2007-11-07, THREE-ARG-BRANCH-BASE, three-arg-branch-merged-2007-11-07, HEAD
Branch point for: THREE-ARG-BRANCH
Changes since 1.20: +1 -1 lines
o Oops.  Fix up a few IN-PACKAGE's for the new package names.

qd-fun.lisp:
o Comment out the old sin/cos routines
o Fix a few mistakes in accurate-sincos-qd
o Rename accurate-sincos-qd to sincos-qd.
1 ;;;; -*- Mode: lisp -*-
2 ;;;;
3 ;;;; Copyright (c) 2007 Raymond Toy
4 ;;;;
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
12 ;;;; conditions:
13 ;;;;
14 ;;;; The above copyright notice and this permission notice shall be
15 ;;;; included in all copies or substantial portions of the Software.
16 ;;;;
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.
25
26
27 (in-package #:octi)
28
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)
34 (qd-0 diff)
35 (cl:/ (qd-0 diff) (abs (qd-0 true))))))
36 (if (zerop (qd-0 diff))
37 t
38 (cl:- (log err 2d0)))))
39
40 (defun print-result (est true)
41 (format t "est: ~/octi::qd-format/~%" est)
42 (format t "tru: ~/octi::qd-format/~%" true)
43 (format t "err: ~A~%" (qd-0 (sub-qd est true)))
44 (format t "bits: ~,1f~%" (bit-accuracy est true)))
45
46 ;; Machin's formula for pi
47 #+nil
48 (defun atan-series (x)
49 (let* ((d 1d0)
50 (eps (make-qd-d (scale-float 1d0 -212)
51 0d0 0d0 0d0))
52 (tmp x)
53 (r (sqr-qd tmp))
54 (s1 (make-qd-dd 0w0 0w0))
55 (k 0)
56 (sign 1))
57 (loop while (qd-> tmp eps) do
58 (incf k)
59 (setf s1
60 (if (minusp sign)
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)))))
63 (incf d 2d0)
64 (setf tmp (mul-qd tmp r))
65 (setf sign (cl:- sign)))
66 s1))
67
68 ;; pi =
69 ;; 3.1415926535897932384626433832795028841971693993751058209749445923078L0
70 (defun test2 (&optional (printp t))
71 ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
72 ;;
73 ;; Arctan is computed using the Taylor series
74 ;;
75 ;; arctan(x) = x - x^3/3 + x^5/5 - x^7/7
76 (flet ((atan-series (x)
77 (let* ((d 1d0)
78 (eps (make-qd-d (scale-float 1d0 -212)))
79 (tmp x)
80 (r (sqr-qd tmp))
81 (s1 (make-qd-d 0d0))
82 (k 0)
83 (sign 1))
84 (loop while (qd-> tmp eps) do
85 (incf k)
86 (setf s1
87 (if (minusp sign)
88 (sub-qd s1 (div-qd tmp (make-qd-d d)))
89 (add-qd s1 (div-qd tmp (make-qd-d d)))))
90 (incf d 2d0)
91 (setf tmp (mul-qd tmp r))
92 (setf sign (cl:- sign)))
93 s1)))
94 (let* ((x1 (div-qd +qd-one+
95 (make-qd-d 5d0)))
96 (s1 (atan-series x1))
97 (x2 (div-qd +qd-one+
98 (make-qd-d 239d0)))
99 (s2 (atan-series x2))
100 (p (mul-qd-d (sub-qd (mul-qd-d s1 4d0)
101 s2)
102 4d0)))
103 (when printp
104 (format t "~2&pi via Machin's atan formula~%")
105 (print-result p +qd-pi+))
106 p)))
107
108 (defun test3 (&optional (printp t))
109 (declare (optimize (speed 3)))
110 ;; Salamin-Brent Quadratic formula for pi
111 (let* ((a +qd-one+)
112 (b (sqrt-qd (make-qd-d 0.5d0)))
113 (s (make-qd-d 0.5d0))
114 (m 1d0)
115 (p (div-qd (mul-qd-d (sqr-qd a) 2d0)
116 s)))
117 (declare (double-float m))
118 (dotimes (k 9)
119 (setf m (cl:* 2 m))
120 (let* ((a-new (mul-qd-d (add-qd a b) .5d0))
121 (b-new (sqrt-qd (mul-qd a b)))
122 (s-new (sub-qd s
123 (mul-qd-d (sub-qd (sqr-qd a-new)
124 (sqr-qd b-new))
125 m))))
126 (setf a a-new)
127 (setf b b-new)
128 (setf s s-new)
129 (setf p (div-qd (mul-qd-d (sqr-qd a) 2d0)
130 s))))
131 (when printp
132 (format t "~2&Salamin-Brent Quadratic formula for pi~%")
133 (print-result p +qd-pi+))
134 p))
135
136 (defun test4 (&optional (printp t))
137 (declare (optimize (speed 3)))
138 ;; Borwein Quartic formula for pi
139 (let* ((a (sub-qd (make-qd-d 6d0)
140 (mul-qd-d (sqrt-qd (make-qd-d 2d0))
141 4d0)))
142 (y (sub-qd-d (sqrt-qd (make-qd-d 2d0))
143 1d0))
144 (m 2d0)
145 (p (div-qd +qd-one+
146 a)))
147 (declare (double-float m))
148 (dotimes (k 9)
149 (setf m (cl:* 4 m))
150 (let ((r (nroot-qd (sub-qd +qd-one+
151 (sqr-qd (sqr-qd y)))
152 4)))
153 (setf y (div-qd (sub-d-qd 1d0
154 r)
155 (add-d-qd 1d0
156 r)))
157 (setf a (sub-qd (mul-qd a
158 (sqr-qd (sqr-qd (add-qd-d y 1d0))))
159 (mul-qd-d (mul-qd y
160 (add-qd-d (add-qd y (sqr-qd y))
161 1d0))
162 m)))
163 (setf p (div-qd +qd-one+
164 a))))
165 (when printp
166 (format t "~2&Borwein's Quartic formula for pi~%")
167 (print-result p +qd-pi+))
168 p))
169
170 ;; e =
171 ;; 2.718281828459045235360287471352662497757247093699959574966967627724L0
172 (defun test5 (&optional (printp t))
173 ;; Taylor series for e
174 (let ((s (make-qd-d 2d0))
175 (tmp +qd-one+)
176 (n 1d0)
177 (delta 0d0)
178 (i 0))
179 (loop while (qd-> tmp (make-qd-d 1d-100)) do
180 (incf i)
181 (incf n)
182 (setf tmp (div-qd tmp
183 (make-qd-d (float n 1d0))))
184 (setf s (add-qd s tmp)))
185 (when printp
186 (format t "~2&e via Taylor series~%")
187 (print-result s +qd-e+))
188 s))
189
190 ;; log(2) =
191 ;; 0.6931471805599453094172321214581765680755001343602552541206800094934L0
192 (defun test6 (&optional (printp t))
193 ;; Taylor series for log 2
194 ;;
195 ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
196 ;;
197 ;; with x = 1/2 to get log(1/2) = -log(2)
198 (let ((s (make-qd-d .5d0))
199 (tt (make-qd-d .5d0))
200 (n 1d0)
201 (i 0))
202 (loop while (qd-> tt (make-qd-d 1d-100)) do
203 (incf i)
204 (incf n)
205 (setf tt (mul-qd-d tt .5d0))
206 (setf s (add-qd s
207 (div-qd tt (make-qd-d (float n 1d0))))))
208 (when printp
209 (format t "~2&log(2) via Taylor series~%")
210 (print-result s +qd-log2+))
211 s))
212
213 (defun test-atan (&optional (fun #'atan-qd) (printp t))
214 ;; Compute atan for known values
215
216 (format t "~2&atan via ~S~%" fun)
217 ;; atan(1/sqrt(3)) = pi/6
218 (let* ((arg (div-qd +qd-one+ (sqrt-qd (make-qd-d 3d0))))
219 (y (div-qd (funcall fun arg) +qd-pi+))
220 (true (div-qd +qd-one+ (make-qd-d 6d0))))
221 (format t "atan(1/sqrt(3))/pi = ~/octi::qd-format/~%" y)
222 (format t "1/6 = ~/octi::qd-format/~%" true)
223 (format t "bits = ~,1f~%"
224 (bit-accuracy y true)))
225 ;; atan(sqrt(3)) = %pi/3
226 (let* ((arg (sqrt-qd (make-qd-d 3d0)))
227 (y (div-qd (funcall fun arg) +qd-pi+))
228 (true (div-qd +qd-one+ (make-qd-d 3d0))))
229 (format t "atan(sqrt(3))/pi = ~/octi::qd-format/~%" y)
230 (format t "1/3 = ~/octi::qd-format/~%" true)
231 (format t "bits = ~,1f~%"
232 (bit-accuracy y true)))
233 ;; atan(1) = %pi/4
234 (let* ((arg +qd-one+)
235 (y (div-qd (funcall fun arg) +qd-pi+))
236 (true (div-qd +qd-one+ (make-qd-d 4d0))))
237 (format t "atan(1)/pi = ~/octi::qd-format/~%" y)
238 (format t "1/4 = ~/octi::qd-format/~%" true)
239 (format t "bits = ~,1f~%"
240 (bit-accuracy y true))))
241
242 (defun test-sin (&optional (printp t))
243 (format t "~2&sin~%")
244 (let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0)))
245 (y (sin-qd arg))
246 (true (make-qd-d 0.5d0)))
247 (format t "sin(pi/6) = ~/octi::qd-format/~%" y)
248 (format t "1/2 = ~/octi::qd-format/~%" true)
249 (format t "bits = ~,1f~%"
250 (bit-accuracy y true)))
251 (let* ((arg (div-qd +qd-pi+ (make-qd-d 4d0)))
252 (y (sin-qd arg))
253 (true (sqrt-qd (make-qd-d 0.5d0))))
254 (format t "sin(pi/4) = ~/octi::qd-format/~%" y)
255 (format t "1/sqrt(2) = ~/octi::qd-format/~%" true)
256 (format t "bits = ~,1f~%"
257 (bit-accuracy y true)))
258 (let* ((arg (div-qd +qd-pi+ (make-qd-d 3d0)))
259 (y (sin-qd arg))
260 (true (div-qd (sqrt-qd (make-qd-d 3d0)) (make-qd-d 2d0))))
261 (format t "sin(pi/3) = ~/octi::qd-format/~%" y)
262 (format t "sqrt(3)/2 = ~/octi::qd-format/~%" true)
263 (format t "bits = ~,1f~%"
264 (bit-accuracy y true)))
265 )
266
267 (defun test-tan (&optional (f #'tan-qd) (printp t))
268 (format t "~2&tan via ~S~%" f)
269 (let* ((arg (div-qd +qd-pi+ (make-qd-d 6d0)))
270 (y (funcall f arg))
271 (true (div-qd +qd-one+ (sqrt-qd (make-qd-d 3d0)))))
272 (format t "tan(pi/6) = ~/octi::qd-format/~%" y)
273 (format t "1/sqrt(3) = ~/octi::qd-format/~%" true)
274 (format t "bits = ~,1f~%"
275 (bit-accuracy y true)))
276 (let* ((arg (div-qd +qd-pi+ (make-qd-d 4d0)))
277 (y (funcall f arg))
278 (true +qd-one+))
279 (format t "tan(pi/4) = ~/octi::qd-format/~%" y)
280 (format t "1 = ~/octi::qd-format/~%" true)
281 (format t "bits = ~,1f~%"
282 (bit-accuracy y true)))
283 (let* ((arg (div-qd +qd-pi+ (make-qd-d 3d0)))
284 (y (funcall f arg))
285 (true (sqrt-qd (make-qd-d 3d0))))
286 (format t "tan(pi/3) = ~/octi::qd-format/~%" y)
287 (format t "sqrt(3) = ~/octi::qd-format/~%" true)
288 (format t "bits = ~,1f~%"
289 (bit-accuracy y true)))
290 )
291
292 (defun test-asin (&optional (printp t))
293 (format t "~2&asin~%")
294 (let* ((arg (make-qd-d 0.5d0))
295 (y (asin-qd arg))
296 (true (div-qd +qd-pi+ (make-qd-d 6d0))))
297 (format t "asin(1/2) = ~/octi::qd-format/~%" y)
298 (format t "pi/6 = ~/octi::qd-format/~%" true)
299 (format t "bits = ~,1f~%"
300 (bit-accuracy y true)))
301 (let* ((arg (sqrt-qd (make-qd-d 0.5d0)))
302 (y (asin-qd arg))
303 (true (div-qd +qd-pi+ (make-qd-d 4d0))))
304 (format t "asin(1/sqrt(2))= ~/octi::qd-format/~%" y)
305 (format t "pi/4 = ~/octi::qd-format/~%" true)
306 (format t "bits = ~,1f~%"
307 (bit-accuracy y true)))
308 (let* ((arg (div-qd (sqrt-qd (make-qd-d 3d0)) (make-qd-d 2d0)))
309 (y (asin-qd arg))
310 (true (div-qd +qd-pi+ (make-qd-d 3d0))))
311 (format t "asin(sqrt(3)/2)= ~/octi::qd-format/~%" y)
312 (format t "pi/3 = ~/octi::qd-format/~%" true)
313 (format t "bits = ~,1f~%"
314 (bit-accuracy y true)))
315 )
316
317 (defun test-log (f &optional (printp t))
318 (format t "~2&Log via ~A~%" f)
319 (let* ((arg (make-qd-d 2d0))
320 (y (funcall f arg))
321 (true +qd-log2+))
322 (format t "log(2) = ~/octi::qd-format/~%" y)
323 (format t "log(2) = ~/octi::qd-format/~%" true)
324 (format t "bits = ~,1f~%"
325 (bit-accuracy y true)))
326 (let* ((arg (make-qd-d 10d0))
327 (y (funcall f arg))
328 (true +qd-log10+))
329 (format t "log(10) = ~/octi::qd-format/~%" y)
330 (format t "log(10) = ~/octi::qd-format/~%" true)
331 (format t "bits = ~,1f~%"
332 (bit-accuracy y true)))
333 (let* ((arg (add-d-qd 1d0 (scale-float-qd (make-qd-d 1d0) -80)))
334 (y (funcall f arg))
335 (true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
336 (format t "log(1+2^-80) = ~/octi::qd-format/~%" y)
337 (format t "log(1+2^-80) = ~/octi::qd-format/~%" true)
338 (format t "bits = ~,1f~%"
339 (bit-accuracy y true)))
340 )
341
342 (defun test-log1p (f &optional (printp t))
343 (format t "~2&Log1p via ~A~%" f)
344 (let* ((arg (make-qd-d 1d0))
345 (y (funcall f arg))
346 (true +qd-log2+))
347 (format t "log1p(1) = ~/octi::qd-format/~%" y)
348 (format t "log(2) = ~/octi::qd-format/~%" true)
349 (format t "bits = ~,1f~%"
350 (bit-accuracy y true)))
351 (let* ((arg (make-qd-d 9d0))
352 (y (funcall f arg))
353 (true (qd-from-string "2.3025850929940456840179914546843642076011014886287729760333279009675726096773525q0")))
354 (format t "log1p(9) = ~/octi::qd-format/~%" y)
355 (format t "log(10) = ~/octi::qd-format/~%" true)
356 (format t "bits = ~,1f~%"
357 (bit-accuracy y true)))
358 (let* ((arg (scale-float-qd (make-qd-d 1d0) -80))
359 (y (funcall f arg))
360 (true (qd-from-string "8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25")))
361 (format t "log1p(2^-80) = ~/octi::qd-format/~%" y)
362 (format t "log(1+2^-80) = ~/octi::qd-format/~%" true)
363 (format t "bits = ~,1f~%"
364 (bit-accuracy y true)))
365 )
366
367 (defun test-exp (f &optional (printp t))
368 (format t "~2&Exp via ~A~%" f)
369 (let* ((arg +qd-zero+)
370 (y (funcall f arg))
371 (true +qd-zero+))
372 (format t "exp(0)-1 = ~/octi::qd-format/~%" y)
373 (format t "0 = ~/octi::qd-format/~%" true)
374 (format t "bits = ~,1f~%"
375 (bit-accuracy y true)))
376 (let* ((arg +qd-one+)
377 (y (funcall f arg))
378 (true (qd-from-string "1.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135966290435729003342952q0")))
379 (format t "exp(1)-1 = ~/octi::qd-format/~%" y)
380 (format t "e-1 = ~/octi::qd-format/~%" true)
381 (format t "bits = ~,1f~%"
382 (bit-accuracy y true)))
383
384 (let* ((arg (scale-float-qd +qd-one+ -100))
385 (y (funcall f arg))
386 (true (qd-from-string "7.888609052210118054117285652830973804370994921943802079729680186943164342372119432861876389514693341738324702996270767390039172777809233288470357147q-31")))
387 (format t "exp(2^-80)-1 = ~/octi::qd-format/~%" y)
388 (format t "exp(2^-80)-1 = ~/octi::qd-format/~%" true)
389 (format t "bits = ~,1f~%"
390 (bit-accuracy y true)))
391
392 )
393
394 (defun all-tests ()
395 (test2)
396 (test3)
397 (test4)
398 (test5)
399 (test6)
400 (test-sin)
401 (test-asin)
402 (dolist (f '(atan-qd/newton atan-qd/cordic atan-qd/duplication))
403 (test-atan f))
404 (dolist (f '(tan-qd/sincos tan-qd/cordic))
405 (test-tan f))
406 (dolist (f '(log-qd/newton log-qd/agm log-qd/agm2 log-qd/agm3 log-qd/halley))
407 (test-log f))
408 (dolist (f '(log1p-qd/duplication))
409 (test-log1p f))
410 (dolist (f (list #'(lambda (x)
411 (sub-qd-d (exp-qd/reduce x) 1d0))
412 #'expm1-qd/series
413 #'expm1-qd/duplication))
414 (test-exp f))
415 )

  ViewVC Help
Powered by ViewVC 1.1.5