# Arctan approximation example

Started by July 30, 2003
```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&#2013265929;
```
```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
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&#2013265929;" <jomfrusti@image.dk> wrote in message
> 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&#2013265929;

```
```"Ren&#2013265929;" <jomfrusti@image.dk> wrote in message
> 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

```
```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&#2013265929;

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

[-Rick-]

```
```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&#2013265929;
>
> 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
>
> [-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.
&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
```
```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&#2013265929;) 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&#2013265929;
```
```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.

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

```
```"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.

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
14400018
GR0    EQU     0               SCRATCH REGISTERS
15300018
GR1    EQU     1
16200018
GR2    EQU     14
17100018
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
46800018
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
54000018
LD      FR2,C6
54900018
DDR     FR2,FR4
55800018
56700018
57600018
LD      FR4,C4
58500018
DDR     FR4,FR2
59400018
60300018
61200018
LD      FR2,C2
62100018
DDR     FR2,FR4
63000018
63900018
MDR     FR0,FR2
64800018
MDR     FR0,FR6
65700018
66600018
SPACE
67500018
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
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

```
```>>>>> "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
```
```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
```