arctangent in Spin

SRLM · 2008-11-30 04:15

How would I go about finding the arctangent of an angle in spin? I've looked at the Float32 object, but I don't want to have to dedicate a cog to it, and I can't extract the assembly code into a form that I can understand. My guess is that I'll have to access the sine table several times, mix those values together some, and see what I get. Any solutions?

Leon · 2008-11-30 05:15

You could use a series to calculate it: x - x^3/3 + x^5/5 - x^7/7 + x^9/9 ...

Leon

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Amateur radio callsign: G1HSM
Suzuki SV1000S motorcycle

SRLM · 2008-11-30 06:34

Hmmm. I never thought about using non trigonometric methods. I guess I should have done this before, but I looked up on wikipedia. There are quite a few ways to solve for arctangent using non trigonometric methods, mostly involving sums of one kind or another. When looking, one caught my eye (printed below), and it was created by Leonhard Euler. Any relation? [noparse]:)[/noparse]

$4a28884d933c74458514c8ee6c6bab1e.png$

In addition to your method posted above, I like this one:

$79a2a1a2beefcbe3c6556528fc976d79.png$

And this one:

$9f4887a0a2ad3a58de53502876cabedd.png$

I'm leaning towards your sum solution, with the logarithmic solution in second. Now, the challenge(for me at least) is to figure out how to do an infinite series on the propeller.

Rayman · 2008-11-30 11:32

I've done atn in assembly before using what I think is the standard technique for microcontrollers, i.e., reducing the range to 0 to pi/12 and then using a very good series approximation...

Unfortunately, my code is somewhat hard to follow, but I found a similar one here:

http://www.dsprelated.com/showmessage/21872/1.php

Look at the comments in this section:

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)))

Rayman · 2008-11-30 11:35

Actually, I see I used a much simpler series, only 3 coefficients and two of them were the same...

Here's my series code:

sciTrigSeries:
·····;store x in RegY
····LD IY,#OP1
····LD HL,#RegY
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
·····;calc x^2
····LD IY,#OP1
····LD HL,#OP2
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
····CAR FMUL
·····;add c3
····LD IY,#C3
····LD HL,#OP2
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
····CAR FADD
·····;store denomenator in RegX
····LD IY,#OP1
····LD HL,#RegX
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
·····;load c2
····LD IY,#C2
····LD HL,#OP1
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
·····;multiply by x twice
····LD IY,#RegY
····LD HL,#OP2
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
····CAR FMUL
····CAR FMUL
·····;add c1
····LD IY,#C1
····LD HL,#OP2
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
····CAR FADD
·····;multiply by x
····LD IY,#RegY
····LD HL,#OP2
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
····CAR FMUL
·····;divide by intermediate result
····LD IY,#RegX
····LD HL,#OP2
····CAR FCOPY···;copy float at [noparse][[/noparse]IY] to [noparse][[/noparse]HL]
····CAR FDIV

And, here's my coefficients:

C1:
····db 03Fh····;1.6867629106
····dw 0D7E8h
C2:
····db 03Dh····;0.4378497304
····dw 0E02Eh
C3:
····db 03Fh····;1.6867633134· ;hmm, same as C1 !
····dw 0D7E8h

SRLM · 2008-11-30 15:39

Wait, what happened to Beau? I thought that he posted here...

Beau Schwabe · 2008-11-30 17:20

SRLM,

I did, but I realized that you were asking about an ARCTAN solution after I had posted an ARCSIN solution, so I removed my post.

How much accuracy and what kind of speed do you need? This will ultimately determine which method you will implement.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Beau Schwabe

IC Layout Engineer
Parallax, Inc.

Post Edited (Beau Schwabe (Parallax)) : 11/30/2008 5:31:25 PM GMT

Beau Schwabe · 2008-11-30 20:41

SRLM,

Here is a polynomial solution that might work... I tested it and it seems to be within 1% or better of the actual ARCTAN radian value.
The largest error is near 0 deg.

x ranges from 0 to 100 (<- 0 to 45 Deg)
Dtheta is in degrees
Rtheta is in radians

