This next phase in my reading of the Go erf
and erfc
implementations requires a more detailed research plan. I have already set my goals on tracing the development of the POSIX mathematics libraries, written down integrals for the error function and its complement, and identified source code in multiple languages for closer inspection.
As a welcome bonus, all the source code deriving from the original 1993 Sun Microsystems implementation preserves the historical commentary. Both the "who" and the "how" have clear documentation.
In particular, the source code for math.Erf
contains a copyright notice identifying the original implementation in the FreeBSD operating system. Likely, the Sun authors copied or modified the code from its predecessor at Berkeley. This same notice, and, what I'm most interested in, the coefficients used for the approximation can be found in MacOS X, in MUSL, in Android Linux, and others.
The comment, which legally must be copied along with the implementation, describes four methods to compute values for erf
. Each method defines a range of input values, computing:
near zero, [0, 0.84375]
near one, [0.8475, 1.25]
near Euler's number e
, [1.25, 2.85714]
below 28, [2.85714, 28]
Upon reflection, these are reasonable constants around which to set the bounds for the approximations. Since the function is anti-symmetric, we can flip the sign for values below zero and continue the calculation accurately, flipping back upon return. Both 1
and e
are so common as to be critical in computing with high precision. And above 28
, the structure of the IEEE 754 Float risks overflow. (Anti)Symmetry also means data may be rescale to fit the computation range.
In reality, the Go implementation cuts off at 6
, returning +/- 1
, respectively. I'm not clear on whether this cutoff is an optimization to reduce compute cycles, as my early experiments with reproducing the approximation functions suggest that an extension of the bound to 28 isn't necessary in practive. I could definitely be wrong! We'll see as I continue digging.
So we see, based on this comment, that the erf
implementation has been developed from four approximations, not just one. The note describes each approximation in turn, but for someone unaccustomed to the mathematics, it'll take a bit of work to see what's happening.
Note: I'm not accustomed to the mathematics. I began learning approximation theory in 2022 to feed my curiosity on this topic. I'm not an expert, but at this point I am a journeyman, of sorts, able to uncover, reproduce, and describe the work.
The theory we need to understand will be:
Rational Approximation, maybe also using Padé tables
Fixed Point Exchange or Remez' Algorithm
Asymptotic Series Expansion
Taylor Series Expansion
Stay awake -- this will be some rewarding work, flexing multiple mathematics and coding skills. Until next time, please subscribe to this series and get new essays in your inbox.