Definitions of sqrt, sin, cos, pow etc. in cmath

问题:

Are there any definitions of functions like sqrt(), sin(), cos(), tan(), log(), exp() (these from math.h/cmath) available ?

I just wanted to know how do they work.

回答1:

This is an interesting question, but reading sources of efficient libraries won't get you very far unless you happen to know the method used.

Here are some pointers to help you understand the classical methods. My information is by no means accurate. The following methods are only the classical ones, particular implementations can use other methods.

• Lookup tables are frequently used
• Trigonometric functions are often implemented via the CORDIC algorithm (either on the CPU or with a library). Note that usually sine and cosine are computed together, I always wondered why the standard C library doesn't provide a sincos function.
• Square roots use Newton's method with some clever implementation tricks: you may find somewhere on the web an extract from the Quake source code with a mind bogging 1 / sqrt(x) implementation.
• Exponential and logarithms use exp(2^n x) = exp(x)^(2^n) and log2(2^n x) = n + log2(x) to have an argument close to zero (to one for log) and use rational function approximation (usually Padé approximants). Note that this exact same trick can get you matrix exponentials and logarithms. According to @Stephen Canon, modern implementations favor Taylor expansion over rational function approximation where division is much slower than multiplication.
• The other functions can be deduced from these ones. Implementations may provide specialized routines.
• pow(x, y) = exp(y * log(x)), so pow is not to be used when y is an integer
• hypot(x, y) = abs(x) sqrt(1 + (y/x)^2) if x > y (hypot(y, x) otherwise) to avoid overflow. atan2 is computed with a call to sincos and a little logic. These functions are the building blocks for complex arithmetic.
• For other transcendental functions (gamma, erf, bessel, ...), please consult the excellent book Numerical Recipes, 3rd edition for some ideas. The good'old Abramowitz & Stegun is also useful. There is a new version at http://dlmf.nist.gov/.
• Techniques like Chebyshev approximation, continued fraction expansion (actually related to Padé approximants) or power series economization are used in more complex functions (if you happen to read source code for erf, bessel or gamma for instance). I doubt they have a real use in bare-metal simple math functions, but who knows. Consult Numerical Recipes for an overview.

回答2:

Every implementation may be different, but you can check out one implementation from glibc's (the GNU C library) source code.

edit: Google Code Search has been taken offline, so the old link I had goes nowhere.

The sources for glibc's math library are located here:

http://sourceware.org/git/?p=glibc.git;a=tree;f=math;h=3d5233a292f12cd9e9b9c67c3a114c64564d72ab;hb=HEAD

回答3:

Have a look at how glibc implements various math functions, full of magic, approximation and assembly.

回答4:

Definitely take a look at the fdlibm sources. They're nice because the fdlibm library is self-contained, each function is well-documented with detailed explanations of the mathematics involved, and the code is immensely clear to read.

回答5:

Having looked a lot at math code, I would advise against looking at glibc - the code is often quite difficult to follow, and depends a lot on glibc magic. The math lib in FreeBSD is much easier to read, if somehow sometimes slower (but not by much).

For complex functions, the main difficulty is border cases - correct nan/inf/0 handling is already difficult for real functions, but it is a nightmare for complex functions. C99 standard defines many corner cases, some functions have easily 10-20 corner cases. You can look at the annex G of the up to date C99 standard document to get an idea. There is also a difficult with long double, because its format is not standardized - in my experience, you should expect quite a few bugs with long double. Hopefully, the upcoming revised version of IEEE754 with extended precision will improve the situation.

回答6:

Most modern hardware include floating point units that implement these functions very efficiently.

回答7:

Those are almost always implemented as system calls. If you want to look at the sources, you'd need access to the OS sources, which means you need to look at an open-source OS like Linux or BSD.