Tracy Allen cordic examples

Peter Verkaik

Hi,
I converted the·basic stamp cordic examples at the top of page
http://www.emesystems.com/BS2mathC.htm
to the javelin and get these results

Test program CORDIC
cos(· 0.0) =· 0.9997·· sin(· 0.0) =· 0.0000
cos(· 0.1) =· 0.9997·· sin(· 0.1) =· 0.0017
cos(· 1.0) =· 0.9996·· sin(· 1.0) =· 0.0175
cos( 30.0) =· 0.8659·· sin( 30.0) =· 0.5000
cos( 45.0) =· 0.7071·· sin( 45.0) =· 0.7071
cos( 60.0) =· 0.5000·· sin( 60.0) =· 0.8659
cos( 89.0) =· 0.0175·· sin( 89.0) =· 0.9996
cos( 89.9) =· 0.0017·· sin( 89.9) =· 0.9997
cos( 90.0) =· 0.0000·· sin( 90.0) =· 0.9997
cos(· 0.1) =· 0.9997·· sin(· 0.1) = -0.0017
cos( -1.0) =· 0.9996·· sin( -1.0) = -0.0175
cos(-30.0) =· 0.8660·· sin(-30.0) = -0.4998
cos(-45.0) =· 0.7071·· sin(-45.0) = -0.7071
cos(-60.0) =· 0.5000·· sin(-60.0) = -0.8659
cos(-89.0) =· 0.0175·· sin(-89.0) = -0.9996
cos(-89.9) =· 0.0018·· sin(-89.9) = -0.9997
cos(-90.0) =· 0.0000·· sin(-90.0) = -0.9997
x =· 0.0001·· y =· 0.0001·· len =· 0.0001·· atn =· 43.3
x =· 0.0010·· y =· 0.0010·· len =· 0.0013·· atn =· 43.3
x =· 0.0100·· y =· 0.0100·· len =· 0.0141·· atn =· 45.0
x =· 0.1000·· y =· 0.1000·· len =· 0.1413·· atn =· 45.0
x =· 1.0000·· y =· 1.0000·· len =· 1.4139·· atn =· 45.0
x =· 0.8660·· y =· 0.5000·· len =· 0.9997·· atn =· 30.0
x = -0.8660·· y =· 0.5000·· len =· 0.9997·· atn =· 120.0
x = -0.8660·· y = -0.5000·· len =· 0.9997·· atn = -120.0
x =· 0.8660·· y = -0.5000·· len =· 0.9997·· atn = -30.0
x =· 0.5000·· y =· 0.8660·· len =· 0.9997·· atn =· 60.0
x = -0.5000·· y =· 0.8660·· len =· 0.9997·· atn =· 150.0
x = -0.5000·· y = -0.8660·· len =· 0.9997·· atn = -150.0
x =· 0.5000·· y = -0.8660·· len =· 0.9997·· atn = -60.0
program finished

All results are ok, except the calculated arctan values +-120.0 and +-150.0
which should be swapped according to the x and y values.
Can anyone confirm this for the basic stamp?
Also, is the calculated arctan value for x=y=0.0001 correct? (should be 45.0)
but as said the error increases for low x and y.
Also, why is the cordic gain (0.607*65536) set to constant 39793 as it should be 39780?

regards peter

·

Peter Verkaik

I found the error in cordic-atn.bse

· IF quadrant23 THEN················· ' wrapper for 2nd and 3rd quadrants
··· IF zs THEN z=z-9000 ELSE z=z+9000
··· ENDIF

should be

· IF quadrant23 THEN················· ' wrapper for 2nd and 3rd quadrants
··· IF zs THEN z=-18000-z ELSE z=18000-z
··· ENDIF
Still want to know about the cordic gain and if the accuracy for low x and y
can be improved.

regards peter
·

Tracy Allen

Hi Peter,

You've been busy!

Thanks for the heads up about the error for the quadrant23 business. I'll correct that in the program (but not today--too much holiday stuff to do!). One of the first steps in the program is to move points that happen to be in quadrants 2 or 3 over into quadrant 1 and 4 by reflection in the y axis (x=-x). So once the angle is computed it makes perfect sense to reflect the computed angles back through the y axis and that is what your correction does. I don't know how the rotation by 90 degrees got in there. There must have been dust in the afterburner.

As to the Cordic gain, when I first wrote the programs, I took 0.6072 as an approximation, and then when I explained it in the text I rounded it off to 0.607. That accounts for the difference between 39793 and 39780. But neither one is right. A closer approximation is 0.607252935, and for that the Stamp ** factor should be 39797. I'll work that into the program too, but that probably won't make much difference.

I did go into the www page and changed the exposition. Always tinkering.

