Commit 00bd409b authored by Raymond Toy's avatar Raymond Toy

Simplify dd-%%sin, dd-%%cos, and dd-%%tan.

These routines did argument reduction, but since we use
__kernel_rem_pio2 to do accurate argument reduction, the argument
reduction in these routines is a waste of time.  This greatly
simplifies the routines to just the polynomial (or rational)
approximations.
parent e2c9ecef
......@@ -995,6 +995,14 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
-1.666666666666666666666666666666666647199w-1
)))
;; Compute sin(x) for |x| < pi/4 (approx).
(defun dd-%%sin (x)
(declare (type (double-double-float -1w0 1w0) x)
(optimize (speed 2) (space 0)
(inhibit-warnings 3)))
(let ((x2 (* x x)))
(+ x (* x (* x2 (poly-eval x2 sincof))))))
;; cos(x) = 1 - .5 x^2 + x^2 (x^2 P(x^2))
;; Theoretical peak relative error = 2.1e-37,
;; relative peak error spread = 1.4e-8
......@@ -1016,101 +1024,17 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
4.166666666666666666666666666666459301466w-2
)))
(defconstant dp1
(scale-float (float #b1100100100001111110110101010001000100001011010001100001000110100110001001100011001100010100010111000000011 1w0) -106))
(defconstant dp2
(scale-float (float #b0111000001110011010001001010010000001001001110000010001000101001100111110011000111010000000010000010111011 1w0) (* 2 -106)))
(defconstant dp3
(scale-float (float #b1110101001100011101100010011100110110010001001010001010010100000100001111001100011100011010000000100110111 1w0) (* 3 -106)))
(defconstant dp4
(scale-float (float #b0111101111100101010001100110110011110011010011101001000011000110110011000000101011000010100110110111110010 1w0) (* 4 -106)))
(defun dd-%%sin (x)
(declare (type double-double-float x)
(optimize (speed 2) (space 0)
(inhibit-warnings 3)))
(when (minusp x)
(return-from dd-%%sin (- (the double-double-float (dd-%%sin (- x))))))
;; y = integer part of x/(pi/4).
(let* ((y (float (floor (/ x dd-pi/4)) 1w0))
(z (scale-float y -4)))
(declare (type double-double-float y z))
(setf z (float (floor z) 1w0)) ; integer part of y/8
(setf z (- y (scale-float z 4))) ; y - 16*(y/16)
(let ((j (truncate z))
(sign 1))
(declare (type (integer -1 1) sign))
(unless (zerop (logand j 1))
(incf j)
(incf y))
(setf j (logand j 7))
(when (> j 3)
(setf sign (- sign))
(decf j 4))
;; Extended precision modular arithmetic
(setf z (- (- (- x (* y dp1))
(* y dp2))
(* y dp3)))
(let ((zz (* z z)))
(if (or (= j 1)
(= j 2))
(setf y (+ (- 1 (scale-float zz -1))
(* zz zz (poly-eval zz coscof))))
(setf y (+ z (* z (* zz (poly-eval zz sincof))))))
(if (< sign 0)
(- y)
y)))))
;; Compue cos(x) for |x| < pi/4 (approx)
(defun dd-%%cos (x)
(declare (type double-double-float x)
(declare (type (double-double-float -1w0 1w0) x)
(optimize (speed 2) (space 0)
(inhibit-warnings 3)))
(when (minusp x)
(return-from dd-%%cos (dd-%%cos (- x))))
;; y = integer part of x/(pi/4).
(let* ((y (float (floor (/ x dd-pi/4)) 1w0))
(z (scale-float y -4)))
(declare (type double-double-float y z))
(setf z (float (floor z) 1w0)) ; integer part of y/8
(setf z (- y (scale-float z 4))) ; y - 16*(y/16)
(let ((i (truncate z))
(j 0)
(sign 1))
(declare (type (integer 0 7) j)
(type (integer -1 1) sign))
(unless (zerop (logand i 1))
(incf i)
(incf y))
(setf j (logand i 7))
(when (> j 3)
(setf sign (- sign))
(decf j 4))
(when (> j 1)
(setf sign (- sign)))
;; Extended precision modular arithmetic. This is basically
;; computing x - y*(pi/4) accurately so that |z| < pi/4.
(setf z (- (- (- x (* y dp1))
(* y dp2))
(* y dp3)))
(let ((zz (* z z)))
(if (or (= j 1)
(= j 2))
(setf y (+ z (* z (* zz (poly-eval zz sincof)))))
(setf y (+ (- 1 (scale-float zz -1))
(* zz (poly-eval zz coscof) zz))))
(if (< sign 0)
(- y)
y)))))
(let ((x2 (* x x)))
(+ (- 1 (scale-float x2 -1))
(* x2 (poly-eval x2 coscof) x2))))
;; Compute tan(x) or cot(x) for |x| < pi/4 (approx). If cotflag is
;; non-nil, cot(x) is returned. Otherwise, return tan(x).
(let ((P (make-array 6 :element-type 'double-double-float
:initial-contents
'(
......@@ -1132,50 +1056,18 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
-4.152206921457208101480801635640958361612w10
8.650244186622719093893836740197250197602w10
))))
(defun dd-tancot (xx cotflag)
(declare (type double-double-float xx)
(optimize (speed 2) (space 0)))
(let ((x 0w0)
(sign 1))
(declare (type double-double-float x)
(type (integer -1 1) sign))
(cond ((minusp xx)
(setf x (- xx))
(setf sign -1))
(t
(setf x xx)))
(let* ((y (float (floor (/ x dd-pi/4)) 1w0))
(z (scale-float y -4))
(j 0))
(declare (type double-double-float y z)
(type fixnum j))
(setf z (float (floor z) 1w0))
(setf z (- y (scale-float z 4)))
(setf j (truncate z))
(unless (zerop (logand j 1))
(incf j)
(incf y))
(setf z (- (- (- x (* y dp1))
(* y dp2))
(* y dp3)))
(let ((zz (* z z)))
(if (> zz 1w-40)
(setf y (+ z
(* z (* zz (/ (poly-eval zz p)
(poly-eval-1 zz q))))))
(setf y z))
(if (not (zerop (logand j 2)))
(if cotflag
(setf y (- y))
(setf y (/ -1 y)))
(if cotflag
(setf y (/ y))))
(if (< sign 0)
(- y)
y))))))
(defun dd-tancot (x cotflag)
(declare (type (double-double-float -1w0 1w0) x)
(optimize (speed 2) (space 0) (inhibit-warnings 3)))
(let* ((xx (* x x))
(y (if (> xx 1w-40)
(+ x
(* x (* xx (/ (poly-eval xx p)
(poly-eval-1 xx q)))))
x)))
(if cotflag
(/ y)
y))))
(defun dd-%%tan (x)
(declare (type double-double-float x))
......@@ -1254,9 +1146,7 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
dd-%sin))
(defun dd-%sin (x)
(declare (double-double-float x))
(cond ((minusp (float-sign x))
(- (dd-%sin (- x))))
((< (abs x) (/ pi 4))
(cond ((< (abs x) (/ pi 4))
(dd-%%sin x))
(t
;; Argument reduction needed
......@@ -1272,9 +1162,7 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
dd-%cos))
(defun dd-%cos (x)
(declare (double-double-float x))
(cond ((minusp x)
(dd-%cos (- x)))
((< (abs x) (/ pi 4))
(cond ((< (abs x) (/ pi 4))
(dd-%%cos x))
(t
;; Argument reduction needed
......@@ -1290,9 +1178,7 @@ pi/4 11001001000011111101101010100010001000010110100011000 010001101001100010
dd-%tan))
(defun dd-%tan (x)
(declare (double-double-float x))
(cond ((minusp (float-sign x))
(- (dd-%tan (- x))))
((< (abs x) (/ pi 4))
(cond ((< (abs x) (/ pi 4))
(dd-%%tan x))
(t
;; Argument reduction needed
......
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