Commit b71b7c85 authored by Raymond Toy's avatar Raymond Toy

Increase accuracy of tan near multiples of pi/4.

code/irrat-dd.lisp
o Make REDUCE-ARG return an extra result.  Since the reduction returns
  3 double's, return a double-double result and the third double
  result for extra accuracy
o Update dd-%tan and dd-%%tan to take the extra arg.
o Add new constant dd-pi/4-lo.
o Increase accuracy of tan by using the relationship
  tan(pi/4-y)=(1-tan(y))/(1+tan(y)).

tests/trig.lisp:
o Update the allowed error threshold for two tests to reflect the
  increased accuracy.
parent 712df0bc
......@@ -56,6 +56,12 @@
0.7853981633974483096156608458198757210492923w0
_N"Pi/4")
;; dd-pi/4-lo is such that dd-pi/4-lo + dd-pi/4 is equal to pi/4 to
;; twice the precision of a double-double-float.
(defconstant dd-pi/4-lo
-7.486924524295848886603985669688687133352026408988280860093283232w-34
)
;; log2-c1 and log-c2 are log(2) arranged in such a way that log2-c1 +
;; log2-c2 is log(2) to an accuracy greater than double-double-float.
(defconstant log2-c1
......@@ -1072,9 +1078,24 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
(/ y)
y))))
(defun dd-%%tan (x)
(declare (type double-double-float x))
(dd-tancot x nil))
(defun dd-%%tan (x extra)
(declare (type double-double-float x)
(double-float extra))
(cond ((>= (abs x) 0.6744)
;; For 0.6744 <= |x| <= pi/4, we want to use the relationship
;;
;; tan(x) = tan(pi/4-y) = (1 - tan(y))/(1 + tan(y))
;; = 1 - 2*(tan(y) - tan(y)^2)/(1+tan(y))
(if (minusp x)
(- (dd-%%tan (- x) (- extra)))
(let* ((z (- dd-pi/4 x))
(w (- dd-pi/4-lo extra))
(tan (dd-tancot (+ z w) nil)))
(- 1
(/ (* 2 (- tan (* tan tan)))
(+ 1 tan))))))
(t
(dd-tancot x nil))))
(declaim (inline %kernel-rem-pi/2))
(alien:def-alien-routine ("__kernel_rem_pio2" %kernel-rem-pi/2) c-call:int
......@@ -1139,10 +1160,9 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
(length parts)
3
(vector-sap two-over-pi))))
(sum (+ (coerce (aref y 2) 'double-double-float)
(coerce (aref y 1) 'double-double-float)
(sum (+ (coerce (aref y 1) 'double-double-float)
(coerce (aref y 0) 'double-double-float))))
(values n sum))))
(values n sum (aref y 2)))))
(declaim (ftype (function (double-double-float) double-double-float)
......@@ -1181,15 +1201,15 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
dd-%tan))
(defun dd-%tan (x)
(declare (double-double-float x))
(cond ((< (abs x) (/ pi 4))
(dd-%%tan x))
(cond ((<= (abs x) (/ pi 4))
(dd-%%tan x 0d0))
(t
;; Argument reduction needed
(multiple-value-bind (n reduced)
(multiple-value-bind (n reduced extra)
(reduce-arg x)
(if (evenp n)
(dd-%%tan reduced)
(- (/ (dd-%%tan reduced))))))))
(dd-%%tan reduced extra)
(- (/ (dd-%%tan reduced extra))))))))
(defun dd-%sincos (x)
(declare (double-double-float x))
......
......@@ -376,12 +376,12 @@
(assert-eq t (rel-or-abs-error
(tan (* 7/4 kernel:dd-pi))
-1.000000000000000000000000000000001844257310064121018312678894979w0
6.467w-33))
3.422w-49))
;; Test for argument reduction with n odd
(assert-eq t (rel-or-abs-error
(tan (* 9/4 kernel:dd-pi))
1.000000000000000000000000000000025802415787810837455445433037983w0
5.773w-33))
0w0))
;; Test for argument reduction, big value
(assert-eq t (rel-or-abs-error
(tan (scale-float 1w0 120))
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment