Quick and dirty four-quadrant integer-only arctangent in Spin

Phil Pilgrim (PhiPi)

For the S2 object I'm writing for the new Scribbler II robot, I needed a quick arctangent function using integer arithmetic, so I could convert (Δx, Δy) target locations to turn angles. I came across a really neat approximation in Abramowitz and Stegun's Handbook of Mathematical Functions, viz:

····arctan(x) = x / (1 + 0.28x²) + ε(x), where |ε(x)| ≤ .005 over -1 ≤ x ≤ 1

I used this approximation to craft the following four-quadrant atan2 function:

[b]PUB[/b] atan2(y, x) | arg, adder, n

'' Four-quadrant arctangent. Y and X are signed integers.
'' FULL_CIRCLE is an integer constant equal to the number of units in a full circle.

  [b]if[/b] ((n := >|(||x <# ||y)) > 21)
    x ~>= n - 21
    y ~>= n - 21 
  [b]if[/b] (||x > ||y)                      
    arg := y << 10 / x 
    adder := constant(FULL_CIRCLE / 2) & (x < 0)
  [b]else[/b]
    arg := -x << 10 / y
    adder := constant(FULL_CIRCLE / 4) + constant(FULL_CIRCLE / 2) & (y < 0)
  [b]result[/b] := arg * constant(FULL_CIRCLE * 14551 / 100000) / (936 + (arg * arg) >> 12) + adder
  [b]result[/b] += FULL_CIRCLE & ([b]result[/b] < 0)
    

Next, I generated 360 data points to test it with, using the following Perl program:

use strict;

foreach (0 .. 359) {
  my $angle = $_ * 3.14159265  / 180;
  my $x = int(cos($angle) * 100000);
  my $y = int(sin($angle) * 100000);
  print "               long      $x,$y\n"
}

By plugging the generated data into the attached Spin program, I was able to generate arctangents with 1/10-degree precision (not the same as accuracy, as we shall see). The results were compared to the correct values, and the differences were plotted against the actual angle:



As you can see the result is accurate to within ± 0.3 degrees, which seems to be good enough for the S2.

-Phil

Addendum: I just noticed that I could have done a much better job of managing truncation error, which could be a factor for small values of FULL_CIRCLE. Tomorrow's job...
_

Post Edited (Phil Pilgrim (PhiPi)) : 5/7/2010 6:58:17 AM GMT

heater

Very cunning.

"if ((n := >|(||x <# ||y)) > 21)" is a statement to be proud of[noparse]:)[/noparse]

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
For me, the past is not over yet.

BradC

heater said...
Very cunning.

"if ((n := >|(||x <# ||y)) > 21)" is a statement to be proud of[noparse]:)[/noparse]

It's a statement only a hardened perl monk would dare dream up...

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
"Are you suggesting coconuts migrate?"

heater

Glad you said that. I had to try very hard not to make a perl comparison.

How about a PASM version?

Might be easier to read[noparse]:)[/noparse]

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
For me, the past is not over yet.

Kye

Mmm, so I thought of a much better way to do arctan using arcsin and table lookup. I have not tested my code yet however to see if it works but the idea is pretty easy and if you could get the code working I'll be able to add a nice and fast integer math object to the obex.



▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Nyamekye,

Phil Pilgrim (PhiPi)

Kye, the reason I didn't use a reverse "table lookup" is that I didn't want the extra compute time required for a binary search. Although your table inversion is undoubtedly more accurate, the quicker approximation is "good enough" for my purposes.

-Phil

BradC

Phil Pilgrim (PhiPi) said...
Kye, the reason I didn't use a reverse "table lookup" is that I didn't want the extra compute time required for a binary search. Although your table inversion is undoubtedly more accurate, the quicker approximation is "good enough" for my purposes.

And likely an imperial crapload smaller (imperial values empirically appear to be bigger than metric as a whole).

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
"Are you suggesting coconuts migrate?"

Phil Pilgrim (PhiPi)

In order to reduce the errors as much as possible for the S2, where FULL_CIRCLE = 960, I've optimized the method a bit for that particular case:

[b]PUB[/b] atan2(y, x) | arg, adder, n

  '' Four-quadrant arctangent. Y and X are signed integers.
  '' FULL_CIRCLE is an integer constant equal to the number of units in a full circle.

  [b]if[/b] ((n := >|(||x <# ||y)) > 21)
    x ~>= n - 21
    y ~>= n - 21 
  [b]if[/b] (||x > ||y)                      
    arg := y << 10 / x 
    adder := constant(FULL_CIRCLE / 2) & (x < 0)
  [b]else[/b]
    arg := -x << 10 / y
    adder := constant(FULL_CIRCLE / 4) + constant(FULL_CIRCLE / 2) & (y < 0)
  [b]result[/b] := (||arg * constant(FULL_CIRCLE * 56841 / 100000) / (914 + (arg * arg) >> 12) + 1) >> 2 - (||arg => 960)
  [b]if[/b] (arg < 0)
    - [b]result[/b] 
  [b]result[/b] += adder
  [b]result[/b] += FULL_CIRCLE & ([b]result[/b] < 0)     

The (||arg => 960)* is a hand optimization based on an observation of the error data. With this addition and the slight amount of rounding that I've added, the errors in the 960-point set of test data have been eliminated entirely. This does not mean, of course, that arguments between those in the test set won't produce results that are a little off; but I'm confident that the method will work in my app.

As a side note, it would be nice if Spin had a sign transfer operator, such that a +- b would return a if b => 0 and -a if b < 0. The negc instruction in PASM makes this easy. There are ways to inline such an operation in Spin now, but they are either excessively wordy and/or involve a multiply.

-Phil

* The 960 here is merely a coincidence, having nothing to do with the value of FULL_CIRCLE.

Jay Kickliter

Forgive my ignorance, but what is ε(x)? Is it the error function? Where is it represented in your code? Thanks.

Phil Pilgrim (PhiPi)

Jay,

It's just an expression of the error as a function of x: IOW, the difference between the real arctangent and the value approximated by the equation. Its only representation in my code is where I attempt to minimize it with the correction factor (||arg => 960).

-Phil