# Patent

## Title
Calculation of trigonometric functions in an integrated circuit device

## Abstract

Circuitry for computing a tangent function of an input value includes first look-up table circuitry that stores pre-calculated tangent values of a limited number of sample values, circuitry for inputting bits of the input value of most significance as inputs to the first look-up table circuitry to look up one of the pre-calculated tangent values as a first intermediate tangent value, circuitry for calculating a second intermediate tangent value from one or more ranges of remaining bits of the input value, and circuitry for combining the first intermediate tangent value and the second intermediate tangent value to yield the tangent function of the input value.

## Background
CROSS REFERENCE TO RELATED APPLICATION
This is a continuation-in-part of copending, commonly-assigned U.S. patent application Ser. No. 12/823,539, filed Jun. 25, 2010, which is hereby incorporated by reference herein in its entirety.


FIELD OF THE INVENTION
This invention relates to calculating trigonometric functions in integrated circuit devices, and particularly in programmable integrated circuit devices such as programmable logic devices (PLDs).
BACKGROUND OF THE INVENTION
Trigonometric functions are generally defined for the relatively small angular range of 0-360°, or 0-2π radians. For angular values above 2π, the values of the trigonometric functions repeat. Indeed, one could restrict the range to 0-π/2, because various trigonometric identities can be used to derive values for trigonometric functions of any angle between π/2 and 2π from trigonometric functions of angles between 0 and π/2.
One method that may be used in integrated circuit devices for calculating trigonometric functions is the CORDIC algorithm, which uses the following three recurrence equations:

xn+1=xn−dnyn2−n 

yn+1=yn+dnxn2−n 

zn+1=zn−dn tan−1(2−n)

For example, to calculate a sine or cosine of an input, the x value is initialized to “1”, the y value is initialized to “0”, and the Z value is initialized to the angle required. Z is then rotated towards zero, which determines the sign of dn, which is ±1—if zn is positive, then so is dn, as the goal is to bring z closer to 0; if zn is negative, then so is dn, for the same reason. x and y represent the x and y components of a unit vector; as z rotates, so does that vector, and when z has reached its final rotation to 0, the values of x and y will have converged to the cosine and sine, respectively, of the input angle.

To account for stretching of the unit vector during rotation, a scaling factor is applied to the initial value of x. The scaling factor is:




  
    
      ∏
      
        n
        =
        0
      
      ∞
    
    ⁢
    
      
    
    ⁢
    
      
        1
        +
        
          2
          
            
              -
              2
            
            ⁢
            n
          
        
      
    
  
  =
  
    1.64676025812106564
    ⁢
    
      
    
    ⁢
    ...
  




The initial x is therefore be set to 1/1.64677 . . .=0.607252935 . . . .

However, CORDIC may become inaccurate as the inputs become small. For example, the actual value of sin(θ) approaches θ as θ approaches 0 (and therefore sin(θ) approaches 0), and the actual value of cos(θ) approaches 1 as θ approaches 0. However, the magnitude of the error between the calculated and actual values increases as e decreases.
Moreover, while CORDIC on initial consideration appears to be easily implemented in integrated circuit devices such as FPGAs, closer analysis shows inefficiencies, at least in part because of multiple, deep arithmetic structures, with each level containing a wide adder. Common FPGA architectures may have 4-6 input functions, followed by a dedicated ripple carry adder, followed by a register. When used for calculating floating point functions, such as the case of single precision sine or cosine functions, the number of hardware resources required to generate an accurate result for smaller input values can become large. Double precision functions can increase the required resources even further.
SUMMARY OF THE INVENTION
According to embodiments of the present invention, different trigonometric functions may be computed using various modified implementations that are based on different trigonometric identities that can be applied.
For sine and cosine functions, a modified CORDIC implementation changes small input angles to larger angles for which the CORDIC results are more accurate. This may be done by using π/2−θ instead of θ for small θ (e.g., for θ<π/4). As discussed above, CORDIC accuracy suffers for smaller angles, but a standard CORDIC implementation may be used for larger angles. A multiplexer can select between the input value θ and the output of a subtractor whose subtrahend and minuend inputs are, respectively, π/2 and the input value θ. A comparison of the input θ to a threshold can be used to control the multiplexer to make the selection. Both sine and cosine are computed by the x and y datapaths of the CORDIC implementation and the desired output path can be selected using another multiplexer, which may be controlled by the same comparison output as the input multiplexer. When π/2−θ has been used as the input, the identities cos(θ)=sin(π/2−θ) and sin(θ)=cos(π/2−θ) can be used to derive the desired result.
For the tangent function, the input angle can be broken up into the sum of different ranges of bits of the input angle, using trigonometric identities for the tangent of a sum of angles. Because some of the component ranges will be small, the identities will be simplified relative to those component ranges. The identities can be implemented in appropriate circuitry, as either a single-precision function or a double-precision function.
For the inverse tangent (i.e., arctan or tan−1) function, the problem is that the potential input range is between negative infinity and positive infinity (unlike, e.g., inverse sine or inverse cosine, where the potential input range is between −1 and +1). In accordance with the invention, trigonometric identities involving the inverse tangent function can be used to break up the input into different ranges, with the most complicated portion of the identity having a contribution below the least significant bit of the result, so that it can be ignored. The identities can be implemented in appropriate circuitry.
For inverse cosine (i.e., arccos or cos−1), the following identity may be used:




  
    arccos
    =
    
      2
      ⁢
      
        arctan
        ⁡
        
          (
          
            
              
                1
                -
                
                  x
                  2
                
              
            
            
              1
              +
              x
            
          
          )
        
      
    
  
  ,




which may be reduced to:





  arccos
  =
  
    2
    ⁢
    
      
        arctan
        ⁡
        
          (
          
            
              1
              -
              x
            
            
              
                1
                -
                
                  x
                  2
                
              
            
          
          )
        
      
      .
    
  




The inverse tangent portion may be calculated as discussed above, simplified because the input range for inverse cosine is limited to between 0 and 1. Known techniques may be used to calculate the inverse square root. For inverse sine (i.e., arcsin or sin−1), which also has an input range limited to between 0 and 1, the inverse cosine can be calculated and then subtracted from π/2, based on the identity arcsin(x)=π/2-arccos(x).

Therefore, in accordance with the present invention there is provided circuitry for computing a tangent function of an input value. That circuitry includes first look-up table circuitry that stores pre-calculated tangent values of a limited number of sample values, circuitry for inputting bits of the input value of most significance as inputs to the first look-up table circuitry to look up one of the pre-calculated tangent values as a first intermediate tangent value, circuitry for calculating a second intermediate tangent value from one or more ranges of remaining bits of the input value, and circuitry for combining the first intermediate tangent value and the second intermediate tangent value to yield the tangent function of the input value.
A corresponding method for configuring an integrated circuit device as such circuitry is also provided. Further, a machine-readable data storage medium encoded with instructions for performing the method of configuring an integrated circuit device is provided.



BRIEF DESCRIPTION OF THE DRAWINGS
Further features of the invention, its nature and various advantages will be apparent upon consideration of the following detailed description, taken in conjunction with the accompanying drawings, in which like reference characters refer to like parts throughout, and in which:
FIG. 1 shows a first embodiment of a CORDIC implementation for calculating sine and/or cosine in accordance with the present invention;
FIG. 2 shows a second embodiment of a CORDIC implementation for calculating sine and/or cosine in accordance with the present invention;
FIG. 3 shows a third embodiment of a CORDIC implementation for calculating sine and/or cosine in accordance with the present invention;
FIG. 4 shows a fourth embodiment of a CORDIC implementation for calculating sine and/or cosine in accordance with the present invention;
FIG. 5 shows an embodiment of an implementation for calculating tangent in accordance with the present invention;
FIG. 6 shows an embodiment of a portion of another implementation for calculating tangent in accordance with the present invention;
FIG. 7 shows a first range of the inverse tangent function;
FIG. 8 shows a second range of the inverse tangent function;
FIG. 9 shows a third range of the inverse tangent function;
FIG. 10 shows an embodiment of an implementation for calculating inverse tangent in accordance with the present invention;
FIG. 11 shows a first portion of an embodiment of an implementation for calculating inverse cosine and/or inverse sine in accordance with the present invention;
FIG. 12 shows a second portion of an embodiment of an implementation for calculating inverse cosine and/or inverse sine in accordance with the present invention;
FIG. 13 is a cross-sectional view of a magnetic data storage medium encoded with a set of machine-executable instructions for performing the method according to the present invention;
FIG. 14 is a cross-sectional view of an optically readable data storage medium encoded with a set of machine executable instructions for performing the method according to the present invention; and
FIG. 15 is a simplified block diagram of an illustrative system employing a programmable logic device incorporating the present invention.



