Method and apparatus for computing an approximation to the reciprocal of a floating point number in IEEE format

Abstract

A method and apparatus allows the quick computation of an estimate of the reciprocal of a floating point number in IEEE format. A table with 2k entries allows the computation of an estimate with 2×k+3 good bits. x is a floating point number in IEEE format for which a reciprocal approximation is to be computed. {circumflex over (x)} be a floating point number in IEEE format derived from x by leaving the sign bit unchanged, complementing the exponent bits, leaving the first k fraction bits unchanged, and complementing the remaining fraction bits. t is another floating point number in IEEE format found by using the first k bits of the fraction of x as an index into a table with 2k entries. The product {circumflex over (x)}×t computed with IEEE floating point arithmetic is an estimate of the reciprocal of x with 2×k+3 good bits (i.e., relative error less than 2−2×k+3).

Description

BACKGROUND OF THE INVENTION

[0001] 1. Field of the Invention

[0002] The present invention generally relates to floating point processors in high speed digital computers and, more particularly, to an improved method to compute an approximation of the reciprocal of a floating point number in the IEEE (Institute of Electrical and Electronics Engineers) format.

[0003] 2. Background Description

[0004] Floating point arithmetic units are commonly used in central processing units (CPUs) of digital computers, especially high performance superscalar processors such as reduced instruction set computers (RISC). An example of a RISC processor is the IBM PowerPC® used in the IBM RISC System/6000 computer. Such processors typically include an input/output (I/O) unit interfacing with an instruction/data cache and a branch processor, one or more fixed point processors and one or more floating point processors.

[0005] The semantics of floating point instructions has not been as clear cut as the semantics of the rest of the instruction set. To address this problem, the computer industry has standardized on the floating point format by IEEE standard 754-1985. The IEEE standard defines 32-bit and 64-bit floating point formats. Each consists of (from left to right) a sign bit, an e-bit exponent and an f-bit fraction. The exponent is assumed to be biased with bias 2(e-1)-1 and it is assumed that an implicit 1 is to be appended to the front (left) of the fraction. The IEEE 32-bit format has e=8 and f=23. The IEEE 64-bit format has ee=11 and f=52. The present invention will work for these and similar formats (with different values for e and f). The IEEE standard also defines the result of arithmetic operations on floating point numbers in these formats.

[0006] Of the four basic arithmetic operations of addition, subtraction, multiplication, and division, division is the most difficult to implement efficiently. One way of performing a floating point division x/y is to compute it as x×(1./y). This can be done by finding a relatively low accuracy initial approximation to 1./y, using an iterative method to refine the low accuracy approximation to the desired accuracy, and then multiplying by x. This method is particularly attractive when it is unnecessary to obtain the exact IEEE specified value for x/y. However, the method requires a method for finding an initial approximation to the reciprocal of the floating point number y. The present invention provides such a method.

[0007] The problem of finding a good easily computed initial approximation to the reciprocal of a floating point number is well known and has been solved in various ways in both hardware and software. The IBM PowerPC® architecture defines a FRES operation which has generally been implemented using a simple lookup table. A table size of 2k will produce an approximation with k+1 good bits. More complicated hardware methods are also known. An initial approximation with 4.3 good bits can be computed with a single integer subtraction. This method has been used in IBM's vector MASS software.

[0008] U.S. Pat. Nos. 5,563,818 to Agarwal et al., 6,163,791 to Schmookler et al., and 6,240,433 to Schmookler et al. include methods for finding an initial approximation to the reciprocal of a floating point number that are similar to the present invention in that they combine table lookup with a floating point multiply-add to produce the initial approximation. However, the Agarwal et al. and Schmookler et al. methods are relatively complex in their implementations. What is needed is a simpler method which simplifies both hardware and software implementations and results in a quicker computation of an estimate of the reciprocal of a floating point number in IEEE format.

SUMMARY OF THE INVENTION

[0009] It is therefore an object of the present invention to provide a method which allows the quick computation of an estimate of the reciprocal of a floating point number in IEEE format.

