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;
}
};