|
Sine
Aug 22, 2019 12:36:55 GMT 1
Post by scalion on Aug 22, 2019 12:36:55 GMT 1
The SinQ dont work very faster. I write this sine function that work 2 times faster.
// Object : Compute a sine more faster // Author : Nicolas Rey (scalion@free.fr) // Created : 22/08/2019 // For this example i allocate 524288 bytes to have 65536 double values // that say 2 * pi / 524288 ~= 0.000012 accuracy ! // It's more than enough to draw circle,ellipse or a spiral in this example... // Unfortunately the more memory we allocate the more the system is slow to respond // to read the data. This is due to the way windows handles memory. // We will have to limit ourselves to a maximum of 65536 values.
Global Double t, tt1 = 0, tt2 = 0, tt3, a, dr = 180 / PI Global Long aw, x, TestsCount = 0, AdrSin, SineValues = 65536, AdrSinSize = SineValues * 8
// Here we write the sines values : AdrSin = mAlloc(AdrSinSize) MemZero AdrSin, AdrSinSize For aw = AdrSin To AdrSin + AdrSinSize - 8 Step 8 DblPoke aw, Sin(2 * PI * (aw - AdrSin) / AdrSinSize) Next aw
FullW 1 : Void SetForegroundWindow(Me.hWnd) : AutoRedraw = 1 PrintScroll = 1 FontName = "lucida console" FontSize = 20 Do Inc TestsCount test1 : tt1 = tt1 + t 'sin test2 : tt2 = tt2 + t 'sinQ test3 : tt3 = tt3 + t 'sinD Print AT(1, 1); " sin = " & tt1 / TestsCount Print "sinQ = " & tt2 / TestsCount Print "sinD = " & tt3 / TestsCount For x = 0 To 10000 t = (1 - Rnd * Rnd) * 20 * PI Plot _X / 2 + SinD(t) * t * 2, _Y / 2 + CosD(t) * t * 2 Next x PeekEvent Loop Until Me Is Nothing Void mFree(AdrSin) Proc test1() Naked t = Timer For x = 0 To 1000000 a = Sin(x) Next x t = Timer - t EndProc Proc test2() Naked t = Timer For x = 0 To 1000000 a = SinQ(x * dr) Next x t = Timer - t EndProc Proc test3() Naked t = Timer For x = 0 To 1000000 a = SinD(x) Next x t = Timer - t EndProc Function SinD(ByVal a#) As Double Naked Static Long p2 = 2 * PI // P2 must be a double, it's a bug...(28/08/2019) If a < 0 a = p2 - Frac(Abs(a) / p2) * p2 Return DblPeek(Add(AdrSin, Shl(Int(Frac(a / p2) * SineValues), 3))) EndFunc Function CosD(ByVal a#) As Double Naked Static Long p2 = 2 * PI // Lol If a < 0 a = p2 - Frac(Abs(a) / p2) * p2 Return DblPeek(Add(AdrSin, Shl(Int(Frac(0.25 + a / p2) * SineValues), 3))) EndFunc
|
|
|
Sine
Aug 28, 2019 20:34:05 GMT 1
Post by dragonjim on Aug 28, 2019 20:34:05 GMT 1
Bonjour Nicolas,
Interesting post. The fact that Sin is now virtually equivalent to SinQ shows the leaps and bounds made in maths co-processors since GFA was first written - I know in early programs I wrote SinQ was very useful.
However, more interesting still is your own custom function which clearly outperform GFA's inbuilt Sin function. However, there seems to be some divergence in values output shown by adding the following code before or after the spiral drawn using SinD:
For x = 0 To 10000 t = (1 - Rnd * Rnd) * 20 * PI Plot _X * 3 / 4 + Sin(t) * t * 2, _Y / 2 + Cos(t) * t * 2 Next x
Please check my code (and maths which is very rusty) but it would be interesting to know why the two spirals differ. Maybe it is possible to tweak your functions to match? For information, GFA uses the output from the fsin IA-32 instruction - again, it is interesting that your code seems to outperform the CPU in terms of speed.
|
|
|
Post by scalion on Aug 28, 2019 21:02:22 GMT 1
Hi DragonJim, Omg, yes it's a serious problem !!! i take a look... i think it's a conceptual pb, not a precision pb, brb. However, I already noticed some differences between sinus calculated by code gfa basic and sine function of gfabasic. Try to code an algorytmh of sin by cutting arc and you will see.
LOL ! I declared Static Long p2 = 2 * PI in Sind and CosD Functions !!!!
P2 = 6 ( must be 6.28........)
Replace Long by Double. And all is good
Thank you DragonJim, BugBuster #1 !
I will post here a demo how to compute a sine.
|
|
|
Sine
Aug 28, 2019 21:36:09 GMT 1
Post by scalion on Aug 28, 2019 21:36:09 GMT 1
Global Double a Global Long OnVaPasYPasserLaJourneeNonPlus FullW 1 Void SetForegroundWindow(Me.hWnd) Me.AutoRedraw = 1 Me.PrintScroll = 1 a = -10 Do a = a + 0.1 Print "Angle = "; a Print " Sin(a)="; Sin(a); " by gfa" Print " Sin(a)="; MySin(a); " by MySin" KeyGet OnVaPasYPasserLaJourneeNonPlus PeekEvent Loop Until Me Is Nothing Function MySin(a#) As Double
Local Double s, sa, s1, s2, x1, y1, x2, y2, bx, by, h, spred Local Long iter = 0, SIGN If a < 0 a = Abs(a) : sa = -1 Else sa = 1 a = Frac(a / (2 * PI))
' 0 à PI/2 > 0 à 0.25 = 0 à 1 ' PI/2 à PI > 0.25 à 0.5 = 1 à 0 ' PI à PI*1.5 > 0.5 à 0.75 = 0 à -1 ' PI*1.5 à 2*PI > 0.75 à 1 = -1 à 0
If a <= 0.25 a = 4 * a SIGN = 1 * sa Else If a <= 0.5 a = 1 - (a - 0.25) * 4 SIGN = 1 * sa Else If a <= 0.75 a = (a - 0.5) * 4 SIGN = -1 * sa Else a = 1 - (a - 0.75) * 4 SIGN = -1 * sa EndIf
s1 = 0 s = 0.5 s2 = 1 x1 = 0 : y1 = 1 : x2 = 1 : y2 = 0 bx = Sqr(0.5) by = Sqr(0.5)
While s <> spred Inc iter spred = s If s > a s2 = s s = (s1 + s) / 2 x2 = bx : y2 = by bx = (x1 + bx) / 2 by = (y1 + by) / 2 Else If s < a s1 = s s = (s2 + s) / 2 x1 = bx : y1 = by bx = (x2 + bx) / 2 by = (y2 + by) / 2 Else Exit Do EndIf h = Sqr(bx * bx + by * by) bx = bx / h by = by / h Wend
Return SIGN * bx
EndFunc
|
|
|
Sine
Aug 29, 2019 14:53:23 GMT 1
Post by scalion on Aug 29, 2019 14:53:23 GMT 1
Improvement :
Sine values are the same on the 4 quarters of circle, just the sign changing.
We can have 4 times more precision :
Just change that :
Global Long aw, AdrSin, SineValues = 65536, AdrSinSize = SineValues * 8 Global Long SineValues4 = SineValues * 4
AdrSin = mAlloc(AdrSinSize) MemZero AdrSin, AdrSinSize For aw = AdrSin To AdrSin + AdrSinSize - 8 Step 8 DblPoke aw, Sin((PI / 2) * (aw - AdrSin) / (AdrSinSize - 8)) Next aw
Function SinD(ByVal a#) As Double Naked Static Double p2 = 2 * PI, ASGN If a < 0 a = 1 - Frac(Abs(a) / p2) Else a = Frac(a / p2) EndIf If a > 0.5 a = a - 0.5 ASGN = -1 Else ASGN = 1 EndIf If a > 0.25 a = 0.5 - a Return ASGN * DblPeek(Add(AdrSin, Shl(Int(a * SineValues4), 3))) EndFunc
Function CosD(ByVal a#) As Double Naked Return SinD(a + PI / 2) EndFunc
|
|
|
Sine
Aug 31, 2019 21:46:10 GMT 1
Post by scalion on Aug 31, 2019 21:46:10 GMT 1
In fact this implementation is slowest than Sin and SinQ.
Try to change data type of X in double... I was really surprised that it could be faster, especially in interpreted mode. Hence the importance of establishing a reliable test bench. What i messed up....
I am both disappointed and reassured.
Here, to make me forgive I give you a program of calculation on large numbers with x significant figures.
|
|
|
Sine
Sept 1, 2019 10:10:25 GMT 1
via mobile
scalion likes this
Post by dragonjim on Sept 1, 2019 10:10:25 GMT 1
Hi scalion,
I am away until Tuesday - hence my lack of response. Once back, I'll have a look at all your code examples and get back to you.
It's a shame if your amended code is slower - it may be possible to make it quicker by using a long static value of 628 then dividing the result of multiplication by 100 - this used to work years ago with old fashioned floating point calculations but haven't tried it recently.
|
|
|
Sine
Sept 1, 2019 18:08:19 GMT 1
Post by scalion on Sept 1, 2019 18:08:19 GMT 1
Hi DragonJim, Ok i write this. Yeah, that's very fastest. Thank's Here the new code (hope this time i dont messed up bench). I will post later a code with more accurate like SinD.
// Object : Compute a sine more faster // Author : Nicolas Rey (scalion@free.fr) // Created : 22/08/2019 // Modified : 01/09/2019 // The first test with SinD was in fact slowest // SinE use instead a long value (thank's to DragonJim) // Now effectively fastest.
Global Double x, t, tt1 = 0, tt2 = 0, tt3 = 0, tt4 = 0, a, dr = 180 / PI Global Long aw, TestsCount = 0, AdrSin, SineValues = 65536 /* must be 65536 for SinE */, AdrSinSize = SineValues * 8 Global Long SineValues4 = SineValues * 4 Global Long AdrSinE, AdrSinSizeE = SineValues * 8 Global Long SineValuesDivBy2Pi = SineValues / (2 * PI)
// Here we write the sines values for sinD : AdrSin = mAlloc(AdrSinSize) MemZero AdrSin, AdrSinSize For aw = AdrSin To AdrSin + AdrSinSize - 8 Step 8 DblPoke aw, Sin((PI / 2) * (aw - AdrSin) / (AdrSinSize - 8)) Next aw
// Sines values for sinE AdrSinE = mAlloc(AdrSinSizeE) MemZero AdrSinE, AdrSinSizeE For aw = AdrSinE To AdrSinE + AdrSinSizeE - 8 Step 8 DblPoke aw, Sin((PI * 2) * (aw - AdrSinE) / (AdrSinSizeE - 8)) Next aw
FullW 1 : Void SetForegroundWindow(Me.hWnd) : AutoRedraw = 1 PrintScroll = 1 FontName = "lucida console" FontSize = 20 Global Double aa = 0.5 Do For x = 0 To 10000 t = (1 - Rnd * Rnd) * 40 * aa * PI Plot _X / 2 + SinE(t) * t , _Y / 2 + CosE(t) * t Next x For x = 0 To 10000 t = (1 - Rnd * Rnd) * 40 * aa * PI Plot _X * 3 / 4 + Sin(t) * t , _Y / 2 + Cos(t) * t Next x Inc TestsCount test1 : tt1 = tt1 + t 'sin test2 : tt2 = tt2 + t 'sinQ test3 : tt3 = tt3 + t 'sinD test4 : tt4 = tt4 + t 'sinE Print AT(1, 1); " sin = " & tt1 / TestsCount Print "sinQ = " & tt2 / TestsCount Print "sinD = " & tt3 / TestsCount Print "sinE = " & tt4 / TestsCount aa = aa + 0.001 PeekEvent Loop Until Me Is Nothing Void mFree(AdrSin) Void mFree(AdrSinE) Proc test1() Naked t = Timer For x = 0 To 1000000 a = Sin(x) Next x t = Timer - t EndProc Proc test2() Naked t = Timer For x = 0 To 1000000 a = SinQ(x * dr) Next x t = Timer - t EndProc Proc test3() Naked t = Timer For x = 0 To 1000000 a = SinD(x) Next x t = Timer - t EndProc Proc test4() Naked t = Timer For x = 0 To 1000000 a = SinE(x) Next x t = Timer - t EndProc Function SinD(ByVal a#) As Double Naked Static Long p2 = 2 * PI, ASGN If a < 0 a = 1 - Frac(Abs(a) / p2) Else a = Frac(a / p2) EndIf If a > 0.5 a = a - 0.5 ASGN = -1 Else ASGN = 1 EndIf If a > 0.25 a = 0.5 - a Return ASGN * DblPeek(Add(AdrSin, Shl(Int(a * SineValues4), 3))) EndFunc Function CosD(ByVal a#) As Double Naked Return SinD(a + PI / 2) EndFunc Function SinE(ByVal a#) As Double Naked Static Long al If Abs(a) > 205894 a=0 // overflow prevent due to the long type of al (if i write return 0 the test is not really true) al = Abs(a) * SineValuesDivBy2Pi SinE = DblPeek(Add(AdrSinE, Shl(And(al, 65535), 3))) // al mod sinevalues = al mod 65536 = and(al,65535) EndFunc Function CosE(ByVal a#) As Double Naked Return SinE(a + PI / 2) EndFunc
|
|
|
Sine
Sept 1, 2019 21:31:14 GMT 1
Post by scalion on Sept 1, 2019 21:31:14 GMT 1
I writted the variation of SinE with a more accurate approximating, but the problem is the type long limit the angle range....
And if i insert a frac(a/p)*p we return to the previous problem and cancel the advantage of long type... Aaaargh...
That mean i have 2 choices :
- Less precision and a wider range of values (the first SinE) - More precision and a smaller range of values (this version)
Well, Here the code of this :
Function SinE(ByVal a#) As Double Naked Static Long al Static Double p = 16384 * PI // not used here... If Abs(a) > 51473 a = 0 ' overflow prevent due to the type long of al (16384*pi) al = Abs(a) * SineValuesDivBy2PiMulBy4 And al, 262143 If al > 196607 Return -DblPeek(Add(AdrSin, Shl(Sub(262143 , al), 3))) If al > 131071 Return -DblPeek(Add(AdrSin, Shl(Sub(al, 131071), 3))) If al > 65535 Return DblPeek(Add(AdrSin, Shl(Sub(131071, al), 3))) SinE = DblPeek(Add(AdrSin, Shl(al, 3))) EndFunc
|
|
|
Sine
Sept 2, 2019 8:10:50 GMT 1
Post by scalion on Sept 2, 2019 8:10:50 GMT 1
Hi JMM, Thank you, your demo use a long type parameter (i%) that can't be test on the same bench of Sin() and SinQ(). Of coure they are faster that is the most intersting. However... in an application who use large ammount of rotating (for example a flight simulation) it can be intersting to use only angle in integer degree, a multiple of 360 like 36000. Then we have just make a modulo to adress the correct sine. But in fact, i think we must use not degree or radian but tr . Where 1 tr = 360 degree = 2*pi radian. Note tr=tour. Then with tr unity we can simply use a 65536 modulo, that's fast with logical And(tr,65535). But at this time there another problem : opengl use degree and not tr.
Then we have to choose between integer multiple of degree or tr with a degree multiplicator = int(65536/360) = 182.
|
|
|
Sine
Sept 2, 2019 14:51:35 GMT 1
Post by scalion on Sept 2, 2019 14:51:35 GMT 1
Yeah, Good Job JMM ! i want to say in french CQFD (perraps QED in inglish..).
|
|
|
Sine
Sept 2, 2019 15:37:47 GMT 1
Post by scalion on Sept 2, 2019 15:37:47 GMT 1
Well now attack the precision question...
|
|
|
Sine
Sept 2, 2019 15:39:05 GMT 1
Post by scalion on Sept 2, 2019 15:39:05 GMT 1
Now after 1.3 billiard (milliard in french) tests the precision i believe i can say :
Computing sines by basic code is 1.67 times precisest of the sin() function for values < 65536*pi and 1.08 times precisest (we can say equivalant) for big values (For the moment with my method).
For my test i used the 16 know values of sine and their correspondant radian, stored in this arrays :
Global Double KnowRadian(1 To 16) Global Double KnowSine(1 To 16)
KnowRadian(1) = 0 : KnowSine(1) = 0 KnowRadian(2) = PI / 6 : KnowSine(2) = 0.5 KnowRadian(3) = PI / 4 : KnowSine(3) = Sqr(0.5) KnowRadian(4) = PI / 3 : KnowSine(4) = Sqr(0.75) KnowRadian(5) = PI / 2 : KnowSine(5) = 1 KnowRadian(6) = 2 * PI / 3 : KnowSine(6) = Sqr(0.75) KnowRadian(7) = 3 * PI / 4 : KnowSine(7) = Sqr(0.5) KnowRadian(8) = 5 * PI / 6 : KnowSine(8) = 0.5 KnowRadian(9) = PI : KnowSine(9) = 0 KnowRadian(10) = 7 * PI / 6 : KnowSine(10) = -0.5 KnowRadian(11) = 5 * PI / 4 : KnowSine(11) = -Sqr(0.5) KnowRadian(12) = 4 * PI / 3 : KnowSine(12) = -Sqr(0.75) KnowRadian(13) = 3 * PI / 2 : KnowSine(13) = -1 KnowRadian(14) = 5 * PI / 3 : KnowSine(14) = -Sqr(0.75) KnowRadian(15) = 7 * PI / 4 : KnowSine(15) = -Sqr(0.5) KnowRadian(16) = 11 * PI / 6 : KnowSine(16) = -0.5 Next, i choose a random know value in this 16, and add a multiple of 2*pi to the know radian value, the sine have not to change in this case.
After calculating the sine of this radian i compare whith the expected sine value.
When my coded sine is nearest to the expected i increase a counter for it, same for the sin() function. When they're equivalant i do nothing. That's what giving the graphical result below : - At left the Sin(), at right the Coded sine SinDA() function. - Lower to greater alpha values, inside to outside. - Pixels positionned by Sin() and SinDA() and their respective Cosinus function.
Red Composant indicate where precision is lowest than other sin function.
Green Composant indicate where precision is better than other sin function.
Blue Composant indicate where precision is low.
You can test yourself with the complete source code :
And now I wonder someone here would be able to make me a sine function that would give a completely green graph! (or almost) <-ninja in ambuch...
|
|
|
Post by dragonjim on Sept 3, 2019 11:38:37 GMT 1
Hi all, I like the way this post is going - there are some intriguing ideas and some interesting maths. scalion : I have finally had a chance to have a quick look through your code examples and the results of your long integer-based SinE are very good - to test I changed the plotting routine to this: For x = 0 To 10000 t = (1 - Rnd * Rnd) * 40 * aa * PI Color RGB(255, 0, 0) ' <<---- Plot using SinE in red Plot _X / 2 + SinE(t) * t , _Y / 2 + CosE(t) * t Next x For x = 0 To 10000 t = (1 - Rnd * Rnd) * 40 * aa * PI Color RGB(0, 255, 0) ' <<---- Plot using Sin in green... Plot _X / 2 + Sin(t) * t , _Y / 2 + Cos(t) * t ' <<---- ...imposed over SinE... Next x Pause 10 ' <<---- ...and pause to see the effect I noticed some red peeking through the green mask, which seemed to show that there were some inaccuracies in the SinE calculations - compared to the GB Sin function and, thus, the IA-32 fsin instruction, at least. Then I tried a control and changed SinE and CosE to Sin and Cos to ensure that every other aspect of the loops were working as they should - the code looked like this: For x = 0 To 10000 t = (1 - Rnd * Rnd) * 40 * aa * PI Color RGB(255, 0, 0) ' <<---- Plot using Sin first in red... Plot _X / 2 + Sin(t) * t , _Y / 2 + Cos(t) * t Next x For x = 0 To 10000 t = (1 - Rnd * Rnd) * 40 * aa * PI Color RGB(0, 255, 0) ' <<---- ...and then in green... Plot _X / 2 + Sin(t) * t , _Y / 2 + Cos(t) * t ' <<---- ...imposed over first Sin spiral... Next x Pause 10 ' <<---- ...and pause to see the effect And I noticed the same red dots peeking through the green mask. I then drew both spirals using SinE and CosE and similar red dots were peeking through the green mask. So I came to the conclusion that most inaccuracies are the result of the Plot command rather than major divergences in your SinE/CosE functions. So, once again, well done on the speed increase - it was slightly faster when compiled but not by any great amount.
|
|
|
Sine
Sept 3, 2019 12:40:16 GMT 1
Post by scalion on Sept 3, 2019 12:40:16 GMT 1
For the F_Sin function of JMM i have good results too. 1.35 times precisest on low values to 1.08 for the big values (same of SinDA). Just i have a question : why using a _Pi variable instead the PI constant ?
|
|