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

Patryk Czachurski picture Patryk Czachurski · Dec 27, 2010 · Viewed 17.7k times · Source

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.

Answer

Alexandre C. picture Alexandre C. · Dec 27, 2010

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.