Welcome to the Parallax Discussion Forums, sign-up to participate.

Heater.
Posts: **21,213**

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

8,348Not 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: EDIT: 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/iteration

final results vary:

Tachyon Forth- compact, fast, forthwright and interactive--->CLICK THE LOGO for more links<---P2+++++TAQOZ INTRO & LINKS+++++P2 SHORTFORM DATASHEETP1+++++Latest binary V5.4 includes EASYFILE+++++Tachyon Forth News BlogBrisbane, Australia8,348Or how about a full circle plot?

Tachyon Forth- compact, fast, forthwright and interactive--->CLICK THE LOGO for more links<---P2+++++TAQOZ INTRO & LINKS+++++P2 SHORTFORM DATASHEETP1+++++Latest binary V5.4 includes EASYFILE+++++Tachyon Forth News BlogBrisbane, Australia21,213I 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.

21,213Not sure which is giving the more accurate result.

The original Minsky is presented using a multiply by number less than one.

51421,21321,213With 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:

259Just 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

21,213Very 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.

21,213Assuming 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.

259Ok 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

21,213Both 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)

21,213The 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:

259A 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

8,348Tachyon Forth- compact, fast, forthwright and interactive--->CLICK THE LOGO for more links<---P2+++++TAQOZ INTRO & LINKS+++++P2 SHORTFORM DATASHEETP1+++++Latest binary V5.4 includes EASYFILE+++++Tachyon Forth News BlogBrisbane, Australia21,213As 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.

21,213You 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.

21,213I 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.

21,213Ah well, here is the code with the magic numbers also:

8,348Using 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.

Tachyon Forth- compact, fast, forthwright and interactive--->CLICK THE LOGO for more links<---P2+++++TAQOZ INTRO & LINKS+++++P2 SHORTFORM DATASHEETP1+++++Latest binary V5.4 includes EASYFILE+++++Tachyon Forth News BlogBrisbane, Australia21,213Hah! 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.

7,986Cannot

Be simply

I think this is a feature of the algorithm (that no temps are needed).

Bean

21,213That 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: 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.

21,2131,981and 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).

21,213I 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.

21,213A 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: 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.

259It 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

21,213You 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.

21,213As 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.