trigonometric functions sin, cos, tan as well as for the exponential function. Keywords: range reduction, elementary functions, computer arithmetic. 1. Introduction.
Journal of VLSI Signal Processing 33, 31–35, 2003 c 2003 Kluwer Academic Publishers. Manufactured in The Netherlands.
On-the-Fly Range Reduction∗ ` VINCENT LEFEVRE INRIA, Projet SPACES, LORIA, Campus Scientifique, B.P. 239, 54506 Vandoeuvre-l`es-Nancy Cedex, France JEAN-MICHEL MULLER CNRS, Projet CNRS/ENS Lyon/INRIA Arenaire, Laboratoire LIP, Ecole Normale Sup´erieure de Lyon, 46 All´ee d’Italie, 69364 Lyon Cedex 07, France Received October 30, 2000; Revised July 26, 2001
Abstract. In several cases, the input argument of an elementary function evaluation is given bit-serially, most significant bit first. We suggest a solution for performing the first step of the evaluation (namely, the range reduction) on the fly: the computation is overlapped with the reception of the input bits. This algorithm can be used for the trigonometric functions sin, cos, tan as well as for the exponential function. Keywords: range reduction, elementary functions, computer arithmetic 1.
Introduction
The algorithms used for evaluating the elementary functions only give a correct result if the argument is within some bounded interval. To evaluate an elementary function f (x) (sine, cosine, exponential, . . .) for any x, one must find some “transformation” that makes it possible to deduce f (x) from some value g(y), where – y, called the reduced argument, is deduced from x; – y belongs to the convergence domain of the algorithm implemented for the evaluation of g. With the usual functions, the only cases for which reduction is not straightforward are the cases where y is equal to x − nC, where n is an integer and C a constant (for instance, for the trigonometric functions, C is a multiple of π /8). Example 1 (Computation of the cosine function). Assume that we want to evaluate cos(x), and that the convergence domain of the algorithm used to evaluate
the sine and cosine of the reduced argument contains [0, +π/4]. We choose C = π/4, and the computation of cos(x) is decomposed in three steps: – Compute y and n such that y ∈ [0, +π/4] and y = x − nπ/4; – Compute g(y, n) cos(y) √ 2 (cos(y) − sin(y)) 2 − sin(y) √ 2 − (cos(y) + sin(y)) 2 = − cos(y) √ 2 (− cos(y) + sin(y)) 2 sin(y) √ 2 (cos(y) + sin(y)) 2
if n mod 8 = 0 if n mod 8 = 1 if n mod 8 = 2 if n mod 8 = 3 if n mod 8 = 4 if n mod 8 = 5 if n mod 8 = 6 if n mod 8 = 7 (1)
∗ This paper is an extended version of a communication to the SPIE’s
45th annual meeting, San Diego, Aug. 2000.
– Obtain cos(x) = g(y, n).
32
Lef`evre and Muller
Example 2 (Computation of the exponential function). Assume that we want to evaluate e x in a radix-2 number system, and that the convergence domain of the algorithm used to evaluate the exponential of the reduced argument contains [0, ln(2)]. We can choose C = ln(2), and the computation of e x is then decomposed in three steps: – Compute y ∈ [0, ln(2)] and n such that y = x − n ln(2); – Compute g(y) = e y ; – Compute e x = 2n g(y).
– 0 ≤ y < C; – n = (x − y)/C is an integer. We also define, for each i, m i (also called 2i mod C) as the unique value between 0 and C such that (2i − m i )/C is an integer. These notations give some constraints on x and C (e.g., C is less than 1, x is less than 2h+1 ). One can easily adapt the algorithms given in the rest of the paper to variables belonging to other domains. We chose these constraints to make the presentation of the algorithms simpler. 3.
Unless multiple-precision arithmetic is used during the intermediate calculations, a straightforward computation of y as x − nC is to be avoided, since this operation will lead to catastrophic cancellations (i.e., to very inaccurate estimates of y) when x is large or close to an integer multiple of C. Many algorithms have been suggested for performing the range reduction accurately [1–5]. Now, there are many cases (on special-purpose systems) where the input argument of a calculation is generated most significant digit first. This happens, for instance, when this argument is the result of a division or a square root obtained through a digit-recurrence algorithm [6, 7], the output of an on-line algorithm [8, 9], or when it is generated by an analog-to-digital converter. In the rest of this paper, we present an adaptation of the Modular Range Reduction Algorithm [3, 10] that accepts such digit serial inputs and performs the range reduction “on the fly”: most of the computation is overlapped with the reception of the input bits, and the reduced argument is produced almost immediately after reception of the last input bit. On-the-fly arithmetic algorithms have already been proposed by Ercegovac and Lang for rounding or converting a number from redundant to non-redundant representation [11, 12].
2.
Notations
In the rest of the paper, x = x h x h−1 · · · x0 .x−1 x−2 · · · x is the input argument, C = 0.C−1 C−2 · · · C− p is the constant of the range reduction (with − p ≤ ), and y = 0.y−1 y−2 · · · y− p is the reduced argument. We assume 1/2 ≤ C < 1 (this assumption is made to simplify the presentation. Modifying the algorithm for a constant larger than 1 or less than 1/2 is straightforward). These values satisfy:
Non-Redundant Algorithm
Algorithm 1 is by far less efficient than the “redundant” algorithm given later. We give it because it is simpler to understand, and because the other algorithm is derived from it. The basic idea is the following: at step i of the algorithm, when we receive input bit x h−i of x, we add x h−i × 2i mod C to an accumulator. If the accumulated value becomes larger than C, we subtract C from it. Let us call Ai+1 the value obtained after this operation. One can easily check that 0 ≤ Ai+1 < C and Ai+1 − x h x h−1 · · · x h−i × 2h−i is an integer multiple of C. Hence the final value stored in the accumulator is equal to the reduced argument y. Algorithm 1 Non-redundant algorithm. A0 = 0 for i = 0 to h − do Ti = Ai + x h−i m h−i if Ti < C then Ai+1 = Ti else Ai+1 = Ti − C y = Ah−+1 A possible variant consists in computing Ui = Ai + x h−i (m h−i − C) in parallel with Ti = Ai + x h−i m h−i , and then choosing Ai+1 equal to Ui if Ui ≥ 0, otherwise Ti . 4.
Redundant Algorithm
Now, to accelerate the reduction, we assume that we perform the accumulations with carry-save additions. The carry-save number system makes it possible to perform very fast, carry-free additions. On the other hand,
On-the-Fly Range Reduction
its intrinsic redundancy makes comparisons somewhat more complex. The accumulator will store the values Ai in carry-save. In the previous algorithm, we needed “exact” comparisons between the Ai ’s and C. Having the Ai ’s stored in carry-save makes these “exact” comparisons difficult. Instead of that, we will perform comparisons based on the examination of the first three carry-save positions of Ai only. This will not allow us to bound the Ai ’s by C. Nevertheless, we will show that the Ai ’s will be upper-bounded by C + 12 (therefore by 3 ), which will suffice for our purpose. We denote: 2 Ai =
(1) (2) (1) (2) ; Ai,−1 , Ai,−1 ; Ai,0 , Ai,0 (1) (1) (2) (2) Ai,−2 , Ai,−2 ; · · · ; Ai,− p , Ai,− p
where Ai,(1)j and Ai,(2)j are in {0, 1} and Ai =
p
Ai,(1)j + Ai,(2)j · 2− j .
j=0
The variable Ti of the non-redundant algorithm is used again, and is also represented in carry-save form: Ti =
(1) (2) (1) (2) Ti,0 , Ti,0 ; Ti,−1 , Ti,−1 ; (1) (1) (2) (2) Ti,−2 , Ti,−2 ; · · · ; Ti,− p , Ti,− p
Algorithm 2 Redundant algorithm. A0 = 0 for i = 0 to h − do Ti = Ai +cs x h−i m h−i (1) (2) (1) (2) (1) (2) , Ti,0 ; Ti,−1 , Ti,−1 ; Ti,−2 , Ti,−2 Tˆi = Ti,0 converted to non-redundant binary using a 3-bit adder if Tˆi < C then Ai+1 = Ti else Ai+1 = Ti −cs C B = Ah−+1 −cs C Convert Ah−+1 and B to non-redundant binary. if B < 0 then y = Ah−+1 else y=B This gives Algorithm 2. In the loop, we do not want to waste time with a full comparison to know whether we need to subtract
33
C from Ti or not. Thus we use a rough approximation Tˆ i to Ti based on the first three digits of Ti . Since
(1) (1) (2) (2) Ti,−3 ; · · · ; Ti,− , Ti,−3 p , Ti,− p ≤ 2 · 2−3 + 2 · 2−4 + · · · + 2 · 2− p
