25
$\begingroup$

Given a real function of real variables, is there software available that can automatically generate numerically-accurate code to calculate the function over all inputs on a machine equipped with IEEE 754 arithmetic?

For example, if the real function to be evaluated were:

f(a, b, c) = \frac{-b - \sqrt{b^2 - 4ac}}{2a}

The software would consider catastrophic cancellation and possibly output table lookups for certain sets of inputs to avoid a loss in computational accuracy.

Alternatively, is there software that can generate a pure table-based lookup routine to calculate a given function to high accuracy?

Improve this question
$\endgroup$
3
  • 5
    $\begingroup$ Hard problem in general. $\endgroup$
    – dmckee --- ex-moderator kitten
    Commented Dec 26, 2011 at 22:31
  • 1
    $\begingroup$ If the problem had specifically been about the root computation (or factorization) of polynoms, there are some C (or C++) libraries out there. $\endgroup$
    – moala
    Commented Dec 27, 2011 at 1:11
  • 2
    $\begingroup$ You might want to check out Richard Harris' excellent series of articles in the ACCU journal Overload about The Floating Point Blues. I indexed them over on Programmers.SX for people who might be interested. $\endgroup$
    – Mark Booth
    Commented Jan 13, 2012 at 11:00

5 Answers 5

Reset to default
25
$\begingroup$

The best solution that I know of is to program the symbolic expressions in Mathematica, Maple, or SymPy; all of the links go directly to the code generation documentation. All of the programs above can generate code in C or Fortran.

None of the programs above mentions accuracy in IEEE 754 arithmetic; in general, it would be difficult to anticipate all sources of catastrophic cancellation, as @dmckee notes. It's hard to replace human expertise in numerical analysis.

To provide a concrete example, consider calculating the trigonometric functions to high precision for arbitrary inputs in $[0, 2\pi]$. There are many strategies for doing so, some even hardware dependent, as see in the Wikipedia article Trigonometric Tables. All of the algorithms require ingenuity and numerical analysis, even algorithms that depend on lookup tables and Taylor series or interpolation (see the Wikipedia article The Table-Maker's Dilemma). For more detail, see the related Stack Overflow question How do Trigonometric Functions work?.

Software that generated code or routines to calculate arbitrary functions to high accuracy would not only need to be aware of cancellation errors, but also series approximants (Taylor, Padé, Chebyshev, rational, etc.) for calculating functions that are not defined in terms of a finite number of additions, subtractions, multiplications, divisions, and bit shifts. (See Approximation Theory.)

Improve this answer
$\endgroup$
3
  • 5
    $\begingroup$ "It's hard to replace human expertise in numerical analysis." - this alone deserves a +1. $\endgroup$
    – J. M.
    Commented Dec 27, 2011 at 1:03
  • $\begingroup$ "It's hard" isn't the same thing as "it's impossible". There are "full employment theorems" for some jobs (e.g. compiler writers). Is there one for numeric analysts? $\endgroup$
    – Pseudonym
    Commented Mar 8, 2013 at 4:19
  • $\begingroup$ Yes. Rice's Theorem. $\endgroup$
    – Geoff Oxberry
    Commented Mar 8, 2013 at 4:44
14
$\begingroup$

If you want an idea of just how far we are away from such a software package, please look at the 2001 LAPACK working note on computing Givens rotations reliably and efficiently. I would expect most non-specialists (and many specialists!) in numerical analysis to be surprised at just how much analysis went into solving such an ostensibly simple problem:

Given $f,g \in \mathbb{C}$, find $c \in \mathbb{R}$ and $s \in \mathbb{C}$ such that

$$R(c,s) \left[\begin{array}{c} f \\ g \end{array}\right] = \left[\begin{array}{cc} c & s \\ -\bar s & c \end{array}\right] \left[\begin{array}{c} f \\ g \end{array}\right] = \left[\begin{array}{c} r \\ 0 \end{array}\right]$$,

where $R(c,s)$ is unitary. Balancing reliability with computational efficiency along with more subtle issues like continuity is highly nontrivial and unlikely to be automated in the forseeable future.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ +1 This is a great example, thank you. I guess if a solution for reals existed, then it might be adapted to complex numbers. $\endgroup$
    – Daniel Trebbien
    Commented Dec 29, 2011 at 20:26
  • $\begingroup$ I should probably mention that the fundamental difficulty is not in the fact that s can be complex, but in avoiding unnecessary overflow and/or underflow. It is related to the hypot function: en.wikipedia.org/wiki/Hypot $\endgroup$
    – Jack Poulson
    Commented Dec 30, 2011 at 23:44
11
$\begingroup$

Code generation and precompilation of mathematical expressions is becoming more popular.

While symbolic packages like SymPy, Mathematica, and Maple may include code generation I'm not confident that any of them also think hard about numerics.

There are a couple of other projects one could look into that are interested both in symbolics and in numerics.

Theano is such a project focused around array operations. They do identify and replace some operations known to be numerically ill-conditioned. I'm not certain that this includes your specific case but it's worth looking into.

Spiral might also be interesting to you. They also precompile an abstract syntax tree and also look out for numeric issues. They are more concerned with scalar operations (like your example). However they are also fairly specialized to a particular domain.

The growth in this field is encouraging however. One can be optimistic that your question will have a better answer in a few years.

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Agreed; perhaps my answer came off as too pessimistic, as there are plenty of domain specific solutions, but the general problem is...hard. $\endgroup$
    – Jack Poulson
    Commented Dec 29, 2011 at 19:18
4
$\begingroup$

Not in general, I can safely say the implementor of the code generator in SymPy didn't even try =P.

Paolo Bientinesi developed a method for generating stability proofs of linear algebra algorithms, which are generated using Robert van de Geijn's FLAME notation.

See this paper, or a longer, working note version.

Improve this answer
$\endgroup$
1
$\begingroup$

Sage lets you express formulas in Cython (a variant of python that generates C code); however, in answer to your more general question: no. Consider Rice's Theorem.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.