/[oct]/oct/tests.lisp
ViewVC logotype

Contents of /oct/tests.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (show annotations)
Mon Oct 15 18:21:47 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.2: +4 -4 lines
o Rename QUAD-DOUBLE-INTERNAL package to OCT-INTERNAL, with nickname
  OCTI instead of QDI.
o Rename OCT package to NET.COMMON-LISP.OCT, with a nickname of OCT
o Remove nickname of QD.  (Conflicts with other packages dealing with
  quad-doubles.)
o Update all uses of QDI: to OCTI:

qd-fun.lisp:
o Add REM-PI/2 to do a simpler computation if the arg is small
  enough.  Otherwise, use the accurate but expensive rem operation.
o Renamed ACCURATE-SIN-QD to SIN-QD, etc.
o Update SIN-QD etc to use REM-PI/2.
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 #:oct)
28
29 (defun bit-accuracy (est true)
30 (let* ((diff (abs (- est true)))
31 (err (float (if (zerop true)
32 diff
33 (/ diff (abs true)))
34 1d0)))
35 (if (zerop diff)
36 t
37 (- (log err 2)))))
38
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)))
44
45 (defconstant +e+
46 (make-instance 'qd-real :value octi::+qd-e+))
47
48 (defconstant +log2+
49 (make-instance 'qd-real :value octi::+qd-log2+))
50
51 (defun test2 ()
52 ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
53 ;;
54 ;; Arctan is computed using the Taylor series
55 ;;
56 ;; arctan(x) = x - x^3/3 + x^5/5 - x^7/7
57 (flet ((atan-series (x)
58 (let* ((d 1d0)
59 (eps (float (scale-float 1d0 -212) #q1))
60 (tmp x)
61 (r (* tmp tmp))
62 (s1 #q0)
63 (k 0)
64 (sign 1))
65 (loop while (> tmp eps) do
66 (incf k)
67 (setf s1
68 (if (minusp sign)
69 (- s1 (/ tmp d))
70 (+ s1 (/ tmp d))))
71 (incf d 2d0)
72 (setf tmp (* tmp r))
73 (setf sign (- sign)))
74 s1)))
75 (let* ((x1 (/ #q1 5))
76 (s1 (atan-series x1))
77 (x2 (/ #q1 239))
78 (s2 (atan-series x2))
79 (p (* (- (* s1 4)
80 s2)
81 4)))
82 (format t "~2&pi via Machin's atan formula~%")
83 (print-result p +pi+)
84 p)))
85
86 (defun test3 ()
87 (declare (optimize (speed 3)))
88 ;; Salamin-Brent Quadratic formula for pi
89 (let* ((a #q1)
90 (b (sqrt #q.5))
91 (s #q.5)
92 (m 1d0)
93 (p (/ (* (* a a)
94 2d0)
95 s)))
96 (declare (double-float m))
97 (dotimes (k 9)
98 (setf m (* 2 m))
99 (let* ((a-new (* (+ a b) .5d0))
100 (b-new (sqrt (* a b)))
101 (s-new (- s
102 (* (- (* a-new a-new)
103 (* b-new b-new))
104 m))))
105 (setf a a-new)
106 (setf b b-new)
107 (setf s s-new)
108 (setf p (/ (* (* a a) 2d0)
109 s))))
110 (format t "~2&Salamin-Brent Quadratic formula for pi~%")
111 (print-result p +pi+)
112 p))
113
114 (defun test4 ()
115 (declare (optimize (speed 3)))
116 ;; Borwein Quartic formula for pi
117 (let* ((a (- 6
118 (* (sqrt #q2)
119 4)))
120 (y (- (sqrt #q2)
121 1))
122 (m 2d0)
123 (p (/ a)))
124 (declare (double-float m))
125 (dotimes (k 9)
126 (setf m (* 4 m))
127 (let ((r (expt (- 1 (expt y 4))
128 1/4)))
129 (setf y (/ (- 1d0 r)
130 (+ 1d0 r)))
131 (setf a (- (* a
132 (expt (+ y 1d0) 4))
133 (* (* y
134 (+ (+ y (expt y 2))
135 1d0))
136 m)))
137 (setf p (/ a))))
138 (format t "~2&Borwein's Quartic formula for pi~%")
139 (print-result p +pi+)
140 p))
141
142 (defun test5 ()
143 ;; Taylor series for e
144 (let ((s #q2)
145 (tmp #q1)
146 (n 1d0)
147 (delta 0d0)
148 (i 0))
149 (loop while (> tmp 1d-100) do
150 (incf i)
151 (incf n)
152 (setf tmp (/ tmp n))
153 (setf s (+ s tmp)))
154 (format t "~2&e via Taylor series~%")
155 (print-result s +e+)
156 s))
157
158 (defun test6 ()
159 ;; Taylor series for log 2
160 ;;
161 ;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
162 ;;
163 ;; with x = 1/2 to get log(1/2) = -log(2)
164 (let ((s #q.5)
165 (tt #q.5)
166 (n 1d0)
167 (i 0))
168 (loop while (> tt 1d-100) do
169 (incf i)
170 (incf n)
171 (setf tt (* tt .5d0))
172 (setf s (+ s
173 (/ tt n))))
174 (format t "~2&log(2) via Taylor series~%")
175 (print-result s +log2+)
176 s))
177
178 (defun test-atan ()
179 (let* ((arg (/ (sqrt #q3)))
180 (y (/ (atan arg) +pi+))
181 (true (/ #q6)))
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+))
188 (true (/ #q3)))
189 (format t "atan(sqrt(3))/pi = 1/3~%")
190 (print-result y true))
191 ;; atan(1) = %pi/4
192 (let* ((arg #q1)
193 (y (/ (atan arg) +pi+))
194 (true (/ #q4)))
195 (format t "atan(1)/pi = 1/4~%")
196 (print-result y true))
197 (let* ((arg #q1q100)
198 (y (/ (atan arg) +pi+))
199 (true #q.5))
200 (format t "atan(1q100)/pi = 1/2~%")
201 (print-result y true))
202 (let* ((arg #q-1q100)
203 (y (/ (atan arg) +pi+))
204 (true #q-.5))
205 (format t "atan(-1q100)/pi = -1/2~%")
206 (print-result y true)))
207
208 (defun test-sin ()
209 (format t "~2&sin for special args~%")
210 (let* ((arg (/ +pi+ 6))
211 (y (sin arg))
212 (true #q.5))
213 (format t "sin(pi/6) = 1/2~%")
214 (print-result y true))
215 (let* ((arg (/ +pi+ 4))
216 (y (sin arg))
217 (true (sqrt #q.5)))
218 (format t "sin(pi/4) = 1/sqrt(2)~%")
219 (print-result y true))
220 (let* ((arg (/ +pi+ 3))
221 (y (sin arg))
222 (true (/ (sqrt #q3) 2)))
223 (format t "sin(pi/3) = sqrt(3)/2~%")
224 (print-result y true)))
225
226 (defun test-tan ()
227 (format t "~2&tan for special args~%")
228 (let* ((arg (/ +pi+ 6))
229 (y (tan arg))
230 (true (/ (sqrt #q3))))
231 (format t"tan(pi/6) = 1/sqrt(3)~%")
232 (print-result y true))
233 (let* ((arg (/ +pi+ 4))
234 (y (tan arg))
235 (true #q1))
236 (format t "tan(pi/4) = 1~%")
237 (print-result y true))
238 (let* ((arg (/ +pi+ 3))
239 (y (tan arg))
240 (true (sqrt #q3)))
241 (format t "tan(pi/3) = sqrt(3)~%")
242 (print-result y true)))
243
244 (defun test-asin ()
245 (format t "~2&asin for special args~%")
246 (let* ((arg #q.5)
247 (y (asin arg))
248 (true (/ +pi+ 6)))
249 (format t "asin(1/2) = pi/6~%")
250 (print-result y true))
251 (let* ((arg (sqrt #q.5))
252 (y (asin arg))
253 (true (/ +pi+ 4)))
254 (format t "asin(1/sqrt(2) = pi/4~%")
255 (print-result y true))
256 (let* ((arg (/ (sqrt #q3) 2))
257 (y (asin arg))
258 (true (/ +pi+ 3)))
259 (format t "asin(sqrt(3)/2) = pi/3~%")
260 (print-result y true)))
261
262 (defun test-log ()
263 (format t "~2&Log for special args~%")
264 (let* ((arg #q2)
265 (y (log arg))
266 (true +log2+))
267 (format t "log(2)~%")
268 (print-result y true))
269 (let* ((arg #q10)
270 (y (log arg))
271 (true (make-instance 'qd-real :value octi::+qd-log10+)))
272 (format t "log(10)~%")
273 (print-result y true))
274 (let* ((arg (+ 1 (scale-float #q1 -80)))
275 (y (log arg))
276 (true #q8.2718061255302767487140834995607996176476940491239977084112840149578911975528492q-25))
277 (format t "log(1+2^-80)~%")
278 (print-result y true)))
279
280 (defun test-sqrt ()
281 (format t "~2&Sqrt for special args~%")
282 (dolist (f '((#q1q200 #q1q100)
283 (#q1q300 #q1q150)
284 (#q1q308 #q1q154)
285 (#q1q-200 #q1q-100)
286 (#q1q-250 #q1q-125)))
287 (destructuring-bind (arg true)
288 f
289 (let ((y (sqrt arg)))
290 (format t "sqrt(~/octi::qd-format/)~%" (qd-value arg))
291 (print-result y true)))))
292
293 (defun all-tests ()
294 (test2)
295 (test3)
296 (test4)
297 (test5)
298 (test6)
299 (test-atan)
300 (test-sin)
301 (test-tan)
302 (test-asin)
303 (test-log)
304 (test-sqrt))

  ViewVC Help
Powered by ViewVC 1.1.5