[oct-cvs] Oct commit: oct branch-test.lisp tests.lisp timing.lisp
rtoy
rtoy at common-lisp.net
Sat Aug 25 18:49:27 UTC 2007
Update of /project/oct/cvsroot/oct
In directory clnet:/tmp/cvs-serv19661
Added Files:
branch-test.lisp tests.lisp timing.lisp
Log Message:
Initial revision.
--- /project/oct/cvsroot/oct/branch-test.lisp 2007/08/25 18:49:27 NONE
+++ /project/oct/cvsroot/oct/branch-test.lisp 2007/08/25 18:49:27 1.1
;;;; -*- Mode: lisp -*-
;;;;
;;;; Copyright (c) 2007 Raymond Toy
;;;;
;;;; Permission is hereby granted, free of charge, to any person
;;;; obtaining a copy of this software and associated documentation
;;;; files (the "Software"), to deal in the Software without
;;;; restriction, including without limitation the rights to use,
;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
;;;; copies of the Software, and to permit persons to whom the
;;;; Software is furnished to do so, subject to the following
;;;; conditions:
;;;;
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;;
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
;;; Some simple tests to see that we're computing the branch cuts
;;; correctly.
;;;
;;; NOTE: the tests assume that the functions for double-float are
;;; computing the values correctly for the branch cuts. We need to
;;; fix this.
(in-package #:qd)
(defun check-signs (fun arg real-sign imag-sign)
(let* ((z (funcall fun arg))
(x (realpart z))
(y (imagpart z)))
(unless (and (= (float-sign x) real-sign)
(= (float-sign y) imag-sign))
(format t "Sign of result doesn't match expected signs~%~
~& fun = ~A~
~& arg = ~A~
~& res = ~A~
~& expected = ~A ~A~%"
fun arg z real-sign imag-sign))))
(defun get-signs (z)
(values (float-sign (realpart z))
(float-sign (imagpart z))))
;; asin branch cut is the real axis |x| > 1. For x < -1, it is
;; continuous with quadrant II; for x > 1, continuous with quadrant
;; IV.
(defun test-asin ()
;; Check x < -1
(multiple-value-bind (tr ti)
(get-signs (asin #c(-2d0 +1d-20)))
(check-signs #'asin -2d0 tr ti)
(check-signs #'asin -2w0 tr ti)
(check-signs #'asin #q-2 tr ti)
(check-signs #'asin #c(-2d0 0d0) tr ti)
(check-signs #'asin #c(-2w0 0w0) tr ti)
(check-signs #'asin #q(-2 0) tr ti)
(check-signs #'asin #c(-2d0 -0d0) tr (- ti))
(check-signs #'asin #c(-2w0 -0w0) tr (- ti))
(check-signs #'asin #q(-2 #q-0q0) tr (- ti))
)
;; Check x > 1
(multiple-value-bind (tr ti)
(get-signs (asin #c(2d0 -1d-20)))
(check-signs #'asin 2d0 tr ti)
(check-signs #'asin 2w0 tr ti)
(check-signs #'asin #q2 tr ti)
(check-signs #'asin #c(2d0 -0d0) tr ti)
(check-signs #'asin #c(2w0 -0w0) tr ti)
(check-signs #'asin #q(2 #q-0q0) tr ti)))
;; acos branch cut is the real axis, |x| > 1. For x < -1, it is
;; continuous with quadrant II; for x > 1, quadrant IV.
(defun test-acos ()
;; Check x < -1
(multiple-value-bind (tr ti)
(get-signs (acos #c(-2d0 +1d-20)))
(check-signs #'acos -2d0 tr ti)
(check-signs #'acos -2w0 tr ti)
(check-signs #'acos #q-2 tr ti))
;; Check x > 1
(multiple-value-bind (tr ti)
(get-signs (acos #c(2d0 -1d-20)))
(check-signs #'acos 2d0 tr ti)
(check-signs #'acos 2w0 tr ti)
(check-signs #'acos #q2 tr ti)))
;; atan branch cut is the imaginary axis, |y| > 1. For y < -1, it is
;; continuous with quadrant IV; for x > 1, quadrant II.
(defun test-atan ()
;; Check y < -1
(multiple-value-bind (tr ti)
(get-signs (atan #c(1d-20 -2d0)))
(check-signs #'atan #c(0d0 -2d0) tr ti)
(check-signs #'atan #c(0w0 -2w0) tr ti)
(check-signs #'atan #q(#q0 #q-2) tr ti))
;; Check y > 1
(multiple-value-bind (tr ti)
(get-signs (atan #c(-1d-20 2d0)))
(check-signs #'atan #c(-0d0 2d0) tr ti)
(check-signs #'atan #c(-0w0 2w0) tr ti)
(check-signs #'atan #q(#q-0 2) tr ti)))
(defun test-atanh ()
;; Check x < -1
(multiple-value-bind (tr ti)
(get-signs (atanh #c(-2d0 -1d-20)))
(check-signs #'atanh -2d0 tr ti)
(check-signs #'atanh -2w0 tr ti)
(check-signs #'atanh #q-2 tr ti))
;; Check x > 1
(multiple-value-bind (tr ti)
(get-signs (atanh #c(2d0 1d-20)))
(check-signs #'atanh 2d0 tr ti)
(check-signs #'atanh 2w0 tr ti)
(check-signs #'atanh #q2 tr ti)))
--- /project/oct/cvsroot/oct/tests.lisp 2007/08/25 18:49:27 NONE
+++ /project/oct/cvsroot/oct/tests.lisp 2007/08/25 18:49:27 1.1
;;;; -*- Mode: lisp -*-
;;;;
;;;; Copyright (c) 2007 Raymond Toy
;;;;
;;;; Permission is hereby granted, free of charge, to any person
;;;; obtaining a copy of this software and associated documentation
;;;; files (the "Software"), to deal in the Software without
;;;; restriction, including without limitation the rights to use,
;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell
;;;; copies of the Software, and to permit persons to whom the
;;;; Software is furnished to do so, subject to the following
;;;; conditions:
;;;;
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;;
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
;;;; OTHER DEALINGS IN THE SOFTWARE.
(in-package #:qd)
(defun bit-accuracy (est true)
(let* ((diff (abs (- est true)))
(err (float (if (zerop true)
diff
(/ diff (abs true)))
1d0)))
(if (zerop diff)
t
(- (log err 2)))))
(defun print-result (est true)
(format t "est: ~A~%" est)
(format t "tru: ~A~%" true)
(format t "err: ~A~%" (float (- est true) 1d0))
(format t "bits: ~,1f~%" (bit-accuracy est true)))
(defconstant +e+
(make-instance 'qd-real :value qdi::+qd-e+))
(defconstant +log2+
(make-instance 'qd-real :value qdi::+qd-log2+))
(defun test2 ()
;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
;;
;; Arctan is computed using the Taylor series
;;
;; arctan(x) = x - x^3/3 + x^5/5 - x^7/7
(flet ((atan-series (x)
(let* ((d 1d0)
(eps (float (scale-float 1d0 -212) #q1))
(tmp x)
(r (* tmp tmp))
(s1 #q0)
(k 0)
(sign 1))
(loop while (> tmp eps) do
(incf k)
(setf s1
(if (minusp sign)
(- s1 (/ tmp d))
(+ s1 (/ tmp d))))
(incf d 2d0)
(setf tmp (* tmp r))
(setf sign (- sign)))
s1)))
(let* ((x1 (/ #q1 5))
(s1 (atan-series x1))
(x2 (/ #q1 239))
(s2 (atan-series x2))
(p (* (- (* s1 4)
s2)
4)))
(format t "~2&pi via Machin's atan formula~%")
(print-result p +pi+)
p)))
(defun test3 ()
(declare (optimize (speed 3)))
;; Salamin-Brent Quadratic formula for pi
(let* ((a #q1)
(b (sqrt #q.5))
(s #q.5)
(m 1d0)
(p (/ (* (* a a)
2d0)
s)))
(declare (double-float m))
(dotimes (k 9)
(setf m (* 2 m))
(let* ((a-new (* (+ a b) .5d0))
(b-new (sqrt (* a b)))
(s-new (- s
(* (- (* a-new a-new)
(* b-new b-new))
m))))
(setf a a-new)
(setf b b-new)
(setf s s-new)
(setf p (/ (* (* a a) 2d0)
s))))
(format t "~2&Salamin-Brent Quadratic formula for pi~%")
(print-result p +pi+)
p))
(defun test4 ()
(declare (optimize (speed 3)))
;; Borwein Quartic formula for pi
(let* ((a (- 6
(* (sqrt #q2)
4)))
(y (- (sqrt #q2)
1))
(m 2d0)
(p (/ a)))
(declare (double-float m))
(dotimes (k 9)
(setf m (* 4 m))
(let ((r (expt (- 1 (expt y 4))
1/4)))
(setf y (/ (- 1d0 r)
(+ 1d0 r)))
(setf a (- (* a
(expt (+ y 1d0) 4))
(* (* y
(+ (+ y (expt y 2))
1d0))
m)))
(setf p (/ a))))
(format t "~2&Borwein's Quartic formula for pi~%")
(print-result p +pi+)
p))
(defun test5 ()
;; Taylor series for e
(let ((s #q2)
(tmp #q1)
(n 1d0)
(delta 0d0)
(i 0))
(loop while (> tmp 1d-100) do
(incf i)
(incf n)
(setf tmp (/ tmp n))
(setf s (+ s tmp)))
(format t "~2&e via Taylor series~%")
(print-result s +e+)
s))
(defun test6 ()
;; Taylor series for log 2
;;
;; -log(1-x) = x + x^2/2 + x^3/3 + x^4/4 + ...
;;
;; with x = 1/2 to get log(1/2) = -log(2)
(let ((s #q.5)
(tt #q.5)
(n 1d0)
(i 0))
(loop while (> tt 1d-100) do
(incf i)
(incf n)
(setf tt (* tt .5d0))
(setf s (+ s
(/ tt n))))
(format t "~2&log(2) via Taylor series~%")
(print-result s +log2+)
s))
(defun test-atan ()
(let* ((arg (/ (sqrt #q3)))
(y (/ (atan arg) +pi+))
(true (/ #q6)))
(format t "~2&atan for special args~%")
(format t "atan(1/sqrt(3))/pi = 1/6~%")
(print-result y true))
;; atan(sqrt(3)) = %pi/3
(let* ((arg (sqrt #q3))
(y (/ (atan arg) +pi+))
(true (/ #q3)))
(format t "atan(sqrt(3))/pi = 1/3~%")
(print-result y true))
;; atan(1) = %pi/4
(let* ((arg #q1)
(y (/ (atan arg) +pi+))
(true (/ #q4)))
(format t "atan(1)/pi = 1/4~%")
(print-result y true))
(let* ((arg #q1q100)
(y (/ (atan arg) +pi+))
(true #q.5))
(format t "atan(1q100)/pi = 1/2~%")
(print-result y true))
(let* ((arg #q-1q100)
(y (/ (atan arg) +pi+))
(true #q-.5))
(format t "atan(-1q100)/pi = -1/2~%")
(print-result y true)))
(defun test-sin ()
(format t "~2&sin for special args~%")
(let* ((arg (/ +pi+ 6))
(y (sin arg))
(true #q.5))
(format t "sin(pi/6) = 1/2~%")
(print-result y true))
(let* ((arg (/ +pi+ 4))
(y (sin arg))
(true (sqrt #q.5)))
(format t "sin(pi/4) = 1/sqrt(2)~%")
(print-result y true))
(let* ((arg (/ +pi+ 3))
(y (sin arg))
(true (/ (sqrt #q3) 2)))
(format t "sin(pi/3) = sqrt(3)/2~%")
(print-result y true)))
(defun test-tan ()
(format t "~2&tan for special args~%")
(let* ((arg (/ +pi+ 6))
(y (tan arg))
(true (/ (sqrt #q3))))
(format t"tan(pi/6) = 1/sqrt(3)~%")
(print-result y true))
(let* ((arg (/ +pi+ 4))
(y (tan arg))
(true #q1))
(format t "tan(pi/4) = 1~%")
(print-result y true))
(let* ((arg (/ +pi+ 3))
(y (tan arg))
(true (sqrt #q3)))
(format t "tan(pi/3) = sqrt(3)~%")
(print-result y true)))
(defun test-asin ()
(format t "~2&asin for special args~%")
(let* ((arg #q.5)
(y (asin arg))
(true (/ +pi+ 6)))
(format t "asin(1/2) = pi/6~%")
(print-result y true))
(let* ((arg (sqrt #q.5))
(y (asin arg))
(true (/ +pi+ 4)))
(format t "asin(1/sqrt(2) = pi/4~%")
(print-result y true))
(let* ((arg (/ (sqrt #q3) 2))
(y (asin arg))
(true (/ +pi+ 3)))
(format t "asin(sqrt(3)/2) = pi/3~%")
(print-result y true)))
(defun test-log ()
(format t "~2&Log for special args~%")
(let* ((arg #q2)
(y (log arg))
(true +log2+))
(format t "log(2)~%")
(print-result y true))
[36 lines skipped]
--- /project/oct/cvsroot/oct/timing.lisp 2007/08/25 18:49:27 NONE
+++ /project/oct/cvsroot/oct/timing.lisp 2007/08/25 18:49:27 1.1
[212 lines skipped]
More information about the oct-cvs
mailing list