DETAILED DESCRIPTION OF THE INVENTION
As discussed above, in a standard CORDIC implementation:

xn+1=xn−dnyn2−n 

yn+1=yn+dnxn2−n 

zn+1=zn−dn tan−1(2−n)

As can be seen from these equations, at the first level (for n=0):

x1=x0−y0 

y1=y0+x0 

z1=z0−tan−1(1)

Similarly, at the second level (for n=1):

x2=x1−d1(y1/2)

y2=y1+d1(x1/2)

z2=z1−d1 tan−1(0.5)

It will be understood that this continues for additional n until zn converges to 0, or as close to 0 as required by a particular implementation. However, as discussed above, the accuracy of a CORDIC implementation decreases as the input angle becomes small.

A logical structure 100 of a first embodiment according to the present invention for implementing CORDIC for sine or cosine is shown in FIG. 1. Structure 100 may be implemented as circuitry. Structure 100 is built around a CORDIC engine 101 which may be any suitable CORDIC engine, including that described in copending, commonly-assigned U.S. patent application Ser. No. 12/722,683, filed Mar. 12, 2010, which is hereby incorporated by reference herein in its entirety. CORDIC engine 101 provides an x output 111 which generally represents the cosine of the input, as well as a y output 112 which generally represents the sine of the input.
Input multiplexer 102 selects either the input variable 103, or difference 104 between π/2 and input variable 103. The control signal 105 for input multiplexer 102 result from a comparison 106 that which indicate(s) whether input variable 103 (θ) is greater than π/4 or less than π/4 (or whatever other threshold may be selected). Because input variable 103 preferably is expressed as a fraction (impliedly multiplied by π/2), the comparison may be performed simply by examining one or more of the most significant bits of input variable 103. Alternatively, a more complex comparison may be made.
In any event, the value passed to CORDIC engine is either θ or π/2−θ, depending on whether θ is greater than (or equal to) or less than π/4.
The same control signal 105 that determines the input to CORDIC engine 101 helps to determine which output 111, 112 is selected by output multiplexer 107. Specifically, cos(θ)=sin(π/2−θ) and sin(θ)=cos(π/2−θ). Additional input 115 may be provided that represents whether the desired output is sine or cosine, and that input is combined at 125 with input 105, to provide a control signal 135, so that if input variable 103 (θ) was used as the CORDIC input directly, then output 111 is selected if the desired function is cosine and output 112 is selected if the desired function is sine. But if difference 104 (π/2−θ) was used as the CORDIC input, then output 112 is selected if the desired function is cosine and output 111 is selected if the desired function is sine.
Known techniques for speeding up CORDIC calculations can be used. For example, once about halfway through the CORDIC calculation of a sine value, the final value of y can be approximated by multiplying the then-current value of x by the then-current value of z, and then subtracting that product from the then-current value of y. Similarly, for cosine, the final value of x can be approximated by multiplying the then-current value of y by the then-current value of z and adding that product to the then-current value of x. This may be referred to as “terminated CORDIC.” FIG. 2 shows a modified CORDIC structure 200 that implements terminated CORDIC in accordance with an embodiment of the invention. Structure 200 may be implemented as circuitry.
The input stage of structure 200 is identical to that of structure 100, including input multiplexer 102, input variable 103, difference 104, comparison circuit 106 and control signal 105. CORDIC engine 201 may be the same as CORDIC engine 101, except that the z datapath is used as an output 211, in addition to outputs 111, 121.
However, the output stage of structure 200 differs from the output stage of structure 100. Instead of one output multiplexer 106, there is a first output multiplexer 206 controlled by signal 135 and a second output multiplexer 216 controlled by the inverse of signal 135. For cosine, this arrangement provides x directly to adder/subtractor 208, and provides y and z to multiplier 207, which provides its output to adder/subtractor 208. For sine, this arrangement provides y directly to adder/subtractor 208, and provides x and z to multiplier 207, which provides its output to adder/subtractor 208.
In some embodiments, structure 100 or 200 may be implemented in a programmable device, such as an FPGA, either in programmable logic, or in a combination of programmable logic and fixed logic (e.g., adders and/or multipliers) if provided. For example, FPGAs in the STRATIX® family of FPGAs available from Altera Corporation, of San Jose, Calif., include digital signal processing blocks having multipliers and adders and programmable interconnect for connecting the multipliers and adders. Such an FPGA may be configured to use the adders and/or multipliers, as well as any programmable logic that may be needed, to implement structure 100 or 200.
The input range can be further limited to a smaller range between π/8 and π/4, which may provide a more accurate CORDIC result than a range between π/4 and π/2. According to such an implementation, which may be carried out using a structure 300 as shown in FIG. 3, which may be implemented as circuitry, when the input is between 0 and π/8, the input is subtracted from π/4, creating a new value between π/8 and π/4. When the input is between 3π/8 and π/2, π/4 is subtracted from the input, also creating a new value between π/8 and π/4. If the input is between π/8 and 3π/8, the input may be passed through unchanged.
As seen in FIG. 3, input multiplexer 302 selects either input 303, or the difference 304 between π/4 and input 303, or the difference 314 between input 303 and π/4. The selected input is used in CORDIC engine 301 which may be substantially identical to CORDIC engine 101.
For the case where difference 304 is used, the outputs of CORDIC engine 301 are processed by output stage 320 in accordance with the following identities:

sin(A−B)=sin(A)cos(B)−sin(B)COS(A)

cos(A−B)=cos(A)cos(B)+sin(A)sin(B)

If A=π/4, and B=π/4−θ, then A−B=π/4−(π/4−θ)=θ. Also, SIN (π/4)=COS (π/4)=2−0.5.

It follows, then, that:




  
    
      
        
          SIN
          ⁡
          
            (
            θ
            )
          
        
        =
        
        ⁢
        
          SIN
          ⁡
          
            (
            
              A
              -
              B
            
            )
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              SIN
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
            ⁢
            
              COS
              ⁡
              
                (
                
                  
                    π
                    /
                    4
                  
                  -
                  θ
                
                )
              
            
          
          -
          
            
              SIN
              ⁡
              
                (
                
                  
                    π
                    /
                    4
                  
                  -
                  θ
                
                )
              
            
            ⁢
            
              COS
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            2
            
              -
              0.5
            
          
          ⁢
          
            
              (
              
                
                  COS
                  ⁡
                  
                    (
                    
                      
                        π
                        /
                        4
                      
                      -
                      θ
                    
                    )
                  
                
                -
                
                  SIN
                  ⁡
                  
                    (
                    
                      
                        π
                        /
                        4
                      
                      -
                      θ
                    
                    )
                  
                
              
              )
            
            .
          
        
      
    
  




Similarly:





  
    
      
        
          COS
          ⁡
          
            (
            θ
            )
          
        
        =
        
        ⁢
        
          COS
          ⁡
          
            (
            
              A
              +
              B
            
            )
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              COS
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
            ⁢
            
              COS
              ⁡
              
                (
                
                  
                    π
                    /
                    4
                  
                  -
                  θ
                
                )
              
            
          
          +
          
            
              SIN
              ⁡
              
                (
                
                  
                    π
                    /
                    4
                  
                  -
                  θ
                
                )
              
            
            ⁢
            
              SIN
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            2
            
              -
              0.5
            
          
          ⁢
          
            
              (
              
                
                  COS
                  ⁡
                  
                    (
                    
                      
                        π
                        /
                        4
                      
                      -
                      θ
                    
                    )
                  
                
                +
                
                  SIN
                  ⁡
                  
                    (
                    
                      
                        π
                        /
                        4
                      
                      -
                      θ
                    
                    )
                  
                
              
              )
            
            .
          
        
      
    
  




