--- ray/src/common/tcos.c	2003/02/25 02:47:22	3.3
+++ ray/src/common/tcos.c	2013/07/12 05:16:02	3.10
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char	RCSid[] = "$Id: tcos.c,v 3.3 2003/02/25 02:47:22 greg Exp $";
+static const char	RCSid[] = "$Id: tcos.c,v 3.10 2013/07/12 05:16:02 greg Exp $";
 #endif
 /*
  * Table-based cosine approximation.
@@ -9,27 +9,23 @@ static const char	RCSid[] = "$Id: tcos.c,v 3.3 2003/02
  *
  * No interpolation in this version.
  *
- * External symbols declared in standard.h
+ * External symbols declared in rtmath.h
  */
 
 #include "copyright.h"
 
 #include <math.h>
 
+#include "rtmath.h"
+
+#ifndef __FAST_MATH__
+
 #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;
@@ -41,8 +37,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]);
@@ -54,4 +50,25 @@ register double	x;
 		return(costab[(4*NCOSENTRY)-i]);
 	}
 	return(0.);		/* should never be reached */
+}
+
+#endif
+
+/* 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);
 }