Minsky Circle Algorithm. Generate sin and cos fast with only add, sub and shift.
Heater.
Posts: 21,230
Sometimes we might want to generate a bunch of points that lie on a circle. Or perhaps a sequence of sin and/or cos values. For a tone generator say. The Minsky Circle Algorithm can do this using only add, subtract and shift operations.
The Prop of course has a sin table in ROM but given the work required to index into that ROM correctly and scale the results there may be occasion when Minsky is less code and faster! The Prop II has no such sin table.
Here is the algorithm in Javascript, I leave it as an exercise to the reader to create that in PASM or Spin if they ever need it, consider it as pseudo code:
It's not a totally accurate circle but if you iterate around that loop forever it will produce stable (repeating) sin and cos values.
Edit: See this post for an optimized version: https://forums.parallax.com/discussion/comment/1342420/#Comment_1342420
The Prop of course has a sin table in ROM but given the work required to index into that ROM correctly and scale the results there may be occasion when Minsky is less code and faster! The Prop II has no such sin table.
Here is the algorithm in Javascript, I leave it as an exercise to the reader to create that in PASM or Spin if they ever need it, consider it as pseudo code:
//
// Minsky Circle Algorithm
//
// Generates circles (sin and cos) using only add, subtract and shift.
//
// From: MIT Programming Hacks: http://www.cl.cam.ac.uk/~am21/hakmemc.html
//
"use strict";
let xOld = 0;
let yOld = -7349; // Amplitude, need not be any special value. Larger the better.
let xNew;
let yNew;
const steps = 202; // This number of steps gives a full circle for the shift size below.
for (let step = 0; step <= steps; step += 1) {
console.log(xOld, yOld); // The output, sin and cos basically.
xNew = xOld - (yOld >> 5); // Use bigger or smaller shifts for different
yNew = yOld + (xNew >> 5); // angular step sizes around the circle.
xOld = xNew;
yOld = yNew;
}
I think I posted this here before some time ago but I stumbled across it again today and thought it was so neat everybody should know about it It's not a totally accurate circle but if you iterate around that loop forever it will produce stable (repeating) sin and cos values.
Edit: See this post for an optimized version: https://forums.parallax.com/discussion/comment/1342420/#Comment_1342420
Comments
Not having a signed shift right although I could create one out of ASR, I just wanted to check it out and just used a signed division by 32 instead. Anyway this is what I get in Forth although I will have to check it properly later.
Copy and Paste straight from my terminal screen:
pub Minsky ( xOld yOld -- xNew yNew ) SWAP OVER 32 / - SWAP ( xNew yOld ) OVER 32 / + ( xNew yNew ) ; ok ok 0 -7349 200 FOR Minsky CR OVER PRINT SPACE DUP PRINT NEXT 229 -7342 458 -7328 687 -7307 915 -7279 1142 -7244 1368 -7202 1593 -7153 1816 -7097 2037 -7034 2256 -6964 2473 -6887 2688 -6803 2900 -6713 3109 -6616 3315 -6513 3518 -6404 3718 -6288 3914 -6166 4106 -6038 4294 -5904 4478 -5765 4658 -5620 4833 -5469 5003 -5313 5169 -5152 5330 -4986 5485 -4815 5635 -4639 5779 -4459 5918 -4275 6051 -4086 6178 -3893 6299 -3697 6414 -3497 6523 -3294 6625 -3087 6721 -2877 6810 -2665 6893 -2450 6969 -2233 7038 -2014 7100 -1793 7156 -1570 7205 -1345 7247 -1119 7281 -892 7308 -664 7328 -435 7341 -206 7347 23 7347 252 7340 481 7325 709 7303 937 7274 1164 7238 1390 7195 1614 7145 1837 7088 2058 7024 2277 6953 2494 6876 2708 6792 2920 6701 3129 6604 3335 6500 3538 6390 3737 6274 3933 6152 4125 6024 4313 5890 4497 5750 4676 5604 4851 5453 5021 5297 5186 5135 5346 4968 5501 4797 5650 4621 5794 4440 5932 4255 6064 4066 6191 3873 6312 3676 6426 3476 6534 3272 6636 3065 6731 2855 6820 2642 6902 2427 6977 2209 7046 1989 7108 1767 7163 1544 7211 1319 7252 1093 7286 866 7313 638 7332 409 7344 180 7349 -49 7348 -278 7340 -507 7325 -735 7303 -963 7273 -1190 7236 -1416 7192 -1640 7141 -1863 7083 -2084 7018 -2303 6947 -2520 6869 -2734 6784 -2946 6692 -3155 6594 -3361 6489 -3563 6378 -3762 6261 -3957 6138 -4148 6009 -4335 5874 -4518 5733 -4697 5587 -4871 5435 -5040 5278 -5204 5116 -5363 4949 -5517 4777 -5666 4600 -5809 4419 -5947 4234 -6079 4045 -6205 3852 -6325 3655 -6439 3454 -6546 3250 -6647 3043 -6742 2833 -6830 2620 -6911 2405 -6986 2187 -7054 1967 -7115 1745 -7169 1521 -7216 1296 -7256 1070 -7289 843 -7315 615 -7334 386 -7346 157 -7350 -72 -7348 -301 -7339 -530 -7323 -758 -7300 -986 -7270 -1213 -7233 -1439 -7189 -1663 -7138 -1886 -7080 -2107 -7015 -2326 -6943 -2542 -6864 -2756 -6778 -2967 -6686 -3175 -6587 -3380 -6482 -3582 -6371 -3781 -6253 -3976 -6129 -4167 -5999 -4354 -5863 -4537 -5722 -4715 -5575 -4889 -5423 -5058 -5265 -5222 -5102 -5381 -4934 -5535 -4762 -5683 -4585 -5826 -4403 -5963 -4217 -6094 -4027 -6219 -3833 -6338 -3635 -6451 -3434 -6558 -3230 -6658 -3022 -6752 -2811 -6839 -2598 -6920 -2382 -6994 -2164 -7061 -1944 -7121 -1722 -7174 -1498 -7220 -1273 -7259 -1047 -7291 -820 -7316 -592 -7334 -363 -7345 okEDIT: Just remembered I did have an opcode defined using the OPCODE directive as SAR so that the definition produces somewhat similar results but only takes 29us/iterationfinal results vary:
pub Minsky ( xOld yOld -- xNew yNew ) SWAP OVER 5 SAR - SWAP ( xNew yOld ) OVER 5 SAR + ( xNew yNew ) ; ok ok 0 -7349 200 FOR Minsky CR OVER 7 SAR 60 + SPACES "*" EMIT NEXT * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * okOr how about a full circle plot?
I first saw this in 1980, done in hardware! All you need is a couple of registers, some adders and a feedback loop. A little board of TTL chips. It was used to generate the rotating wave forms needed in those old vector graphic radar displays. It had probably been designed 10 years earlier before we had microprocessors.
Thanks for the comment on the FFT commentary. Good to hear it made a bit of sense. It's not really close to explaining the FFT yet. That "weird" rule I have there for multiplying values gives rise to rotations (sin/cos) in much the same way as this algorithm. I have been puzzling how to write that FFT without complex numbers document since the "Fourier Transform For Dummies" document was posted here 5 years ago. I might get there eventually.
Not sure which is giving the more accurate result.
The original Minsky is presented using a multiply by number less than one.
With the addition of a three more operations, a compare, three ANDs, and an ADD, we can correct the rounding errors caused by using shifts instead of proper division. Here is the modified version. I changed x and y to sin and cos in case that is any clearer. Also removed the redundant cosNew:
"use strict"; let sin = 0; let cos = -7349; // Amplitude, need not be any special value. Larger the better. let sinNew; const steps = 202; // This number of steps gives a full circle for the shift size below. for (let step = 0; step <= steps; step += 1) { console.log(sin, cos); // The output sinNew = sin - div32(cos); // Use bigger or smaller shifts for different cos = cos + div32(sinNew); // angular step sizes around the circle. sin = sinNew; } // Divide by 32 using shift. You would in line this in assembler. function div32(x) { let result = x >> 5; // Use shift for division if ((x & 0x80000000) && (x & 0x1F)) result += 1; // Optional line: Correct for rounding errors compared to division return(result); }Just trying to understand what is happening here
if ((x < 0) && (x & 0x1F)) result += 1;
My poor understanding leads me to think it means-
0 = false anything else = true
both halves (x<0) (x & 0x1F) must be true (ie not zero) to increment x ie-
IF x is negative AND any of the 5 least significant bits of x are set
THEN add 1 to x
Now if that is so then why? especially why &0x1F
I can see that you are checking five bits after doing five shifts but....why?
Dave
Very good questions. I'm not sure if I can answer very clearly.
Firstly you have to be aware that this "pseudo code" is actually Javascript. JS only has one number type, a 64 bit floating point number. However that allows for 53 bit integer values. But then the logical operators like &, |, <<, >>, will always produce 32 bit signed integer results. Turns out his fact can be used to write JS that will execute nearly as fast as C !
But never mind all that Javascript specific weirdness.
We find that when we do a division by a power of 2, like 2, 4, 8, 16, 32..., using a shift rather than an actual divide we get rounding errors. This is true no matter what language you are using.
As far as I can tell the worst case of this is -1. Which is 0xffffffff. Shifting that just gives you 0xffffffff back again, which is -1, instead of the zero you want!
So far we see the problem occurs for negative numbers, hence the test for "(x < 0)". Or in my recent edit (x & 0x80000000) which is just checking the sign bit of a 32 bit integer.
Then we find the correction for rounding error need not be applied when the number we are dividing is a multiple of 32 (The amount we are dividing by). Hence the (x & 0x1F) checking that the low 5 bits are zero. Which is the same as a test with the modulus operator (-x % 32) == 0. Except of course that an & is a lot faster than a modulus which is actually a division.
In short the checking of the low five bits is a check to see if the number is exactly divisible by 32, which it will be if the low 5 bits are zeros.
Now, this is possibly not quite right. I can't really justify it to myself yet! But testing so far against actual division has always produced the correct result.
Perhaps someone wiser than be can explain why.
Assuming we can iterate this at a constant "n" times per second, we have an output frequency that depends on the little fractions we are using 1/m, in this case 1/32. So the output frequency must be proportional to n and some how related to m.
If I understood that we could actually specify n and m to get a frequency we want rather than just guessing and tweaking.
Ok Thanks for clearing that up, however I cant resist trying out new algorithms- even (especially) when I don't understand them, so I compared this with one I found by Adafruit (Bressingham?) using propbasic and plotting one circle on top of the other using 320 X 240 4 colour vga driver.
It wasn't straight forward because the number of steps, shift number and amplitude (Radius?) all interact. In order to plot a 100 pixel radius circle I used an amplitude of -6400, a shift of 7 (not 5) and divided the result by 64 (by shifting 6 right).
Plotting each in its own colour the results were very close, but applying the rounding error had the effect of apparently distorting (squashing) the circle by a pixel or so with perfect overlap at NE and SW and misalignment at SE and NW (compass)
Dave
Both minski's are a little bit distorted and it arguable whether the rounding error correction makes it better or worse. (I forget which minsky is which on that plot now)
The down side is that yon need to iterate of 64 times as many steps to get all the way around the circle. Partly mitigated by having made the loop a bit faster.
We end up with:
"use strict"; let sin = 0; let cos = 0x10000000; // Amplitude, need not be any special value. Larger the better. // I find it blows up after 0x7fffffe6 using 32 bit math. let sinNew; const steps = 6464; // This number of steps gives a full circle for the shift size used below. for (let step = 0; step <= steps; step += 1) { console.log(sin, cos); // The output sinNew = sin - (cos >> 10); // Use bigger or smaller shifts for different cos = cos + (sinNew >> 10); // angular step sizes around the circle. sin = sinNew; }A few more observations just changing the numbers around.
with a ystart of -3200 and steps of 18000 and shift of 10 the circle starts to look like a polygon and the 18000 steps gives 2 circles on top of each other which are not perfectly aligned.
I realise that shifting by 10 = divide by 1024 which coarsens up the calculations - a shift of 11 actually gives a nice 8 sided polygon with beautifully straight sides, the vertical and horizontal sides being twice the length of the corners.
I'm sure there are perfectly good mathematical explanations for this, and those sins and coses take up some strange values, but its so far over my head as to be rather frustratingly amusing.
Dave
As you radius gets bigger, more bits are left over after the shift and more directions of travel are possible. See attachment.
Of course if you are working with such small radii, and don't need so much accuracy then smaller shifts can be used. Which gets you your circle points faster.
You don't need any fancy maths to get an idea how this works.
Imagine a graphs of a sine wave and a cosine wave. We can think of the "rate of change" of these curves. Rate of change is simply how much the line goes up or down as we move along the x axis.
Clearly the sine and cosine have a zero rate of change at the top and bottom. They are neither moving up or down at those points. Also we see that the rate of change is maximum at the point they cross the x axis.
Now, magically, it turns out that the rate of change of sine is actually the value of cosine, the bigger cosine is the faster sine is changing. Similarly the rate of change of cosine is minus the value of sine, the bigger sine is the faster cosine is changing. The size of one is the speed of the other, as it were. You can sort of see that from looking at the graphs for a few seconds.
What this algorithm does is, every time around the loop it adds small part of sine to cosine and subtracts a small part of cosine from sine. That means the bigger one is the more the other one is being changed. The rate of change of one is proportional to the size of the other. Which is exactly what we want to draw out the sine and cosine waves.
I found that if you set the radius to 8192 and have a shift of 10. Minsky comes back to exactly it starting point after 6487 iterations. Like so:
-512 7680
-520 7679
-527 7678
-534 7677
-541 7676
-548 7675
-555 7674
-562 7673
-569 7672
-576 7671
...
-440 7680
-448 7680
-456 7680
-464 7680
-472 7680
-480 7680
-488 7680
-496 7680
-504 7680
-512 7680
Further, subtracting 512 from the every output centres our circle about 0,0 nicely and when drawn against a "proper" circle it shows no distortion. See attachment. We end up with a circle of radius 7680. No doubt there are other powers of 2 will produce nice results as well.
Ah well, here is the code with the magic numbers also:
"use strict"; let sin = 0; let cos = 8192; // Amplitude, need not be any special value. Larger the better. // I find it blows up after 0x7fffffe6 using 32 bit math. const steps = 6487; // This number of steps gives a full circle for the shift size used below. let sinNew; for (let step = 0; step < steps; step += 1) { console.log(sin - 512, cos - 512); // The output sinNew = sin - (cos >> 10); // Use bigger or smaller shifts for different cos = cos + (sinNew >> 10); // angular step sizes around the circle. sin = sinNew; }Using your latest settings and scaling this somewhat still:
@Heater: BTW, if you missed inserting the image into your post you can pretend to start a new post and insert the image into that as you would but grab the image link, save the current post as a draft, edit the previous post and paste in that link. Works like a charm and you can remove the link from the draft to keep it clean for next time although I haven't yet tried deleting the draft to see what happens.
Hah! That must be the highest resolution ASCII art I have ever seen. Well it would be if I could make out the characters.
Other Hah! That is some awesomely convoluted way of attaching to a previous post.
Cannot
sinNew = sin - (cos >> 10); // Use bigger or smaller shifts for different cos = cos + (sinNew >> 10); // angular step sizes around the circle. sin = sinNew;Be simply
sin = sin - (cos >> 10); // Use bigger or smaller shifts for different cos = cos + (sin >> 10); // angular step sizes around the circle.I think this is a feature of the algorithm (that no temps are needed).
Bean
That is a very good point and one that has been bugging be for some time.
Firstly the algorithm, as described to me in 1980 or so, when it was implemented in TTL hardware was, as you say, without those temporary variables.
But then I found the "Minsky" version which uses those temps. As seen here : http://www.cl.cam.ac.uk/~am21/hakmemc.html. It even emphasises the point by writing: A year or so ago I tried to write this thing as you suggest without the temporary variables because I was sure they were not needed. It did not work, my circles were all kind of weird ellipses.
Now you suggest the idea of removing the temps and I try it again. It works perfectly! As follows:
// // Minsky Circle Algorithm // // Generates circles (sin and cos) using only add, subtract and shift. // // From: MIT Programming Hacks: http://www.cl.cam.ac.uk/~am21/hakmemc.html // // Note: Does not use the temporary "NEW" variables of the code in that link. // "use strict"; let sin = 0; let cos = 8192; // Amplitude, need not be any special value. Larger the better. // I find it blows up after 0x7fffffe6 using 32 bit math. const steps = 6487; // This number of steps gives a full circle for the shift size used below. for (let step = 0; step < steps; step += 1) { console.log(sin - 512, cos - 512); // The output sin = sin - (cos >> 10); // Use bigger or smaller shifts for different cos = cos + (sin >> 10); // angular step sizes around the circle. }And the resulting plot of the generated circle and a "real" circle. See attachment.Now I have no idea why I had such weird results before. Probably a silly error.
Thanks for the push there.
and flexible approach and doesn't risk accumulating errors.
For circles you use r^2 - x^2-y^2 as the value you keep as close to zero as you can
while varying one of the variables (ie x or y). You use the fact that as x is incremented
x^2 increases by 2x+1 to lose all the multiplies. However many steps you take the fact
that r^2 - x^2-y^2 is minimized means no accumulation of error, and besides all the
operations are exact integer ops if you constrain the centre of the circle to be on the
grid (though you can relax this with fixpoint ops).
I do agree that there is not much call for this now a days.
However, can you produce the circle above, in less operations with Bresenham's algorithm? Let's talk generating the highest possible frequency of sine wave, with about the accuracy of the above, on a Propeller.
Up for the challenge? I was thinking I should do this in PASM anyway.
A shift right is not the same as integer divide. That is x >> 10 is not the same as x / 1024. The rounding comes out differently. Turns out you can make the rounding come out the same if you do a little work on the negative results. Like so:
result = x >> 10; if ((x < 0) && (-x & 0x3FF)) result += 1;I thought that getting the same rounding behaviour when using shift instead of divide would probably produce a more accurate circle in the minsky algorithm. Turns out this is not the case. It's much worse! See attachment. The "circle" plot there is done with floating point, it is overlaid with minsky using shift, minsky using integer divide or corrected shift is way out.This is very weird, we have the following result:
Minsky with floating point divide or plain shift produces very accurate circles,
Minsky with integer divide or the corrected shift produces horrible circles!
Some how the rounding errors of plain shift compensate for it not being floating point.
Just to nerd out on numbers I get a maximum error in radius of 18 and an RMS error of 8. So 0.1%. Going up to a radius of a billion the absolute errors were about the same so 0.000001%. Not bad.
It might just be me, but it seems that the green circle, and yes even the blue one at closer inspection, is actually a series of joined up straight lines, or many sided polygon. Even in my low resolution screen attempts I have noticed this, apart from the obvious 8 sided already mentioned.
Dave
You are right. When you are working with a low number of bits you are going to end up with straight lines and "blocky" results. We are working with a division of 1024 here so you soon eat bits.
As I showed here: https://forums.parallax.com/discussion/download/115079/plot.png
Like I said above, start with a radius of a billion or so and your error is very small.
As seen here in my earlier post https://forums.parallax.com/discussion/comment/1342422/#Comment_1342422
Here we are mashing around with multiples of 2, divide by 1024, radius 4096, and we get a multiple of ten popping up. Weird.
Not what I would have guessed given that a radius of 1024 gives a square and 2048 is a sort of octagon.
How many "sides" does that radius of 8196 have?
Yeah, I know, I'm getting obsessed with this silly thing.