For the case where difference 314 is used, the outputs of CORDIC engine 301 are processed by output stage 320 in accordance with the following identities:

sin(A+B)=sin(A)cos(B)+sin(B)COS(A)

cos(A+B)=cos(A)cos(B)−sin(A)sin(B)

If A=π/4, and B=θ−π/4, then A+B=π/4+(θ−π/4)=θ. Again, SIN(π/4)=COS(π/4)=2−0.5.

It follows, then, that:




  
    
      
        
          SIN
          ⁡
          
            (
            θ
            )
          
        
        =
        
        ⁢
        
          SIN
          ⁡
          
            (
            
              A
              +
              B
            
            )
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              SIN
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
            ⁢
            
              COS
              ⁡
              
                (
                
                  θ
                  -
                  
                    π
                    /
                    4
                  
                
                )
              
            
          
          +
          
            
              SIN
              ⁡
              
                (
                
                  θ
                  -
                  
                    π
                    /
                    4
                  
                
                )
              
            
            ⁢
            
              COS
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            2
            
              -
              0.5
            
          
          ⁢
          
            
              (
              
                
                  COS
                  ⁡
                  
                    (
                    
                      θ
                      -
                      
                        π
                        /
                        4
                      
                    
                    )
                  
                
                +
                
                  SIN
                  ⁡
                  
                    (
                    
                      θ
                      -
                      
                        π
                        /
                        4
                      
                    
                    )
                  
                
              
              )
            
            .
          
        
      
    
  




Similarly:





  
    
      
        
          COS
          ⁡
          
            (
            θ
            )
          
        
        =
        
        ⁢
        
          COS
          ⁡
          
            (
            
              A
              +
              B
            
            )
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              COS
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
            ⁢
            
              COS
              ⁡
              
                (
                
                  θ
                  -
                  
                    π
                    /
                    4
                  
                
                )
              
            
          
          -
          
            
              SIN
              ⁡
              
                (
                
                  θ
                  -
                  
                    π
                    /
                    4
                  
                
                )
              
            
            ⁢
            
              SIN
              ⁡
              
                (
                
                  π
                  /
                  4
                
                )
              
            
          
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            2
            
              -
              0.5
            
          
          ⁢
          
            
              (
              
                
                  COS
                  ⁡
                  
                    (
                    
                      θ
                      -
                      
                        π
                        /
                        4
                      
                    
                    )
                  
                
                -
                
                  SIN
                  ⁡
                  
                    (
                    
                      θ
                      -
                      
                        π
                        /
                        4
                      
                    
                    )
                  
                
              
              )
            
            .
          
        
      
    
  




This is implemented in output stage 320 by adder/subtractor 321 which adds the y (SIN) datapath to, or subtracts it from, the x (COS) datapath, and by multiplier 322 which multiplies that difference by 2−0.5 (indicated in FIG. 3 as SIN(π/4)). Output stage 320 also may include a multiplexing circuit 323 (similar to the combination of multiplexers 206, 216 in FIG. 2) for implementing the pass-through of the correct datapath (SIN or COS, depending on the desired function) for a case where the input 303 is between π/8 and 3π/8 and was passed through input multiplexer 302 unchanged.
Thus, as compared to embodiment 100 of FIG. 1, only one additional subtractor and one additional constant multiplication are needed for increased accuracy.
Alternatively, if the input is between π/4 and 3π/8, it also falls under embodiment 100 of FIG. 1, and it can be reflected around π/4 by subtracting the input from π/2, switching SIN and COS results to get the desired output. This also may be useful for implementing other types of algorithms to calculate SIN and COS values where a small input range can be used to improve the convergence rate.
A “terminated CORDIC” implementation similar to embodiment 200 of FIG. 2 can be used with embodiment 300 of FIG. 3. Such an implementation 400 is shown in FIG. 4.
In a case where signal 105 selects the direct input 408, then, as in embodiment 200, only one of the SIN/COS datapaths 401, 402 is multiplied at 411, 412 by the z datapath 403 and then, at 421, is added to or subtracted from the other of the SIN/COS datapaths 401, 402, and that result is multiplied at 422 by 2−0.5. Whether datapath 401 or 402 is multiplied by datapath 403 is determined by signal(s) 435, output by logic 425 which, based on comparison signal 105 and signal 115 which indicates whether sine or cosine is desired, causes one of multiplexers 413, 423 to select datapath 403 for input to a respective one of multipliers 411, 412, and other of multiplexers 413, 423 to select the value ‘1’ for input to the other respective one of multipliers 411, 412. In this case, signal(s) 435 also determines whether adder/subtractor 421 adds or subtracts, and causes multiplexer 443 to select the value ‘0’ for addition at 453 to sum/difference 421. Sum 453 is then multiplied by sin(π/4) (i.e., 2−0.5) at 422.
In a case where signal 105 selects difference input 404 (π/4−θ), because input 408 is less than π/8, then, depending on whether sine or cosine is desired, the following relationships, as discussed above, will apply:

SIN(θ)=2−0.5(COS(π/4−θ)−SIN(π/4−θ))

COS(θ)=2−0.5(COS(π/4−θ)+SIN(π/4−θ))

Similarly, in a case where signal 105 selects difference input 444 (θ−π/4), because input 408 is greater than 3π/8, then, depending on whether sine or cosine is desired, the following relationships, as discussed above, will apply:

SIN(θ)=2−0.5(COS(θ−π/4)+SIN(θ−π/4))

COS(θ)=2−0.5(COs(θ−π/4)−SIN(θ−π/4))

In a terminated CORDIC implementation,

COS(φ)=x+yz 

SIN(gyp)=y−xz 