Dtheta = 619550 * x - 1659 * x * x - 300970
Dtheta = Dtheta / 1000000
Rtheta = ( Dtheta * Pi ) / 180

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Beau Schwabe

IC Layout Engineer
Parallax, Inc.

Rayman · 2008-11-30 22:57

I think you can do much better with the way I outlined above... You can do the math using the "FloatMath" object, which doesn't use a cog.
You can get accuracy to 1 ppm this way...

But, I guess it depends on hou accurate you want to be and whether you want to do an all integer calculation or not...

Andrew E Mileski · 2008-12-01 03:30

From a C math library for fixed point S31.32 I wrote long ago:

The standard series for atan converges too slowly, doing this instead is extremely fast:

1st iteration: atan(x) = approximation1 + something1
2nd iteration: atan(x) = approximation1 + approximation2 + something2
3nd iteration: atan(x) = approximation1 + approximation2 + approximation3 + something3

Where an approximation is 2-terms of the standard series: approximation = x - x³ / 3

The trick is in figuring out how to calculate the "something". If "something" is expressed as an arctan itself, then it can also be approximated:

1st iteration: atan(x) = approximation of atan(x) + atan(A)
2nd iteration: atan(x) = approximation of atan(x) + approximation of atan(A) + atan(B)
3nd iteration: atan(x) = approximation of atan(x) + approximation of atan(A) + approximation of atan(B) + atan(C)

The final formula is:

atan(x) = x - x³ / 3 + atan((tan(x) - t) / (1 + x * t)), where t = tan(x - x³ / 3)

To further speed things up, x is scaled to be centered around pi /2

/*
 * atan x = SUM[noparse][[/noparse]n=0,) (-1)^n * x^(2n + 1)/(2n + 1), |x| < 1
 *
 * Two term approximation
 *    a = x - x^3/3
 * Gives us
 *    atan x = a + ??
 * Let ?? = arctan ?
 *    atan x = a + arctan ?
 * Rearrange
 *    atan x - a = arctan ?
 * Apply tan to both sides
 *    tan (atan x - a) = tan arctan ?
 *    tan (atan x - a) = ?
 * Applying the standard formula
 *    tan (u - v) = (tan u - tan v) / (1 + tan u * tan v)
 * Gives us
 *    tan (atan x - a) = (tan atan x - tan a) / (1 + tan arctan x * tan a)
 * Let t = tan a
 *    tan (atan x - a) = (x - t) / (1 + x * t)
 * So finally
 *    arctan x = a + arctan ((tan x - t) / (1 + x * t))
 * And the typical worst case is x = 1.0 which converges in 3 iterations.
 */
sll _sllatan(sll x)
{
    sll a, t, retval;

    /* First iteration */
    a = sllmul(x, _sllsub(CONST_1, sllmul(x, sllmul(x, CONST_1_3))));
    retval = a;

    /* Second iteration */
    t = slldiv(_sllsin(a), _sllcos(a));
    x = slldiv(_sllsub(x, t), _slladd(CONST_1, sllmul(t, x)));
    a = sllmul(x, _sllsub(CONST_1, sllmul(x, sllmul(x, CONST_1_3))));
    retval = _slladd(retval, a);

    /* Third  iteration */
    t = slldiv(_sllsin(a), _sllcos(a));
    x = slldiv(_sllsub(x, t), _slladd(CONST_1, sllmul(t, x)));
    a = sllmul(x, _sllsub(CONST_1, sllmul(x, sllmul(x, CONST_1_3))));
    return _slladd(retval, a);
}

sll sllatan(sll x)
{
    sll retval;

    if (x < sllneg(CONST_1))
        /* if x < -1 then atan x = PI/2 + atan 1/x */
        retval = sllneg(CONST_PI_2);
    else if (x > CONST_1)
        /* if x > 1 then atan x = PI/2 - atan 1/x */
        retval = CONST_PI_2;
    else
        return _sllatan(x);
    return _sllsub(retval, _sllatan(sllinv(x)));
}

Post Edited (Andrew E Mileski) : 12/1/2008 4:46:09 AM GMT

arctangent in Spin

Comments