The reason I introduced my Loopy version was to treat and hopefully cure Albert of infinite loops. My Loopy version is just CORDIC, pretty much the same technique as used by the Intel FPU for the FSINCOS assembly instruction.
I think the aim of this topic is to find an approximation for Sin and Cos that is significantly faster, though less accurate, than the internal FB Sin and Cos functions. Since the application was for drawing circular arcs on the screen the aim should be to generate the fastest code that will give a maximum error less than pixel size, about 1 in 2000. If the code is too accurate then speed it up. All scatter points must be within the inner circle with mag = 100. That is a 1 in 1000 accuracy absolute limit.
@Muttonhead. If I remember rightly, the Bezier curve is a parametric curve that generates x,y pairs as you step through a parameter t. Now the thing with a parametric curve is that the path taken in x,y is well defined, but the velocity of the point on the graph as you step through parameter t is not. I believe that velocity sets the phase angle error while the accuracy of the x,y path sets the radius error. The diagonal line of Bezier scatter points shows that your phase is accurate but amplitude is more variable. This is counter intuitive and surprises me. It seems that with the Bezier curve, phase error may not be as big a problem as radius error.
The use of a "5th-order-bezier" is too accurate and therefore slow. There are too many references to variables. Rather than continuous improvement of accuracy, what is the fastest Bezier code that meets the 1 in 2000 specification. Can it be done with a cubic?
Circles
@Richard
Take a look at this thread...
http://www.freebasic.net/forum/viewtopic.php?t=14783
I posted (and MichaelW optimized) a pretty fast SQRT function that can be used to find SIN ( ie: FSQRT(1-x*x) ).
-Vince
Take a look at this thread...
http://www.freebasic.net/forum/viewtopic.php?t=14783
I posted (and MichaelW optimized) a pretty fast SQRT function that can be used to find SIN ( ie: FSQRT(1-x*x) ).
-Vince
Thanks Richard the s=y/r thingy was what i was trying to get to but couldn't figure it out.
I was trying to get the difference and split it to the x and y ratios.
but the sqr() makes a bigger number rather than a smaller number so i had to use the raw diff and put it into a loop.
I tested it to 1,000,000,000 mag and its still on the circle.
now the problem is trying to get the degrees right:
the way the degree lines cross the hypotenuse are not equal divisions so my formula to simply mul the grad by a set division doesn't exactly work.
Any formulas for that Richard? ,or anyone else.
I was trying to get the difference and split it to the x and y ratios.
but the sqr() makes a bigger number rather than a smaller number so i had to use the raw diff and put it into a loop.
I tested it to 1,000,000,000 mag and its still on the circle.
now the problem is trying to get the degrees right:
the way the degree lines cross the hypotenuse are not equal divisions so my formula to simply mul the grad by a set division doesn't exactly work.
Any formulas for that Richard? ,or anyone else.
My revised formula, thanks to Richard.
You can set the mag to 100,000,000,000 and the yellow dots are still perfectly on the circle.
Now I have to fix the degree errors.
Richard; If you just use x and y the formula makes a perfect octogon and the degrees seem evenly spaced at .0142......
What i was trying to accomplish was to finish out the octogon to the radius.
using the s=y/r destroys the primary octogon and messes up the degree spacing.
Is there a way to just finish the octogon to the radius?????
You can set the mag to 100,000,000,000 and the yellow dots are still perfectly on the circle.
Now I have to fix the degree errors.
Code: Select all
'================================================================
' Albert's(yellow), and Richard's(pink)
' approximations measured relative to the reference circle ... blue
Const As Double mag = 100 ' magnification of errors
'================================================================
Dim As Double t, tt, grad, q, s, c, dx, dy, x, y, r = 1.1
Screen 19
Dim As Integer w, h, depth
Screeninfo w, h, depth
Window (-r * w / h, -r) - (r * w / h, r) ' maintain pixel aspect
'----------------------------------------------------------------
Line(-c, 0)-(c, 0), 1 ' x axis
Line(0, -c)-(0, c), 1 ' y axis
Const As Double Pion2 = 2 * Atn(1)
For t = 0 To 4 * Pion2 Step Pion2 / 1000
x = Cos(t)
y = Sin(t)
Pset(x, y), 1 ' reference circle
Pset(x*.1, y*.1), 1 ' one tenth of circle at centre
Next t
'================================================================
' Albert's formula in yellow
Dim As Double sqr_root
For grad = 0 To 100 Step 1
'----------------------------------------
If grad <= 50 Then
x= ((50-grad) * (1 - Sqr(.5)) / 50) + Sqr(.5)
y= grad * (Sqr(.5) / 50)
Else
x= (100-grad) * (Sqr(.5) / 50)
y= ((grad-50) * (1 - Sqr(.5)) / 50) + Sqr(.5)
End If
sqr_root = Sqr(x*x + y*y)
s = y / sqr_root
c = x / sqr_root
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 14
Pset (+x, +y), 14
Pset (-x, -y), 14
Pset (-x, +y), 14
t = grad * Pion2 / 100
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 14 ' plot the scatter
Next grad
'================================================================
' Richard's constant parameter version using radians
Const As Double a = -0.166 ' given a, then b is closure parameter
Const As Double b = .00759 ' b = Pion2^-5 - Pion2^-4 - a*Pion2^-2
For t = 0 To Pion2 Step pion2 / 100
'----------------------------------------
tt = t * t
s = t * ( 1 + tt * (tt * b + a))
c = Pion2 - t
tt = c * c
c = c * ( 1 + tt * (tt * b + a))
'----------------------------------------
' operation count 8* 4+ 1- 5c 16v
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 13
Pset (+x, +y), 13
Pset (-x, -y), 13
Pset (-x, +y), 13
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 13 ' plot the scatter
Next t
'================================================================
Sleep
'================================================================
What i was trying to accomplish was to finish out the octogon to the radius.
using the s=y/r destroys the primary octogon and messes up the degree spacing.
Is there a way to just finish the octogon to the radius?????
Richard; Here are the octogon and the circle side by side.
The octogon doesn't have the gaps that the circle has.
Is there a way to correct the circle formula to get the same spacing on the degrees as the octogon?
The octogon doesn't have the gaps that the circle has.
Is there a way to correct the circle formula to get the same spacing on the degrees as the octogon?
Code: Select all
dim as double xctr,yctr,radius
screen 19
xctr=205
yctr=305
radius = 200
circle(xctr,yctr),radius,12
dim as double x,y,grad,s,c,sqr_root
For grad = 0 To 100 step 1
'----------------------------------------
If grad <= 50 Then
x= ((50-grad) * (1 - Sqr(.5))/ 50) + Sqr(.5)
y= grad * (Sqr(.5) / 50)
Else
x= (100-grad) * (Sqr(.5) / 50)
y= ((grad-50) * (1 - Sqr(.5)) / 50) + Sqr(.5)
End If
sqr_root = Sqr(x*x + y*y)
s = (y / sqr_root)
c = (x / sqr_root)
x = c * radius
y = s * radius
pset(xctr+x,yctr+-y),14
pset(xctr+x,yctr+ y),14
pset(xctr-x,yctr+-y),14
pset(xctr-x,yctr+ y),14
line(xctr,yctr)-(xctr + x, yctr + -y),14
line(xctr,yctr)-(xctr + x, yctr + y),14
line(xctr,yctr)-(xctr - x, yctr + -y),14
line(xctr,yctr)-(xctr - x, yctr + y),14
'print deg;" - ";s,c
'print atan2(s,c)*100/(atn(1)*2)
next grad
xctr=605
yctr=305
radius = 200
circle(xctr,yctr),radius,12
'dim as double x,y,grad,s,c,sqr_root
For grad = 0 To 100 step 1
'----------------------------------------
If grad <= 50 Then
x= ((50-grad) * (1 - Sqr(.5))/ 50) + Sqr(.5)
y= grad * (Sqr(.5) / 50)
Else
x= (100-grad) * (Sqr(.5) / 50)
y= ((grad-50) * (1 - Sqr(.5)) / 50) + Sqr(.5)
End If
x = x * radius
y = y * radius
pset(xctr+x,yctr+-y),14
pset(xctr+x,yctr+ y),14
pset(xctr-x,yctr+-y),14
pset(xctr-x,yctr+ y),14
line(xctr,yctr)-(xctr + x, yctr + -y),14
line(xctr,yctr)-(xctr + x, yctr + y),14
line(xctr,yctr)-(xctr - x, yctr + -y),14
line(xctr,yctr)-(xctr - x, yctr + y),14
'print deg;" - ";s,c
'print atan2(s,c)*100/(atn(1)*2)
next grad
sleep
@Albert. Your first approximations are straight lines that touch the circle at 0, 45 and 90 degrees, sorry that should be 0, 50 and 100 grad. They are straight because there is only one grad term in each of the equations, which limits them to straight lines.
If grad <= 50 Then
x= ((50-grad) * (1 - Sqr(.5))/ 50) + Sqr(.5)
y= grad * (Sqr(.5) / 50)
Else
x= (100-grad) * (Sqr(.5) / 50)
y= ((grad-50) * (1 - Sqr(.5)) / 50) + Sqr(.5)
End If
The problem with phase angle is that your angle is not being measured around a circle but along a straight chord. Your first approximation needs to be curved similar to the sine or cosine function in order to give a more linear initial phase angle. The cosine approximation will have the form of 1 – a*x*x.
If grad <= 50 Then
x= ((50-grad) * (1 - Sqr(.5))/ 50) + Sqr(.5)
y= grad * (Sqr(.5) / 50)
Else
x= (100-grad) * (Sqr(.5) / 50)
y= ((grad-50) * (1 - Sqr(.5)) / 50) + Sqr(.5)
End If
The problem with phase angle is that your angle is not being measured around a circle but along a straight chord. Your first approximation needs to be curved similar to the sine or cosine function in order to give a more linear initial phase angle. The cosine approximation will have the form of 1 – a*x*x.
Thanks for all your input and help Richard you've been a great help.
My newest formula (LOOK AT THIS PERFECT CIRCLE!!!!!).
I think sines and cosines are going to be obsoleteed soon!
My newest formula (LOOK AT THIS PERFECT CIRCLE!!!!!).
I think sines and cosines are going to be obsoleteed soon!
Code: Select all
dim as double xctr,yctr,radius
screen 19
xctr=305
yctr=305
radius = 270
circle(xctr,yctr),radius,12
dim as double x,y,grad,s,c,sqr_root
For grad = 0 To 100 step 1
'----------------------------------------
If grad <= 50 Then
'x= ((50-grad) * (1 - Sqr(.5))/ 50) + Sqr(.5)
'y= grad * (Sqr(.5) / 50)
x=(100-grad+(.5*grad))/100
y=(grad+(.5*grad))/100
Else
'x= (100-grad) * (Sqr(.5) / 50)
'y= ((grad-50) * (1 - Sqr(.5)) / 50) + Sqr(.5)
x=(100-grad+(.5*(100-grad)))/100
y=(grad+(.5*(100-grad)))/100
End If
sqr_root = Sqr(x*x + y*y)
s = (y / sqr_root)
c = (x / sqr_root)
x = c * radius
y = s * radius
pset(xctr+x,yctr+-y),14
pset(xctr+x,yctr+ y),14
pset(xctr-x,yctr+-y),14
pset(xctr-x,yctr+ y),14
line(xctr,yctr)-(xctr + x, yctr + -y),14
line(xctr,yctr)-(xctr + x, yctr + y),14
line(xctr,yctr)-(xctr - x, yctr + -y),14
line(xctr,yctr)-(xctr - x, yctr + y),14
'print deg;" - ";s,c
'print atan2(s,c)*100/(atn(1)*2)
next grad
sleep
@Albert. It does not matter how round your circle looks if the points are in the wrong place on the circumference. You cannot tell how small the phase error is by looking at a circle of pixels on the screen. It is an improvement, but you still have significant phase errors. Look at a scatter plot of your approximation with mag = 50 and you will see why sine and cosine are going to be around for quite some time yet.
Richard the s=y/r ; thingy works with yours too, here is your with 10,000 mag.
Can you dial it in a little,to get the degrees closer?
Can you dial it in a little,to get the degrees closer?
Code: Select all
'================================================================
' Richard's(pink)
' approximations measured relative to the reference circle ... blue
Const As Double mag = 10000 ' magnification of errors
'================================================================
Dim As Double t, tt, grad, q, s, c, dx, dy, x, y, r = 1.1
Screen 19
Dim As Integer w, h, depth
Screeninfo w, h, depth
Window (-r * w / h, -r) - (r * w / h, r) ' maintain pixel aspect
'----------------------------------------------------------------
Line(-c, 0)-(c, 0), 1 ' x axis
Line(0, -c)-(0, c), 1 ' y axis
Const As Double Pion2 = 2 * Atn(1)
For t = 0 To 4 * Pion2 Step Pion2 / 1000
x = Cos(t)
y = Sin(t)
Pset(x, y), 1 ' reference circle
Pset(x*.1, y*.1), 1 ' one tenth of circle at centre
Next t
'================================================================
' Richard's constant parameter version using radians
dim as double sqr_root
Const As Double a = -0.166 ' given a, then b is closure parameter
Const As Double b = .00759 ' b = Pion2^-5 - Pion2^-4 - a*Pion2^-2
For t = 0 To Pion2 Step pion2 / 100
'----------------------------------------
tt = t * t
s = t * ( 1 + tt * (tt * b + a))
c = Pion2 - t
tt = c * c
c = c * ( 1 + tt * (tt * b + a))
sqr_root = sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
' operation count 8* 4+ 1- 5c 16v
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 13
Pset (+x, +y), 13
Pset (-x, -y), 13
Pset (-x, +y), 13
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 13 ' plot the scatter
Next t
locate 1,1:print "Magnification = ";mag
'================================================================
Sleep
'================================================================
I got mine inside the circle at 100 mag. finally!
Its smaller than the little circle.
I applied the s-c=s-c / sqr(c*c+s*s) to both Richards and mine.
Its smaller than the little circle.
I applied the s-c=s-c / sqr(c*c+s*s) to both Richards and mine.
Code: Select all
'================================================================
' Albert's(yellow), and Richard's(pink)
' approximations measured relative to the reference circle ... blue
Const As Double mag = 100 ' magnification of errors
'================================================================
Dim As Double t, tt, grad, q, s, c, dx, dy, x, y, r = 1.1
Dim As Double sqr_root
Screen 19
Dim As Integer w, h, depth
Screeninfo w, h, depth
Window (-r * w / h, -r) - (r * w / h, r) ' maintain pixel aspect
'----------------------------------------------------------------
Line(-c, 0)-(c, 0), 1 ' x axis
Line(0, -c)-(0, c), 1 ' y axis
Const As Double Pion2 = 2 * Atn(1)
For t = 0 To 4 * Pion2 Step Pion2 / 1000
x = Cos(t)
y = Sin(t)
Pset(x, y), 1 ' reference circle
Pset(x*.1, y*.1), 1 ' one tenth of circle at centre
Next t
'================================================================
' Albert's formula in yellow
For grad = 0 To 100 Step 1
'----------------------------------------
if grad <= 50 then s = grad *(sqr(.5)/50)
if grad > 50 then s = ((grad-50)*((1-sqr(.5))/50))+sqr(.5)
if (100-grad) > 50 then c=((100-grad-50)*((1-sqr(.5))/50))+ sqr(.5)
if (100-grad) <= 50 then c= (100-grad)*(sqr(.5)/50)
'makes octogon , we have to apply an arc to 0-50 and 50-100
if grad > 50 then s=s + ((50-(grad-50))/12.5)*((grad-50)/100) *.099
if grad > 50 then c=c + ((50-(grad-50))/12.5)*((grad-50)/100) *.04099
if grad <= 50 then s=s + ((50-grad)/12.5)*(grad/100) *.04099
if grad <= 50 then c=c + ((50-grad)/12.5)*(grad/100) *.099
if grad <= 25 then s=s - ((25-grad)/6.25)*(grad/25) *.0008
if grad <= 25 then c=c - ((25-grad)/6.25)*(grad/25) *.0017
if grad > 25 and grad<= 50 then s=s + ((25-(grad-25))/6.25)*((grad-25)/25) *.001
if grad > 25 and grad<= 50 then c=c + ((25-(grad-25))/6.25)*((grad-25)/25) *.001
if grad > 50 and grad<= 75 then s=s + ((25-(grad-50))/6.25)*((grad-50)/25) *.001
if grad > 50 and grad<= 75 then c=c + ((25-(grad-50))/6.25)*((grad-50)/25) *.001
if grad > 75 and grad<= 100 then s=s - ((25-(grad-75))/6.25)*((grad-75)/25) *.0017
if grad > 75 and grad<= 100 then c=c - ((25-(grad-75))/6.25)*((grad-75)/25) *.0008
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 14
Pset (+x, +y), 14
Pset (-x, -y), 14
Pset (-x, +y), 14
t = grad * Pion2 / 100
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 14 ' plot the scatter
Next grad
'================================================================
' Richard's constant parameter version using radians
Const As Double a = -0.166 ' given a, then b is closure parameter
Const As Double b = .00759 ' b = Pion2^-5 - Pion2^-4 - a*Pion2^-2
For t = 0 To Pion2 Step pion2 / 100
'----------------------------------------
tt = t * t
s = t * ( 1 + tt * (tt * b + a))
c = Pion2 - t
tt = c * c
c = c * ( 1 + tt * (tt * b + a))
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
' operation count 8* 4+ 1- 5c 16v
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 13
Pset (+x, +y), 13
Pset (-x, -y), 13
Pset (-x, +y), 13
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 13 ' plot the scatter
Next t
'================================================================
Sleep
'================================================================
Last edited by albert on Nov 22, 2009 23:23, edited 2 times in total.
@Albert. I only need to be better than 1 part in 2000 radially and angularly. There is no point in wasting time using a very slow square root to correct the radius if the angle error is bigger. That is why I do not use normalisation. Your fixation on the radius error is distracting you from the real problem that is your angular error.
Take the following as approximate times for double precision floating point operations....
Sin or Cos 100, Square root 60, Divide 30, Multiply 5, Add or Subtract 3.
Variable fetch or store 2, Constant reference 1.
Try to avoid square roots and division, they slow you down too much.
Take the following as approximate times for double precision floating point operations....
Sin or Cos 100, Square root 60, Divide 30, Multiply 5, Add or Subtract 3.
Variable fetch or store 2, Constant reference 1.
Try to avoid square roots and division, they slow you down too much.
Richard the point is i'm writing a calculator and want to add sine,cosine,and tangent functions to it.
It doesn't matter how long it takes to do an opperation (Well within a minute or so.)
What i'm looking for is something to replace the power series within reason, my calc can extract a square root in less than a minute.
with your help i can now plot perfect circles,I just have to get the degree angles correct or as close as possible.
I'm trying to get everything to under a minute for an answer.
In the above post by me I got the errors to all fall well within the limits of the smaller circle at 100 mag. (back to the drawing board!)
It doesn't matter how long it takes to do an opperation (Well within a minute or so.)
What i'm looking for is something to replace the power series within reason, my calc can extract a square root in less than a minute.
with your help i can now plot perfect circles,I just have to get the degree angles correct or as close as possible.
I'm trying to get everything to under a minute for an answer.
In the above post by me I got the errors to all fall well within the limits of the smaller circle at 100 mag. (back to the drawing board!)
@Albert. Your calculator is multi-precision so it requires results that are precise to however many digits the operator chooses. No single approximation will solve that problem. You will not discover a quick accurate fix for the transcendental functions. I would be happy to bet a million to one that Sine and Cosine will still be around long after you are gone, but I could not collect because you would be dead.
For the multi-precision situation there are really only two ways to generate accurate Sines and Cosines. Those two techniques are Power series and CORDIC. I showed that CORDIC was not applicable for the variable precision case because of the size of the precomputed table. http://www.freebasic.net/forum/viewtopi ... 170#126170
That left power series, but you need to avoid factorial and raising to powers. For that reason I gave you a version that avoided those. http://www.freebasic.net/forum/viewtopi ... 780#125780
If you define your big number format as a UDT and rewrite your calculator as an Overlay of the FB arithmetic operators and functions then you will be able to program your calculator in FB to implement those functions. If you don't know how to do it, then take a look at the example in the body of the my dims.bi file in http://www.freebasic.net/forum/viewtopi ... 177#127177
For the multi-precision situation there are really only two ways to generate accurate Sines and Cosines. Those two techniques are Power series and CORDIC. I showed that CORDIC was not applicable for the variable precision case because of the size of the precomputed table. http://www.freebasic.net/forum/viewtopi ... 170#126170
That left power series, but you need to avoid factorial and raising to powers. For that reason I gave you a version that avoided those. http://www.freebasic.net/forum/viewtopi ... 780#125780
If you define your big number format as a UDT and rewrite your calculator as an Overlay of the FB arithmetic operators and functions then you will be able to program your calculator in FB to implement those functions. If you don't know how to do it, then take a look at the example in the body of the my dims.bi file in http://www.freebasic.net/forum/viewtopi ... 177#127177
Thanks Richard!
I put your loopy version to the test and at 55 loops
(doesn't correct any better at higher loops.)
mag = 100,000,000,000,000 the errors are still within the small circle.
I put the s-c = s-c / sqr(x*x+y*y) in to it also. two times one after the other. The second time corrects degree errors within the small circle.
but were getting processor noise, the circle is jagged.
Richard how do you set your loopy version to dial in a specfic degree?
Can you walk me thru it i don't understand whats going on in the looping.
THANKS!
I put your loopy version to the test and at 55 loops
(doesn't correct any better at higher loops.)
mag = 100,000,000,000,000 the errors are still within the small circle.
I put the s-c = s-c / sqr(x*x+y*y) in to it also. two times one after the other. The second time corrects degree errors within the small circle.
but were getting processor noise, the circle is jagged.
Code: Select all
'================================================================
' Albert's(yellow),Albert's .0009 in green, and Richard's(pink)
' approximations measured relative to the reference circle ... blue
Const As Double mag = 100000000000000 ' magnification of errors
'================================================================
Dim As Double t, tt, grad, q, s, c, dx, dy, x, y, r = 1.1
Dim As Double sqr_root
Screen 19
Dim As Integer w, h, depth
Screeninfo w, h, depth
Window (-r * w / h, -r) - (r * w / h, r) ' maintain pixel aspect
'----------------------------------------------------------------
Line(-c, 0)-(c, 0), 1 ' x axis
Line(0, -c)-(0, c), 1 ' y axis
Const As Double Pion2 = 2 * Atn(1)
For t = 0 To 4 * Pion2 Step Pion2 / 1000
x = Cos(t)
y = Sin(t)
Pset(x, y), 1 ' reference circle
Pset(x*.1, y*.1), 1 ' one tenth of circle at centre
Next t
'================================================================
'Albert's .000945 formula for a 200 degree pi
dim as double sarc,carc,sgrad2,cgrad2
dim as double deg = atn(1)/100
For grad = 0 To 100 Step 1
sgrad2 = grad * 2
cgrad2 = 200 - (grad*2)
sarc = deg * (grad*2)
carc = deg * (200-(grad*2))
s = sarc * ((100- ((sgrad2 * .000945 )*sgrad2 )) / 100)
c = carc * ((100- ((cgrad2 * .000945 )*cgrad2 )) / 100)
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 10
Pset (+x, +y), 10
Pset (-x, -y), 10
Pset (-x, +y), 10
t = grad * Pion2 / 100
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 10 ' plot the scatter
Next grad
'================================================================
' Albert's formula in yellow
For grad = 0 To 100 Step 1
'----------------------------------------
if grad <= 50 then
s = grad *(sqr(.5)/50)
c=((100-grad-50)*((1-sqr(.5))/50))+ sqr(.5)
else
s = ((grad-50)*((1-sqr(.5))/50))+sqr(.5)
c= (100-grad)*(sqr(.5)/50)
end if
'makes octogon , we have to apply an arc to 0-50 and 50-100
if grad > 50 then s=s + ((50-(grad-50))/12.5)*((grad-50)/100) *.099
if grad > 50 then c=c + ((50-(grad-50))/12.5)*((grad-50)/100) *.04099
if grad <= 50 then s=s + ((50-grad)/12.5)*(grad/100) *.04099
if grad <= 50 then c=c + ((50-grad)/12.5)*(grad/100) *.099
if grad <= 25 then s=s - ((25-grad)/6.25)*(grad/25) *.0008
if grad <= 25 then c=c - ((25-grad)/6.25)*(grad/25) *.0017
if grad > 25 and grad<= 50 then s=s + ((25-(grad-25))/6.25)*((grad-25)/25) *.001
if grad > 25 and grad<= 50 then c=c + ((25-(grad-25))/6.25)*((grad-25)/25) *.001
if grad > 50 and grad<= 75 then s=s + ((25-(grad-50))/6.25)*((grad-50)/25) *.001
if grad > 50 and grad<= 75 then c=c + ((25-(grad-50))/6.25)*((grad-50)/25) *.001
if grad > 75 and grad<= 100 then s=s - ((25-(grad-75))/6.25)*((grad-75)/25) *.0017
if grad > 75 and grad<= 100 then c=c - ((25-(grad-75))/6.25)*((grad-75)/25) *.0008
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 14
Pset (+x, +y), 14
Pset (-x, -y), 14
Pset (-x, +y), 14
t = grad * Pion2 / 100
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 14 ' plot the scatter
Next grad
'================================================================
' Richard's constant parameter version using radians
Const As Double a = -0.166 ' given a, then b is closure parameter
Const As Double b = .00759 ' b = Pion2^-5 - Pion2^-4 - a*Pion2^-2
For t = 0 To Pion2 Step pion2 / 100
'----------------------------------------
tt = t * t
s = t * ( 1 + tt * (tt * b + a))
c = Pion2 - t
tt = c * c
c = c * ( 1 + tt * (tt * b + a))
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
' operation count 8* 4+ 1- 5c 16v
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 13
Pset (+x, +y), 13
Pset (-x, -y), 13
Pset (-x, +y), 13
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 13 ' plot the scatter
Next t
'================================================================
' Richard's Loopy version
'================================================================
Const As Integer n = 55 ' number of loop iterations
Dim As Double table(0 To n) ' first set it all up once
Dim As Double aa = 1, k = 1, z, dz
For i As Integer = 0 To n
table(i) = Atn(aa)
k = k / Sqr(1 + aa*aa)
aa = aa * 0.5
Next i
For t = 0 To Pion2 Step pion2 / 100
'----------------------------------------
c = k
s = 0
z = t
aa = 1
For i As Integer = 0 To n
dx = s * aa
dy = c * aa
dz = table(i)
If z >= 0 Then
c = c - dx
s = s + dy
z = z - dz
Else
c = c + dx
s = s - dy
z = z + dz
End If
aa = aa * 0.5
Next i
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
sqr_root=sqr(c*c+s*s)
s=s/sqr_root
c=c/sqr_root
'----------------------------------------
r = 1 + mag * (Sqr(c*c + s*s) - 1) ' magnify radius error
x = c * r
y = s * r
Pset (+x, -y), 15
Pset (+x, +y), 15
Pset (-x, -y), 15
Pset (-x, +y), 15
'line(0,0)-(+x,-y),12 'the lines are just about perfect also.
'line(0,0)-(+x, y),12
'line(0,0)-(-x,-y),12
'line(0,0)-(-x, y),12
dx = c - Cos(t)
dy = s - Sin(t)
Circle(mag * dx, mag * dy), 0.01, 15 ' plot the scatter
Next t
'================================================================
Sleep
'================================================================
Can you walk me thru it i don't understand whats going on in the looping.
THANKS!
There is no point going beyond n = 52 with the Loopy version. The double precision data format only uses 52 bits in the mantissa. The Loopy version effectively does a binary search for the solution one bit at the time. It uses a table of Atan(1^-i) calculated only once. That table is also restricted to 52 bit accuracy. 2^52 = 10^15.5 = 15.5 digits, just like double precision. With mag = 1e15 you are looking at the noise floor of the floating point processor.
There is no point applying normalization to the loopy version as it has been corrected for radius internally and the normalization will just add noise, there is certainly a disadvantage in applying normalization twice as it slows it down twice as much and adds twice the noise.
The Loopy version is one sixth of Walther's Unified CORDIC, the best place to start is http://en.wikipedia.org/wiki/CORDIC You will not be able to understand fully how the loopy version works, it needs 2'nd year University applied maths.
There is no point applying normalization to the loopy version as it has been corrected for radius internally and the normalization will just add noise, there is certainly a disadvantage in applying normalization twice as it slows it down twice as much and adds twice the noise.
The Loopy version is one sixth of Walther's Unified CORDIC, the best place to start is http://en.wikipedia.org/wiki/CORDIC You will not be able to understand fully how the loopy version works, it needs 2'nd year University applied maths.