From 00bd409b8d9de4f5f6223bf4d017e6f1c0826e48 Mon Sep 17 00:00:00 2001
From: Raymond Toy <toy.raymond@gmail.com>
Date: Thu, 12 Dec 2013 18:48:41 -0800
Subject: [PATCH] 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.
---
 src/code/irrat-dd.lisp | 174 +++++++----------------------------------
 1 file changed, 30 insertions(+), 144 deletions(-)

diff --git a/src/code/irrat-dd.lisp b/src/code/irrat-dd.lisp
index 4c57165ae..170bc06a9 100644
--- a/src/code/irrat-dd.lisp
+++ b/src/code/irrat-dd.lisp
@@ -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
-- 
GitLab