[0010] The method according to the invention produces a better combination of speed, accuracy and simplicity than other known methods. The invention produces twice the number of correct bits as a simple table lookup (for the same size table) with the additional cost of a floating point multiply. Since it takes two floating point multiply-add instructions to double the accuracy of the initial approximation (using an iterative method), this is a significant improvement over the prior art. The advantages of the invention are particularly apparent for software implementations where simplicity is very important as it reduces the number of instructions required to generate the initial approximation.

[0011] In the practice of the invention, a table with 2k entries allows the computation of an estimate with 2×k+3 good bits (as compared with k+1 good bits for a simple table lookup). By way of illustrative example, let x be a floating point number in IEEE format for which a reciprocal approximation is to be computed. Let {circumflex over (x)} be a floating point number in IEEE format derived from x as follows: Leave the sign bit unchanged, complement the exponent bits, leave the first k fraction bits unchanged, complement the remaining fraction bits. Note that {circumflex over (x)} can be computed from x with a single exclusive OR operation. Let t be another floating point number in IEEE format found by using the first k bits of the fraction of x as an index into a table with 2k entries. Then, if the table is chosen correctly, the product {circumflex over (x)}×t computed with IEEE floating point arithmetic will be an estimate of the reciprocal of x with 2×k+3 good bits (i.e., relative error less than 2−2×k+3).

BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The foregoing and other objects, aspects and advantages will be better understood from the following detailed description of a preferred embodiment of the invention with reference to the drawings, in which:

[0013] FIG. 1 is a block diagram illustrating the principle components and data flow of the preferred embodiment of a hardware implementation of the invention; and

[0014] FIG. 2 is a flow diagram illustrating the logic of a software implementation of the preferred embodiment of a software implementation of the invention.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION

[0015] Referring now to the drawings, and more particularly to FIG. 1, there is shown a register 11 which receives the input x. The input x is a binary floating point number in the IEEE format. This consists of (from left to right) a 1-bit sign, an e-bit exponent, and an f-bit fraction. The exponent is assumed to be biased with bias 2e-1-1 and it is assumed that an implicit 1 is to be added at the front of the fraction. The IEEE 32-bit format has e=8 and f=23. The IEEE 64-bit format has e=11 and f=52.

[0016] The value x is input to function block 12 which extracts an index i which, in turn, is used to address a lookup table 13. The lookup table 13 has 2k entries, where k is an adjustable parameter. A larger value of k gives a better approximation but requires a larger lookup table. The index i consists of the leftmost k fraction bits of x. It can assume values of 0 through 2k-1 inclusive.

[0017] A register 14 holds a mask m which depends on k. The mask m consists of (from left to right) one zero bit corresponding to the sign bit of x, e one bits corresponding to the exponent bits of x, k zero bits corresponding to the left most k fraction bits of x, and f-k one bits corresponding to the remaining fraction bits of x. The value x from register 11 and the mask m from the register 14 are combined in a bitwise exclusive OR 15 to produce {circumflex over (x)}, the bitwise exclusive OR of x and m which consists of (from left to right) the sign bit of x, the complement of the exponent bits of x, the first k fraction bits of x and the complement of the remaining f-k fraction bits of x.

[0018] The lookup table 13 contains 2k entries in the same format as x. Let δ=2.−k, ε=2.−f and z=1+i×δ+0.5(δ−ε), where i is the index used to address the lookup table. Then, the ith entry, T(i), in the lookup table is 1$\frac{4}{\left(8z^2-{\left(\delta -\epsilon \right)}^2-{\epsilon }^2\right)}.$

[0019] This is output to IEEE floating point multiplier 16, which also receives {circumflex over (x)} from exclusive OR 15. The multiplier output at register 17 is an estimate of the reciprocal of x which is r=T(i)×{circumflex over (x)}.

