<p>A common way to generate polynomial minimax approximation is to use the Remez exchange algorithm published by Russian mathematician Evgeny Remez in 1934. This is a numerical procedure that often involves ill-conditioned systems of equations. As a consequence, it is is usually implemented with the help of an arbitrary-precision library. For example, in the implementation of the Remez algorithm that I use I configured the library for 1024-bit precision.</p>
<p>For reasonably well-behaved functions, various variants of the Remez algorithm can find approximations very close to the mathematical minimax polynomial. The problem, as noted in the question, is what happens when one moves the generated coefficients of the polynomial to a finite-precision floating-point computation. One often finds the minimax property of the approximation impaired, sometimes significantly so. There are two sources of error at play. First, the generated coefficients cannot be represented accurately in the finite precision floating-point format. Second, the evaluation of the polynomial uses finite-precision operations instead of mathematical operations with infinite precision.</p>
<p>The first problem is the easier one to address. As one can see from some quick experiments, simply rounding the coefficients to the finite-precision format doesn't accomplish the desired near minimax result. By using a finite-precision format, we basically transform from an N-dimensional continuous space to an N-dimensional discrete lattice, and to do this properly, we need to find the closest lattice points. This is a solvable but hard problem, which is usually made easier through the use of heuristics. Relevant literature:</p>
<p>N. Brisebarre, J.-M. Muller, and A. Tisserand, &quot;Computing machine-efficient polynomial approximations&quot;. <em>ACM Transactions on Mathematical Software</em>, Vol. 32. No. 2, June 2006, pp. 236-256.
(<a href="http://perso.ens-lyon.fr/jean-michel.muller/TruncToms.pdf" rel="nofollow noreferrer">online</a>)</p>
<p>Nicolas Brisebarre and Sylvain Chevillard, &quot;Efficient polynomial L<sup>∞</sup>-approximations&quot;, In <em>18th IEEE Symposium on Computer Arithmetic</em>, June 2007, pp. 169-176
(<a href="https://hal.inria.fr/inria-00119513/document" rel="nofollow noreferrer">online</a>)</p>
<p>Florent de Dinechin and Christoph Lauter, &quot;Optimizing polynomials for floating-point implementation&quot;, ArXiv preprint 2008
(<a href="https://arxiv.org/pdf/0803.0439" rel="nofollow noreferrer">online</a>)</p>
<p>The <a href="http://sollya.gforge.inria.fr/" rel="nofollow noreferrer">Sollya tool</a> uses these techniques from the literature for its <code>fpminimax</code> <a href="http://sollya.gforge.inria.fr/sollya-7.0/help.php?name=fpminimax" rel="nofollow noreferrer">command</a>. Worth checking out in addition to Maple's and Mathematica's facilities for generating minimax polynomial approximations, as it often results in superior approximations in my experience.</p>
<p>The second problem, how to account for evaluation with finite-precision floating-point computation and how to adjust the coefficients of a polynomial approximation accordingly, are still subject to research. Some initial results:</p>
<p>Tor Myklebust, &quot;Computing accurate Horner form approximations to special functions in finite precision arithmetic&quot;, ArXiv manuscript 2015
(<a href="https://arxiv.org/pdf/1508.03211" rel="nofollow noreferrer">online</a>)</p>
<p>Denis Arzelier, Florent Bréhard, Mioara Joldes, &quot;Exchange algorithm for evaluation and approximation error-optimized polynomials&quot;, In <em>26th IEEE Symposium on Computer Arithmetic</em>, June 2019, pp. 30-37
(<a href="https://hal.archives-ouvertes.fr/hal-02006606/document" rel="nofollow noreferrer">online</a>)</p>
<p>Note that the first publication came about due to a <a href="https://stackoverflow.com/questions/26692859/best-machine-optimized-polynomial-minimax-approximation-to-arctangent-on-1-1">question</a> I asked on Stackoverflow.</p>
<p>For my own use I am using a heuristic search for finding approximations optimized to account for both representational error in the coefficients and evaluation error in the polynomial evaluation. it can be loosely described as a form of simulated annealing. I have also checked into the use of genetic programming, but the preliminary results did not look promising, so I stopped pursuing this approach.</p>
Transfer of approximation from double precision to single precision
<p>I would like to transfer a double precision approximation of a given function into a single precision C implementation (target device provides single precision ALU only).</p>
<p>It's not too complicated to generate a high precision (e.g. max error 0.1e-12) approximation using double precision. I have used maple minimax function, but I have as well found some implementations using double precision <a href="http://www.ganssle.com/approx.htm" rel="nofollow noreferrer">example</a>.</p>
<p>But as soon as it comes to transferring this approximation to a single precision method, I face a loss of accuracy, when I simply convert the coefficients to float. I'm aiming an approximation (single precision) which is precise for about +/-5 ulp. Simply converting the coefficients to float doesnt seem to do the job. I already learned to split up constants like pi/2 into a rounded part and an error part, and I think there is some kind of trick to transfer the coefficients (the core computation of the approximations are usually polynomials, i want to focus on them in this question), that I don't know yet.</p>
<p>I am grateful for every hint, paper concerning the <strong>transfer</strong> aiming an implementation. I have already studied some papers about <a href="https://randomascii.wordpress.com/2012/03/08/float-precisionfrom-zero-to-100-digits-2/" rel="nofollow noreferrer">float precision</a>, but for the last two weeks did not make to much progress.</p>
<p>Thanks in advance!</p>