Ticket #2109: sincos.patch

source/maths/Fixed.cpp

@@ -184 +184 @@
     return CFixed_15_16(205887); // = pi << 16
 }
 
+// Based on sup-norm optimal approximations of sin+cos, sin-cos on [0,pi/2] (Remez algorithm)
+// 2^15*(sin(pi/4*x)-cos(pi/4*x)) ~ 36361*(x-1) - 3600*(x-1)^3
+// 2^15*(sin(pi/4*x)+cos(pi/4*x)) ~ 46340 - 14281*(x-1)^2 + 709 * (x-1)^4
+// The idea is to use these symmetric terms in order to get both sin and cos
+// Errors: e ~ 1.1*10^-4 in both cases.
+// TODO : Check if this breaks parts of the program that rely on special cases.
 void sincos_approx(CFixed_15_16 a, CFixed_15_16& sin_out, CFixed_15_16& cos_out)
 {
-    // Based on http://www.coranac.com/2009/07/sines/
-
-    // TODO: this could be made a bit more precise by being careful about scaling
-
     typedef CFixed_15_16 fixed;
 
-    fixed c2_pi;
-    c2_pi.SetInternalValue(41721); // = 2/pi << 16
+    fixed c4_pi;
+    c4_pi.SetInternalValue(83443); // = 4/pi << 16
 
-    // Map radians onto the range [0, 4)
-    fixed z = a.Multiply(c2_pi) % fixed::FromInt(4);
+    // Map radians onto the range [0, 8)
+    fixed z = a.Multiply(c4_pi) % fixed::FromInt(8);
 
-    // Map z onto the range [-1, +1] for sin, and the same with z+1 to compute cos
-    fixed sz, cz;
-    if (z >= fixed::FromInt(3)) // [3, 4)
-    {
-        sz = z - fixed::FromInt(4);
-        cz = z - fixed::FromInt(3);
-    }
-    else if (z >= fixed::FromInt(2)) // [2, 3)
-    {
-        sz = fixed::FromInt(2) - z;
-        cz = z - fixed::FromInt(3);
-    }
-    else if (z >= fixed::FromInt(1)) // [1, 2)
+    // Map z onto the range [0, 2], compute (sin+cos)/2 and (sin-cos)/2
+    fixed c1,c2,c3,c4,c5;
+    c1.SetInternalValue(36361);
+    c2.SetInternalValue(3600);
+    c3.SetInternalValue(46340);
+    c4.SetInternalValue(14281);
+    c5.SetInternalValue(709);
+
+    if (z >= fixed::FromInt(4)) // [4,8)
     {
-        sz = fixed::FromInt(2) - z;
-        cz = fixed::FromInt(1) - z;
+        if (z >= fixed::FromInt(6)) // [6,8)
+        {
+            // sin(3*pi/2+x) = -cos(x), cos(3*pi/2+x) = sin(x)
+            fixed zn = z - fixed::FromInt(7);
+            z  = zn.Multiply(zn);
+            fixed scdif = zn.Multiply(c1 - z.Multiply(c2));
+            fixed scsum = c3 - z.Multiply(c4 - z.Multiply(c5));
+            cos_out = scsum + scdif;
+            sin_out = scdif - scsum;
+        }
+        else // [4,6)
+        {
+            // sin(pi+x) = -sin(x), cos(pi+x) = -cos(x)
+            fixed zn = z - fixed::FromInt(5);
+            z  = zn.Multiply(zn);
+            fixed scdif= zn.Multiply(c1 - z.Multiply(c2));
+            fixed scsum = c3 - z.Multiply(c4 - z.Multiply(c5));
+            sin_out = -scsum - scdif;
+            cos_out = scdif - scsum;
+        }
     }
-    else // [0, 1)
+    else // [0,4)
     {
-        sz = z;
-        cz = fixed::FromInt(1) - z;
+        if (z >= fixed::FromInt(2)) // [2,4)
+        {
+            // sin (pi/2+x) = cos(x), cos(pi/2+x) = -sin(x)
+            fixed zn = z - fixed::FromInt(3);
+            z  = zn.Multiply(zn);
+            fixed scdif = zn.Multiply(c1 - z.Multiply(c2));
+            fixed scsum = c3 - z.Multiply(c4 - z.Multiply(c5));
+            cos_out = -scsum - scdif;
+            sin_out =  scsum - scdif;
+        }
+        else // [0,2)
+        {
+            fixed zn = z - fixed::FromInt(1);
+            z  = zn.Multiply(zn);
+            fixed scdif = zn.Multiply(c1 - z.Multiply(c2));
+            fixed scsum = c3 - z.Multiply(c4 - z.Multiply(c5));
+            sin_out = scsum + scdif;
+            cos_out = scsum - scdif;
+        }
     }
-
-    // Third-order (max absolute error ~0.02)
-
-//  sin_out = (sz / 2).Multiply(fixed::FromInt(3) - sz.Multiply(sz));
-//  cos_out = (cz / 2).Multiply(fixed::FromInt(3) - cz.Multiply(cz));
-
-    // Fifth-order (max absolute error ~0.0005)
-
-    fixed sz2 = sz.Multiply(sz);
-    sin_out = sz.Multiply(fixed::Pi() - sz2.Multiply(fixed::Pi()*2 - fixed::FromInt(5) - sz2.Multiply(fixed::Pi() - fixed::FromInt(3)))) / 2;
-
-    fixed cz2 = cz.Multiply(cz);
-    cos_out = cz.Multiply(fixed::Pi() - cz2.Multiply(fixed::Pi()*2 - fixed::FromInt(5) - cz2.Multiply(fixed::Pi() - fixed::FromInt(3)))) / 2;
-}
+}
+ No newline at end of file