From c63d920e9f9f65b36dfd2e2a96657cf56cb2ae6a Mon Sep 17 00:00:00 2001
From: rtoy <rtoy>
Date: Wed, 19 Jul 2006 02:54:31 +0000
Subject: [PATCH] Add the trig argument reduction routines from Sun's fdlibm so
 we can accurately reduce the arg and therefore compute the value of trig
 functions accurately.

lisp/Config.linux_gencgc:
o Compile e_rem_pio2.c and k_rem_pio2.c

code/irrat.lisp:
o Disable %sin, %cos, %tan functions.
o Implement %sin, %cos, and %tan to call the fdlibm routine
  __ieee754_rem_pio2 to do argument reduction before calling the sin,
  cos, tan vops.

compiler/x86/float.lisp:
o Disable the vops for %sin, %cos, and %tan, so the Lisp code in
  irrat.lisp is used.
---
 code/irrat.lisp          | 48 +++++++++++++++++++++++++++++++++++++++-
 compiler/x86/float.lisp  |  5 ++++-
 lisp/Config.linux_gencgc |  2 +-
 3 files changed, 52 insertions(+), 3 deletions(-)

diff --git a/code/irrat.lisp b/code/irrat.lisp
index 25087903e..66d0b08e1 100644
--- a/code/irrat.lisp
+++ b/code/irrat.lisp
@@ -5,7 +5,7 @@
 ;;; Carnegie Mellon University, and has been placed in the public domain.
 ;;;
 (ext:file-comment
-  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.46 2006/06/30 18:41:22 rtoy Exp $")
+  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/irrat.lisp,v 1.47 2006/07/19 02:54:30 rtoy Exp $")
 ;;;
 ;;; **********************************************************************
 ;;;
@@ -84,6 +84,7 @@
 
 #+x86 ;; These are needed for use by byte-compiled files.
 (progn
+  #+nil
   (defun %sin (x)
     (declare (double-float x)
 	     (values double-float))
@@ -92,6 +93,7 @@
     (declare (double-float x)
 	     (values double-float))
     (%sin-quick x))
+  #+nil
   (defun %cos (x)
     (declare (double-float x)
 	     (values double-float))
@@ -100,6 +102,7 @@
     (declare (double-float x)
 	     (values double-float))
     (%cos-quick x))
+  #+nil
   (defun %tan (x)
     (declare (double-float x)
 	     (values double-float))
@@ -168,6 +171,49 @@
     (%sqrt x))
   )
 
+#+x86
+(progn
+(alien:def-alien-routine ("__ieee754_rem_pio2" %ieee754-rem-pi/2) c-call:int
+  (x double-float)
+  (y (* double-float)))
+
+(let ((y (make-array 2 :element-type 'double-float)))
+  (defun %sin (x)
+    (declare (double-float x))
+    (if (< (abs x) (/ pi 4))
+	(%sin-quick x)
+	;; Argument reduction needed
+	(let* ((n (%ieee754-rem-pi/2 x (vector-sap y)))
+	       (reduced (+ (aref y 0) (aref y 1))))
+	  (case (logand n 3)
+	    (0 (%sin-quick reduced))
+	    (1 (%cos-quick reduced))
+	    (2 (- (%sin-quick reduced)))
+	    (3 (- (%cos-quick reduced)))))))
+  (defun %cos (x)
+    (declare (double-float x))
+    (if (< (abs x) (/ pi 4))
+	(%cos-quick x)
+	;; Argument reduction needed
+	(let* ((n (%ieee754-rem-pi/2 x (vector-sap y)))
+	       (reduced (+ (aref y 0) (aref y 1))))
+	  (case (logand n 3)
+	    (0 (%cos-quick reduced))
+	    (1 (- (%sin-quick reduced)))
+	    (2 (- (%cos-quick reduced)))
+	    (3 (%sin-quick reduced))))))
+  (defun %tan (x)
+    (declare (double-float x))
+    (if (< (abs x) (/ pi 4))
+	(%tan-quick x)
+	;; Argument reduction needed
+	(let* ((n (%ieee754-rem-pi/2 x (vector-sap y)))
+	       (reduced (+ (aref y 0) (aref y 1))))
+	  (if (evenp n)
+	      (%tan-quick reduced)
+	      (- (/ (%tan-quick reduced))))))))
+)
+
 
 ;;;; Power functions.
 
diff --git a/compiler/x86/float.lisp b/compiler/x86/float.lisp
index 3f129c896..b8f710416 100644
--- a/compiler/x86/float.lisp
+++ b/compiler/x86/float.lisp
@@ -7,7 +7,7 @@
 ;;; Scott Fahlman or slisp-group@cs.cmu.edu.
 ;;;
 (ext:file-comment
-  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/x86/float.lisp,v 1.43 2006/07/01 13:52:48 rtoy Exp $")
+  "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/compiler/x86/float.lisp,v 1.44 2006/07/19 02:54:31 rtoy Exp $")
 ;;;
 ;;; **********************************************************************
 ;;;
@@ -2511,6 +2511,7 @@
 ;;; These versions of fsin, fcos, and ftan simply load a 0.0 result if
 ;;; the argument is out of range 2^63 and would thus be hopelessly
 ;;; inaccurate.
+#+nil
 (macrolet ((frob (func trans op)
 	     `(define-vop (,func)
 		(:translate ,trans)
@@ -2546,6 +2547,7 @@
 	  (frob fsin  %sin fsin)
 	  (frob fcos  %cos fcos))
 	     
+#+nil
 (define-vop (ftan)
   (:translate %tan)
   (:args (x :scs (double-reg) :target fr0))
@@ -3656,6 +3658,7 @@
 ;;; These versions of fsin, fcos, and ftan simply load a 0.0 result if
 ;;; the argument is out of range 2^63 and would thus be hopelessly
 ;;; inaccurate.
+#+nil
 (macrolet ((frob (func trans op)
 	     `(define-vop (,func)
 		(:translate ,trans)
diff --git a/lisp/Config.linux_gencgc b/lisp/Config.linux_gencgc
index 22c67f944..c64900cb1 100644
--- a/lisp/Config.linux_gencgc
+++ b/lisp/Config.linux_gencgc
@@ -27,7 +27,7 @@ NM = $(PATH1)/linux-nm
 UNDEFSYMPATTERN = -Xlinker -u -Xlinker &
 ASSEM_SRC = x86-assem.S linux-stubs.S
 ARCH_SRC = x86-arch.c
-OS_SRC = Linux-os.c os-common.c elf.c
+OS_SRC = Linux-os.c os-common.c elf.c e_rem_pio2.c k_rem_pio2.c
 OS_LINK_FLAGS = -rdynamic -Xlinker --export-dynamic -Xlinker -Map -Xlinker foo
 OS_LIBS = -ldl
 #GC_SRC = gencgc.c
-- 
GitLab