Is it possible to get a fast fix point Arctan routine written in C using the principle with a lookup table or perhaps another principle? Regards, Ren�
Arctan approximation example
Started by ●July 30, 2003
Reply by ●July 30, 20032003-07-30
Rene, What is "fast" often depends on the low level instruction set of the DSP/computer. You could take a look at the Cordic algorithm for an atan2() or atan(). The underlying idea is that you are going to rotate an (x,y) vector to align with the positive x-axis through a series of known shifts that are computationally efficient to apply. The final result will be the negative of the sum of shifts applied to get as close to the x-axis as you wanted. Generally, each additional phase bit output costs a sign check, 2 shifts and 2 accumulates (conditionally add or subtract), plus any addressing overhead. There is a known scaling required (multiply) that can be at the end, beginning, or could be distributed. There is a discussion of it in an appendix in Frerking. Besides a typo, he needs to extend his index range to one lower to be able to compute atan2 on any angle in the [0,90] degree range he claims. You need to rotate the vector into the [0,90] degree range using x,y exchanges and sign adjustments, accumulating the corresponding phase rotations in the process, before you launch into his algorithm. Dirk Dirk A. Bell DSP Consultant "Ren�" <jomfrusti@image.dk> wrote in message news:c9a8aff3.0307300402.5dab2c11@posting.google.com...> Is it possible to get a fast fix point Arctan routine written in C > using the principle with a lookup table or perhaps another principle? > > Regards, > Ren�
Reply by ●July 30, 20032003-07-30
"Ren�" <jomfrusti@image.dk> wrote in message news:c9a8aff3.0307300402.5dab2c11@posting.google.com...> Is it possible to get a fast fix point Arctan routine written in C > using the principle with a lookup table or perhaps another principle?The lookup table algorithms that I have seen use one table for a coarse (in domain) approximation, and another to do linear interpolation on that. You don't say how much accuracy you need. -- glen
Reply by ●July 30, 20032003-07-30
On 30 Jul 2003 05:02:11 -0700, jomfrusti@image.dk (=?ISO-8859-1?Q?Ren=E9?=) wrote:>Is it possible to get a fast fix point Arctan routine written in C >using the principle with a lookup table or perhaps another principle? > >Regards, >Ren�Hi, Look-up table methods are the fastest, but require much memory. (Table look-up isn't computation-free though. You still have to compute the table address value.) Cordic takes many iterations (10-12) to get 0.1 degree accuracy. Chebyshev polynomials are popular. The problem is, arctan is a very non-linear function, and resistent to quick polynomial evaluation with reasonable accuracy. You didn't say what accuracy you need nor did you say how many multiplies and adds you have available to compute your arctan. [-Rick-]
Reply by ●July 30, 20032003-07-30
Rick Lyons wrote:> > On 30 Jul 2003 05:02:11 -0700, jomfrusti@image.dk > (=?ISO-8859-1?Q?Ren=E9?=) wrote: > > >Is it possible to get a fast fix point Arctan routine written in C > >using the principle with a lookup table or perhaps another principle? > > > >Regards, > >Ren� > > Hi, > > Look-up table methods are the fastest, but require much > memory. (Table look-up isn't computation-free though. > You still have to compute the table address value.) > > Cordic takes many iterations (10-12) to get 0.1 degree > accuracy. Chebyshev polynomials are popular. > > The problem is, arctan is a very non-linear function, > and resistent to quick polynomial evaluation with > reasonable accuracy. > > You didn't say what accuracy you need nor did you > say how many multiplies and adds you have available > to compute your arctan. > > [-Rick-] >One way to make very non-linear functions more amenable to polynomial evaluation is to break them into shorter, more tractable sections. Deciding where to cut is an art, and compromise is needed to make all sections the same length. For integer calculations, I usually divide the range into sections that correspond to the most significant N bits, making N small enough to meet my accuracy needs. There's an example at http://users.erols.com/jyavins/typek.htm with code in Forth. A thermocouple needs only quadratics, but higher-order polynomials can be used. One can trade off polynomial order for segment count. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●July 30, 20032003-07-30
A hardware multiplier is mighty powerful so a technique that utilzes the multiplier, such as polynomial approximation, tends to work well in dsp software. Matt jomfrusti@image.dk (Ren�) wrote in message news:<c9a8aff3.0307300402.5dab2c11@posting.google.com>...> Is it possible to get a fast fix point Arctan routine written in C > using the principle with a lookup table or perhaps another principle? > > Regards, > Ren�
Reply by ●July 31, 20032003-07-31
In article b90ff073.0307301726.5b92dfe5@posting.google.com, Matt Boytim at maboytim@yahoo.com wrote on 07/30/2003 21:26:> A hardware multiplier is mighty powerful so a technique that utilzes > the multiplier, such as polynomial approximation, tends to work well > in dsp software.yeah, but arctan() is pretty funky. you can approximate finite segments of the function pretty well, but not the whole thing. at least i had trouble. i tried a polynomial (that had odd symmetry) that was post-processed through something like pi/2*tanh(x) (tanh() is not hard to compute with polynomials because you can build it with exp() and a division). tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x)) = 2/(exp(-2*x) + 1) - 1 both arctan() and tanh() are odd symmetry functions that have asymptotic behavior at +/- inf. that pre-polynomial had to approximate: p(x) = log( (1 - 2/pi*arctan(x)) / (1 + 2/pi*arctan(x)) ) so arctan(x) = pi/2 * tanh( -p(x)/2 ) or arctan(x) = pi/(exp(p(x)) + 1) - pi/2 that was the best i could do to try to get arctan(x) to work in the middle and at the edges. i wonder how a standard C or Fortran (blech) math library does arctan()? has anyone got some code for that? do they just split it into |x|<1 and |x|>1 regions and run one of the two converging series on that? r b-j
Reply by ●July 31, 20032003-07-31
"robert bristow-johnson" <rbj@surfglobal.net> wrote in message news:BB4E0ED6.28F8%rbj@surfglobal.net... (snip)> i wonder how a standard C or Fortran (blech) math library does arctan()? > has anyone got some code for that? do they just split it into |x|<1 and > |x|>1 regions and run one of the two converging series on that?Here is the IBM OS/360 Fortran library DATAN function. The comments explain argument reduction. LATN TITLE 'ARCTANGENT FUNCTION (LONG)' 00900018 IHCLATAN CSECT 01800018 * ARCTANGENT FUNCTION (LONG) 02700018 * 1. REDUCE THE CASE TO THE 1ST OCTANT BY USING 03600018 * ATAN(-X) = -ATAN(X), ATAN(1/X) = PI/2-ATAN(X). 04500018 * 2. REDUCE FURTHER TO THE CASE /X/ LESS THAN TAN(PI/12) 05400018 * BY ATAN(X) PI/6+ATAN((X*SQRT3-1)/(X+SQRT3)). 06300018 * 3. FOR THE BASIC RANGE (X LESS THAN TAN(PI/12)), USE 07200018 * A FRACTIONAL APPROXIMATION. 08100018 SPACE 09000018 ENTRY DATAN 09900018 SPACE 10800018 GRA EQU 1 ARGUMENT POINTER 11700018 GRS EQU 13 SAVE AREA POINTER 12600018 GRR EQU 14 RETURN REGISTER 13500018 GRL EQU 15 LINK REGISTER 14400018 GR0 EQU 0 SCRATCH REGISTERS 15300018 GR1 EQU 1 16200018 GR2 EQU 14 17100018 FR0 EQU 0 ANSWER REGISTER 18000018 FR2 EQU 2 SCRATCH REGISTERS 18900018 FR4 EQU 4 19800018 FR6 EQU 6 20700018 SPACE 21600018 USING *,GRL 22500018 DATAN BC 15,LATAN 23400018 DC AL1(5) 24300018 DC CL5'DATAN' 25200018 SPACE 26100018 LATAN STM GRR,GRL,12(GRS) SAVE REGISTERS 27000018 L GR1,0(GRA) 27900018 LD FR0,0(GR1) OBTAIN ARGUMENT 28800018 L GR0,0(GR1) SAVE ARG FOR SIGN CONTROL 29700018 LPER FR0,FR0 AND SET SIGN POSITIVE 30600018 LD FR4,ONE 31500018 SR GR1,GR1 GR1, GR2 FOR DISTINGUISHING CASES 32400018 LA GR2,ZERO 33300018 CER FR0,FR4 34200018 BC 12,SKIP1 35100018 LDR FR2,FR4 IF X GREATER THAN 1, TAKE INVERSE 36000018 DDR FR2,FR0 AND INCREMEMENT GR1 BY 16 36900018 LDR FR0,FR2 37800018 LA GR1,16 38700018 SPACE 39600018 SKIP1 CE FR0,SMALL IF ARG LESS THAN 16**-7, ANS=ARG. 40500018 BC 12,READY THIS AVOIDS UNDERFLOW EXCEPTION 41400018 CE FR0,TAN15 42300018 BC 12,SKIP2 43200018 LDR FR2,FR0 IF X GREATER THAN TAN(PI/12), 44100018 MD FR0,RT3M1 REDUCE X TO (X*SQRT3-1)/(X+SQRT3) 45000018 SDR FR0,FR4 COMPUTE X*SQRT3-1 AS 45900018 ADR FR0,FR2 X*(SQRT3-1)-1+X 46800018 AD FR2,RT3 TO GAIN ACCURACY 47700018 DDR FR0,FR2 48600018 LA GR2,8(GR2) INCREMENT GR2 BY 8 49500018 SPACE 50400018 SKIP2 LDR FR6,FR0 COMPUTE ATAN OF REDUCED ARGUMENT BY 51300018 MDR FR0,FR0 ATAN(X) = X+X*X**2*F, WHERE 52200018 LD FR4,C7 F = C1+C2/(X**2+C3+C4/ 53100018 ADR FR4,FR0 (X**2+C5+C6/(X**2+C7))) 54000018 LD FR2,C6 54900018 DDR FR2,FR4 55800018 AD FR2,C5 56700018 ADR FR2,FR0 57600018 LD FR4,C4 58500018 DDR FR4,FR2 59400018 AD FR4,C3 60300018 ADR FR4,FR0 61200018 LD FR2,C2 62100018 DDR FR2,FR4 63000018 AD FR2,C1 63900018 MDR FR0,FR2 64800018 MDR FR0,FR6 65700018 ADR FR0,FR6 66600018 SPACE 67500018 READY AD FR0,0(GR1,GR2) DEPENDING ON THE CASE, 68400018 LNR GR1,GR1 EITHER ADD 0 OR PI/6, OR 69300018 SD FR0,ZERO(GR1) SUBTRACT FROM PI/3 OR PI/2 70200018 LPER FR0,FR0 71100018 LTR GR0,GR0 DO THE LATTER IN TWO STEPS 72000018 BC 10,*+6 IF SIGN WAS NEGATIVE, 72900018 LNER FR0,FR0 ANSWER IS NEGATIVE 73800018 SPACE 74700018 EXIT L GRR,12(GRS) 75600018 MVI 12(GRS),X'FF' 76500018 BCR 15,GRR RETURN 77400018 SPACE 78300018 DS 0F 79200018 SMALL DC X'3A100000' 80100018 DS 0D 81000018 C1 DC X'BF1E31FF1784B965' -0.7371899082768562E-2 81900018 C2 DC X'C0ACDB34C0D1B35D' -0.6752198191404210 82800018 C3 DC X'412B7CE45AF5C165' 0.2717991214096480E+1 83700018 C4 DC X'C11A8F923B178C78' -0.1660051565960002E+1 84600018 C5 DC X'412AB4FD5D433FF6' 0.2669186939532663E+1 85500018 C6 DC X'C02298BB68CFD869' -0.1351430064094942 86400018 C7 DC X'41154CEE8B70CA99' 0.1331282181443987E+1 87300018 RT3M1 DC X'40BB67AE8584CAA8' SQRT(3)-1 88200018 ONE DC X'4110000000000000' THESE 89100018 RT3 DC X'411BB67AE8584CAB' SQRT(3) SIX 90000018 ZERO DC D'0' 0 CONSTANTS 90900018 DC X'40860A91C16B9B2C' PI/6 MUST 91800018 DC X'C0921FB54442D184' -PI/2+1 BE 92700018 DC X'BFC152382D736574' -(PI/3-F)+1 CONSECUTIVE 93600018 TAN15 DC X'40449851' 94500018 END 95400018
Reply by ●July 31, 20032003-07-31
>>>>> "robert" == robert bristow-johnson <rbj@surfglobal.net> writes:robert> i wonder how a standard C or Fortran (blech) math library does arctan()? robert> has anyone got some code for that? do they just split it into |x|<1 and There's all kinds of C and Fortran code for this at www.netlib.org. robert> |x|>1 regions and run one of the two converging series on that? I'd use the identity atan(x) = pi/2 - atan(1/x) for x > 1. Then only one series is needed. Hart et. al. has tables for smaller ranges and explains how to do compute atan this way. But to do this well takes some care because you want monotonicity and everything, just like the real function. Ray
Reply by ●July 31, 20032003-07-31
robert bristow-johnson <rbj@surfglobal.net> wrote in message news:<BB4E0ED6.28F8%rbj@surfglobal.net>...> In article b90ff073.0307301726.5b92dfe5@posting.google.com, Matt Boytim at > maboytim@yahoo.com wrote on 07/30/2003 21:26: > > > A hardware multiplier is mighty powerful so a technique that utilzes > > the multiplier, such as polynomial approximation, tends to work well > > in dsp software. > > yeah, but arctan() is pretty funky. you can approximate finite segments of > the function pretty well, but not the whole thing. at least i had trouble.I guess I glossed over the divide. Since atn(x)=+-pi/2-atn(1/x) you only need atn(x) for |x|<1. I don't remember any particular difficulties with the polynomial approximation in this range, and of coures it's odd which allows calculation in x^2 instead of x. Similarly, with atn2 you can use symmetry arguments to get it withn +-pi/4 and compute the same poly in y/x. So in either case you (may) have a divide. Logic, branching, etc. tends to be pretty costly on dsp. Often it is better to just add a few orders to the poly and launch into the calculation as soon as it is managable than taking further advantage of symmetry, etc. This is based on fixed point dsp experience; maybe not so for floating point. In any event, an algorithm which is mac centric such as polynomial approximation is usually a good choice for dsp's. Matt> > i tried a polynomial (that had odd symmetry) that was post-processed through > something like pi/2*tanh(x) (tanh() is not hard to compute with polynomials > because you can build it with exp() and a division). > > tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x)) > > = 2/(exp(-2*x) + 1) - 1 > > both arctan() and tanh() are odd symmetry functions that have asymptotic > behavior at +/- inf. > > that pre-polynomial had to approximate: > > p(x) = log( (1 - 2/pi*arctan(x)) / (1 + 2/pi*arctan(x)) ) > > so > > arctan(x) = pi/2 * tanh( -p(x)/2 ) > > or > > arctan(x) = pi/(exp(p(x)) + 1) - pi/2 > > > that was the best i could do to try to get arctan(x) to work in the middle > and at the edges. > > i wonder how a standard C or Fortran (blech) math library does arctan()? > has anyone got some code for that? do they just split it into |x|<1 and > |x|>1 regions and run one of the two converging series on that? > > r b-j
