Changeset 900


Ignore:
Timestamp:
Sep 4, 2011, 2:28:56 PM (8 years ago)
Author:
sam
Message:

core: improve tan() accuracy by tweaking higher order Taylor coefficients.

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/trig.cpp

    r899 r900  
    4444static const double SC[] =
    4545{
    46     -1.6449340668482264364724e-0, // pi^2/3!
    47     +8.1174242528335364363700e-1, // pi^4/5!
    48     -1.9075182412208421369647e-1, // pi^6/7!
    49     +2.6147847817654800504653e-2, // pi^8/9!
    50     -2.3460810354558236375089e-3, // pi^10/11!
    51     +1.4842879303107100368487e-4, // pi^12/13!
    52     -6.9758736616563804745344e-6, // pi^14/15!
    53     +2.5312174041370276513517e-7, // pi^16/17!
     46    -1.6449340668482264364724e-0, // π^2/3!
     47    +8.1174242528335364363700e-1, // π^4/5!
     48    -1.9075182412208421369647e-1, // π^6/7!
     49    +2.6147847817654800504653e-2, // π^8/9!
     50    -2.3460810354558236375089e-3, // π^10/11!
     51    +1.4842879303107100368487e-4, // π^12/13!
     52    -6.9758736616563804745344e-6, // π^14/15!
     53    +2.5312174041370276513517e-7, // π^16/17!
    5454};
    5555
     56/* Note: the last value should be -1.3878952462213772114468e-7 (ie.
     57 * π^18/18!) but we tweak it in order to get the better average precision
     58 * required for tan() computations when close to π/2+kπ values. */
    5659static const double CC[] =
    5760{
    58     -4.9348022005446793094172e-0, // pi^2/2!
    59     +4.0587121264167682181850e-0, // pi^4/4!
    60     -1.3352627688545894958753e-0, // pi^6/6!
    61     +2.3533063035889320454188e-1, // pi^8/8!
    62     -2.5806891390014060012598e-2, // pi^10/10!
    63     +1.9295743094039230479033e-3, // pi^12/12!
    64     -1.0463810492484570711802e-4, // pi^14/14!
    65     +4.3030695870329470072978e-6, // pi^16/16!
     61    -4.9348022005446793094172e-0, // π^2/2!
     62    +4.0587121264167682181850e-0, // π^4/4!
     63    -1.3352627688545894958753e-0, // π^6/6!
     64    +2.3533063035889320454188e-1, // π^8/8!
     65    -2.5806891390014060012598e-2, // π^10/10!
     66    +1.9295743094039230479033e-3, // π^12/12!
     67    -1.0463810492484570711802e-4, // π^14/14!
     68    +4.3030695870329470072978e-6, // π^16/16!
     69    -1.3777e-7,
    6670};
    6771
    6872/* These coefficients use Sloane’s http://oeis.org/A002430 and
    6973 * http://oeis.org/A036279 sequences for the Taylor series of tan().
    70  * Note: the last value should be 443861162*pi^18/1856156927625, ie.
    71  * 2.12485922978838540352881e5, but we tweak it in order to get
    72  * sub 1e-11 precision in a larger range. */
     74 * Note: the last value should be 2.12485922978838540352881e5 (ie.
     75 * 443861162*π^18/1856156927625), but we tweak it in order to get
     76 * sub 1e-11 average precision in a larger range. */
    7377static const double TC[] =
    7478{
    75     3.28986813369645287294483e0, // pi^2/3
    76     1.29878788045336582981920e1, // 2*pi^4/15
    77     5.18844961612069061254404e1, // 17*pi^6/315
    78     2.07509320280908496804928e2, // 62*pi^8/2835
    79     8.30024701695986756361561e2, // 1382*pi^10/155925
    80     3.32009324029001216460018e3, // 21844*pi^12/6081075
    81     1.32803704909665483598490e4, // 929569*pi^14/638512875
    82     5.31214808666037709352112e4, // 6404582*pi^16/10854718875
    83     2.373e5, // XXX: last value tweaked to improve precision
    84     //2.12485922978838540352881e5, // 443861162*pi^18/1856156927625
     79    3.28986813369645287294483e0, // π^2/3
     80    1.29878788045336582981920e1, // 2*π^4/15
     81    5.18844961612069061254404e1, // 17*π^6/315
     82    2.07509320280908496804928e2, // 62*π^8/2835
     83    8.30024701695986756361561e2, // 1382*π^10/155925
     84    3.32009324029001216460018e3, // 21844*π^12/6081075
     85    1.32803704909665483598490e4, // 929569*π^14/638512875
     86    5.31214808666037709352112e4, // 6404582*π^16/10854718875
     87    2.373e5,
    8588};
    8689
     
    344347        double x4 = x2 * x2;
    345348
    346         double subs1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
     349        /* Computing the Taylor series to the 11th order is enough to get
     350         * x * 1e-11 precision, but we push it to the 13th order so that
     351         * tan() has a better precision. */
     352        double subs1 = ((SC[5] * x4 + SC[3]) * x4 + SC[1]) * x4 + ONE;
    347353        double subs2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
    348354        double taylors = subs2 * x2 + subs1;
     
    407413    double x4 = x2 * x2;
    408414#if defined LOL_FEATURE_VERY_CHEAP_BRANCHES
    409     double subs1 = (SC[3] * x4 + SC[1]) * x4 + ONE;
     415    double subs1 = ((CC[5] * x4 + SC[3]) * x4 + SC[1]) * x4 + ONE;
    410416    double subs2 = (SC[4] * x4 + SC[2]) * x4 + SC[0];
    411417    double subc1 = ((CC[5] * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
     
    414420    double subs1 = (((SC[7] * x4 + SC[5]) * x4 + SC[3]) * x4 + SC[1]) * x4 + ONE;
    415421    double subs2 = ((SC[6] * x4 + SC[4]) * x4 + SC[2]) * x4 + SC[0];
     422    /* Push Taylor series to the 19th order to enhance tan() accuracy. */
    416423    double subc1 = (((CC[7] * x4 + CC[5]) * x4 + CC[3]) * x4 + CC[1]) * x4 + ONE;
    417     double subc2 = ((CC[6] * x4 + CC[4]) * x4 + CC[2]) * x4 + CC[0];
     424    double subc2 = (((CC[8] * x4 + CC[6]) * x4 + CC[4]) * x4 + CC[2]) * x4 + CC[0];
    418425#endif
    419426    double taylors = subs2 * x2 + subs1;
  • trunk/test/trig.cpp

    r899 r900  
    9494        }
    9595
    96         for (int i = -10000; i < 10000; i++)
     96        for (int i = -100000; i < 100000; i++)
    9797        {
    98             double f = (double)i * (1.0 / 1000.0);
     98            double f = (double)i * (1.0 / 10000.0);
    9999            double a = __builtin_tan(f);
    100100            double b = lol_tan(f);
    101             CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
     101            if (fabs(a) > 1e4)
     102               CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * fabs(a) * 1e-11);
     103            else if (fabs(a) > 1.0)
     104               CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
     105            else
     106               CPPUNIT_ASSERT(fabs(a - b) <= fabs(f) * 1e-11);
    102107        }
    103108
     
    107112            double a = __builtin_tan(f);
    108113            double b = lol_tan(f);
    109             CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
     114            if (fabs(a) > 1e4)
     115               CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * fabs(a) * 1e-11);
     116            else if (fabs(a) > 1.0)
     117               CPPUNIT_ASSERT(fabs(a - b) <= fabs(a) * 1e-11);
     118            else
     119               CPPUNIT_ASSERT(fabs(a - b) <= fabs(f) * 1e-11);
    110120        }
    111121    }
Note: See TracChangeset for help on using the changeset viewer.