Basic operations edit

In a math library, it is frequently useful to use type punning to interpret a floating-point number as an unsigned integer of the same size. This allows for faster inspection of certain numeral properties (positive or not) and number comparison. In more advanced cases, bit twiddling may be used to modify a number in a specific way.

For more exact operation, a double double or even triple double format may be used. In this case, a high-precision number is expressed as the sum of two or three floating-point numbers.^[5]

Transcendental functions edit

Transcendental functions such as log, exponential, and trig functions make up the backbone of any math library. These functions are generally implemented by a polynomial fit, usually a Taylor polynomial or a Chebyshev polynomial derived by the Remez algorithm (having the benefit of an improved error bound), but the pre-processing steps are equally important.^[5]

Trignometry edit

Range reduction (also argument reduction, domain-spltting) is the first step for any function, after checks for unusual values (infinity and NaN) are performed. The goal here is to reduce the domain of the argument for the polynomial to process, using the function's symmetry and periodicity (if any), setting flags to indicate e.g. whether to negate the result in the end (if needed). It is worth noting that periodic functions require higher-than-input precision when reducing, with the prototypical method being the Payne–Hanek–Corbett algorithm.^[6] After range reduction, near-zero values may be subject to a "fast path": for example, a tiny input x can be returned directly in sin, while 1 − |x| may be used for cos.

The next step is the evaluation of the polynomial, with a conventional Horner's method used. After that, the sign of the result can be flipped according to the information from the range-reduction routine before being returned.

Logarithm and exponential edit

Logarithm in base 2 is relatively straightforward, as the integer part k is already in the floating-point exponent; a preliminary range reduction is accordingly performed, yielding k. The mantissa x (where log2(x) is between -1/2 and 1/2) is then compared to a table and intervals for further reduction into a z with known log2 and an in-range x/z, along with polynomial coefficients used for the interval x/z is in. The result is then log(z) + log(x/z) + k.^[7] Logarithm in other bases follow a similar approach, with the differences being (a) a different table is used; (b) k needs to be multiplied by log2(new-base).^[7]

Exponential in base 2 is similarly the "base case" due to floating-point structure. The procedure is simply a combination of range reduction (usually by lookup) and a polynomial over the remaining mantissa.^[7] The natural exponent may be implemented with a separate table for precision, while exp10 can simply be exp(x × log2(10)) when in-range.^[7] Finally, the any-base exponent function pow() is built around exp() and log() in the general case.^[7]

• Intrinsic function
• List of numerical libraries