Ticket #2109: 2109rc1.patch

source/maths/Fixed.cpp

@@ -1 @@
-/* Copyright (C) 2010 Wildfire Games.
+/* Copyright (C) 2016 Wildfire Games.
  * This file is part of 0 A.D.
  *
  * 0 A.D. is free software: you can redistribute it and/or modify
@@ -143 +143 @@
     return r.str();
 }
 
-// Based on http://www.dspguru.com/dsp/tricks/fixed-point-atan2-with-self-normalization
 CFixed_15_16 atan2_approx(CFixed_15_16 y, CFixed_15_16 x)
 {
-    CFixed_15_16 zero;
+    // Overflow for length of (x,y) > 19910.8.
+    // Error <= arctan(1/16)/2 ~ 1.788 degrees.
 
-    // Special case to avoid division-by-zero
-    if (x.IsZero() && y.IsZero())
-        return zero;
+    // Binary search using the following angles:
+    // 45, 26.565, 14.036, 7.125, 3.576 degrees (= arctan(1/2^n))
+    // Idea: [1, 2^(-n); -2^(-n), 1] are multiples of the rotation matrices for these angles.
+    i32 tan_angle_table_rad[5] = {51472, 30386, 16055, 8150, 4091};
 
-    CFixed_15_16 c1;
-    c1.SetInternalValue(51472); // pi/4 << 16
+    i32 angleint = 2045; // = 1/2*arctan(1/16) << 16
+    i32 xint = x.GetInternalValue(), yint = y.GetInternalValue(), tmpint;
 
-    CFixed_15_16 c2;
-    c2.SetInternalValue(154415); // 3*pi/4 << 16
+    // Handle the case (0,0).
+    if (xint == 0 && yint == 0)
+        return CFixed_15_16::FromInt(0);
 
-    CFixed_15_16 abs_y = y.Absolute();
-
-    CFixed_15_16 angle;
-    if (x >= zero)
+    // Transform the problem into the case x,y >= 0, i.e. 0 <= angle <= pi/2
+    if (yint >= 0) // (0, pi)
     {
-        CFixed_15_16 r = (x - abs_y) / (x + abs_y);
-        angle = c1 - c1.Multiply(r);
+        if (xint < 0) // (pi/2, pi)
+        {
+            // (x, y) -> (y, -x)
+            tmpint = xint;
+            xint = yint;
+            yint = -tmpint;
+            angleint = 102944+2045; // pi/2 << 16
+        }
     }
-    else
+    else // (-pi, 0)
     {
-        CFixed_15_16 r = (x + abs_y) / (abs_y - x);
-        angle = c2 - c1.Multiply(r);
+        if (xint < 0) // (-pi, -pi/2)
+        {
+            // (x, y) -> (-x, -y)
+            xint *= -1;
+            yint *= -1;
+            angleint = -205887+2045; // -pi << 16
+        }
+        else // (-pi/2, 0)
+        {
+            // (x, y) -> (-y, x)
+            tmpint = xint;
+            xint = -yint;
+            yint = tmpint;
+            angleint = -102944+2045; // -pi/2 << 16
+        }
     }
 
-    if (y < zero)
-        return -angle;
-    else
-        return angle;
-}
+    // The following will have massive rounding errors for small x, y. If needed, scale them up. E.g.:
+    // if (xint < 128 && yint < 128)
+    // {
+    //  xint << 7;
+    //  yint << 7;
+    // }
+    for (u8 loop_var = 0; loop_var < 5; ++loop_var)
+    {
+        tmpint = xint >> loop_var;
+        // If remaining angle > arctan(1/2^n) "rotate" counterclockwise.
+        if (yint >= tmpint)
+        {
+            // (x,y) -> (x,y) + 1/2^n * (-y,x)
+            xint += yint >> loop_var;
+            yint -= tmpint;
+            angleint += tan_angle_table_rad[loop_var];
+        }
+    }
+    // Remaining angle is in [0,arctan(1/16)], arctan(1/16)/2 is the best possible approximation.
 
-template<>
-CFixed_15_16 CFixed_15_16::Pi()
-{
-    return CFixed_15_16(205887); // = pi << 16
+    // If needed, get the value into the desired interval of length 2*pi.
+
+    CFixed_15_16 angle;
+    angle.SetInternalValue(angleint);
+    return angle;
 }
 
 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/
+    // Based on roughly sup-norm optimal approximations of sin+cos, sin-cos on [0,pi/2]
+    // 2^16*(sin(pi/4*x)-cos(pi/4*x)) ~ 72730*(x-1) - 7210*(x-1)^3
+    // 2^16*(sin(pi/4*x)+cos(pi/4*x)) ~ 92682 - 28562*(x-1)^2 + 1416 * (x-1)^4
+    // The idea is to use these symmetric terms in order to get both sin and cos
+    // Errors: e ~ 1.3*10^-4 in both cases.
 
-    // TODO: this could be made a bit more precise by being careful about scaling
+    // The error could be decreased to around 1.1*10^-4 with further tweaking.
+    // TODO : Check if this breaks parts of the program that rely on specifics.
 
     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 and sin-cos
+    // Actually compute half of the above.
+    fixed c1,c2,c3,c4,c5;
+    c1.SetInternalValue(36365); // ~ 0.5549
+    c2.SetInternalValue(3605); // ~ 0.055
+    c3.SetInternalValue(46341); // ~ 0.7071
+    c4.SetInternalValue(14281); // ~ 0.2179
+    c5.SetInternalValue(708); // ~ 0.0108
+
+    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;
 }

source/maths/Fixed.h

@@ -126 +126 @@
 
     static CFixed Zero() { return CFixed(0); }
     static CFixed Epsilon() { return CFixed(1); }
-    static CFixed Pi();
+    static CFixed Pi() { return CFixed(205887); } // = pi << 16
+
 
     T GetInternalValue() const { return value; }
     void SetInternalValue(T n) { value = n; }

source/maths/tests/test_Fixed.h

@@ -281 +281 @@
         fixed s, c;
 
         sincos_approx(fixed::FromInt(0), s, c);
-        TS_ASSERT_EQUALS(s, fixed::FromInt(0));
-        TS_ASSERT_EQUALS(c, fixed::FromInt(1));
+        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.00015);
+        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.00015);
 
         sincos_approx(fixed::Pi(), s, c);
-        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.00005);
-        TS_ASSERT_EQUALS(c, fixed::FromInt(-1));
+        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.00015);
+        TS_ASSERT_DELTA(c.ToDouble(),-1.0, 0.00015);
 
         sincos_approx(fixed::FromDouble(M_PI*2.0), s, c);
-        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.0001);
-        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.0001);
+        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.00015);
+        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.00015);
 
         sincos_approx(fixed::FromDouble(M_PI*100.0), s, c);
-        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.004);
-        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.004);
+        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.0004);
+        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.0004);
 
         sincos_approx(fixed::FromDouble(M_PI*-100.0), s, c);
-        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.004);
-        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.004);
+        TS_ASSERT_DELTA(s.ToDouble(), 0.0, 0.0004);
+        TS_ASSERT_DELTA(c.ToDouble(), 1.0, 0.0004);
 
 /*
         for (double a = 0.0; a < 6.28; a += 0.1)
@@ -317 +317 @@
         }
 //      printf("### Max error %f = %f/2^16\n", err, err*65536.0);
 
-        TS_ASSERT_LESS_THAN(err, 0.00046);
+        TS_ASSERT_LESS_THAN(err, 0.00015);
     }
 };