gistfile1.cs

        private static double InterpolantRootfinder(double wt1, double wt2, double x)
        {
            double t = 0.5;

            if (wt1 == 1 / 3.0d && wt2 == 1 / 3.0d) return x;

            double wt2S = 1 - wt2;

            int i = 10;

            // it is about 2.5x faster to do this inline that it is to call an abstract function with callbacks
            while (i-- > 0)
            {
                double t2 = 1 - t;
                double fg = 3 * t2 * t2 * t * wt1 + 3 * t2 * t * t * wt2S + t * t * t - x;
                if (Math.Abs(fg) < 2 * EPS)
                    break;

                // third order householder method
                double fpg = 3 * t2 * t2 * wt1 + 6 * t2 * t * (wt2S - wt1) + 3 * t * t * (1 - wt2S);
                double fppg = 6 * t2 * (wt2S - 2 * wt1) + 6 * t * (1 - 2 * wt2S + wt1);
                double fpppg = 18 * wt1 - 18 * wt2S + 6;

                t -= (6 * fg * fpg * fpg - 3 * fg * fg * fppg) / (6 * fpg * fpg * fpg - 6 * fg * fpg * fppg + fg * fg * fpppg);
            }

            return t;
        }