Implementing the Exponential Function(pseudorandom.com) |
Implementing the Exponential Function(pseudorandom.com) |
The Sinclair Scientific calculator was notable for cramming transcendental functions into a chip designed as a 4-function calculator, so they took some severe shortcuts. Exponentials were based on computing 0.99^n through repeated multiplication. Multiplying by 0.99 is super-cheap in BCD since you shift by two digits and subtract. To compute 10^x, it computes .99^(-229.15*x). Slow and inaccurate, but easy to implement.
http://files.righto.com/calculator/sinclair_scientific_simul...
And then I noticed who I'm replying to...
Despite its flaws, this calculator was quite an achievement, doing everything it did with a ROM that holds only 320 instructions.
There are a number of facts/folklore about Chebychev polynomials that seem to go stated without proof in a lot of papers, so a few years ago I wrote up some notes [3,4] with the most common claims. Maybe someone will find them useful!
[1] (p35) http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradi...
[2] http://blog.mrtz.org/2013/09/07/the-zen-of-gradient-descent....
[3] https://benrbray.com/static/notes/chebychev-polynomials_dec1...
[4] https://benrbray.com/static/notes/minimax-approx_nov16.pdf
I will say that the article didn't really touch on techniques that minimize IEEE 754 subnormal[2] performance impact, which is a very interesting problem in itself. A lot of real-life implementations of something like e^x will have various clever ways to avoid subnormals.
I used the Remez algorithm, modified to minimize equivalent input error. This algorithm lets you find a polynomial with the smallest maximum error. You start with a set of X coordinates, and create a polynomial which oscillates up and down around the target function, so p(x0) = f(x0) + E, p(x1) = f(x1) - E, p(x2) = f(x2) + E, etc.
You then iterate by replacing the set of x values with the x values which have maximum error. This will converge to a polynomial with smallest maximum error.
From my experiments, you can use a fairly low order approximation for musical applications. Used to calculate frequencies of musical notes, 2nd order is within 3.0 cents, and 3rd order is within 0.13 cents.
> Any function with this property can be uniquely represented as a Taylor series expansion “centered” at a
It's fairly easy to construct tame looking functions where this isn't true.
Here is an open access link to just one of many papers on the subject: https://projecteuclid.org/euclid.em/1120145574
It includes a link to the actual software they used.
[EDIT] Obviously applying the technique above to dx/dt = x reproduces the Taylor Series in the article!
Some other resources you may find interesting
https://fpbench.org/index.html - I think I saw a talk at FPTalks that felt related to your topic "Creating correctly rounded math libraries for real number approximations"
https://fpbench.org/community.html
https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804 - Correctly Rounded libm
https://www.mpfr.org/ - GNU Multi precision floating point
http://sollya.gforge.inria.fr/ - A tool for floating point code development
https://github.com/sipa/minisketch/blob/master/src/false_pos...
I was interested in the exponential function too, but approached it from the standpoint of arbitrary-precision arithmetic. My explanation is an order of magnitude shorter, but of course my code is a lot slower than a typical floating-point algorithm. I offer a way to compute correctly rounded exponentials without resorting to a big math package like Mathematica/Maple/etc. https://www.nayuki.io/page/approximating-eulers-number-corre...
I knew this in general but this specific statistic was eye-opening to me. Systems like OpenGL normalize coordinates so they only use values in this [-1,1] range. That means they effectively lose half the expressive power of the 32 bits. Are there other machine representations of real numbers which would be better for that use case? i.e. restricted to that range and evenly distributed?
So there's very little loss in using single-precision floating point, and a lot of gains in smoothly handling larger numbers that arise from addition et al.
edit: I'm half wrong. The bit in mantissa isn't wasted, the only bit that's wasted in the sign bit in the exponent, so the efficiency loss is ~1.5%
I don’t recommend it.
That said I liked that there’s a specific exp-1 function (otherwise you drop most precision).
x87 has pretty much all of them, including two remainder implementations \o/
$$ \log\left(\frac{x e}{n}\right) \leq \frac{1}{n}\log(\eta) $$
$$ \frac{xe}{n} \leq \eta^{-n} $$
$$ \frac{xe}{n} \leq \eta^{1/n} $$
Certainly understandable flub, though. Working with logarithms, it's easy to mix up minus sign vs. 1/n conversions.
First, `x + x + x + x + x == 5x`. This is true for values up to 5, but is not true for 6. Proving this is a fun exercise, but is a little bit painful and has a lot of cases.
Second, SMT solvers can actually prove stuff like this relatively easy! In Z3py, one can prove this in 4 lines of code.
from z3 import x = FP('x', FPSort(8, 24)) set_default_rounding_mode(RNE()) solve(5*x != x + x+ x+ x + x)
I only wish Trefethen had written it with general C/C++ instead of his own custom MATLAB library, chebfun. One of these days I'd like to get around to implementing chebfun in C++.
* The article looks beautiful, but am I the only one that hates left-centered content? Wide monitors are the norm these days, and reading left centered content literally hurts my neck (need to resize manually to center the content, approximately.
* Why would someone make the effort to write such detailed analysis and conceal their identity? I just don't get what's the objective of not providing any sort of provenance... just a gripe of mine.
It is funny you mention the Padé approximation; I was debating whether or not to include that in this article. Originally I planned to stop after Minimax (using Remez), but if I include rational approximation schemes anyway (e.g. Carathéodory-Fejer), it probably makes sense to implement and analyze that as well.
That's not a deep exploration of the numerical analysis (especially compared with OP), but it does make one point: traditional "numerical analysis" techniques tend to assume the traditional four mathematical operators as canonical, but if your goal is to get something computed on actual computer hardware quickly and reasonably accurately, you have a wider palette to play with. In particular, inverse square root is a primitive about as efficient as reciprocal on modern hardware (or if using Carmack's wtf trick).
http://mpfr.loria.fr/mpfr-current/mpfr.html#index-mpfr_005fs...
It is way to deep in the grubby details for my pay grade, but who knows?
[1] https://github.com/dvx/lofi/commit/285bf80ff6d7f0784c14270de...
Instead, it's better to exploit the definition of IEEE754 as an exponent and a mantissa. You calculate the exponent directly (because e^x = 2^(x/ln2)), and use the Remez algorithm to find a polynomial that fits just the mantissa.
You have to explicitly choose a range [x1,x2] for the Remez algorithm. It minimizes the maximum error over that range.
I'm ordinarily happy to share the work I'm doing but this comment reminds me of the reasons I don't share much with HN.
If anyone knows some nice articles about this topic I would love to read them
For a domain like [0,1], the Remez algorithm will do a great job -- because the exponent is barely changing, just the mantissa. For a range like [-10, 10], you won't find any polynomial of reasonable order that has acceptable error.
But, there's a really neat trick that does let the Remez algorithm work over the entire domain of floats, and get excellent performance with just a fourth or fifth order polynomial, or even third order if you're pressed for time, and minimal extra computation.
That is, with just a little simplification, you multiply x by (1/ln(2)), floor it, and stuff the resulting integer into the exponent bits. Then, you take the remainder (what is left over after taking the floor), and make that your input to the polynomial, and stuff the result into the mantissa bits.
There's a little more to it; signs and NaNs need to be handled correctly, subnormal numbers can be treated as a special case with a different Remez polynomial, and a few other corner cases and exact values might be important (e^0, e^1, e^-1). But the result works wonders and it's basically just using Remez but transforming the input to a domain that behaves more like a polynomial, leveraging the magic of IEEE-754.
That reminds me, I'm using the domain [-1, 1] for range reduction. But by symmetry of e^x, I can actually use [0, 1]. Thank you!
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.9.4...
I honestly don't think this has anything to do with your gender or the gender of the person you're replying to.