--- ray/src/common/tcos.c	1998/12/16 18:15:07	3.1
+++ ray/src/common/tcos.c	2013/02/08 16:39:01	3.9
@@ -1,9 +1,6 @@
-/* Copyright (c) 1998 Silicon Graphics, Inc. */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ SGI";
+static const char	RCSid[] = "$Id: tcos.c,v 3.9 2013/02/08 16:39:01 greg Exp $";
 #endif
-
 /*
  * Table-based cosine approximation.
  *
@@ -11,24 +8,23 @@ static char SCCSid[] = "$SunId$ SGI";
  * to avoid conversion and guarantee that tsin(x)^2 + tcos(x)^2 == 1.
  *
  * No interpolation in this version.
+ *
+ * External symbols declared in rtmath.h
  */
 
+#include "copyright.h"
+
 #include <math.h>
 
+#include "rtmath.h"
+
 #ifndef NCOSENTRY
-#define NCOSENTRY	256
+#define NCOSENTRY	1024
 #endif
 
-#ifdef M_PI
-#define PI		((double)M_PI)
-#else
-#define PI		3.14159265358979323846
-#endif
 
-
 double
-tcos(x)				/* approximate cosine */
-register double	x;
+tcos(double x)				/* approximate cosine */
 {
 	static double	costab[NCOSENTRY+1];
 	register int	i;
@@ -40,8 +36,8 @@ register double	x;
 	if (x < 0.)
 		x = -x;
 	i = (NCOSENTRY*2./PI) * x  +  0.5;
-	if (i >= 4*NCOSENTRY)
-		i %= 4*NCOSENTRY;
+	while (i >= 4*NCOSENTRY)
+		i -= 4*NCOSENTRY;
 	switch (i / NCOSENTRY) {
 	case 0:
 		return(costab[i]);
@@ -53,4 +49,24 @@ register double	x;
 		return(costab[(4*NCOSENTRY)-i]);
 	}
 	return(0.);		/* should never be reached */
+}
+
+
+/* Fast arctangent approximation due to Rajan et al. 2006 */
+double
+atan2a(double y, double x)
+{
+	double	ratio, aratio, val;
+
+	if (x == 0)
+		return (y > 0) ? PI/2. : 3./2.*PI;
+
+	aratio = (ratio = y/x) >= 0 ? ratio : -ratio;
+
+	if (aratio > 1.01)
+		return PI/2. - atan2a(x, y);
+
+	val = PI/4.*ratio - ratio*(aratio - 1.)*(0.2447 + 0.0663*aratio);
+
+	return val + PI*(x < 0);
 }