Code Poetry
and Text Adventures

by catid posted (>30 days ago) 11:10pm Sat. Jul 2nd 2016 PDT
Just a little toy project to share:

/*

Generate a tight polynomial fit from example data
and calculate error level in the exact solution

(1) Use scaled input data.

Based on the included benchmark:
   For unscaled data, average error level = 2.19262e-06
   For (0..1) scaled data, average error level = 2.54577e-07
The result is that scaled data interpolates 10x better on average.

(2) Using doubles hurts accuracy.

I also determined that using doubles instead of floats does not help at all
and can actually be slightly worse for accuracy due to conversion.

Technical contact: Christopher.Taylor@oculus.com
*/
struct LowErrorPolynomialFit
{

// p(t) = alpha + beta * t + gamma * t * t
float Alpha, Beta, Gamma;

float Evaluate(float t)
{
   // ~10% less error with Horner's rule as compared to calculating t * t
   return Alpha + t * (Beta + t * Gamma);
}

// Returns error level at perfect fit
// If the points are too close together the error will be high
float Init(float t0, float y0, float t1, float y1, float t2, float y2)
{
   const float sub_y1y0 = y1 - y0, sub_t2t0 = t2 - t0, sub_t1t0 = t1 - t0;

   // Terms have been shifted around to reduce numbers of divisions
   Gamma = sub_t1t0 * (y2 - y0) - sub_y1y0 * sub_t2t0;
   Gamma /= sub_t1t0 * sub_t2t0 * (t2 - t1);

   Beta = sub_y1y0 - sub_t1t0 * (t1 + t0) * Gamma;
   Beta /= sub_t1t0;

   // Apply Horner's rule to cut error in half here rather than doing t0 * t0
   Alpha = y0 - t0 * (Beta + Gamma * t0);

   return std::abs(Evaluate(t0) - y0)
      + std::abs(Evaluate(t1) - y1)
      + std::abs(Evaluate(t2) - y2);
}

// Returns average error level over a million contrived trials
// The input points are scaled to (0, 1) to make the best use of mantissa bits
// The effect of scaling is a 10x increase in accuracy, so it is worthwhile
static float Benchmark()
{
   LowErrorPolynomialFit fitf;
   srand(0);
   float accumErrorFitF = 0.f;
   static const int Trials = 1000000; // 1 meeellion trials
   for (int i = 0; i < Trials; ++i)
   {
      float t0 = 600.f + (10.f + (rand() / (float)RAND_MAX) * 200.f);
      float t1 = t0 + (10.f + (rand() / (float)RAND_MAX) * 200.f);
      float t2 = t1 + (10.f + (rand() / (float)RAND_MAX) * 200.f);

      static const float Low = 600.f;
      static const float High = 600.f + 210.f * 3 + 1.f;

      accumErrorFitF += fitf.Init(
           (t0 - Low) / (High - Low), 0.f,
           (t1 - Low) / (High - Low), 0.5f,
           (t2 - Low) / (High - Low), 1.f);
   }
   return accumErrorFitF / Trials;
}
};