[0020] FIG. 2 is a flow diagram showing the logic of a software implementation of the invention. It will be understood from the following description that the software implements the process performed by the hardware shown in FIG. 1. The process is initiated in function blocks 21 and 22 where the input x is loaded into a first fixed point register and the mask m is loaded into a second fixed point register, respectively. Then, in function block 23, the bitwise exclusive OR of x and m is computed, i.e., {circumflex over (x)}=Λm, and {circumflex over (x)} is loaded into a third fixed point register. In function block 24, {circumflex over (x)} is moved to a floating point register. Then, in function block 25, the index i, consisting of the leftmost k bits of the fractional part of x, is isolated and moved to a fixed point register. In function block 26, the index i is used to load the ith entry of a table, T(i), into a floating point register. Now, the reciprocal r is computed as r={circumflex over (x)}×T(i) and the result is loaded in a floating point register in function block 27. Finally, in function block 28, the reciprocal r is output as an approximation to 1./x.

[0021] While the invention has been described in terms of preferred hardware and software embodiments, those skilled in the art will recognize that the invention can be practiced with modification within the spirit and scope of the appended claims.

Claims

1. A method of computing an approximation to the reciprocal of a floating point number x comprising the steps of: generating a floating point number {circumflex over (x)} by leaving a sign bit of x unchanged, complementing exponent bits of x, leaving a first k fraction bits of x unchanged, and complementing remaining fraction bits of x; generating another floating point number t found by using the first k bits of the fraction of x as an index into a table with 2k entries; and computing the product {circumflex over (x)}×t computed with floating point arithmetic as an estimate of the reciprocal of x with 2×k+3 good bits.

2. The method of computing an approximation to the reciprocal of a floating point number x as recited in claim 1, wherein the floating point number x is in IEEE format, the floating point number {circumflex over (x)} is in IEEE format, the floating point number t is in IEEE format, and the product {circumflex over (x)}×t is in IEEE format.

3. The method of computing an approximation to the reciprocal of a floating point number x as recited in claim 1, wherein {circumflex over (x)} is computed from x with a single bitwise exclusive OR operation with a mask m.

4. The method of computing an approximation to the reciprocal of a floating point number x as recited in claim 3, wherein the mask m consists of one zero bit corresponding to the sign bit of x, e one bits corresponding to the exponent bits of x, k zero bits corresponding to the first k fraction bits of x and f-k one bits corresponding to the remaining f-k fraction bits of x.

5. The method of computing an approximation to the reciprocal of a floating point number x as recited in claim 1, wherein the table a lookup table containing 2k entries in the same format as x, the ith entry, T(i), in the lookup table is

6. An apparatus for computing an approximation to the reciprocal of a floating point number x comprising: a first register storing the number x comprising a 1-bit sign, an e-bit exponent, and an f-bit fraction; a second register storing a mask m mask comprised of one zero bit corresponding to the sign bit of x, e one bits corresponding to the exponent bits of x, k zero bits corresponding to the first k fraction bits of x and f-k one bits corresponding to the remaining f-k fraction bits of x; a bitwise exclusive OR function combining the number x from the first register and the mask m from the second register to produce a bitwise exclusive OR output {circumflex over (x)}=xΛm; an index extractor extracting the first k bits of the fraction of x as an index i; a lookup table with 2k entries responding to the index i and outputting an entry T(i) as the ith entry of the table; and a multiplier computing the product {circumflex over (x)}×T(i) with floating point arithmetic an output r=T(i)×{circumflex over (x)} as an estimate of the reciprocal of x with 2×k+3 good bits.

7. The apparatus for computing an approximation to the reciprocal of a floating point number x as recited in claim 1, wherein the floating point number x is in IEEE format, the floating point number {circumflex over (x)} is in IEEE format, the floating point number T(i) is in IEEE format, and the product {circumflex over (x)}×T(i) is in IEEE format.

8. The apparatus for computing an approximation to the reciprocal of a floating point number x recited in claim 6, wherein the lookup table contains 2k entries in the same format as x, the ith entry, T(i), in the lookup table is