Where to find algorithms for standard math functions?

Question

I'm looking to submit a patch to the D programming language standard library that will allow much of std.math to be evaluated at compile time using the compile-time function evaluation facilities of the language. Compile-time function evaluation has several limitations, the most important ones being:

1. You can't use assembly language.
2. You can't call C code or code for which the source is otherwise unavailable.

Several std.math functions violate these and compile-time versions need to be written. Where can I get information on good algorithms for computing things such as logarithms, exponents, powers, and trig functions? I prefer just high level descriptions of algorithms to actual code, for two reasons:

1. To avoid legal ambiguity and the need to make my code look "different enough" from the source to make sure I own the copyright.

2. I want simple, portable algorithms. I don't care about micro-optimization as long as they're at least asymptotically efficient.

Edit: D's compile time function evaluation model allows floating point results computed at compile time to differ from those computed at runtime anyhow, so I don't care if my compile-time algorithms don't give exactly the same result as the runtime version as long as they aren't less accurate to a practically significant extent.

Answer 1

See the book "Elementary Functions: Algorithms and Implementation" by Jean-Michel Muller.

Answer 2

Maybe this helps you (at least for some of the functions): http://en.wikipedia.org/wiki/CORDIC

Answer 3

several sources, including:

Abramowitz and Stegun, "Handbook of Mathematical Functions" (available online!)

Hart, "Computer Approximations" (out of print but good)

also see the several other SO questions about trigonometry, including "How do Trigonometric Functions Work?" and "Trigonometric Functions on Embedded Systems".

Answer 4

The source that I recommend is Numerical Methods for Scientists and Engineers by R. W. Hamming.

This book is published by Dover Press, and is an inexpensive paperback.

Answer 5

Jean-Michel Muller's book is an excellent recommendation, as is Hart.

Is it actually necessary for you to own the copyright? It's usually a bad idea to get into the business of writing math library functions if you can avoid it (and I say that as someone who does so professionally). I don't know whether or not D can take in BSD-licenced code, but there are several good open-source implementations that may prove helpful. You may want to look at Sun's FDLIBM, for example. Stephen Moshier's Cephes would also be a possibility, though its licensing situation is a little bit odd, but I believe that he's been willing to let people redistribute his code under other licenses in the past.

On a side-note, unless you're supporting arbitrary-precision floating-point (I don't think that D does), there isn't usually a notion of "asymptotic efficiency" for libm functions.

Answer 6

John Hart Computer Approximations 1968 by John Wiley & Sons.

The calculations ideally should match precisely what they would if done at runtime. That can be tricky. For many functions, no series converges quickly over the full domain, so algoritms paste together various methods.

Also, there are various floating point formats. Most platforms (I think) now use IEEE 754. When I wrote a compiler ca. 1985, I had to deal with cross-platform floating point formats. It was very tedious to get it right, because you have to piece the numbers together bit by bit, being sure that you get precisely the value that would be calculated on the target machine. I don't know if you have to deal with that.

Answer 7

See Stand-alone code for numerical computing for links to code for a few special functions and for random number generation. All the code there is public domain. The code is implemented in C++ and Python, but it's easy to translate to any other language.

Answer 8

As you'd expect, similar issues arise in other languages:

http://java.sun.com/j2se/1.5.0/docs/api/java/lang/StrictMath.html

To help ensure portability of Java programs, the definitions of some of the numeric functions in this package require that they produce the same results as certain published algorithms. These algorithms are available from the well-known network library netlib as the package "Freely Distributable Math Library," fdlibm. These algorithms, which are written in the C programming language, are then to be understood as executed with all floating-point operations following the rules of Java floating-point arithmetic.

I don't know what D's rules are for runtime calculation of math functions, but you may be able to pull a similar trick - re-interpret the C source of fdlibm as D. If D calls platform-specific C libraries, then you have the problem that it may be impossible at compile time to predict the runtime value.

I think fdlibm's license is very permissive, you'd have to check for yourself whether it's suitable for redistribution in D. One version I've seen requires a copyright notice to be preserved, and that's it.