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

Question

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.

Solution 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.

Solution 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

Solution 3

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

Solution 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.

Solution 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.

View more solutions

Comments

[cHao]

sin as a syscall? That'd be quite a waste of kernel space, unless the OS in question is extremely graphical -- and even then, it'd be slower than a RTL implementation due to context switching. There are instructions that can calculate sin, cos, etc, but syscalls? I doubt it.

[wkl]

+1 for the glibc sourceware link, but wow the site is slow right now. (edited)

[AntonyFalegk]

Afaik these are the slow versions with faster implementations in arch. specific subfolders.

[wnoise]

Almost never. The functions are normal C calls, provided by the compiler orC library. For systems without FPUs,the instructions these functions use may be be trapped by the OS and then emulated, but this is a rare case.

[KriptSkitty]

Math functions implemented as system calls? Really??

[AntonyFalegk]

LOL, internet is always slow here can't see a difference ;)

[user502515]

No. Whatever can be executed from user mode will be done from there. It's pointless to context-switch to kernel mode for something that is not a privileged operation.

[Ben Jackson]

+1 for actually explaining the math. I felt a lot better when I realized trig functions were just truncated Taylor series expansions. Otherwise the approximations look like serious magic!

[Mateen Ulhaq]

@Thomas Especially if the OS is not Windows, Mac, Linux, and any other OS out there. (Do calculators count as OSes?)

[Stephen Canon]

@Ben: good libraries typically do not use truncated taylor series -- other polynomial approximations (Minimax, Chebyshev, Padé) have much more desirable error characteristics and allow one to get the same accuracy with fewer arithmetic operations.

[Stephen Canon]

@Alexandre: rational approximation has mostly gone out of fashion for exp and log because multiplication and addition are so much faster than division on modern hardware. It is still used for stiffer functions, however.

[Janus]

Abramowitz & Stegun is the colloquial name of the NIST publication "Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables" from 1964. Please link to the 2010 edition "NIST Digital Library of Mathematical Functions" freely available online at dlmf.nist.gov

[Dr. belisarius]

@Janus "The" Abramowitz is such a classic ... I almost shed a few tears for founding it cited in SO ...

[Janus]

@belisarius If you haven't already, go and check out the new version I link to. Very nice to always have the authoritative reference available, although I must admit I haven't checked in detail how much is changed.

[Dr. belisarius]

@Janus Of course I did! My comment was just a nostalgic one :)

[Alexandre C.]

@Janus: The NIST dlmf references Abramowitz & Stegun for the numeric tables, or numerical coefficients that you sometimes crucially need (or eg. are too lazy to devise a good scheme to compute them inductively)

[Piskvor left the building]

@muntoo: Probably not in most cases. There is some sort of limited OS for the high-end calculators like TI-89, but I'd guess that 1) the question was about general-purpose computing devices and 2) calculators could implement the math functions directly in the HW (which is what the low-end ones are doing).

[Alexandre C.]

Good point about the corner cases. They can easily become bottlenecks in some cases (the implementation of ldexp on MSVC renders the function pretty much useless for instance)

[Alexandre C.]

@Piskvor: IIRC on the TI89 there were no FPU and everything was provided by the OS (as binary coded decimal !)

[Sandeep Kumar]

Welcome to SO, please refer to this for better presenting your question/answer: meta.stackoverflow.com/a/251362/3519504