Therefore, in a terminated CORDIC case where signal 105 selects difference input 404 (π/4−θ),




  
    
      
        
          
            SIN
            ⁡
            
              (
              θ
              )
            
          
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                +
                yz
                -
                
                  (
                  
                    y
                    -
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
    
      
        
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                -
                y
                +
                
                  (
                  
                    yz
                    +
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
  




  
    
      
        
          
            COS
            ⁡
            
              (
              θ
              )
            
          
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                +
                yz
                +
                
                  (
                  
                    y
                    -
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
    
      
        
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                +
                y
                +
                
                  (
                  
                    yz
                    -
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
  



Similarly, in a terminated CORDIC case where signal 105 selects difference input 444 (θ−π/4),




  
    
      
        
          
            SIN
            ⁡
            
              (
              θ
              )
            
          
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                +
                yz
                +
                
                  (
                  
                    y
                    -
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
    
      
        
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                +
                y
                +
                
                  (
                  
                    yz
                    -
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
  




  
    
      
        
          
            COS
            ⁡
            
              (
              θ
              )
            
          
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                +
                yz
                -
                
                  (
                  
                    y
                    -
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
    
      
        
          =
          
          ⁢
          
            
              2
              
                -
                0.5
              
            
            ⁢
            
              (
              
                x
                -
                y
                +
                
                  (
                  
                    yz
                    +
                    xz
                  
                  )
                
              
              )
            
          
        
      
    
  



For these implementations, signal(s) 435 would cause both of multiplexers 413, 423 to select datapath 403 for input to multipliers 411, 412. Signal(s) 435 also would determine whether adders/subtractors 421 and 431 add or subtract, respectively, and would causes multiplexer 443 to select the sum/difference 431 for addition at 453 to sum/difference 421. Sum 453 is then multiplied by sin(π/4) (i.e., 2−0.5) at 422.
When implementing embodiment 400 in one of the aforementioned STRATIX® FPGAs, one of the aforementioned digital signal processing blocks, having multipliers and adders, can be used at 410 to provide multipliers 411, 412 and adder/subtractor 421. The digital signal processing blocks of such FPGAs, for example, are well-suited for performing two 36-bit-by-18-bit multiplications which may be used for this purpose. One such implementation is described at page 5-21 of the Stratix III Device Handbook, Volume 1 (ver. 2, March 2010), published by Altera Corporation, which is hereby incorporated by reference herein.
Other identity-based approaches can be used to simplify the calculations of other trigonometric functions.
For example, for tan(θ), the normalized or range-reduced input is −π/2≦θ≦π/2, while the output is between negative infinity and positive infinity. The following identity holds true for the tangent function:




  
    tan
    ⁡
    
      (
      
        a
        +
        b
      
      )
    
  
  =
  
    
      
        tan
        ⁡
        
          (
          a
          )
        
      
      +
      
        tan
        ⁡
        
          (
          b
          )
        
      
    
    
      1
      -
      
        
          tan
          ⁡
          
            (
            a
            )
          
        
        ×
        
          tan
          ⁡
          
            (
            b
            )
          
        
      
    
  




By judiciously breaking up θ, one can break up the problem of calculating tan(θ) into easily calculable pieces. In this case, it may be advantageous to break up tan(θ) into three pieces. This case be done by substituting c for b, and a+b for a, above. Thus:





  
    tan
    ⁡
    
      (
      
        a
        +
        b
        +
        c
      
      )
    
  
  =
  
    
      
        tan
        ⁡
        
          (
          
            a
            +
            b
          
          )
        
      
      +
      
        tan
        ⁡
        
          (
          c
          )
        
      
    
    
      1
      -
      
        
          tan
          ⁡
          
            (
            
              a
              +
              b
            
            )
          
        
        ×
        
          tan
          ⁡
          
            (
            c
            )
          
        
      
    
  




Expanding further:





  
    tan
    ⁡
    
      (
      
        a
        +
        b
        +
        c
      
      )
    
  
  =
  
    
      
        
          
            tan
            ⁡
            
              (
              a
              )
            
          
          +
          
            tan
            ⁡
            
              (
              b
              )
            
          
        
        
          1
          -
          
            
              tan
              ⁡
              
                (
                a
                )
              
            
            ×
            
              tan
              ⁡
              
                (
                b
                )
              
            
          
        
      
      +
      
        tan
        ⁡
        
          (
          c
          )
        
      
    
    
      1
      -
      
        
          (
          
            
              
                tan
                ⁡
                
                  (
                  a
                  )
                
              
              ×
              
                tan
                ⁡
                
                  (
                  b
                  )
                
              
            
            
              1
              -
              
                
                  tan
                  ⁡
                  
                    (
                    a
                    )
                  
                
                ×
                
                  tan
                  ⁡
                  
                    (
                    b
                    )
                  
                
              
            
          
          )
        
        ×
        
          tan
          ⁡
          
            (
            c
            )
          
        
      
    
  




Although this looks much more complex than the original identity, the properties of the tangent function, and the precision of single precision arithmetic, can be used to greatly simplify the calculation.
As with many trigonometric functions, tan(θ)≈θ for small θ. The input range for the tangent function is defined as −π/2≦θ≦π/2. In single precision floating point arithmetic—e.g., under the IEEE754-1985 standard—the exponent is offset by 127 (i.e., 2° becomes 2127). If the input exponent is 115 or less (i.e., a true exponent of −12 or less), the error between tan(x) and x is below the precision of the number format, therefore below that value, tan(θ) can be considered equal to θ.
The tangent function is therefore defined for a relatively narrow exponent range, between 115 and 127, or 12 bits of dynamic range. For IEEE754-1985 arithmetic, the precision is 24 bits. The input number can therefore be represented accurately as a 36-bit fixed point number (24 bits precision +12 bits range).
Such a 36-bit fixed point number can then be split into three components. If θ=a+b+c as indicated above, the upper 9 bits can be designated the c component, the next 8 bits can be designated the a component, and 19 least significant bits may be designated the b component.
As discussed above, tan(θ)=θ for any value of θ that is smaller than 2−12. As the 19 least significant bits of a 36-bit number, b is smaller than 2−17. Therefore, tan(b)=b and we can write:




  
    tan
    ⁡
    
      (
      
        a
        +
        b
        +
        c
      
      )
    
  
  =
  
    
      
        
          
            tan
            ⁡
            
              (
              a
              )
            
          
          +
          b
        
        
          1
          -
          
            
              tan
              ⁡
              
                (
                a
                )
              
            
            ×
            b
          
        
      
      +
      
        tan
        ⁡
        
          (
          c
          )
        
      
    
    
      1
      -
      
        
          (
          
            
              
                tan
                ⁡
                
                  (
                  a
                  )
                
              
              +
              b
            
            
              1
              -
              
                
                  tan
                  ⁡
                  
                    (
                    a
                    )
                  
                
                ×
                b
              
            
          
          )
        
        ×
        
          tan
          ⁡
          
            (
            c
            )
          
        
      
    
  




The tangent of a is relatively small. The maximum value of a is slightly less than 0.0039062510 (tan(a)=0.0039062710) and the maximum value of b is 0.0000152610, which is also its tangent. Therefore, the maximum value of tan(a)×b is 5.96×10−8, therefore the minimum value 1−tan(a)×b=0.9999999410. The maximum value of tan(a)+b is 0.0039215286610. The difference between that value, and that value divided by the minimum value of 1−tan(a)×b is 2.35×10−10, or about 32 bits, which is nearly the entire width of a, b and c combined. In the worst case, where c is zero, this error would not be in the precision of the result either, which is only 24 bits.
The foregoing equation therefore can be further simplified to:




  
    tan
    ⁡
    
      (
      
        a
        +
        b
        +
        c
      
      )
    
  
  =
  
    
      
        tan
        ⁡
        
          (
          a
          )
        
      
      +
      b
      +
      
        tan
        ⁡
        
          (
          c
          )
        
      
    
    
      1
      -
      
        
          (
          
            
              tan
              ⁡
              
                (
                a
                )
              
            
            +
            b
          
          )
        
        ×
        
          tan
          ⁡
          
            (
            c
            )
          
        
      
    
  




Insofar as a and c are 8 and 9 bits respectively, the tangents for all possible bit combinations can be stored in a table with 36-bit data. Therefore the problem is reduced to a 36-bit fixed point multiplication, a 36-bit fixed point division, and a fixed point subtraction, although the additions are floating point additions as described below.

An embodiment 500 of this tangent calculation is shown in FIG. 5 and may be implemented in circuitry.
The input value 501 (θ) is first converted to a 36-bit fixed-point number by shifting at 502 by the difference 503 between its exponent and 127 (the IEEE754-1985 exponent offset value). The converted fixed-point input 504 is then split into three numbers: c—bits [36:28], a—bits [27:20], and b—bits [19:1]. Tan(c) is determined in lookup table 505 and tan(a) is determined in lookup table 506. Tan(a) and b, which are both in fixed-point format, are added at 507. That sum must then be normalized to the exponent of c, which can range from 0 to 19 (127 to 146 in single precision offset equivalent). The ‘tan(a)+b’ sum has a maximum exponent of −8 (119 in single precision offset equivalent), and is normalized at 508 for multiplying by tan(c) at 509 (for the denominator) and adding to tan(c) at 510 (for the numerator).
The numerator is normalized at 511 and now exists as a floating-point number. The local exponent (‘15’−‘c exponent’+‘a+b exponent’) is a number that is relative to ‘1.0’, and is used to denormalize the product to a fixed point number again at 512.
The denominator product is subtracted from ‘1’ at 513. The difference is normalized at 514. The difference is then inverted at 515 to form the denominator.
The denominator is multiplied by the numerator at 516. Before that multiplication, the numerator exponent is normalized at 517. The exponent is ‘119’ (which is the minimum value of c, or the maximum value of b—the reference point to which the internal exponents are normalized) plus the numerator exponent plus the denominator exponent. The denominator exponent is the shift value from the final denominator normalization 514—normally this would be considered a negative relative exponent, and subtracted from any final exponent. However, because the denominator is arithmetically inverted immediately following the normalization, the exponent is converted from negative to positive, and can therefore be added at this stage.
The result is rounded at 518 and is ready for use. However, if the exponent of the original input value 501 is less than 115, the output and the input are considered the same—i.e., tan(θ)=θ. This is implemented with the multiplexer 519, which selects as the final output either the rounded calculation result 520, or input 501, based on control signal 521 which is determined (not shown) by the size of the exponent of input 501.
A double-precision tangent calculation may also be implemented.
As noted above, the following identity holds true for the tangent function:




  
    tan
    ⁡
    
      (
      
        a
        +
        b
      
      )
    
  
  =
  
    
      
        tan
        ⁡
        
          (
          a
          )
        
      
      +
      
        tan
        ⁡
        
          (
          b
          )
        
      
    
    
      1
      -
      
        
          tan
          ⁡
          
            (
            a
            )
          
        
        ×
        
          tan
          ⁡
          
            (
            b
            )
          
        
      
    
  




By judiciously breaking up θ, again as above one can break up the problem of calculating tan(θ) into easily calculable pieces. Once again, it may be advantageous to break up tan(θ) into three pieces by substituting b+c for b, above. Thus:





  
    tan
    ⁡
    
      (
      
        a
        +
        b
        +
        c
      
      )
    
  
  =
  
    
      
        tan
        ⁡
        
          (
          a
          )
        
      
      +
      
        tan
        ⁡
        
          (
          
            b
            +
            c
          
          )
        
      
    
    
      1
      -
      
        
          tan
          ⁡
          
            (
            a
            )
          
        
        ×
        
          tan
          ⁡
          
            (
            
              b
              +
              c
            
            )
          
        
      
    
  




One can define a to be a 9-bit value, for the range of 1.11111111 down to 0.0000000—i.e., the implied leading ‘1’ for the mantissa when the number of radians is 1 or greater, and the most significant 8 bits of the mantissa to the right of the decimal point. b can be defined as the next 9 bits down, while c value will be the least significant bits of the mantissa. For smaller arguments, where the exponent is less than 1023 (i.e., for double-precision values less than 1), c can still have 53 significant bits (one implied bit followed by 52 mantissa bits), but the invention applies regardless of the value of the exponent.
One can start with the less significant bits, b and c:




  
    tan
    ⁡
    
      (
      
        b
        +
        c
      
      )
    
  
  =
  
    
      
        tan
        ⁡
        
          (
          b
          )
        
      
      +
      
        tan
        ⁡
        
          (
          c
          )
        
      
    
    
      1
      -
      
        
          tan
          ⁡
          
            (
            b
            )
          
        
        ×
        
          tan
          ⁡
          
            (
            c
            )
          
        
      
    
  




For c as defined above, the maximum value of c will be just below 2−17, and therefore, because for small α, tan(α)≈α, the maximum value of tan(c) will be around 2−17 as well. If a series expansion of tan(α) is:





  
    
      tan
      ⁡
      
        (
        α
        )
      
    
    =
    
      α
      +
      
        
          α
          3
        
        3
      
      +
      
        
          2
          ⁢
          
            α
            5
          
        
        15
      
      +
      ...
    
  
  ⁢
  
    
  
  ,




then when α=c, tan(c) is at most around 2−17, then the second (c3) term will be about 35 binary positions to the right of the input argument, and the third (c5) term will be almost 70 binary positions to the right of the input argument. The second term will therefore contribute up to 52−35=17 bits to the mantissa. The third term may contribute a rounding input very occasionally, and can be ignored for double precision calculations.

The second (c3) term can be reasonably accurately be calculated with three 18-bit fixed-point multipliers, in a tree format. One multiplier can be used to multiply c by itself to yield c2 (which is then right shifted by 18 bits). A second multiplier can be used to multiply c by 174763 to scale it by 1⁄3 to account for the use of fixed-point (and that result is then right shifted by 18 bits). Those two results are multiplied by the third multiplier and that result, right-shifted by 35 bits, is then added to c. At each step, the upper 18 bits of c are used for the multiplier inputs. For higher accuracy, 27-bit multipliers can also be used, operating on the upper 27 bits of c. This can be implemented using three 18×18 multipliers and one 55-bit adder, which, in a programmable device, may be implemented in programmable logic.
As described above, b is a 9-bit value and therefore, if non-zero, is in the range of around 2−17 to around 2−8, which means that tan(b) is about the same. The 9-bit input range corresponds to 512 values, which can easily be pre-computed in stored in a look-up table. 512 64-bit values would occupy 32 Kb, which can easily be accommodated, e.g., in two M20K memories of a programmable logic device sold by Altera Corporation under the name STRATIX® V.
To calculate tan(b+c), the look-up table result for tan(b) is added to the calculation of tan(c). Because tan(b) has only a narrow range, this addition is smaller than a full floating-point addition, so can be efficiently implemented as a wider fixed-point addition, followed by a normalization.
The next step is to divide by (1−tan(b)tan(c)). This can be done using the binomial theorem, according to which:




  
    N
    
      1
      -
      x
    
  
  =
  
    
      
        N
        ⁡
        
          (
          
            1
            +
            x
          
          )
        
      
      ⁢
      
        (
        
          1
          +
          
            x
            2
          
        
        )
      
    
    
      1
      -
      
        x
        4
      
    
  




Here, we can set N=tan(b)+tan(c) and x=tan(b)tan(c). The maximum value of x would be about 2−25, so x4 would be around 2−100 which can be ignored in this case. If b or c were smaller, then x would be proportionally smaller as well.

The (1−x) term can be constructed using only fixed-point arithmetic, multiplying the 27-bit fixed-point values of b and c. The multiplication of the numerator N by (1+x)—i.e., (1+tan(b)tan(c))—can be implemented in a 54×54 multiplier. The (1+x2) term does not have to be multiplied out, but can be accurately estimated.
Specifically, because the value of x2 is very small (about 50 bits to the right of the ‘1’), it is easier to add a correction factor to the least significant portion of the mantissa of N(1+x), rather than multiplying N(1+x) by (1+x2). First, the value of x2 can be estimated by squaring some of the most significant bits (MSBs) of the tan(b)tan(c) product (i.e., some of the MSBs of x). Only a few bits (e.g., 3-5 bits) are needed because the contribution to the (1+x2) total is so small. Therefore, this multiplication could be implemented in either hard multipliers or programmable logic (if implemented in a programmable logic device).
Remembering that N=tan(b)+tan(c) and that the correction factor is meant to approximate the contribution of (1+x2) to N(1+x)(1+x2), in view of the magnitudes of b and c, only the most significant bits of the tan(b) term would have an impact on the result, so the correction factor can be accurately estimated by multiplying a few MSBs (again, e.g., 3-5 bits) of tan(b) by a few MSBs of the square of tan(b)tan(c). This result is then added to the least significant bits of the output of the upper 54 bits of the 54×54 multiplier (i.e., to N(1+x)), resulting in an accurate approximation of tan(b+c). This is achieved using one double-precision normalization, one 27×27 multiplier and one 54×54 multiplier.
The final identity, combining tan(a)—which, like tan(b), can be derived from a look-up table—with tan(b+c), can be calculated as a floating-point operation, using a floating-point multiplier, two floating-point adders, and a floating-point divider. The range of tan(b+c) is small compared to that of tan(a). tan(b+c) has a range of 0 to 2−8, while tan(a) has a range of 2−8 to about 28, ensuring that the product tan(a)tan(b+c) will never exceed 1.
An embodiment 550 of a portion of this double-precision tangent calculation, and particularly the tan(b)+tan(c) portion, is shown in FIG. 6 and may be implemented in circuitry.
As noted above, values for tan(b) may be pre-calculated and stored in a look-up table 565, while tan(c) may be calculated as a series expansion as discussed above. These may be 18 bits or 27 bits as discussed above. tan(b) is input from look-up table 565 at 551 and tan(c) (the least significant bits of θ) is input at 552, and the two values are added at 553 to provide the N term (above) and multiplied at 554, which may be a 54×54 multiplier, to provide the x term (above).
The x term is right-shifted at 555 by an amount equal to the different between the exponent of tan(b), which is output from the look-up table, and the exponent of tan(c), which is determined in the calculation described above. The x term is then explicitly added to ‘1’ or inverted to produce a one's-complement value (neither option shown) to generate the (1+x) term (above), which is multiplied by the N term at 556. Because x is much smaller than ‘1’, the addition of ‘1’ to x is a simple bit insertion to the left of the decimal point.
The aforementioned correction factor to replace the (1+x2) term may be derived as follows:

N(1+x)(1+x2)=(N+Nx)(1+x2)=N+Nx+Nx2+Nx3 

Nx3 will be smaller than the least significant bit of the result and so can be ignored, leaving N+Nx+Nx2=N(1+x)+Nx2. Nx2 may be approximated by multiplying the upper bits of x (e.g., the four uppermost bits in FIG. 5) by themselves at 557, and then multiplying that product by the upper bits (e.g., the eight uppermost bits in FIG. 5) of tan(b) at 558. Again, only the upper bits are kept (e.g., the eight uppermost bits in FIG. 5) and right-shifted at 559. The amount of the shift is derived from the exponents as for the 1+x term, but is doubled to account for the squaring. This result is then added at 560 to the least significant bits of the output of the upper bits (e.g., the 54 uppermost bits in FIG. 5) of multiplier 554 (i.e., to N(1+x)) to provide an accurate approximation of tan(b+c).

The remainder of the double-precision tangent calculation, combining tan(a) with tan(b+c), is not shown, but as discussed above can be calculated as a floating-point operation.
It should be noted that for either the single-precision tangent calculation described above or the double-precision tangent calculation described above, any one or more of the ranges a, b and c can be further broken down into subranges and a tangent value can be determined for each subrange by either a look-up table or by a calculation such as those described above.
For inverse tangent (i.e., arctan(x) or tan−1(x)), the situation is reversed in that the input is between negative infinity and positive infinity, while the normalized or range-reduced output is −π/2≦θ≦+π/2. Once again, the problem can be broken down into input ranges. Thus, as shown in FIG. 7, for −1≦x≦+1, arctan(x) is relatively linear, and has an output in the range −π/4≦θ≦+π/4. Viewed on an intermediate scale in FIG. 8, for −1010≦x≦+1010, arctan(x) shows inflection points past ±1, flattening out to an output in the range −π/2≦θ≦+π/2. Viewed on a larger scale in FIG. 9, for −10010≦x≦+10010, arctan(x) remains in an output range of −π/2≦θ≦+π/2. In IEEE754-1985 single-precision arithmetic, with 24 bits of precision and an exponent up to 127, the output must be exact for a set of 833860810 (223) points in 127 segments along the input curve.
The calculation can be simplified by separately handling inputs of magnitude less than 1 and inputs of magnitude greater than 1.
Considering first inputs of magnitude less than 1, the following identity may be applied:




  
    arctan
    ⁡
    
      (
      a
      )
    
  
  =
  
    
      arctan
      ⁡
      
        (
        b
        )
      
    
    +
    
      arctan
      ⁡
      
        (
        
          
            a
            -
            b
          
          
            1
            +
            ab
          
        
        )
      
    
  




If b is close enough to a, then c=(a−b)/(1+ab) is very small and arctan(c)≈c. If a has a maximum value of 1, and b is made equal to the 8 or 10 most significant bits of a, then c will have a maximum value of 1/256 for 8 bits, or 1/1024 for 10 bits. The inverse tangent 1/256 (0.0039062510) is 0.0039062310. The error is almost at the floor of the precision of the input (23 bits). The inverse tangent of 1/1024 (0.000976562510) is 0.00097652210, which has an error below the least significant bit of the input range.

The subrange b can easily be separated from a by truncation. The (a−b) term is made up of the truncated bits. ab will always be less than 1, so the term 1+ab can be calculated without an adder, by directly concatenating a ‘1’ to the ab term. c can therefore be calculated easily. Values for arctan(b) where b<1 can be stored in a lookup table. Therefore, arctan(a) can be determined by looking up arctan(b) in the table and adding c to the lookup result.
If a is greater than (or equal to) 1, then the following identity can be used:




  
    arctan
    ⁡
    
      (
      a
      )
    
  
  =
  
    
      π
      2
    
    -
    
      arctan
      ⁡
      
        (
        
          1
          a
        
        )
      
    
  




The inverse tangent of the inverse of the input may be determined as described above for inputs less than 1 and then the desired result is obtained by subtracting from π/2.

The input mantissa is converted to a 36-bit number, by appending a ‘1’ to the left of the most significant bit, and appending a number of ‘0’ bits to the right of the least significant bit. The exponent of this number is 127−(input_exponent−127)−1=253−input exponent. Once the correct floating point input has been selected, the number is converted to a fixed point equivalent by right shifting it from the reference point of 1. This should not reduce accuracy of the result given that 36-bit numbers are used. If the number has an exponent of 117 (for a right shift of 10 bits), there will still be 26 valid bits in the 36-bit magnitude, and only 24 bits are needed for single-precision floating point representation.
An embodiment 900 of this inverse tangent calculation is shown in FIG. 10 and may be implemented in circuitry.
The mantissa of the input number (with leading ‘1’ and trailing zero(es) appended as discussed above) is input at 901, while the exponent of the input number is input at 902. Multiplexers 903, 904 select the unaltered input mantissa and exponent if the input value is less than 1.0 as determined at 905. Otherwise, multiplexers 903, 904 select inverse 906 and the new exponent (original exponent subtracted from ‘253’). Inverse 906 can be computed using any suitable inverse calculation module 907. The input is then normalized at 908 to a fixed-point representation 909.
The uppermost bits 910 are input to lookup table 911 and are also input to multiplier 912 along with all bits 913—this is the ab calculation discussed above. a-b is the remaining bits 914.
The ab product 915 is added to ‘1’ at 916 to form the 1+ab sum 917 which is inverted, again using any suitable inversion module 918 to form 1/(1+ab) quotient 919, which is multiplied at 920 by a-b term 914, forming the c term (a-b)/(1+ab). c term 921 is added at 922 to arctan(b) as determined in lookup table 911.
If the original input was less than ‘1’, then sum 922 is the result, which is selected by control signal 905 at multiplexer 923 following normalization at 924. Sum 922 is also subtracted from π/2 at 925 and normalized at 926, and if the original input was greater than (or equal to) ‘1’, then difference 925, as normalized at 926, is the result, which is selected by control signal 905 at multiplexer 923. Any necessary rounding, exception handling, etc., is performed at 927 to provide result 928.
In the input range of exponents 115-120 (as an example), some inaccuracies in the output (still limited to a small number of least significant bits) may occur. One way to solve this is to use a second, smaller lookup table (not shown) for a limited subset of most significant valid bits—e.g., 6 bits. The b value would be the upper 6 bits of the subrange, and the a−b value would be the lower 20 bits. To maintain the maximum amount of precision, the table could contain results that are normalized to the subrange—for example if the largest exponent in the subrange were 120, then 1.999910×2120 would be a fully normalized number, with all other table entries relative to that one. The c value would have to be left-shifted so that it would have the correct magnitude in relation to the b table output. One way (and possibly the most accurate way) to implement this would be to take the a−b value from before the fixed point shifter. That is, instead of fixed-point representation 909, the output of multiplexer 903 could be used directly. An additional multiplexer (not shown) could be provided to select between the output of multiplexer 903 for smaller exponents and the output of normalizer 908 for larger exponents.
Once inverse tangent can be calculated, inverse cosine and inverse sine can easily be calculated based on:



  
    
      arccos
      ⁡
      
        (
        x
        )
      
    
    =
    
      2
      ⁢
      
        arctan
        ⁡
        
          (
          
            
              1
              -
              x
            
            
              
                1
                -
                
                  x
                  2
                
              
            
          
          )
        
      
    
  




  and




  
    
      arcsin
      ⁡
      
        (
        x
        )
      
    
    =
    
      
        π
        /
        2
      
      -
      
        
          arccos
          ⁡
          
            (
            x
            )
          
        
        .
      
    
  



An embodiment 1000 of the inverse cosine calculation is shown in FIG. 11 and may be implemented in circuitry.
The input argument x is input as sign 1001, mantissa 1002 and exponent 1003 to preprocessing module 1100, shown in more detail in FIG. 12. Preprocessing module 1100 prepares the 1−x numerator term 1004 and the 1−x2 denominator term 1005 from x. The inverse square root of denominator term 1005 is taken at inverse square root module 1006, which be any suitable inverse square root module.
The mantissa 1007 of the inverse square root is multiplied at 1008 by the mantissa 1009 of the numerator term, while the exponent 1010 of the inverse square root is added at 1011 to the exponent 1012 of the numerator term, and ‘127’ is subtracted from the exponent at 1013. The result is input to an inverse tangent module 1014 which may be inverse tangent module 900, above.
The inverse tangent module 1014 outputs a 36-bit fixed point value between 0 and π/2. If the input number is positive, the inverse cosine of that input must lie in the first (or fourth) quadrant, and the output of inverse tangent module 1014 is used directly. This is implemented by exclusive-OR gate 1024 and AND-gate 1034. Sign bit 1001 will be a ‘0’, meaning that XOR-gate 1024 will pass the output of inverse tangent module 1014 without change, and there will be no contribution at adder 1044 from AND-gate 1034. If the input number is negative, the inverse cosine of that input must lie in the second (or third) quadrant. In that case, the inverse cosine value can be calculated by subtracting the output of inverse tangent module 1014 from π. Sign bit 1001 will be a ‘1’, meaning that XOR-gate 1024 will pass the 1's-complement negative of the output of inverse tangent module 1014, while AND-gate 1034 will pass the value π. The sign bit is also used as a carry input (not shown) to adder 1044, converting the 1's-complement number to 2's-complement format, and adder 1044 outputs the difference between π and the output of inverse tangent module 1014.
Output 1015 is the inverse cosine. By subtracting inverse cosine output 1015 from π/2 at 1016, inverse sine 1017 can be determined. However, for inputs having real exponents less than −12 (IEEE754-1985 exponents less than 115), it would be more accurate to rely on arcsin(x)≈x than to rely on the calculated value.
FIG. 12 shows preprocessing module 1100, including numerator portion 1101 and denominator portion 1102. The input mantissa (with leading ‘1’ and trailing zero(es)) is input to numerator portion 1101 at 1103 and to denominator portion 1102 at 1104, while the input exponent is input to numerator portion 1101 at 1105 and to denominator portion 1102 at 1106.
On the numerator side, ‘127’ is subtracted from the exponent at 1107 to determine the “real” exponent, which is then used in shifter 1108 to turn the input mantissa into a fixed-point number 1109, which is subtracted at 1110 from ‘1’ to yield the numerator 1−x. The number of leading zeroes in the result are counted at count-leading-zeroes module 1111 and used at shifter 1112 to normalize the numerator mantissa and at subtractor 1113 to determine the IEEE754-1985 exponent by subtracting from ‘127’.
On the denominator side, the input mantissa is multiplied by itself at multiplier 1114 to determine x2 value 1115. The input exponent is left-shifted by one place at 1116 and subtracted from ‘253’ at 1117 to determine how far to right-shift value 1115 at 1118 to yield a fixed-point representation of x2. The fixed-point representation of x2 is subtracted from ‘1’ at 1119 to yield denominator 1−x2 at 1120. The number of leading zeroes in value 1120 are counted at count-leading-zeroes module 1121 and used at shifter 1122 to normalize the denominator mantissa and at subtractor 1123 to determine the IEEE754-1985 exponent by subtracting from ‘127’.
The trigonometric function calculating structures described above can be implemented as dedicated circuitry or can be programmed into programmable integrated circuit devices such as FPGAs. As discussed, in FPGA implementations, certain portions of the circuitry, particularly involving multiplications and combinations of multiplications as indicated, can be carried out in specialized processing blocks of the FPGA, such as a DSP block, if provided in the FPGA.
Instructions for carrying out a method according to this invention for programming a programmable device to implement the trigonometric function calculating structures described above may be encoded on a machine-readable medium, to be executed by a suitable computer or similar device to implement the method of the invention for programming or configuring PLDs or other programmable devices to perform addition and subtraction operations as described above. For example, a personal computer may be equipped with an interface to which a PLD can be connected, and the personal computer can be used by a user to program the PLD using a suitable software tool, such as the QUARTUS® II software available from Altera Corporation, of San Jose, Calif.
FIG. 13 presents a cross section of a magnetic data storage medium 1200 which can be encoded with a machine executable program that can be carried out by systems such as the aforementioned personal computer, or other computer or similar device. Medium 1200 can be a floppy diskette or hard disk, or magnetic tape, having a suitable substrate 1201, which may be conventional, and a suitable coating 1202, which may be conventional, on one or both sides, containing magnetic domains (not visible) whose polarity or orientation can be altered magnetically. Except in the case where it is magnetic tape, medium 1200 may also have an opening (not shown) for receiving the spindle of a disk drive or other data storage device.
The magnetic domains of coating 1202 of medium 1200 are polarized or oriented so as to encode, in manner which may be conventional, a machine-executable program, for execution by a programming system such as a personal computer or other computer or similar system, having a socket or peripheral attachment into which the PLD to be programmed may be inserted, to configure appropriate portions of the PLD, including its specialized processing blocks, if any, in accordance with the invention.
FIG. 14 shows a cross section of an optically-readable data storage medium 1210 which also can be encoded with such a machine-executable program, which can be carried out by systems such as the aforementioned personal computer, or other computer or similar device. Medium 1210 can be a conventional compact disk read-only memory (CD-ROM) or digital video disk read-only memory (DVD-ROM) or a rewriteable medium such as a CD-R, CD-RW, DVD-R, DVD-RW, DVD+R, DVD+RW, or DVD-RAM or a magneto-optical disk which is optically readable and magneto-optically rewriteable. Medium 1210 preferably has a suitable substrate 1211, which may be conventional, and a suitable coating 1212, which may be conventional, usually on one or both sides of substrate 1211.
In the case of a CD-based or DVD-based medium, as is well known, coating 1212 is reflective and is impressed with a plurality of pits 1213, arranged on one or more layers, to encode the machine-executable program. The arrangement of pits is read by reflecting laser light off the surface of coating 1212. A protective coating 1214, which preferably is substantially transparent, is provided on top of coating 1212.
In the case of magneto-optical disk, as is well known, coating 1212 has no pits 1213, but has a plurality of magnetic domains whose polarity or orientation can be changed magnetically when heated above a certain temperature, as by a laser (not shown). The orientation of the domains can be read by measuring the polarization of laser light reflected from coating 1212. The arrangement of the domains encodes the program as described above.
A PLD 140 programmed according to the present invention may be used in many kinds of electronic devices. One possible use is in a data processing system 1400 shown in FIG. 15. Data processing system 1400 may include one or more of the following components: a processor 1401; memory 1402; I/O circuitry 1403; and peripheral devices 1404. These components are coupled together by a system bus 1405 and are populated on a circuit board 1406 which is contained in an end-user system 1407.
System 1400 can be used in a wide variety of applications, such as computer networking, data networking, instrumentation, video processing, digital signal processing, or any other application where the advantage of using programmable or reprogrammable logic is desirable. PLD 140 can be used to perform a variety of different logic functions. For example, PLD 140 can be configured as a processor or controller that works in cooperation with processor 1401. PLD 140 may also be used as an arbiter for arbitrating access to a shared resources in system 1400. In yet another example, PLD 140 can be configured as an interface between processor 1401 and one of the other components in system 1400. It should be noted that system 1400 is only exemplary, and that the true scope and spirit of the invention should be indicated by the following claims.
Various technologies can be used to implement PLDs 140 as described above and incorporating this invention.
It will be understood that the foregoing is only illustrative of the principles of the invention, and that various modifications can be made by those skilled in the art without departing from the scope and spirit of the invention. For example, the various elements of this invention can be provided on a PLD in any desired number and/or arrangement. One skilled in the art will appreciate that the present invention can be practiced by other than the described embodiments, which are presented for purposes of illustration and not of limitation, and the present invention is limited only by the claims that follow.

## Claims

1. Circuitry for computing a tangent function of an input value, said circuitry comprising:
first look-up table circuitry that stores pre-calculated tangent values of a limited number of sample values;
circuitry for inputting bits of said input value of most significance as inputs to said first look-up table circuitry to look up one of said pre-calculated tangent values as a first intermediate tangent value;
circuitry for calculating a second intermediate tangent value from one or more ranges of remaining bits of said input value; and
circuitry for combining said first intermediate tangent value and said second intermediate tangent value to yield said tangent function of said input value.

2. The circuitry of claim 1 wherein:
said circuitry for calculating a second intermediate tangent value from one or more ranges comprises, for each respective one of said one or more ranges:
respective additional look-up table circuitry for outputting a respective third intermediate tangent value from a respective first subrange of said respective one of said one or more ranges;
respective additional circuitry for calculating a respective fourth intermediate tangent value from a respective second subrange of said respective one of said one or more ranges; and
respective additional circuitry for combining said respective third intermediate tangent value and said respective fourth intermediate tangent value to yield said second intermediate tangent value.

3. The circuitry of claim 2 further comprising, for each of one or more of said subranges, respective further look-up table circuitry, respective further circuitry for calculating, and respective further circuitry for combining.

4. The circuitry of claim 1 wherein circuitry for calculating one of said intermediate tangent values for one of said ranges of remaining bits of least significance approximates said one of said intermediate tangent values by setting said one of said intermediate tangent values equal to said one of said ranges of remaining bits of least significance.

5. The circuitry of claim 1 wherein circuitry for calculating one of said intermediate tangent values for one of said ranges of remaining bits of least significance calculates said one of said intermediate tangent values as a power series.

6. The circuitry of claim 1 wherein said circuitry for combining includes circuitry for dividing.

7. The circuitry of claim 6 wherein said circuitry for dividing performs division by approximation.

8. The circuitry of claim 7 wherein said circuitry for dividing performs division as a power series.

9. The circuitry of claim 8 wherein said circuitry for dividing calculates higher-order terms of said power series as an additive approximation.

10. The circuitry of claim 1 wherein:
said circuitry for calculating a second intermediate value comprises circuitry for computing a series expansion of a tangent of a sum of bits of said input value of intermediate significance and bits of said input value of least significance.

11. The circuitry of claim 10 wherein said circuitry for computing a series expansion comprises:
circuitry for adding said first intermediate tangent value and said bits of said input value of least significance to provide a first interim value;
circuitry for multiplying said first intermediate tangent value and said bits of said input value of least significance to provide a second interim value;
circuitry for adding said second interim value to ‘1’ to provide a third interim value;
circuitry for multiplying said first interim value by said third interim value to provide a fourth interim value;
circuitry for approximating a sum of ‘1’ and a square of said second interim value; and
circuitry for adding said approximated sum to said fourth interim value.

12. The circuitry of claim 11 wherein said circuitry for approximating a sum of ‘1’ and a square of said second interim value comprises:
circuitry for multiplying upper bits of said second interim value by themselves to provide an approximation of said square of said second interim value; and
circuitry for multiplying said approximation of said square of said second interim value by said first intermediate tangent value.

13. The circuitry of claim 12 wherein:
said upper bits of said second interim value comprise between three bits and five bits of said second interim value.

14. The circuitry of claim 12 wherein:
circuitry for adding said approximated sum to said fourth interim value adds, to said fourth interim value, upper bits of output of said circuitry for multiplying said approximation of said square of said second interim value by said first intermediate tangent value.

15. A method of configuring a programmable integrated circuit device as circuitry for computing a tangent function of an input value, said method comprising:
configuring first look-up table circuitry that stores pre-calculated tangent values of a limited number of sample values;
configuring circuitry for inputting bits of said input value of most significance as inputs to said first look-up table circuitry to look up one of said pre-calculated tangent values as a first intermediate tangent value;
configuring circuitry for calculating a second intermediate tangent value from one or more ranges of remaining bits of said input value; and
configuring circuitry for combining said first intermediate tangent value and said second intermediate tangent value to yield said tangent function of said input value.

16. The method of claim 15 wherein:
said configuring circuitry for calculating a second intermediate tangent value from one or more ranges comprises configuring, for each respective one of said one or more ranges:
respective additional look-up table circuitry for outputting a respective third intermediate tangent value from a respective first subrange of said respective one of said one or more ranges;
respective additional circuitry for calculating a respective fourth intermediate tangent value from a respective second subrange of said respective one of said one or more ranges; and
respective additional circuitry for combining said respective third intermediate tangent value and said respective fourth intermediate tangent value to yield said second intermediate tangent value.

17. The method of claim 16 further comprising configuring, for each of one or more of said subranges, respective further look-up table circuitry, respective further circuitry for calculating, and respective further circuitry for combining.

18. The method of claim 15 wherein said configuring circuitry for calculating one of said intermediate tangent values for one of said ranges of remaining bits of least significance comprises configuring circuitry that approximates said one of said intermediate tangent values by setting said one of said intermediate tangent values equal to said one of said ranges of remaining bits of least significance.

19. The method of claim 15 wherein said configuring circuitry for calculating one of said intermediate tangent values for one of said ranges of remaining bits of least significance comprises configuring circuitry that calculates said one of said intermediate tangent values as a power series.

20. The method of claim 15 wherein said configuring circuitry for combining includes configuring circuitry for dividing.

21. The method of claim 20 wherein said configuring circuitry for dividing comprises configuring circuitry that performs division by approximation.

22. The method of claim 21 wherein said configuring circuitry for dividing comprises configuring circuitry that performs division as a power series.

23. The method of claim 22 wherein said configuring circuitry for dividing comprises configuring circuitry that calculates higher-order terms of said power series as an additive approximation.

24. The method of claim 15 wherein:
said configuring circuitry for calculating a second intermediate value comprises configuring circuitry for computing a series expansion of a tangent of a sum of bits of said input value of intermediate significance and bits of said input value of least significance.

25. The method of claim 24 wherein said configuring circuitry for computing a series expansion comprises:
configuring circuitry for adding said first intermediate tangent value and said bits of said input value of least significance to provide a first interim value;
configuring circuitry for multiplying said first intermediate tangent value and said bits of said input value of least significance to provide a second interim value;
configuring circuitry for adding said second interim value to ‘1’ to provide a third interim value;
configuring circuitry for multiplying said first interim value by said third interim value to provide a fourth interim value;
configuring circuitry for approximating a sum of ‘1’ and a square of said second interim value; and
configuring circuitry for adding said approximated sum to said fourth interim value.

26. The method of claim 25 wherein said configuring circuitry for approximating a sum of ‘1’ and a square of said second interim value comprises:
configuring circuitry for multiplying upper bits of said second interim value by themselves to provide an approximation of said square of said second interim value; and
configuring circuitry for multiplying said approximation of said square of said second interim value by said first intermediate tangent value.

27. The method of claim 26 wherein:
said upper bits of said second interim value comprise between three bits and five bits of said second interim value.

28. The method of claim 26 wherein:
said configuring circuitry for adding said approximated sum to said fourth interim value comprises configuring circuitry that adds, to said fourth interim value, upper bits of output of said circuitry for multiplying said approximation of said square of said second interim value by said first intermediate tangent value.

29. A non-transitory machine-readable data storage medium encoded with non-transitory machine-executable instructions for configuring a programmable integrated circuit device as circuitry for computing a tangent function of an input value, said instructions comprising:
instructions to configure first look-up table circuitry that stores pre-calculated tangent values of a limited number of sample values;
instructions to configure circuitry for inputting bits of said input value of most significance as inputs to said first look-up table circuitry to look up one of said pre-calculated tangent values as a first intermediate tangent value;
instructions to configure circuitry for calculating a second intermediate tangent value from one or more ranges of remaining bits of said input value; and
instructions to configure circuitry for combining said first intermediate tangent value and said second intermediate tangent value to yield said tangent function of said input value.