Gingerbreadman Map

[DDRM]

I came across this in a rather roundabout way while revising, so for a bit of fun I implemented it. I've given a link in the header to the (very brief) wikipedia page - that in turn links to a slightly less brief page on Wolfram MathWorld.

The process takes each point in the plane, and uses the given functions to iterate, until it comes back to its start point (or not!). The time it takes to do that determines the colour of the plot.The code runs happily in either BB4W or BBC-SDL - slightly quicker in the latter, I think, but not much in it.

D

Code: Select all

      REM Gingerbreadman map
      REM https://en.wikipedia.org/wiki/Gingerbreadman_map
      REM Each point in the plane is taken, and iterated using the following rules:
      REM x(n+1) = 1 - y(n) + ABS(x(n))
      REM y(n+1) = x(n)
      REM We count how many iterations it takes to get back to the start (capped at some depth limit)
      REM Points are coloured by their cycle length.

      REM Uncomment whichever of these you fancy... or choose your own values
      min = -10:max = 10: depth = 1000: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -12:max = 12: depth = 1000: step = 0.0625:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -15:max = 15: depth = 250: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -20:max = 30: depth = 500: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2

      REM Define a custom mode into which the plot fits
      VDU 23,22,sf*(max-min)/2 + 2;sf*(max-min)/2 + 2;16,16,16,0
      GCOL 1
      FOR x = min TO max STEP step
        FOR y = min TO max STEP step
          tx = x
          ty = y
          n = 0
          c = 0
          REPEAT
            tx2 = 1 - ty + ABS(tx)
            ty = tx
            tx=tx2
            n += 1
            IF ABS(tx - x)<0.01 AND  ABS(ty - y) < 0.01  THEN c = n
          UNTIL c OR n = depth
          REM Define a custom colour that covers a wide range of cycle lengths
          COLOUR 1, (c MOD 16) * 16,(c MOD 64)*4,c MOD 255
          PLOT 69,(x-min)*sf,(y-min)*sf
        NEXT y
      NEXT x
      END

[DDRM]

As noted above, the basic plot is coloured according to the number of iterations it takes to return to the start point. I was also interested in the "path" it traces out in getting back there - which can take hundreds, or even thousands, of steps.

The code here allows you to explore those paths interactively. As you move the mouse around the space in the window, the path taken by the iterations from that point (whose coordinates are shown at the top) are shown in real time.

Note that the axes for choosing the point and plotting it are different !

Code: Select all

      MODE 21:min = -12:max = 12: depth = 500: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      *REFRESH OFF
      GCOL 1
      cx=0
      cy=0
      REPEAT
        MOUSE x,y,z
        x = x/1600 * (max-min) + min
        y = y/1200 *(max-min) + min
        IF x <> cx AND y <> cy THEN
          cx = x
          cy = y
          c = 0
          n=0
          tx = x
          ty = y
          CLS
          REPEAT
            tx2 = 1 - ty + ABS(tx)
            ty = tx
            tx=tx2
            n += 1
            PLOT 69,(tx-min)*sf,(ty-min)*sf
            IF ABS(tx - x)<0.01 AND  ABS(ty - y) < 0.01  THEN c = n
          UNTIL c OR n = depth
          PRINT "X=";INT(x*100)/100;" Y=";INT(y*100)/100;" Cycle length= ";c
          *REFRESH
        ENDIF
      UNTIL FALSE
      *REFRESH ON
      END

[Hated Moron]

The code runs happily in either BB4W or BBC-SDL - slightly quicker in the latter, I think, but not much in it.

There are a few speed-ups possible. Firstly by replacing the REPEAT loop with a FOR loop, secondly by noticing that you can reduce the first three statements in the inner loop to two statements by using SWAP, and thirdly by using IF rather than AND to avoid an unnecessary y calculation if the x test has already failed. Here's the code with all those changes:

