Abstract:
Circuitry for computing a trigonometric function of an input includes circuitry for relating the input to another value to generate an intermediate value, circuitry for selecting one of the input and the intermediate value as a trigonometric input value, circuitry for determining respective initial values of a plurality of trigonometric functions for the trigonometric input value, and circuitry for deriving, based at least in part on a trigonometric identity, a final value of the first trigonometric function from the respective initial values of the plurality of trigonometric functions. The trigonometric function may be any of sine, cosine and tangent and their inverse functions. The trigonometric identities used allow a computation of a trigonometric function to be broken into pieces that either are easier to perform or can be performed more accurately.

Description:
BACKGROUND 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). 
     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:
 
 x   n+1   =x   n   −d   n   y   n 2 −n  
 
 y   n+1   =y   n   +d   n   x   n 2 −n  
 
 z   n+1   =z   n   −d   n  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 d n , which is ±1—if z n  is positive, then so is d n , as the goal is to bring z closer to 0; if z n  is negative, then so is d n , 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 θ 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. 
     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 θ&lt;π/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. 
     For the inverse tangent (i.e., arc tan 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., arc cos 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 arc sin(x)=π/2−arc cos(x).
 
     Therefore, in accordance with the present invention there is provided circuitry for computing a trigonometric function of an input. The circuitry includes circuitry for relating the input to another value to generate an intermediate value, circuitry for selecting one of the input and the intermediate value as a trigonometric input value, circuitry for determining respective initial values of a plurality of trigonometric functions for the trigonometric input value, and circuitry for deriving, based at least in part on a trigonometric identity, a final value of the first trigonometric function from the respective initial values of the plurality of trigonometric functions. 
     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 a first range of the inverse tangent function; 
         FIG. 7  shows a second range of the inverse tangent function; 
         FIG. 8  shows a third range of the inverse tangent function; 
         FIG. 9  shows an embodiment of an implementation for calculating inverse tangent in accordance with the present invention; 
         FIG. 10  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. 11  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. 12  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. 13  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. 14  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:
 
 x   n+1   =x   n   −d   n   y   n 2 −n  
 
 y   n+1   =y   n   +d   n   x   n 2 −n  
 
 z   n+1   =z   n   −d   n  tan −1 (2 −n )
 
As can be seen from these equations, at the first level (for n=0):
 
 x   1   =x   0   −y   0  
 
 y   1   =y   0   +x   0  
 
 z   1   =z   0 −tan −1 (1)
 
Similarly, at the second level (for n=1):
 
 x   2   =x   1   −d   1 ( y   1 /2)
 
 y   2   =y   1   +d   1 ( x   1 /2)
 
 z   2   =z   1   −d   1  tan −1 (0.5)
 
It will be understood that this continues for additional n until z n  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 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 ,  121  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(•)= 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 
                         - 
                         yz 
                         + 
                         
                           ( 
                           
                             y 
                             + 
                             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  404 ,  405  and adder/subtractor  406 . 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 0  becomes 2 127 ). 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.00390625 10  (tan(a)=0.00390627 10 ) 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.99999994 10 . The maximum value of tan(a)+b is 0.00392152866 10 . 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 . 
     For inverse tangent (i.e., arc tan(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. 6 , for −1≦x≦+1, arc tan(x) is relatively linear, and has an output in the range −π/4≦θ≦+π/4. Viewed on an intermediate scale in  FIG. 7 , for −10 10 ≦x≦+10 10 , arc tan(x) shows inflection points past ±1, flattening out to an output in the range −π/2≦θ≦+π/2. Viewed on a larger scale in  FIG. 8 , for −100 10 ≦x≦+100 10 , 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 (2 23 ) 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 arc tan(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.00390625 10 ) is 0.00390623 10 . The error is almost at the floor of the precision of the input (23 bits). The inverse tangent of 1/1024 (0.0009765625 10 ) is 0.000976522 10 , 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 arc tan(b) where b&lt;1 can be stored in a lookup table. Therefore, arc tan(a) can be determined by looking up arc tan(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. 9  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 arc tan(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.9999 10 ×2 120  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
 
arc sin( x )=π/2−arc cos( x ).
 
     An embodiment  1000  of the inverse cosine calculation is shown in  FIG. 10  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. 11 . Preprocessing module  1100  prepares the 1−x numerator term  1004  and the 1−x 2  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&#39;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&#39;s-complement number to 2&#39;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 arc sin(x)≈x than to rely on the calculated value. 
       FIG. 11  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 x 2  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 x 2 . The fixed-point representation of x 2  is subtracted from ‘1’ at  1119  to yield denominator 1−x 2  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. 12  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. 13  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. 13 . 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.