Fast Trigonometric Functions using the RLIBM Approach

TL;DR

Problem: delivering correctly rounded trigonometric functions across multiple FP representations and rounding modes is difficult because range reduction hinges on precise handling of $2^t/\pi$. Approach: RLibm formulates the task as approximating the correctly rounded result and solves a linear program to enforce rounding constraints, extended to multiple representations via an $(n+2)$-bit representation with round-to-odd. Key contributions include FP- and integer-based range reduction strategies that preserve sufficient bits of $2^t/\pi$, plus a successful implementation for sin, cos, and tan that is fast and correctly rounded from 10- to 32-bit representations. Significance: the work yields portable, high-precision trig functions competitive with mainstream libraries, with practical impact for multi-representation computing and future GPU support.

Abstract

This paper describes our experience developing polynomial approximations for trigonometric functions that produce correctly rounded results for multiple representations and rounding modes using the RLIBM approach. A key challenge with trigonometric functions concerns range reduction with "pi", which reduces a given input in the domain of a 32-bit float to a small domain. Any rounding error in the value of "pi" is amplified during range reduction, which can result in wrong results. We describe our experience implementing fast range reduction techniques that maintain a large number of bits of "pi" both with floating-point and integer computations. The resulting implementations for trigonometric functions are fast and produce correctly rounded results for all inputs for multiple representations up to 32-bits with a single implementation.

Figures (7)

• Figure 1: (A) The rounding interval of a correctly rounded result v2. (B) The linear inequality generated with the RLibm approach to produce the correctly rounded result for an input $x$ given the rounding interval $[l, h]$.
• Figure 1: (a) High level range reduction producing a reduced input $x' \in [-\frac{\pi}{512}, \frac{\pi}{512}]$ and the index $k'$ for pre-computed tables. Here, $[\frac{256x}{\pi}]$ computes the integer nearest to $\frac{256x}{\pi}$. (b) A FP-based implementation of Algorithm \ref{['alg:mysin']} that produces reduced inputs for small inputs (i.e., $|x| < 2^{30}$). The values $xp$ and $kp$ denote $x'$ and $k'$ respectively.
• Figure 2: (A) Illustration of round-to-odd rounding mode with a 34-bit FP representation. (B) Intuition on why double rounding is innocuous with the round-to-odd mode. Rounding the round-to-odd result to a target representation has the same truncated value, rounding bit, and sticky bit as rounding that real value directly to the given target representation.
• Figure 5: (a) The FP-based implementation of Algorithm \ref{['alg:mysin']} for large inputs. We denote the reduced input and the lower 9-bits of $k$ as $xp$ and $kp$, respectively. We store 7 28-bit pieces of $\frac{256}{\pi}$ in a double array named $\_256\_over\_pi\_28$. We maintain the exponent of the last bit of each piece in the array $\_256\_over\_pi\_28\_exp$. Here, $k$ and $k\_hi\_int$ are signed 64-bit integers, $idx$ is a signed 32-bit integer, and the rest are FP double values. The function $bits(z)$ provides the IEEE-754 FP bit-pattern of the input as an integer. (b) An alternative FP-based strategy for large inputs. We store 4 53-bit pieces of $\frac{256}{\pi}$ in a double array named $\_256\_over\_pi\_53$. The fused-multiply-add ($fma$) operation that performs both FP multiplication and addition with a single instance of rounding.
• Figure 6: Range reduction for small inputs when $\frac{\pi}{128} \leq |x| < 2^{30}$ and $t=8$. Here, we are creating a reduced input $x' \in [-\frac{\pi}{512}, \frac{\pi}{512})$. (A) Computation with 80 bits of x * 256_Over_Pi with 40 bits each in two 64-bit integers P1 and P0. (B) Computing $k$ and $r$ with integer and bitwise operations. (c) Intuition for finding the bit whose exponent is 0 when we perform computation with integers and do not explicitly maintain exponents.