Code: Select all

      REM Gingerbreadman map
      REM https://en.wikipedia.org/wiki/Gingerbreadman_map
      REM Each point in the plane is taken, and iterated using the following rules:
      REM x(n+1) = 1 - y(n) + ABS(x(n))
      REM y(n+1) = x(n)
      REM We count how many iterations it takes to get back to the start (capped at some depth limit)
      REM Points are coloured by their cycle length.

      REM Uncomment whichever of these you fancy... or choose your own values
      min = -10:max = 10: depth = 1000: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -12:max = 12: depth = 1000: step = 0.0625:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -15:max = 15: depth = 250: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -20:max = 30: depth = 500: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2

      REM Define a custom mode into which the plot fits
      VDU 23,22,sf*(max-min)/2 + 2;sf*(max-min)/2 + 2;16,16,16,0
      GCOL 1
      FOR x = min TO max STEP step
        FOR y = min TO max STEP step
          tx = x
          ty = y
          C% = 0
          FOR N% = 0 TO depth
            SWAP tx,ty
            tx = 1 - tx + ABS(ty)
            IF ABS(tx - x)<0.01 IF ABS(ty - y)<0.01 THEN C% = N% : EXIT FOR
          NEXT
          REM Define a custom colour that covers a wide range of cycle lengths
          COLOUR 1, (C% MOD 16) * 16,(C% MOD 64)*4,C% MOD 255
          PLOT 69,(x-min)*sf,(y-min)*sf
        NEXT y
      NEXT x
      END

Click here to run it in your browser.

[Hated Moron]

There are a few speed-ups possible.

In addition (although you might think this is cheating) you can assume the symmetry:

Code: Select all

      REM Gingerbreadman map
      REM https://en.wikipedia.org/wiki/Gingerbreadman_map
      REM Each point in the plane is taken, and iterated using the following rules:
      REM x(n+1) = 1 - y(n) + ABS(x(n))
      REM y(n+1) = x(n)
      REM We count how many iterations it takes to get back to the start (capped at some depth limit)
      REM Points are coloured by their cycle length.

      REM Uncomment whichever of these you fancy... or choose your own values
      min = -10:max = 10: depth = 1000: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -12:max = 12: depth = 1000: step = 0.0625:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -15:max = 15: depth = 250: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2
      REM min = -20:max = 30: depth = 500: step = 0.125:sf = 4/step   :REM keep step as a reciprocal of a power of 2

      REM Define a custom mode into which the plot fits
      VDU 23,22,sf*(max-min)/2 + 2;sf*(max-min)/2 + 2;16,16,16,0
      GCOL 1
      FOR x = min TO max STEP step
        FOR y = min TO x STEP step
          tx = x
          ty = y
          C% = 0
          FOR N% = 0 TO depth
            SWAP tx,ty
            tx = 1 - tx + ABS(ty)
            IF ABS(tx - x)<0.01 IF ABS(ty - y)<0.01 THEN C% = N% : EXIT FOR
          NEXT
          REM Define a custom colour that covers a wide range of cycle lengths
          COLOUR 1, (C% MOD 16) * 16,(C% MOD 64)*4,C% MOD 255
          PLOT 69,(x-min)*sf,(y-min)*sf
          PLOT 69,(y-min)*sf,(x-min)*sf
        NEXT y
      NEXT x
      END

[Hated Moron]

Of course if you want real speed, you can use this alternative approach (runs in BB4W or BBCSDL). Arguably if the original algorithm is expressed as C code this is an even more straightforward conversion, as well as being hundreds of times faster. Whether you consider it BASIC is another matter, but it's no worse than incorporating assembly-language code:

