Hi to all who need more COGs
cessnapilot
Posts: 182
A small and complete Fixed-point arithmetic package (http://obex.parallax.com/objects/501/) with the necessary Trig in SPIN has just appeared in··Obex.
Enjoy,
Istvan····
Enjoy,
Istvan····
Comments
One comment, it would be worth adding a comment to the deg_atan2 routine that the parameters must be < 256. That caught me for a few minutes - you do a square of the 2 parameters and they overflow. And/or do a range check on the input parameters.
Thanks for the bug report. You are right. The ugly head of Fixed-point overflow rose again. The range is checked before the int atan2 for the product of x*y, but is not checked before the (later appended) Newton-Raphson refinement for x*x and y*y, separately. My fault, I am sorry.
In your code, as a quick treatment of the symptoms, you can Qdiv x and y with 100.0, for example, before the atan2, and take ten times the calculated distance (if needed).
The disease will be cured in the driver by downscaling·· x and y parallel for both ABS(x) and ABS(y) to be less that 128.0 (Qvalue). This will· affect the precision of the result only a very little. The debugged code will be uploaded within a week.
Numerically much better solution will be a QRadius(x, y) procedure, where, like in the Qmuldiv, 64-bit arithmetic is used for the (x*x+y*y) sum, and a 64-bit SQR is done. However, I do not know yet how to make 64-bit SQR in SPIN quickly. Maybe someone out there can help us.
Cheers,
Istvan
Post Edited (cessnapilot) : 9/1/2009 6:05:34 AM GMT
It was written by a mathematician while locked up
in a concentration camp. It is supposed to be a method of
doing advanced math quickly and all in your head.
I have been wondering if those methods might be useful
translated into pasm for the prop...I never actually read the
book, I have it here somewhere in all this mess.
The guys name who wrote that book was Trachtenberg.
Might be worth looking into.
▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔▔
- Some mornings I wake up cranky.....but usually I just let him sleep -
The code I am converting runs ~ 2x faster from floatmath to fixed int (i.e. the spin version of float), that contains 2 atan2, 10 mul and 6 div. With the 2x atan2 taking as much as the mul and divs.
Post Edited (Timmoore) : 9/1/2009 6:31:06 AM GMT
www.geocities.com/cnowlen/Cathy/Emat4680/Squareroot.htm
www.homeschoolmath.net/teaching/square-root-algorithm.php
humanoido
Place the 64-bit result of the operations
·· (X*X) + (Y*Y)
··
in the [noparse][[/noparse]Hi] and [noparse][[/noparse]Lo] 32-bit registers:
IF [noparse][[/noparse]Hi] is zero
·· RETURN :=·Qsqr([noparse][[/noparse]Lo])·· 'As a valid and precise Qvalue result of SQR
ELSE
· ·Take 32-bit SPIN SQR of [noparse][[/noparse]Hi] to get a 1st approximation of the integer part of SQR
·· [noparse][[/noparse]1st_Approx] := ^^[noparse][[/noparse]Hi]
···1st_Approx contains at least 1 (or usually more) decimal digit of the correct result , I guess. This is· an integer approximation of SQR(something x 2^32), so shift it to be a Qvalue Fixed-point number
·· [noparse][[/noparse]1st_Approx] := [noparse][[/noparse]1st_Approx << 16]
··
Unrolled REPEAT 3 follows
···Double the number of the correct digits by a· 64-bit / 32-bit division (similar to that in Qmuldiv) + averaging,
·· [noparse]/noparse]2nd_Approx] := ([noparse][[/noparse][noparse][[/noparse]Hi][noparse][[/noparse]Lo / [noparse][[/noparse]1st_Approx] + [noparse][[/noparse]1st_Approx]) / 2
·· Now we have at least 2 (or usually more) correct decimal digits. Double the number of the correct decimal digits again
·· [noparse]/noparse]3rd_Approx] := ([noparse][[/noparse][noparse][[/noparse]Hi][noparse][[/noparse]Lo / [noparse][[/noparse]2st_Approx] + [noparse][[/noparse]2st_Approx]) / 2
·· And again(to be sure of =>4 digit precision for both the integer and the fractional part of··the RESULT)
Precise ·64-bit Fixed-point SQR done. (As a first trial.)
To avoid unnecessary 64-bit divisions· new approximations will be tested to be close enough to the previous ones. If so Quit REPEAT loop
Cheers,
Istvan
Post Edited (cessnapilot) : 9/1/2009 4:39:40 PM GMT
1. They return 90degree different. I believe I compared mine with the float32full implementation and got the same results though it was a while ago and I could be mistaken. You might want to check.
2. Looking at the timing the other maths stuff mul/div I have used from your fixed implementation is 3-5 times faster than using floatmath. The atan2 implementations are the same speed both take ~20ms to run (on 80Mhz system). I am attaching my float implementation in case there is anything interesting
i added
if s
-result
at end of qmul
Couple of other comments, you can speed up atan2 (> 2x faster) with 2 changes
1.
The above code does the same thing with 2xmultiplies and 1div rather than 3x divides
2.
The RadtoDeg call in atan2 is only working with small degrees so you can use a qmul rather than the 64bit multiply in radtodeg since you will not get overflow in this case.
Post Edited (Timmoore) : 9/2/2009 4:45:28 AM GMT
Thanks for good tip for that Newton-Raphson step. Your version is not only faster, but more accurate, as well. SPIN_TrigPack is running with it now.
The missing 'sign restore' in Qmul is not a bug, or at least a harmless one, as we do not need sign check/restore at all with SPIN * and ** operators. I have left an unnecessary check before the operations. It's only the laziness of the coder (me). The results were good, so I didn't care enough to clean up. So, the Qmul now looks like
·The Qradius(X,Y) procedure goes as follows
PUB Qradius(qX, qY) : qR | hx, lx, hy, ly, h, l, cf, ap '------------------------------------------------------------------------- '--------------------------------┌─────────┐------------------------------ '--------------------------------│ QRadius │------------------------------ '--------------------------------└─────────┘------------------------------ '------------------------------------------------------------------------- '' Action: Calculates Distance of point (X,Y) from the origo '' Parameters: X, Y Descartes coordinates in Qs15_16 Fixed-point format '' Result: Distance from the origo in Qs15_16 Fixed-point format ''+Reads/Uses: None '' +Writes: None '' Calls: Add64, Div64 '' Note: QAngle(qX, qY) is encrypted as Deg_ATAN2(qX, qY) '' ToDo: Swap ap, qR for clarity '------------------------------------------------------------------------- hx := qX ** qX lx := qX * qX hy := qY ** qY ly := qY * qY Add64(hx, lx, hy, ly, @dQval) 'LONG dQval[noparse][[/noparse]2] global array hx := dQval[noparse][[/noparse]0] lx := dQval[noparse][[/noparse]1] 'Check for zero hx IF (hx == 0) RETURN (Qsqr(lx) >> 8) ELSE 'Prepare SQR iteration loop ap := (^^hx) << 16 'Do iteration 3 times REPEAT 3 h := hx l := lx 'Perform 64-bit division qR := Div64(hx, lx, ap) 'Calculate next approximation as ap ap := (qR + ap) >> 1 IF (qR == ap) QUIT RETURN ap '------------------------------------------------------------------------- PRI Add64(ah, al, bh, bl, res_) | la, lb '------------------------------------------------------------------------- '---------------------------------┌───────┐------------------------------- '---------------------------------│ Add64 │------------------------------- '---------------------------------└───────┘------------------------------- '------------------------------------------------------------------------- ' Action: Adds 64-bit numbers ' Parameters: - Hi and Lo part of Arg1 ' - Hi and Lo part of Arg2 ' - Addr of 64-bit dQval ' Result: None '+Reads/Uses: None ' +Writes: 64-bit result into LONG[noparse][[/noparse]address] ' Calls: None ' Note: Result is passed by reference '------------------------------------------------------------------------- LONG[noparse][[/noparse]res_][noparse][[/noparse]0] := ah + bh LONG[noparse][[/noparse]res_][noparse][[/noparse]1] := al + bl IF ((al < 0) AND (bl < 0)) LONG[noparse][[/noparse]res_][noparse][[/noparse]0]++ RETURN ELSEIF ((al => 0) AND (bl => 0)) RETURN ELSE la := al & $7FFF_FFFF lb := bl & $7FFF_FFFF IF ((la + lb) < 0) LONG[noparse][[/noparse]res_][noparse][[/noparse]0]++ RETURN '------------------------------------------------------------------------- PRI Div64(dh, dl, dr) : qQ | cf '------------------------------------------------------------------------- '---------------------------------┌───────┐------------------------------- '---------------------------------│ Div64 │------------------------------- '---------------------------------└───────┘------------------------------- '------------------------------------------------------------------------- ' Action: Divides a 64-bit dQvalue with q 32-bit Qvalue ' Parameters: - Hi and Lo part of dividend ' - Divisor ' Result: Quotient in Qs15_16 Fixed-point Qvalue format '+Reads/Uses: None ' +Writes: None ' Calls: None ' Note: None ' ToDo: Some optimization can be done for speed here '------------------------------------------------------------------------- qQ~ REPEAT 32 cf := dh < 0 dh := (dh << 1) + (dl >> 31) dl <<= 1 qQ <<= 1 IF ((dh > dr) OR cf) ++qQ dh -= dr RETURN qQ '-------------------------------------------------------------------------Qradius·gives for arguments
X = 10000.0
Y = 10000.0
>R = 14142.1356
Note the 9 digits of precision.
If someone has a better way to do 64-bit addition in SPIN, let us know. I feel my· method a little bit clumsy.
SPIN_TrigPack upgrade· on OBEX < 3 days.
Cheers,
Istvan
Post Edited (Timmoore) : 9/2/2009 7:56:03 PM GMT
Cheers,
Istvan
SPIN_TrigPack upgrade < 2 days
s is not set except if qV is -ve, if it happened to be 1 in memory then the result is -ve.
SPIN_TrigPack upgrade has been placed on OBEX. You can try the first Fixed-point package there. As a bonus this object now contains the first True Random Generator with the Propeller microcontroller using only SPIN. Sequences will repeat themselves only after the end of Times. Thanks to all who helped in the development.
Cheers,
Istvan