Minsky Circle Algorithm. Generate sin and cos fast with only add, sub and shift.

Heater.Heater. Posts: 19,254
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:
//
// 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
«1

Comments

  • 33 Comments sorted by Date Added Votes
  • Peter JakackiPeter Jakacki Posts: 6,458
    edited August 2015 Vote Up0Vote Down
    Thanks Heater, that looks neat. I love your simplified explanation of an FFT too!
    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 ok
    
    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:
    -1163 -7277
    -935 -7307
    -706 -7330
    -476 -7345
    -246 -7353 ok
    
    Tachyon Forth - compact, fast, forthwright and interactive
    useforthlogo-s.png
    Brisbane, Australia
  • Peter JakackiPeter Jakacki Posts: 6,458
    edited August 2015 Vote Up0Vote Down
    Couldn't help but ASCII plot the X axis to see how it looked and since I can't seem to add this to the first comment I have had to post again:
    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 
                                                                 *
                                                                   *
                                                                     *
                                                                       *
                                                                        *
                                                                          *
                                                                            *
                                                                              *
                                                                               *
                                                                                 *
                                                                                   *
                                                                                     *
                                                                                      *
                                                                                        *
                                                                                          *
                                                                                           *
                                                                                             *
                                                                                              *
                                                                                                *
                                                                                                 *
                                                                                                   *
                                                                                                    *
                                                                                                     *
                                                                                                       *
                                                                                                        *
                                                                                                         *
                                                                                                           *
                                                                                                            *
                                                                                                             *
                                                                                                              *
                                                                                                               *
                                                                                                                *
                                                                                                                 *
                                                                                                                  *
                                                                                                                   *
                                                                                                                    *
                                                                                                                    *
                                                                                                                     *
                                                                                                                      *
                                                                                                                      *
                                                                                                                       *
                                                                                                                       *
                                                                                                                        *
                                                                                                                        *
                                                                                                                        *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                         *
                                                                                                                        *
                                                                                                                        *
                                                                                                                        *
                                                                                                                       *
                                                                                                                       *
                                                                                                                      *
                                                                                                                     *
                                                                                                                     *
                                                                                                                    *
                                                                                                                   *
                                                                                                                  *
                                                                                                                  *
                                                                                                                 *
                                                                                                                *
                                                                                                               *
                                                                                                              *
                                                                                                            *
                                                                                                           *
                                                                                                          *
                                                                                                         *
                                                                                                        *
                                                                                                      *
                                                                                                     *
                                                                                                    *
                                                                                                  *
                                                                                                 *
                                                                                               *
                                                                                              *
                                                                                            *
                                                                                           *
                                                                                         *
                                                                                       *
                                                                                      *
                                                                                    *
                                                                                  *
                                                                                 *
                                                                               *
                                                                             *
                                                                           *
                                                                          *
                                                                        *
                                                                      *
                                                                    *
                                                                   *
                                                                 *
                                                               *
                                                             *
                                                           *
                                                          *
                                                        *
                                                      *
                                                    *
                                                  *
                                                 *
                                               *
                                             *
                                            *
                                          *
                                        *
                                       *
                                     *
                                   *
                                  *
                                *
                               *
                             *
                            *
                           *
                         *
                        *
                       *
                     *
                    *
                   *
                  *
                 *
                *
               *
              *
             *
            *
           *
           *
          *
         *
         *
        *
        *
       *
       *
       *
      *
      *
      *
      *
      *
      *
      *
      *
       *
       *
       *
       *
        *
        *
         *
          *
          *
           *
            *
            *
             *
              *
               *
                *
                 *
                  *
                   *
                    *
                      *
                       *
                        *
                          *
                           *
                            *
                              *
                               *
                                 *
                                  *
                                    *
                                     *
                                       *
                                         *
                                          *
                                            *
                                              *
                                               *
                                                 *
                                                   *
                                                     *
                                                      *
                                                        *
                                                          *
                                                            *
                                                              * ok
    

    Or how about a full circle plot?
    Screenshot%20from%202015-08-23%2013%3A45%3A08.png
    Tachyon Forth - compact, fast, forthwright and interactive
    useforthlogo-s.png
    Brisbane, Australia
  • Of course the "Minsky" algorithm is simply making use of the fact that the differential of sin(x) is cos(x) and the differential of cos(x) = -sin(x). That is to say the rate of change of one is proportional to the value of the other. So you just add and subtract a small fraction of one from the other and you get sin and cos.

    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.
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    Tuns out that using shift is not the same as divide. As you have noticed. The rounding errors are not the same.
    Not sure which is giving the more accurate result.

    The original Minsky is presented using a multiply by number less than one.
  • Thanks for sharing this! I assume that this will also calculate sine values, right? :)
  • Sure does. It produces a sin and cos value on each iteration. That's how we get a circle out of it. See the comments in the code in the first post.
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    Edit: The technique presented here for making the rounding behavior of a shift right the same as an integer divide does work as intended. BUT it using it in the Minsky algorithm seriously degrades its accuracy.Counter intuitive. It's better not to use it. See discussion here: http://forums.parallax.com/discussion/comment/1342503/#Comment_1342503



    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);
    }
    
  • Hi
    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
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    tritonium,

    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.

  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    What I would really like is if someone smarter than me can tell us what frequency this runs at.

    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.
  • Hi
    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
  • Here is a plot of minsky circle, with and without the rounding error correction, and an actual circle computed from cos and sin using floating point.

    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)
    816 x 768 - 40K
  • Of course we can start with a bigger amplitude/radius and use bigger shifts, say $10000000 and 10. Then we get a circle that is a lot more accurate. I can't see any difference from a proper circle when I plot them here. Then the rounding errors are so small the rounding correction is not needed and the code gets simpler and faster.

    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;
    }
    

  • Hi
    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
  • Peter JakackiPeter Jakacki Posts: 6,458
    edited August 2015 Vote Up0Vote Down
    That's strange, I tried a shift of 11 and I still got a circle.
    pub Minsky ( xOld yOld step -- xNew yNew )
    	  >L SWAP OVER IX SAR - SWAP ( xNew yOld )
    	  OVER L> SAR + ( xNew yNew )
    	;
    
    0 -10000 7 11 << FOR 11 Minsky OVER 7 SAR 100 + OVER 8 SAR 44 + XY "*" EMIT NEXT HOME
    
    Screenshot%20from%202015-08-27%2022%3A49%3A10.png
    Tachyon Forth - compact, fast, forthwright and interactive
    useforthlogo-s.png
    Brisbane, Australia
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    If you are shifting by 10, a division by 1024, and your radius is only 4096 then most of your bits disappear in the shift. There is only 2 bits left. That means there can only be 3 possible values added/subtracted to sin and cos in the loop. 0, 1, -1. So it can only draw horizontal, vertical and lines at 45 degrees. Sometimes it amazes me that does actually come out as a nice octagon.

    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.
    878 x 804 - 46K
  • tritonium,

    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.


  • Some magic numbers.

    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.
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    Oops, had to put the attachment here. Stupid forum.

    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;
    }
    
    878 x 804 - 44K
  • Peter JakackiPeter Jakacki Posts: 6,458
    edited August 2015 Vote Up1Vote Down
    I've turned my terminal screen into a "graphics" screen by setting the font size to 0.5 to give me plenty of resolution so each character is a pixel.
    Using your latest settings and scaling this somewhat still:
    CLS !SP BOLD white PEN 0 8192 6487 FOR 10 Minsky OVER 15 / 600 + OVER 32 / 240 + XY "O" EMIT NEXT HOME
    
    Screenshot%20from%202015-08-28%2001%3A02%3A18.png

    @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
    useforthlogo-s.png
    Brisbane, Australia
  • Peter,

    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.
  • BeanBean Posts: 7,805
    edited August 2015 Vote Up0Vote Down
    Heater,
    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
    logo.png?503
    Esterline Research & Design
    thitt@esterlineresearch.com

    We offer consulting on the following areas of expertise:
    Frequency Control - Micro-Controller/Processor Projects
    Test and Automation - General Programming and Coding
    Circuit Design - Board Layouts
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    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:
    NEW X = OLD X - epsilon * OLD Y
    NEW Y = OLD Y + epsilon * NEW(!) X
    
    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.
    1086 x 1009 - 40K
  • Blow me down, works perfectly for generating squares, sort of octagons, whateveragons, all the way to nice circles.

    1077 x 999 - 48K
  • I think this is now a footnote in the textbooks, Bressenhams (ie DDA) is a more general
    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).
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    Mark_T,

    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.
  • Heater.Heater. Posts: 19,254
    edited August 2015 Vote Up0Vote Down
    This little algorithm has some weird and puzzling behaviours. Perhaps someone can explain this:

    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.
    878 x 804 - 47K
  • Hi
    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
  • tritonium,

    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.
  • How come, when I set the radius to 4096 we get a 20 sided polygon? An "icosagon" by the way.

    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.
Sign In or Register to comment.