As to the errors to calculate arctangent when x=y=0.0001. It is due to integer roundoff. The algorithm needs digits to work with internally to carry through the precision. One solution would be to put another wrapper on it, so that if the input numbers are "small", they will be multiplied times a constant at the outset and then rounded down or handled appropriately at the end. It would be kind of ad-hoc floating point, so that the algorithm always works at maximum precision. That way one could represent the angle in that particular case correctly as 45 degrees and also the length of the vector at 0.00007071. There is no way to get the 0.00007071 without scaling it at the outset. Another way to do it is to go to double precision, 32 bits, but call it 24 bit precision and keep 8 throwaway bits in reserve to slosh around internally. On the demo page, I was more interested in just getting it working on the Stamp--without getting into the complications of double precision.

Here is the sequence of values of the product of the {cos(atn(1/2^i)}, to show how fast it converges.
0.707106781186547
0.632455532033676
0.613571991077896
0.608833912517752
0.607648256256168
0.607351770141296
0.607277644093526
0.607259112298893
0.607254479332562
0.607253321089875
0.607253031529134
0.607252959138945
0.607252941041397
0.607252936517010
0.607252935385914
0.607252935103139
0.607252935032446 <-- 16th product on Stamp **39797
0.607252935014772
0.607252935010354
0.607252935009249
0.607252935008973
0.607252935008904
0.607252935008887
0.607252935008883
0.607252935008882
0.607252935008881
0.607252935008881
0.607252935008881
0.607252935008881

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Peter Verkaik

Thanks Tracy for the explanantion.
I found the factor 1/1.64676 = 0.607253 in other references, and this yields 39797.
When I use·that my calculation for sin(90)=0.9998, sin(45)=0.7071
Increasing to 39798 yields sin(90)=0.9999, sin(45)=0.7072
So by using 39798 I have·0.0001 errors overall, instead of 0.0000 and 0.0002 errors.
Tweaking this constant appears to allow spreeading the errors evenly.
Changing the constant did not change the calculated atn angles.
Here are the results for 39798

Test program CORDIC

cos(· 0.0) =· 0.9999·· sin(· 0.0) =· 0.0000
cos(· 0.1) =· 0.9999·· sin(· 0.1) =· 0.0017
cos(· 0.4) =· 0.9999·· sin(· 0.4) =· 0.0069
cos(· 0.5) =· 0.9999·· sin(· 0.5) =· 0.0087
cos(· 1.0) =· 0.9997·· sin(· 1.0) =· 0.0176
cos( 15.0) =· 0.9660·· sin( 15.0) =· 0.2588
cos( 30.0) =· 0.8660·· sin( 30.0) =· 0.5000
cos( 45.0) =· 0.7072·· sin( 45.0) =· 0.7072
cos( 60.0) =· 0.5000·· sin( 60.0) =· 0.8660
cos( 75.0) =· 0.2588·· sin( 75.0) =· 0.9660
cos( 89.0) =· 0.0176·· sin( 89.0) =· 0.9997
cos( 89.4) =· 0.0104·· sin( 89.4) =· 0.9999
cos( 89.5) =· 0.0086·· sin( 89.5) =· 0.9999
cos( 89.9) =· 0.0017·· sin( 89.9) =· 0.9999
cos( 90.0) =· 0.0000·· sin( 90.0) =· 0.9999
cos(· 0.1) =· 0.9999·· sin(· 0.1) = -0.0017
cos(· 0.4) =· 0.9999·· sin(· 0.4) = -0.0069
cos(· 0.5) =· 0.9999·· sin(· 0.5) = -0.0086
cos( -1.0) =· 0.9997·· sin( -1.0) = -0.0176
cos(-15.0) =· 0.9660·· sin(-15.0) = -0.2588
cos(-30.0) =· 0.8660·· sin(-30.0) = -0.5000
cos(-45.0) =· 0.7072·· sin(-45.0) = -0.7072
cos(-60.0) =· 0.5000·· sin(-60.0) = -0.8660
cos(-75.0) =· 0.2588·· sin(-75.0) = -0.9660
cos(-89.0) =· 0.0176·· sin(-89.0) = -0.9997
cos(-89.4) =· 0.0104·· sin(-89.4) = -0.9999
cos(-89.5) =· 0.0087·· sin(-89.5) = -0.9999
cos(-89.9) =· 0.0017·· sin(-89.9) = -0.9999
cos(-90.0) =· 0.0000·· sin(-90.0) = -0.9999
x =· 0.0001·· y =· 0.0001·· len =· 0.0001·· atn =· 43.3
x =· 0.0010·· y =· 0.0010·· len =· 0.0013·· atn =· 43.3
x =· 0.0100·· y =· 0.0100·· len =· 0.0141·· atn =· 45.0
x =· 0.1000·· y =· 0.1000·· len =· 0.1413·· atn =· 45.0
x =· 1.0000·· y =· 1.0000·· len =· 1.4141·· atn =· 45.0
x =· 0.8660·· y =· 0.5000·· len =· 0.9998·· atn =· 30.0
x = -0.8660·· y =· 0.5000·· len =· 0.9998·· atn =· 150.0
x = -0.8660·· y = -0.5000·· len =· 0.9998·· atn = -150.0
x =· 0.8660·· y = -0.5000·· len =· 0.9998·· atn = -30.0
x =· 0.5000·· y =· 0.8660·· len =· 0.9998·· atn =· 60.0
x = -0.5000·· y =· 0.8660·· len =· 0.9998·· atn =· 120.0
x = -0.5000·· y = -0.8660·· len =· 0.9998·· atn = -120.0
x =· 0.5000·· y = -0.8660·· len =· 0.9998·· atn = -60.0
program finished

I think I will try the extra scaling for small input numbers.

regards peter

Peter Verkaik

Tracy, I applied the scaling such that while abs(x) and abs(y) are both <=5000,
I double both x and y (10000 is my unit vector length), and at the end I halve
len for each scale step. That really gives great results (error <= 0.0002 in len)
The results
cos(· 0.0) =· 0.9999·· sin(· 0.0) =· 0.0000
cos(· 0.1) =· 0.9999·· sin(· 0.1) =· 0.0017
cos(· 0.4) =· 0.9999·· sin(· 0.4) =· 0.0069
cos(· 0.5) =· 0.9999·· sin(· 0.5) =· 0.0087
cos(· 1.0) =· 0.9997·· sin(· 1.0) =· 0.0176
cos( 15.0) =· 0.9660·· sin( 15.0) =· 0.2588
cos( 30.0) =· 0.8660·· sin( 30.0) =· 0.5000
cos( 45.0) =· 0.7072·· sin( 45.0) =· 0.7072
cos( 60.0) =· 0.5000·· sin( 60.0) =· 0.8660
cos( 75.0) =· 0.2588·· sin( 75.0) =· 0.9660
cos( 89.0) =· 0.0176·· sin( 89.0) =· 0.9997
cos( 89.4) =· 0.0104·· sin( 89.4) =· 0.9999
cos( 89.5) =· 0.0086·· sin( 89.5) =· 0.9999
cos( 89.9) =· 0.0017·· sin( 89.9) =· 0.9999
cos( 90.0) =· 0.0000·· sin( 90.0) =· 0.9999
cos(- 0.1) =· 0.9999·· sin(- 0.1) = -0.0017
cos(- 0.4) =· 0.9999·· sin(- 0.4) = -0.0069
cos(- 0.5) =· 0.9999·· sin(- 0.5) = -0.0086
cos(- 1.0) =· 0.9997·· sin(- 1.0) = -0.0176
cos(-15.0) =· 0.9660·· sin(-15.0) = -0.2588
cos(-30.0) =· 0.8660·· sin(-30.0) = -0.5000
cos(-45.0) =· 0.7072·· sin(-45.0) = -0.7072
cos(-60.0) =· 0.5000·· sin(-60.0) = -0.8660
cos(-75.0) =· 0.2588·· sin(-75.0) = -0.9660
cos(-89.0) =· 0.0176·· sin(-89.0) = -0.9997
cos(-89.4) =· 0.0104·· sin(-89.4) = -0.9999
cos(-89.5) =· 0.0087·· sin(-89.5) = -0.9999
cos(-89.9) =· 0.0017·· sin(-89.9) = -0.9999
cos(-90.0) =· 0.0000·· sin(-90.0) = -0.9999
x =· 0.0001·· y =· 0.0001·· len =· 0.0002·· atn =· 45.0
x =· 0.0010·· y =· 0.0010·· len =· 0.0015·· atn =· 45.0
x =· 0.0100·· y =· 0.0100·· len =· 0.0142·· atn =· 45.0
x =· 0.1000·· y =· 0.1000·· len =· 0.1415·· atn =· 45.0
x =· 1.0000·· y =· 1.0000·· len =· 1.4141·· atn =· 45.0
x =· 0.8660·· y =· 0.5000·· len =· 0.9998·· atn =· 30.0
x = -0.8660·· y =· 0.5000·· len =· 0.9998·· atn =· 150.0
x = -0.8660·· y = -0.5000·· len =· 0.9998·· atn = -150.0
x =· 0.8660·· y = -0.5000·· len =· 0.9998·· atn = -30.0
x =· 0.5000·· y =· 0.8660·· len =· 0.9998·· atn =· 60.0
x = -0.5000·· y =· 0.8660·· len =· 0.9998·· atn =· 120.0
x = -0.5000·· y = -0.8660·· len =· 0.9998·· atn = -120.0
x =· 0.5000·· y = -0.8660·· len =· 0.9998·· atn = -60.0

regards peter

Peter Verkaik

Tracy, I found that scaling is also necessary for big input values.
If either abs(x) or abs(y) > 19898 than the calculated atn changes
sign or decreases. This makes sense because x and y are scaled by the cordic gain
so 32767/1.64676 = 19898 is highest allowable input value
In my calculation up to 19900 was allowed before showing errors.
So for big input values both x and y must be scaled down.
I have applied that too and the calculated atn is correct, even for big values.
x=100 and y=20000 returns atn =89.7·which is correct.
regards peter
·

Tracy Allen

Yeah, great, that makes sense. I'll put that wrapper on the Stamp code too. I'm really glad you're taking an interest in this.

By the way, someone contacted me recently with an updated URL for the Knaust article, and I found a new link to Ingo Cyliax' BASIC Stamp code in the onliine reader for a course at the U. Alberta. Those updates are in the page at
www.emesystems.com/BS2mathC.htm

Have you posted you Javelin code (class?). I know Andy Lindsey was doing something with CORDIC a while ago for the SX or SX/B, too. I'll make sure he sees this.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Peter Verkaik

Tracy, I have not yet posted it.
I attached my test class. I will create a special cordic class with entries
sin(), cos(), tan(), atan2() etc.
I have attached my testclass for you to see how I adapted your code.
My code allows for angle input in 0.1 degree units, atn results
are also in 0.1 degree units (your code specified 0.01 degree units
for atn results but I think that is not necessary)

regards peter

Peter Verkaik

Tracy, the atan2(y,x) function really behaves weird when x=0
x =· 0.0000·· y =· 0.0000·· atan2(y,x) = - 80.1
x =· 0.0000·· y =· 0.0001·· atan2(y,x) = - 90.0
x =· 0.0000·· y =· 0.0010·· atan2(y,x) =·· 90.0
x =· 0.0000·· y =· 0.0100·· atan2(y,x) =·· 90.0
x =· 0.0000·· y =· 0.1000·· atan2(y,x) = - 90.0
x =· 0.0000·· y =· 1.0000·· atan2(y,x) =·· 90.0
x =· 0.0000·· y = -0.0001·· atan2(y,x) =·· 90.0
x =· 0.0000·· y = -0.0010·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -0.0100·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -0.1000·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -1.0000·· atan2(y,x) = - 90.0
x =· 0.0001·· y =· 0.0001·· atan2(y,x) =·· 45.0
x =· 0.0001·· y =· 0.0010·· atan2(y,x) =·· 84.3
x =· 0.0001·· y =· 0.0100·· atan2(y,x) =·· 89.4
x =· 0.0001·· y =· 0.1000·· atan2(y,x) =·· 89.9
x =· 0.0001·· y =· 1.0000·· atan2(y,x) =·· 90.0
x =· 0.0001·· y = -0.0001·· atan2(y,x) = - 45.0
x =· 0.0001·· y = -0.0010·· atan2(y,x) = - 84.3
x =· 0.0001·· y = -0.0100·· atan2(y,x) = - 89.4
x =· 0.0001·· y = -0.1000·· atan2(y,x) = - 89.9
x =· 0.0001·· y = -1.0000·· atan2(y,x) = - 90.0
x = -0.0001·· y =· 0.0001·· atan2(y,x) =· 135.0
x = -0.0001·· y =· 0.0010·· atan2(y,x) =·· 95.7
x = -0.0001·· y =· 0.0100·· atan2(y,x) =·· 90.6
x = -0.0001·· y =· 0.1000·· atan2(y,x) =·· 90.1
x = -0.0001·· y =· 1.0000·· atan2(y,x) =·· 90.0
x = -0.0001·· y = -0.0001·· atan2(y,x) = -135.0
x = -0.0001·· y = -0.0010·· atan2(y,x) = - 95.7
x = -0.0001·· y = -0.0100·· atan2(y,x) = - 90.6
x = -0.0001·· y = -0.1000·· atan2(y,x) = - 90.1
x = -0.0001·· y = -1.0000·· atan2(y,x) = - 90.0

When x=0 and y=0 it really is a discontinuity and requires special treatment.
When x=0 and y<>0 the result sometimes have correct sign, sometimes incorrect sign.
The absolute value of the result is correct.
So I test for x=0 as a special case and return angles based on sign of y. Then the results are
x =· 0.0000·· y =· 0.0000·· atan2(y,x) =·· 45.0
x =· 0.0000·· y =· 0.0001·· atan2(y,x) =·· 90.0
x =· 0.0000·· y =· 0.0010·· atan2(y,x) =·· 90.0
x =· 0.0000·· y =· 0.0100·· atan2(y,x) =·· 90.0
x =· 0.0000·· y =· 0.1000·· atan2(y,x) =·· 90.0
x =· 0.0000·· y =· 1.0000·· atan2(y,x) =·· 90.0
x =· 0.0000·· y = -0.0001·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -0.0010·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -0.0100·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -0.1000·· atan2(y,x) = - 90.0
x =· 0.0000·· y = -1.0000·· atan2(y,x) = - 90.0
x =· 0.0001·· y =· 0.0001·· atan2(y,x) =·· 45.0
x =· 0.0001·· y =· 0.0010·· atan2(y,x) =·· 84.3
x =· 0.0001·· y =· 0.0100·· atan2(y,x) =·· 89.4
x =· 0.0001·· y =· 0.1000·· atan2(y,x) =·· 89.9
x =· 0.0001·· y =· 1.0000·· atan2(y,x) =·· 90.0
x =· 0.0001·· y = -0.0001·· atan2(y,x) = - 45.0
x =· 0.0001·· y = -0.0010·· atan2(y,x) = - 84.3
x =· 0.0001·· y = -0.0100·· atan2(y,x) = - 89.4
x =· 0.0001·· y = -0.1000·· atan2(y,x) = - 89.9
x =· 0.0001·· y = -1.0000·· atan2(y,x) = - 90.0
x = -0.0001·· y =· 0.0001·· atan2(y,x) =· 135.0
x = -0.0001·· y =· 0.0010·· atan2(y,x) =·· 95.7
x = -0.0001·· y =· 0.0100·· atan2(y,x) =·· 90.6
x = -0.0001·· y =· 0.1000·· atan2(y,x) =·· 90.1
x = -0.0001·· y =· 1.0000·· atan2(y,x) =·· 90.0
x = -0.0001·· y = -0.0001·· atan2(y,x) = -135.0
x = -0.0001·· y = -0.0010·· atan2(y,x) = - 95.7
x = -0.0001·· y = -0.0100·· atan2(y,x) = - 90.6
x = -0.0001·· y = -0.1000·· atan2(y,x) = - 90.1
x = -0.0001·· y = -1.0000·· atan2(y,x) = - 90.0

I treat the case x=0 and y=0·just like·y=x and so 45.0 degrees is returned.
I will post my class shortly.

regards peter
·

Peter Verkaik

Tracy, I have uploaded my class here:
http://groups.yahoo.com/group/JavelinCode/files/Javelin%20Stamp%20IDE/lib/stamp/math/

Final results:
x=cos(· 0.0) =· 0.9999·· y=sin(· 0.0) =· 0.0000·· tan(· 0.0) =··· 0.00·· atan2(y,x) =··· 0.0
x=cos(· 0.1) =· 0.9999·· y=sin(· 0.1) =· 0.0017·· tan(· 0.1) =··· 0.00·· atan2(y,x) =··· 0.1
x=cos(· 0.4) =· 0.9999·· y=sin(· 0.4) =· 0.0069·· tan(· 0.4) =··· 0.01·· atan2(y,x) =··· 0.4
x=cos(· 0.5) =· 0.9999·· y=sin(· 0.5) =· 0.0087·· tan(· 0.5) =··· 0.01·· atan2(y,x) =··· 0.5
x=cos(· 1.0) =· 0.9997·· y=sin(· 1.0) =· 0.0176·· tan(· 1.0) =··· 0.02·· atan2(y,x) =··· 1.0
x=cos( 15.0) =· 0.9660·· y=sin( 15.0) =· 0.2588·· tan( 15.0) =··· 0.27·· atan2(y,x) =·· 15.0
x=cos( 30.0) =· 0.8660·· y=sin( 30.0) =· 0.5000·· tan( 30.0) =··· 0.58·· atan2(y,x) =·· 30.0
x=cos( 45.0) =· 0.7072·· y=sin( 45.0) =· 0.7072·· tan( 45.0) =··· 1.00·· atan2(y,x) =·· 45.0
x=cos( 60.0) =· 0.5000·· y=sin( 60.0) =· 0.8660·· tan( 60.0) =··· 1.73·· atan2(y,x) =·· 60.0
x=cos( 75.0) =· 0.2588·· y=sin( 75.0) =· 0.9660·· tan( 75.0) =··· 3.73·· atan2(y,x) =·· 75.0
x=cos( 89.0) =· 0.0176·· y=sin( 89.0) =· 0.9997·· tan( 89.0) =·· 56.80·· atan2(y,x) =·· 89.0
x=cos( 89.4) =· 0.0104·· y=sin( 89.4) =· 0.9999·· tan( 89.4) =·· 96.14·· atan2(y,x) =·· 89.4
x=cos( 89.5) =· 0.0086·· y=sin( 89.5) =· 0.9999·· tan( 89.5) =· 114.93·· atan2(y,x) =·· 89.5
x=cos( 89.8) =· 0.0034·· y=sin( 89.8) =· 0.9999·· tan( 89.8) =· 294.09·· atan2(y,x) =·· 89.8
x=cos( 89.9) =· 0.0017·· y=sin( 89.9) =· 0.9999·· tan( 89.9) =· 327.67·· atan2(y,x) =·· 89.9
x=cos( 90.0) =· 0.0000·· y=sin( 90.0) =· 0.9999·· tan( 90.0) =· 327.67·· atan2(y,x) =·· 90.0
x=cos( 90.1) = -0.0017·· y=sin( 90.1) =· 0.9999·· tan( 90.1) = -327.67·· atan2(y,x) =·· 90.1
x=cos( 90.2) = -0.0034·· y=sin( 90.2) =· 0.9999·· tan( 90.2) = -294.08·· atan2(y,x) =·· 90.2
x=cos(- 0.1) =· 0.9999·· y=sin(- 0.1) = -0.0017·· tan(- 0.1) =··· 0.00·· atan2(y,x) = -· 0.1
x=cos(- 0.4) =· 0.9999·· y=sin(- 0.4) = -0.0069·· tan(- 0.4) = -· 0.01·· atan2(y,x) = -· 0.4
x=cos(- 0.5) =· 0.9999·· y=sin(- 0.5) = -0.0086·· tan(- 0.5) = -· 0.01·· atan2(y,x) = -· 0.5
x=cos(- 1.0) =· 0.9997·· y=sin(- 1.0) = -0.0176·· tan(- 1.0) = -· 0.02·· atan2(y,x) = -· 1.0
x=cos(-15.0) =· 0.9660·· y=sin(-15.0) = -0.2588·· tan(-15.0) = -· 0.27·· atan2(y,x) = - 15.0
x=cos(-30.0) =· 0.8660·· y=sin(-30.0) = -0.5000·· tan(-30.0) = -· 0.58·· atan2(y,x) = - 30.0
x=cos(-45.0) =· 0.7072·· y=sin(-45.0) = -0.7072·· tan(-45.0) = -· 1.00·· atan2(y,x) = - 45.0
x=cos(-60.0) =· 0.5000·· y=sin(-60.0) = -0.8660·· tan(-60.0) = -· 1.73·· atan2(y,x) = - 60.0
x=cos(-75.0) =· 0.2588·· y=sin(-75.0) = -0.9660·· tan(-75.0) = -· 3.73·· atan2(y,x) = - 75.0
x=cos(-89.0) =· 0.0176·· y=sin(-89.0) = -0.9997·· tan(-89.0) = - 57.13·· atan2(y,x) = - 89.0
x=cos(-89.4) =· 0.0104·· y=sin(-89.4) = -0.9999·· tan(-89.4) = - 96.14·· atan2(y,x) = - 89.4
x=cos(-89.5) =· 0.0087·· y=sin(-89.5) = -0.9999·· tan(-89.5) = -114.93·· atan2(y,x) = - 89.5
x=cos(-89.8) =· 0.0034·· y=sin(-89.8) = -0.9999·· tan(-89.8) = -294.08·· atan2(y,x) = - 89.8
x=cos(-89.9) =· 0.0017·· y=sin(-89.9) = -0.9999·· tan(-89.9) = -327.67·· atan2(y,x) = - 89.9
x=cos(-90.0) =· 0.0000·· y=sin(-90.0) = -0.9999·· tan(-90.0) = -327.67·· atan2(y,x) = - 90.0
x=cos(-90.1) = -0.0017·· y=sin(-90.1) = -0.9999·· tan(-90.1) =· 327.67·· atan2(y,x) = - 90.1
x=cos(-90.2) = -0.0034·· y=sin(-90.2) = -0.9999·· tan(-90.2) =· 294.08·· atan2(y,x) = - 90.2

regards peter

Newzed

Peter, cos(1) = .9998 and sin(1) = .0174.· Are these minor differences important?

Sid

Peter Verkaik

Newzed, the absolute error is <= 0.0002 which I find very good considering
all values are calculated and not read from a table. Also note that the atan2
calculates the original angle given the x=cos and y=sin values (which have
an error in the last decimal as you remarked).

regards peter
·

Tracy Allen

Hi Sid,
SIN(1) = 0.0174524
COS(1) = 0.9998477

so the rounding is less than 0.0001 in those cases. The algorithm tends to round down.


Peter, I noticed on the Javelin forum that you also have some trig class code there, but I haven't looked at the Yahoo file site. Is that done with standard polynomials? I wonder, how do you find those compare in terms of speed and accuracy with the CORDIC?

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Peter Verkaik

Hi Tracy,
The Trig class uses a 90 words table with sine values for whole degrees. In between
whole degrees the sine and cosine are approximated using
lineair interpolation. The accuray is of the same order
as cordic I guess, but I haven't compare the cordic results
with the table results yet.
regards peter

Peter Verkaik

Tracy,
Here are the results using the 90 words table
x=cos(· 0.0) =· 1.0000·· y=sin(· 0.0) =· 0.0000·· tan(· 0.0) =··· 0.00·· atan2(y,x) =··· 0.0
x=cos(· 0.1) =· 0.9999·· y=sin(· 0.1) =· 0.0017·· tan(· 0.1) =··· 0.00·· atan2(y,x) =··· 0.1
x=cos(· 0.4) =· 0.9999·· y=sin(· 0.4) =· 0.0070·· tan(· 0.4) =··· 0.01·· atan2(y,x) =··· 0.4
x=cos(· 0.5) =· 0.9999·· y=sin(· 0.5) =· 0.0087·· tan(· 0.5) =··· 0.01·· atan2(y,x) =··· 0.5
x=cos(· 1.0) =· 0.9998·· y=sin(· 1.0) =· 0.0175·· tan(· 1.0) =··· 0.02·· atan2(y,x) =··· 1.0
x=cos( 15.0) =· 0.9659·· y=sin( 15.0) =· 0.2588·· tan( 15.0) =··· 0.27·· atan2(y,x) =·· 15.1
x=cos( 30.0) =· 0.8660·· y=sin( 30.0) =· 0.5000·· tan( 30.0) =··· 0.58·· atan2(y,x) =·· 30.2
x=cos( 45.0) =· 0.7071·· y=sin( 45.0) =· 0.7071·· tan( 45.0) =··· 1.00·· atan2(y,x) =·· 45.0
x=cos( 60.0) =· 0.5000·· y=sin( 60.0) =· 0.8660·· tan( 60.0) =··· 1.73·· atan2(y,x) =·· 59.8
x=cos( 75.0) =· 0.2588·· y=sin( 75.0) =· 0.9659·· tan( 75.0) =··· 3.73·· atan2(y,x) =·· 74.9
x=cos( 89.0) =· 0.0175·· y=sin( 89.0) =· 0.9998·· tan( 89.0) =·· 57.13·· atan2(y,x) =·· 89.0
x=cos( 89.4) =· 0.0105·· y=sin( 89.4) =· 0.9998·· tan( 89.4) =·· 95.22·· atan2(y,x) =·· 89.4
x=cos( 89.5) =· 0.0087·· y=sin( 89.5) =· 0.9999·· tan( 89.5) =· 114.93·· atan2(y,x) =·· 89.5
x=cos( 89.8) =· 0.0035·· y=sin( 89.8) =· 0.9999·· tan( 89.8) =· 285.68·· atan2(y,x) =·· 89.8
x=cos( 89.9) =· 0.0017·· y=sin( 89.9) =· 0.9999·· tan( 89.9) =· 327.67·· atan2(y,x) =·· 89.9
x=cos( 90.0) =· 0.0000·· y=sin( 90.0) =· 1.0000·· tan( 90.0) =· 327.67·· atan2(y,x) =·· 90.0
x=cos( 90.1) = -0.0017·· y=sin( 90.1) =· 0.9999·· tan( 90.1) = -327.67·· atan2(y,x) =·· 90.1
x=cos( 90.2) = -0.0035·· y=sin( 90.2) =· 0.9999·· tan( 90.2) = -285.68·· atan2(y,x) =·· 90.2
x=cos(- 0.1) =· 0.9999·· y=sin(- 0.1) = -0.0017·· tan(- 0.1) =··· 0.00·· atan2(y,x) = -· 0.1
x=cos(- 0.4) =· 0.9999·· y=sin(- 0.4) = -0.0070·· tan(- 0.4) = -· 0.01·· atan2(y,x) = -· 0.4
x=cos(- 0.5) =· 0.9999·· y=sin(- 0.5) = -0.0087·· tan(- 0.5) = -· 0.01·· atan2(y,x) = -· 0.5
x=cos(- 1.0) =· 0.9998·· y=sin(- 1.0) = -0.0175·· tan(- 1.0) = -· 0.02·· atan2(y,x) = -· 1.0
x=cos(-15.0) =· 0.9659·· y=sin(-15.0) = -0.2588·· tan(-15.0) = -· 0.27·· atan2(y,x) = - 15.1
x=cos(-30.0) =· 0.8660·· y=sin(-30.0) = -0.5000·· tan(-30.0) = -· 0.58·· atan2(y,x) = - 30.2
x=cos(-45.0) =· 0.7071·· y=sin(-45.0) = -0.7071·· tan(-45.0) = -· 1.00·· atan2(y,x) = - 45.0
x=cos(-60.0) =· 0.5000·· y=sin(-60.0) = -0.8660·· tan(-60.0) = -· 1.73·· atan2(y,x) = - 59.8
x=cos(-75.0) =· 0.2588·· y=sin(-75.0) = -0.9659·· tan(-75.0) = -· 3.73·· atan2(y,x) = - 74.9
x=cos(-89.0) =· 0.0175·· y=sin(-89.0) = -0.9998·· tan(-89.0) = - 57.13·· atan2(y,x) = - 89.0
x=cos(-89.4) =· 0.0105·· y=sin(-89.4) = -0.9998·· tan(-89.4) = - 95.22·· atan2(y,x) = - 89.4
x=cos(-89.5) =· 0.0087·· y=sin(-89.5) = -0.9999·· tan(-89.5) = -114.93·· atan2(y,x) = - 89.5
x=cos(-89.8) =· 0.0035·· y=sin(-89.8) = -0.9999·· tan(-89.8) = -285.68·· atan2(y,x) = - 89.8
x=cos(-89.9) =· 0.0017·· y=sin(-89.9) = -0.9999·· tan(-89.9) = -327.67·· atan2(y,x) = - 89.9
x=cos(-90.0) =· 0.0000·· y=sin(-90.0) = -1.0000·· tan(-90.0) = -327.67·· atan2(y,x) = - 90.0
x=cos(-90.1) = -0.0017·· y=sin(-90.1) = -0.9999·· tan(-90.1) =· 327.67·· atan2(y,x) = - 90.1
x=cos(-90.2) = -0.0035·· y=sin(-90.2) = -0.9999·· tan(-90.2) =· 285.68·· atan2(y,x) = - 90.2

It shows the sin and cos are slightly more accurate (and thus tan also), but atan2 is slightly worse.
This atan2 is calculated using the approximation

··················· (y/x)··········· 1800
··· atn =

---

* ----················ 0.1 deg units, for |y/x|<1
········· 1 + 0.280872*(y/x)*(y/x)··· pi

So I favor cordic for atan2 calculation (and that also calculates sqrt(x*x+y*y) ).

regards peter


·

Peter Verkaik

Tracy, here are some timing results
Trig uses tables and interpolation, cordic uses calculation

Test speed of trig functions
calculating Trig.sin(x) 100 times takes 107 milliseconds
calculating cordic.sin(x) 100 times takes 4662 milliseconds
calculating Trig.atan2(y,x) 100 times takes 4383 milliseconds
calculating cordic.atan2(y,x) 100 times takes 3561 milliseconds
program finished

So the table approach takes far less time for sin and cos, the atan2
takes less time for cordic calculation (on Javelin). I guess the
basic stamp does a bit better due to its builtin ** operator.

regards peter
·

Tracy Allen

Hi Peter,

That makes sense. Table lookup, even with interpolation, will always be faster than iterative calculation.

On the other hand, CORDIC has an edge speedwise on polynomials. CORDIC always improves precision on the order of one bit for each iteration, while polynomials (in general) require an increasing number of iterations for each digit of precision.

I'm not sure how SIN, COS, ATN and HYP functions are implemented in the Stamp, (table or CORDIC). They are limited precision, but knowing Chip Gracey, who has a deep appreciation of this stuff, I'd guess they are CORDIC.

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Tracy Allen

Hi Peter,

I updated the cordic-atn.bpe program for the BASIC Stamp to correct the bug you found with reflection into quandrant 2 & 3.
www.emesystems.com/BS2mathC.htm
Also
-- Added a wrapper for scaling small values up to the range of 8192 to 16384 (maintaining the ratio y:x) for better accuracy in the calculation, using Stamp NCD operator.
-- Changed the value used to represent 45 degrees from 16384 to 16200. That avoids the numerical overflow error that occurs in the angle accumulator if it passes the 32767 mark. Along with the scaling, that resolves the strange results you found when x=0.
-- Applied the CORDIC gain before going into the iteration. Applying the gain first scales down large numbers so it can work with larger numbers. Now it works okay for numbers up to at least x=23172 and y=23172. The cutoff is that the length of the vector hits 32768 and thus wraps into negative territory in twos complement.

** with regards

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Peter Verkaik

Thanks Tracy,
I noticed the loop executes for idx=0 to 15 but there are
17 entries in the atans table. Is that a typo or must the
loop execute for idx=0 to 16?

regards peter
·

Peter Verkaik

For those interested,
I found this article that deals with approximation of functions.
There are lots of examples and explanations.
regards peter

Tracy Allen

Right, there need only be 16 entries in the table. I pulled the values out of an excel spreadsheet and copied one extra, which is zero anyway, so it would not affect the result (?) even if the loop extended FOR 0 TO 16.

A lot of the residual error has to do with roundoff errors in the table entries for the binary arctangents. It occurs to me that a little more precision could be had by applying one more wrapper, to map all the angles into within +/- 45 degrees, and then use a larger number to represent 45 degrees. The wrapper would test x and y and exchange them if abs y > abs x, and after the iteration, adjust the angle by reflection in the 45 degree line and then reflection in the y axis, if the respective reflection flags are set. That in lieu of going to double precision.

I had trouble opening the "numeric analysis" file. Neither Acrobat nor Preview (mac) recognized it and the contents page opened as a graphic in bbedit. It looks interesting. Do you have another URL?

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Peter Verkaik

Tracy,

The attachement opens perfectly in adobe reader 7.0.5.
Good point about the 45 degrees. The polynome approximation
uses that also in case |y/x|>1.

regards peter

Tracy Allen

Okay, it does open fine with Acrobat reader 7. I was trying earlier with Acrobat 5 and with Preview.

It seems to be chapter 3 from
Peter Turner, _Guide to Numerical Analysis_
www.clarkson.edu/mcs/faculty/turner/
Did you find the whole PDF online somewhere?

▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
Tracy Allen
www.emesystems.com

Peter Verkaik

I found this after googling. Unfortenately I lost the exact url (it was turner),
but I am sure this was all. I found this the best resource I looked at because
of the eexamples and the explanations about different approaches.
I don't know if the book is still available.

regards peter