New ATAN approximation in Code snippet...

AllyCat

Senior Member
Hi,

Slow, but accuracy better than 0,001° :
Code:
....
[color=Green]' Y=ATAN(X)   45°=45000  [b]ATN(45°)=10000[/b][/color]
[color=Black]ATAN1:
      [/color][color=Purple]I[/color][color=DarkCyan]=[/color][color=Purple]X[/color][color=DarkCyan]/[/color][color=Navy]100[/color][color=DarkCyan]*[/color][color=Navy]2
      [/color][color=Blue]read [/color][color=Purple]I[/color][color=Black],[/color][color=Blue]WORD [/color][color=Purple]T1
      i[/color][color=DarkCyan]=[/color][color=Purple]i[/color][color=DarkCyan]+[/color][color=Navy]2
      [/color][color=Blue]read [/color][color=Purple]I[/color][color=Black],[/color][color=Blue]WORD [/color][color=Purple]T2
      T2[/color][color=DarkCyan]=[/color][color=Purple]T2[/color][color=DarkCyan]-[/color][color=Purple]T1
      Y[/color][color=DarkCyan]=[/color][color=Purple]X[/color][color=DarkCyan]//[/color][color=Navy]100[/color][color=DarkCyan]*[/color][color=Purple]T2[/color][color=DarkCyan]/[/color][color=Navy]100[/color][color=DarkCyan]+[/color][color=Purple]T1[/color][color=DarkCyan]+[/color][color=Navy]1[/color]
[color=Blue]RETURN[/color]
I think you mean TAN(45°) = 10000. The program code seems to use a simple 100+ points (200+ bytes) lookup table with 100 point linear interpolation between each pair of adjacent reference points.

So I thought it would be interesting to compare your results with the similar interpolation routine for ATAN2 that I posted in #18 of this fairly recent thread. My "General Purpose" interpolation routine is itself documented in this code snippet and uses about 60 bytes of program code compared with about 50 for the above

The initial results from my code were much less "accurate" for various reasons, particularly because I used only 9 reference points (8 lines) for the interpolation, and a result (Y) range of 360 degrees (for full ATAN2). So I increased my lookup table to 33 points (66 bytes) to give 8192 "input" values (256 x 32) compared with the 10000 for the above (100 x 100). However, when calculating the lookup table values I noticed that even 32 points are rather "overkill" because the TAN function (from 0 to 45 degrees) is not very curved, so the "worst fit" (at the centre and ends of each interpolation line) is much less than 1/2 of the Least Significant Bit (input).

But my results are still less "accurate" than yours in the final digit, but is the accuracy really 0.001 (degree)? Yes, if you use only exact integer input (X) values, but an "analogue" input value of (e.g.) 0.00015 must be rounded to either 0.0001 or 0.0002, which have ATAN values of 0.006 and 0.012 degrees, so the "error" must be at least 0.003 degree ? Thus it is mainly the scaling (or mapping) of my 8192 points to 10000 (i.e. 1.0000) which introduces the main "errors", together with incomplete rounding of my intermediate integer values. But in a PICaxe application is the input (X) value not more likely to be scaled from a binary value (e.g. READADC10) than a decimal value?

For comparison, here is my adapted test code:

Code:
#simspeed 10
symbol X=w9
symbol Y=w10
symbol T1=w11
symbol T2=w12
symbol I=b11

symbol tempw = w1		; Used explicitly
symbol wx = w2
symbol wy = w3
symbol atanbase = 24		; 0 for 8 segment table, 24 for 32 segments

data 0,(0,0,17,5,251,9,158,14,228,18,194,22,56,26,74,29,255,31)      ; 8 segments, 8192 = unity Y
data 24,(0,0,23,5,44,10,60,15,68,20,67,25,53,30,25,35)		      ; 32 segments, 32768 = unity Y
data(237,39,175,44,93,49,246,53,121,58,228,62,55,67,112,71)
data(144,75,150,79,130,83,83,87,10,91,166,94,41,98,145,101)
data(224,104,22,108,52,111,57,114,39,117,254,119,191,122,106,125,255,127)

table (0,0,60,2,121,4,182,6,242,8,46,11,105,13,164,15,221,17,22,20,78,22,133,24,186,26,238,28,33,31,82,33,130,35,176,37,219,39,5,42,45,44,83,46,119,48,152,50,183,52,212,54,238,56,5,59,26,61,44,63,59,65,71,67,80,69,86,71,90,73,90,75,86,77,80,79,70,81,57,83,41,85,21,87,254,88,227,90,197,92,163,94,126,96,85,98,41,100,248,101,197,103,141,105,82,107,19,109,209,110,138,112,64,114,243,115,161,117,76,119,243,120,151,122,54,124,210,125,107,127,255,128,144,130,30,132,167,133,45,135,176,136,46,138,169,139,33,141,149,142,5,144,114,145,220,146,66,148,164,149,3,151,95,152,183,153,12,155,94,156,172,157,247,158,63,160,131,161,197,162,3,164,62,165,118,166,170,167,220,168,11,170,54,171,95,172,133,173,168,174,200,175,229,176)