Code: Select all

      REM Gingerbreadman map
      REM https://en.wikipedia.org/wiki/Gingerbreadman_map
      REM Each point in the plane is taken, and iterated using the following rules:
      REM x(n+1) = 1 - y(n) + ABS(x(n))
      REM y(n+1) = x(n)
      REM We count how many iterations it takes to get back to the start (capped at some depth limit)
      REM Points are coloured by their cycle length.

      REM Define a custom mode into which the plot fits
      VDU 23,22,640;640;16,16,16,0

      REM Install library:
      INSTALL @lib$ + "shaderlib"

      REM Create arrays and structures:
      DIM Vertex$(10), Fragment$(999), Float{0%,1%}

      REM Fill vertex and fragment shader arrays from DATA statements:
      PROC_readshader(Vertex$())
      PROC_readshader(Fragment$())

      REM Ensure cleanup on exit:
      ON CLOSE PROCcleanup : QUIT
      ON ERROR PROCcleanup : MODE 7 : PRINT REPORT$ : END

      REM Initialise shader library:
      PROC_shaderinit(oVertex%, oFragment%)

      REM Compile shaders:
      PROC_compileshader(oVertex%, Vertex$(), "Vertex")
      PROC_compileshader(oFragment%, Fragment$(), "Fragment")

      REM Link and use shaders:
      oProgram% = FN_useshaders(oVertex%, oFragment%)

      REM Render loop:
      REPEAT
        PROC_render
        WAIT 10
      UNTIL FALSE
      END

      DEF PROCcleanup
      ON ERROR OFF
      oProgram% += 0 : oVertex% += 0 : oFragment% += 0
      PROC_shaderexit(oProgram%, oVertex%, oFragment%)
      *REFRESH ON
      ENDPROC

      REM Minimal vertex shader:
      DATA "attribute vec4 vPosition;"
      DATA "void main()"
      DATA "{"
      DATA "gl_Position = vPosition ;"
      DATA "}"
      DATA ""

      REM Gingerbread Man map:
      DATA "const int DEPTH = 1000;"
      DATA "const float MIN = -10.0;"
      DATA "const float MAX =  10.0;"

      DATA "void main()"
      DATA "{"
      DATA "float x = gl_FragCoord.x * (MAX - MIN) / 640.0 + MIN;"
      DATA "float y = gl_FragCoord.y * (MAX - MIN) / 640.0 + MIN;"
      DATA "float tx = x;"
      DATA "float ty = y;"
      DATA "float c = 0.0;"
      DATA "for( int n=0; n<=DEPTH; n++ )"
      DATA " {"
      DATA "  float tx2 = 1. - ty + abs(tx);"
      DATA "  ty = tx;"
      DATA "  tx = tx2;"
      DATA "  if ((abs(tx - x) < 0.01) && (abs(ty - y) < 0.01)) {c = float(n); break;}"
      DATA " }"
      DATA "float r = mod(c, 16.0) * 16.0 / 255.0;"
      DATA "float g = mod(c, 64.0) * 4.0 / 255.0;"
      DATA "float b = mod(c, 255.0) / 255.0;"
      DATA "gl_FragColor = vec4( r, g, b, 1.0 );"
      DATA "}"
      DATA ""

[DDRM]

Hi Richard,

Thanks: the swap of order to avoid the reassignment is elegant! I also like the nested IFs instead of the AND conditional. Symmetry: yes, since I'm assuming that the range applies equally to x and y, that's a perfectly reasonable speed-up!

I hadn't realised that a FOR loop, even with a break out of it, would be significantly faster than a REPEAT...

...but the shader version is shockingly fast! Presumably a combination of (a) compiling the HLSL code, and (b) parallel processing by the graphics card. I'll obviously need to revisit HLSL in the context of BBC BASIC programming.

Is there documentation of SHADERLIB somewhere? I can find the D3D11 tutorials on the wiki, but as far as I can see SHADERLIB is not mentioned there or in the main help. Presumably it's a version developed for OPENGL as part of the BBC-SDL project?

I see the Conway demo - I'll investigate that, too...

Best wishes,

D

[Hated Moron]

Is there documentation of SHADERLIB somewhere?

Sadly no, but really there's very little to know apart from the boilerplate code that you can simply copy from any of the example programs that use it (fluid.bbc, mandel.bbc, raytrace.bbc, seascape.bbc, slitscan.bbc and voronoi.bbc).

The key part (and the only part that varies significantly between these examples) is the fragment shader (here 'fragment' really means 'pixel') which is the last piece of code in every case. All that does - and the Gingerbread Man example illustrates this very well - is to take the coordinates of a pixel (they are in gl_FragCoord) and from them calculate the RGB value of that pixel (which is returned in gl_FragColor).

So long as you can derive the latter from the former, writing it as shader code is pretty straightforward (do test it both on a desktop platform and in a browser, because there are slight differences that can catch you out).

conway.bbc is a more complicated example, because the RGB value of a pixel depends on its coordinates and the output from the previous frame.