Formula Translation: The Error Function

Not At All Elementary, Watson

Previously:

Differentiating f(x)^g(x)
The Differentiator
Graphing Derivatives
The Complexities Of The Simplifier

Project Differentiator is growing rapidly (!) into a customizable calculus teaching tool. Such resources are sometimes more valuable than any generally available mathematics code library, utility, or even online suite like Wolfram Alpha or Mathematica. Just ask Grant Sanderson, whose magnificent YouTube channel 3blue1brown - in addition to its inestimable educational payload - delivers and promotes beautiful mathematics as an art form, rolling along atop its very own powerful, handcrafted, dynamic, mathematical graphics engine. When you want to explain things in great detailed graphics, and especially in animations, nothing beats running your own source code.

Waiting for the dust to settle on the latest UI iteration, I thought I'd mention this short function here: a three-line translation, based on some old 1992 code found in a Fortran77 book, of the error function, y=erf(x), as defined in the integral equation alongside. That FORmula TRANslation provenance remains visible in my choice of variable names in the code below: well of course n is the integer loop variable.

Incidentally, you should never use the Fortran-recommended I as your loop variable. If a single pin breaks in your dot matrix printer, the resultant printout on your music ruled paper might be mistaken for the digit 1, crashing your lander on the moon.

The interesting thing about erf(x) is that it's non-elementary, making it actually one of the uninteresting ones, globally speaking.

What's Special About Error?

    double t = 1 / (1 + Math.Abs(x) / 2), u = 1, v = 0;
    for (var n = 0; n < 10; v += erf_a[n] * u, n++, u *= t) ;
    return (1 - t * Math.Exp(v - x * x)) * Math.Sign(x);

Yeah, absolutely nothing. Non-elementary functions might appear curious to engineers learning calculus, but the proportion of all possible continuous functions that are elementary is zero (assuming that's possible for the ratio of two uncountable infinities). An elementary function, as originally introduced by Joseph Liouville starting in 1833, is one built from certain basic functions (constant, algebraic, exponential and logarithmic, along with their inverses if they exist), by recursive composition and application of arithmetic operators, including root extraction.

On the other hand, despite its non-elementary nature, erf(x) is in great demand, mostly because it's literally the area under the normal distribution graph in statistics. And when you hear there's no way to integrate y=exp(x²), well, a constant multiple of erf(x) is the answer to that too. These reasons, together with its eminent differentiability (it's defined as an integral, after all), are why I wanted it in the Lego™ box.

Speaking of which, the list of elementary functions available in Sid, as the project is now dubbed, has more than doubled since you last looked. And if you ever needed an inverse trigonometric or hyperbolic function, but discovered it missing from the supplied System.Math set, you'll find the entire lot implemented here. Anyway, the erf code is the bottom third of this file:

namespace FormulaBuilder
{
    using System;

    /// <summary>
    /// Several math functions are either absent from System.Math (Asec, Sech, Asech,
    /// Step) or unsuitable (Sign, which has an integer output), so they are provided
    /// in this class. Might as well include the whole System.Math lot here too; then
    /// we never need perform two reflection operations to find any function by name.
    /// </summary>
    public static class Functions
    {
        #region Elementary Functions

        public static double Abs(double x) => Math.Abs(x);
        public static double Acos(double x) => Math.Acos(x);
        public static double Acosh(double x) => Math.Log(x + Math.Sqrt(Math.Pow(x, 2) - 1));
        public static double Acot(double x) => Math.Atan(1 / x);
        public static double Acoth(double x) => Math.Log((x + 1) / (x - 1)) / 2;
        public static double Acsc(double x) => Math.Asin(1 / x);
        public static double Acsch(double x) => Math.Log((1 + Math.Sqrt(1 + Math.Pow(x, 2))) / x);
        public static double Asec(double x) => Math.Acos(1 / x);
        public static double Asech(double x) => Math.Log((1 + Math.Sqrt(1 - Math.Pow(x, 2))) / x);
        public static double Asin(double x) => Math.Asin(x);
        public static double Asinh(double x) => Math.Log(x + Math.Sqrt(Math.Pow(x, 2) + 1));
        public static double Atan(double x) => Math.Atan(x);
        public static double Atanh(double x) => Math.Log((1 + x) / (1 - x)) / 2;
        public static double Ceiling(double x) => Math.Ceiling(x);
        public static double Cos(double x) => Math.Cos(x);
        public static double Cosh(double x) => Math.Cosh(x);
        public static double Cot(double x) => 1 / Math.Tan(x);
        public static double Coth(double x) => 1 / Math.Tanh(x);
        public static double Csc(double x) => 1 / Math.Sin(x);
        public static double Csch(double x) => 1 / Math.Sinh(x);
        public static double Exp(double x) => Math.Exp(x);
        public static double Floor(double x) => Math.Floor(x);
        public static double Ln(double x) => Math.Log(x);
        public static double Log10(double x) => Math.Log10(x);
        public static double Round(double x) => Math.Round(x);
        public static double Sec(double x) => 1 / Math.Cos(x);
        public static double Sech(double x) => 1 / Math.Cosh(x);
        public static double Sign(double x) => Math.Sign(x);
        public static double Sin(double x) => Math.Sin(x);
        public static double Sinh(double x) => Math.Sinh(x);
        public static double Sqrt(double x) => Math.Sqrt(x);
        public static double Step(double x) => x < 0 ? 0 : 1;
        public static double Tan(double x) => Math.Tan(x);
        public static double Tanh(double x) => Math.Tanh(x);

        #endregion

        #region Non-elementary Functions

        private static readonly double[] erf_a =
        {
            -1.26551223, +1.00002368, +0.37409196, +0.09678418, -0.18628806,
            +0.27886807, -1.13520398, +1.48851587, -0.82215223, +0.17087277
        };

        /// <summary>
        /// An approximation, with a worst case error of 1.2e-7, from the book
        /// "Numerical Recipes in Fortran 77: The Art of Scientific Computing"
        /// (ISBN 0-521-43064-X), 1992, page 214, Cambridge University Press.
        /// </summary>
        /// <param name="x">The input value to the function Erf(x).</param>
        /// <returns>An approximation to the value of erf(x).</returns>
        public static double Erf(double x)
        {
            double t = 1 / (1 + Math.Abs(x) / 2), u = 1, v = 0;
            for (var n = 0; n < 10; v += erf_a[n] * u, n++, u *= t) ;
            return (1 - t * Math.Exp(v - x * x)) * Math.Sign(x);
        }

        #endregion
    }
}

I'll confess to considering the elision of that redundant final u *= t in the for loop, before remembering the pipeline exception would certainly waste more picoseconds than a blast through the floating point silicon. Old optimisation rabbits die hard.

The algorithm first sets t=1/(1+|x|/2). Then in the for loop, it forms a ninth degree polynomial v(t), using ten supplied coefficients erf_a[0..9]. Finally, according to the original sign of x, your answer is either 1-t*e^(v(t)-x²), or its negative.

Next time: let's get that parser up and running...

Image credits: https://en.wikipedia.org/wiki/Error_function.