sertxd("----",13,10)
for x=0 to 100 step 10   	' ATAN from 0.0000 to 0.0100
	gosub ATAN1
	sertxd ("ATAN(",#X,")=",#Y," ",#w1,13,10)
next X

sertxd(".....",13,10)
for x=4060 to 4070   		' ATAN from 0.4060 to 0.4070
	gosub ATAN1
	sertxd ("ATAN(",#X,")=",#Y," ",#w1,13,10)
next X

sertxd(".....",13,10)
for x=9980 to 10001   		' ATAN from 0.9980 to 1.0000
	gosub ATAN1
	sertxd ("ATAN(",#X,")=",#Y," ",#w1,13,10)
next X
end 

' Y=ATAN(X)   45°=45000=Y  TAN(45°)=10000=X
ATAN1:							; [50 bytes]
      I=X/100*2
      readtable I,WORD T1 
      i=i+2
      readtable I,WORD T2
      T2=T2-T1
      Y=X//100*T2/100+T1+1
; RETURN                              ; Fall Into ATAN2
ATAN2:  						; [74 bytes with scaling]
	w1 = x ** 53687				; X scale 10000 -> 8192 
interpXtoY:	 					; [61 bytes]
	wx = w1					; Frees up b2 and b3
	b2 = atanbase				; Address of start (word) of table
	b2 = b2 + b3 + b3			; Add 2* wx/256 (b3 still contains hi byte of w1) = base of cell
	b3 = b2 + 2					; Point to top of cell
	read b3,WORD wy
	b3 = b2 + 1					; Point to high byte
	read b2,b2					; Read low byte
	read b3,b3					; Read high byte
	wy = wy - w1 				; Calculate height of cell (wy - wx)
	wx = wx / 256 * 65280 + wx		; Horizontal distance across cell (65280 = -1 * 256)
	w1 = wx * 256 ** wy + w1		; Height inside cell (0<wx<256) * gradient + height of cell base
;	w1 = w1 * 2 ** 45000			; Y scale up to 45 degrees (simpler formula)
	w1 = w1 ** 24464 + w1		; 32768 -> 45000 degrees [24464 = ((45000/32768)-1)*65536
return
Cheers, Alan.
 

BESQUEUT

Senior Member
I think you mean TAN(45°) = 10000. The program code seems to use a simple 100+ points (200+ bytes) lookup table with 100 point linear interpolation between each pair of adjacent reference points....
Cheers, Alan.
Thank you for your interest !
I must apologize : of course TAN(45°) = 10000
I will look at all that next week, but at first reading your are true for all...
Thanks
 

AllyCat

Senior Member
Hi,

I've now updated my Code Snippet thread (Warning, it's LONG) that I linked above, to show functions for SIN, COS, TAN, ARCSIN, ARCCOS, ARCTAN and ATAN2. ;)

The functions can be used with any (recent) PICaxe and should give at least 100 times better "accuracy" (resolution) than the embedded X2 functions. They work for all four quadrants (0 - 360 degrees), internally using full 16-bit signed (twos complement) words (i.e. 15 bits + sign), with a choice of overall accuracy determined by choosing the number of points (joined by line segments) in each lookup table.

Internally, angles are represented as "binary fractions" of a full circle (360 degrees or 2 * PI Radians) but can be considered as the proportion of a full circle multiplied by 65536. Used with the ** operator, they can be easily converted to degrees {or hundredths} by ** 360{00} or to Radians*10,000 (as used by many maths programs such as Excel), by **62832 (i.e. 2 * PI * 10,000), or even Compass directions (N, S, E, W and all valid permutations) directly from the Most Significant Bits of the word (e.g. bit15 = N/S, bit14 = E/W, etc.).

Similarly, the internal numbers are handled as "signed binary fractions" (i.e. 32768 is equivalent to Unity) so can be converted to (e.g.) four decimal places by ** 20000 , etc.. Generally, the results will be accurate to about the 4th significant digit, but one must be careful with non-linear expressions such as trigonometric functions.


In passing, I did notice a few interesting features of the embedded (X2) trig functions. Angles are in degrees, but the maximum value is specified as 65535, i.e. approximately (but not exactly) 182 revolutions. Of course that simply means that the algorithm does a (modulus) // 360 before the main calculation. The results are in hundredths of unity (i.e. to two decimal places), and have a "sign" flag, but are NOT twos complement, so any further calculations need to be performed with care. Also, the resolution appears to be in increments of two degrees (and are not rounded to the nearest value) so "odd" angles (such as 15 degrees) can be in error by 2 digits.

Fundamentally, the TAN function is "nasty" because it becomes infinite (and negatively infinite as well) at 90 degrees, due to division by zero. So the X2 ATAN function makes a fair stab with 5729 (57.28) giving 89 degrees (actually TAN 89 degrees = 57.28996), but it gives only 88 degrees for 5728. Also, I was amused and surprised that my "Demo" code for TAN 90 degrees (in #4 of the Code Snippet) gets as close to infinity as it can with "32767.9999". ;)


Finally, an example of how careful one needs to be with non-linear functions. I was quite happy that the ARCSIN function for "+1" produced a value of 89.98 degrees (it should be 90 degrees of course) when using only a 16-segment table, but surprised that the 64-segment calculation produces a "worse" value around 89.8 degrees. The reason is that any SINE above about 89.7 degrees should be reported as "1.0000" in a 16-bit signed calculation, if rounded to the nearest integer. The value would be even lower (at about 89.4 degrees) for a calculation rounded to 4 decimal places. Therefore, the "reverse search" for ARCSIN first finds a value of around 89.8 and reports that.

It might be possible to produce a "cosmetically" more pleasing result of 90 degrees by searching in the opposite direction (i.e. through a "mirrored" table starting at 90 degrees) and/or by a more rigorous handling of bit rounding and maximum values, etc.. But this is all largely of "academic interest" because any real benefits will probably be lost in any previous or subsequent PICaxe Basic calculations, even if performed to full 16-bit resolution.

Cheers, Alan.
 
Top