UxMand: Multi-colour Mandlebrot for the Unexpanded VIC-20!

[Snial]

Hi folks,

For probably my last VIC-20 program of the year, I thought I'd write a colour Mandlebrot program for the unexpanded VIC-20, as a bit of a foray into multi-colour bitmapped graphics. The code is actually based on the modifications I submitted to Ken Shirrif's Mandlebrot generator for the Xerox Alto: https://www.righto.com/2017/06/.

It uses a similar technique as for my UxHiLife program I wrote about a year or two ago, where the .prg calls some machine code which remaps the video to $1000 and then the character set can occupy the last 3200 bytes of RAM up to $1fff.

At 5:49 for a full screen, it's not particularly optimised, I think I could improve it by about 50% and increase the accuracy by 1 bit. Maybe I could add an interactive zoom feature (but with 16-bit arithmetic, that won't go very far ).



But it was a bit of fun! The source code (assembled with mac2c64 and tested with VICE's xVIC emulator is attached along with the assembled .prg.

[Mike]

Nice! Sub 10 minutes for a (near) full screen Mandelbrot fractal on the VIC-20 is quite a feat.

You should check the bail-out condition, though. It seems you're using a lower value than 2 in |z| >= 2, which introduces visible artifacts (see the red-yellow stripe in the large blue area, repeated via self-similarity into the deeper iteration bands).

Some years ago, I did a Mandelbrot in B/W, with a display algorithm that highlights the filaments around the set. It uses the arithmetic routines in the BASIC interpreter (with two shortcuts) and achieves roughly 120 iterations per second.

BTW, it's the Mandelbrot fractal, not "Mandlebrot".

Greetings,

Michael

[Snial]

You should check the bail-out condition, though. It seems you're using a lower value than 2 in |z| >= 2, which introduces visible artifacts (see the red-yellow stripe in the large blue area, repeated via self-similarity into the deeper iteration bands).

You’re right, I thought there might be some math errors remaining. This was the first version that basically looked like a Mandelbrot. Also, I probably keep getting the ‘el’ the wrong way around, because it looks so much like a French word and I tend to expect French words to have a ‘le’ order instead.

I’ll check the math and post an update. I believe I bail out at >=2.0 though.

[Mike]

The circle around the set with radius 2 (i.e., |z|=2) is the smallest one that can be used as bail-out condition. Any smaller circle truncates the 'antenna' at the 'head' of the Mandelbrot set and adds those artifact bands. Conversely, bigger values for comparison do not change the shape of the set itself but the appearance of the iteration bands only (they become more 'dense') and, mostly, this slows the calculation (it takes longer for the values of z to bail out - if they ever bail out, that is).

From the source I see you add the two squares z.Re² and z.Im² to calculate the magnitude squared of z, i.e. |z|². That's o.k. (I use the same method in my own program), but that value should then be compared with 4, actually - here's the inner loop of my program before I translated it to machine code:

17 FOR I=1 TO N : X2 = X*X : Y2 = Y*Y : IF X2+Y2 < 4 THEN XY = X*Y : X = X2-Y2 + R : Y = XY+XY + J : NEXT

(the eventually dangling NEXT gets popped by a named NEXT S later on, so the two FOR loops in S and T work without issues)

[Snial]

From the source I see you add the two squares z.Re² and z.Im² to calculate the magnitude squared of z, i.e. |z|². That's o.k. (I use the same method in my own program), but that value should then be compared with 4, actually

Yes, it should be doing that, because in 4:12 fixed-point arithmetic the value 4 is 0100:0000 0000 0000 which is 16384. So, ANDing the top byte with $c0 will trap: 0100 to 1111, i.e. anything from 4 up.

Code: Select all

	clc
	lda x2sp
	adc y2sp
	lda x2sp+1
	adc y2sp+1
	and #$c0 ;>=16K?
	bne MandXy5 ;done.

We don't need any STAs in that sequence, because we don't need to store the result. z.Re²=x2sp and z.Im²=y2sp, so in fact that sequence should add them.

Having another look.

[Snial]

I found out the real problem. There are overflow errors in the original calculation. As it stands, 4:12 arithmetic doesn't provide enough integer bits, because some values for X can be ≤0xc000 or ≥0x4000. Consequently, when X2 is generated and occasionally when X2+Y2 are generated, there's an overflow. I resolved these in an ad-hoc fashion by testing for X and Y before squaring and I also tested for overflow when X2+Y2 is calculated.

With those changes I get a much cleaner image:



In addition I took the opportunity to speed up the multiplication routine by combining the shifted Multiplicand with the 32-bit result (which gets shifted on each pass). This means the entire image is generated in 4:09mins instead of 5:49mins.

The new code is attached, but it's also at:

https://sites.google.com/site/libby8dev/home/vic20

[Mike]

There are overflow errors in the original calculation. [...]

I'm glad to see you sorted it out.

Do you have any ambitions to port the program over to use MINIGRAFIK? Granted, that requires at least a +8K RAM expansion then, but you get a slightly higher resolution of 80x192 multi-colour pixels and you have enough room for a table-based implementation of the multiply routine - using a·b = ((a+b)²-(a-b)²)/4 - that should at least give twice the speed.

[Snial]

I'm glad to see you sorted it out.

Thanks!

Do you have any ambitions to port the program over to use MINIGRAFIK?

I've come across MINIGRAFIK a few times and had a look at the API and it is admirable. It's more that I'm interested in programs that fit on an unexpanded machine. Hence the Conway's Life I wrote a couple of years back: UxLife, then my foray into UxForth (a partly ItsyForth inspired attempt to squeeze a Forth onto a standard VIC-20); then Blastar for the VIC-20 too and now UxMand.

table-based implementation of the multiply routine - using a·b = ((a+b)²-(a-b)²)/4

Yes, it'd be nice to do that!

[srowe]

Hence the Conway's Life I wrote a couple of years back: UxLife, then my foray into UxForth (a partly ItsyForth inspired attempt to squeeze a Forth onto a standard VIC-20);

You might be interested in my Forth implementation

http://www.sleepingelephant.com/ipw-web ... f=2&t=7557

[MrSterlingBS]

Hello,

here is a little faster version of the program.
New Version V1.2 is 3:39 minutes.
That is 26 seconds faster and the code is 280 bytes shorter. I moved the screen at $0200.

[JonBrawn]

The set is symmetric about the X-axis, so you can double your speed by calculating the +Y values and then plotting both +Y and -Y.

[Mike]

The set is symmetric about the X-axis, so you can double your speed by calculating the +Y values and then plotting both +Y and -Y.

As soon as you are going to zoom in on interesting parts, i.e. on the border of the Mandelbrot set, you will very likely be handling with views that are only located fully on either side of the x-axis. *) In practice, in my own versions in particular, I never bothered to implement that 'optimization', especially not within the inner calculation loop. The act of merely testing for that condition slows down the rendering process in 99%+ of all cases!

If anything, I would put that test on the top level, and in case the user inputs two y-coordinates y_min, y_max with different signs, let the program calculate either the [y_min;0] or [0;y_max] spans, whichever is bigger; then mirror the missing part as graphics operation on the bitmap. That would however also only be correct, when y=0 actually is one of the calculated spans.


*) the exception being zooms into the 'antenna' to the left of the main Mandelbrot set.