jnw@mcnc.UUCP (John White) (09/06/85)
[The above line was for the line-eater]
This is a quickly hacked together version 0.0 of a program to
plot the number of iterations of (z <- z*z+b) required for |z|>=2 starting
at z=0. The plot is done on the complex plane of b. See the august 1985
issue of "Scientific American" for more info.
The DeSmet C compiler is used and in-line assembler is used for
the graphics (everything is done through BIOS) and for the critical
inner loop. Note that integer arithmetic is used in the inner loop so
no floating point processor is needed.
Currently the limit of iterations is 255.
If |z|<2 after that, then b is assumed to be in the Mandelbrot set which
is plotted as black. Other colors are used when |z|>=2 for a smaller count.
Alternation of colors is used to display the contors of increasing count.
Since I don't have a floating point processor, I had to use integer
arithmetic and fudge things a little. Because of this, the scale is
specified by an integer. A value of 180 is used to display the whole set.
The maximum magnification is with a value of 1. This is the scale on the
real axis, the program multiplies this by 1.15 to get the scale on the
imaginary axis but you can override this. The center of the plot can also
be specified. Note, the displayed area is real=320*scale/16384,
imaginary=200*scale/16384.
The real and imaginary parts of z are represented with 16 bit integers
with the value 2 being represented as 32k and -2 as -32k. Because of this,
specifing a center for the plot with part of the real or imaginary axis
less than -2 or greater than 2 will cause a repeated image to occur. This
is an artifact. Note that with small values of scale the aspect ratio will
be off.
A typical plot might take ten minutes. Time taken is strongly dependant
on the fraction of the points that are in the Mandelbrot set as these require
the full number of iterations.
This program was originaly written for a Tandy 2000 with high res color
graphics (640 by 400 by 8 colors) and the plots on the Tandy are far better
than those on the IBM. But the IBM plots are still interesting.
Areas you might want to zoom in on (use scale=2):
whorls, at -.72, .25
hooks, at -.915, .27
front, at -1.41, 0
miniature Mandelbrot, at -.16, 1.03
star, at -.55, .64
connection cusp, at -.76, .11
- John N. White <mcnc!jnw>
-------------------------------------------------------
/* Program to compute the Mandelbrot set.
* (convergence of z = z*z + b using complex numbers)
* Arithmetic is done in ordinary integers with 32k corresponding to 2
* and -32k corresponding to -2.
* Compile with the DeSmet compiler.
* Originaly for a Tandy 2000 with color graphics (640 by 400 by 8 color).
* Converted to run on the IBM in 320 by 200 by 4 color mode (mode 4).
* Written by John N. White, 9/2/85
*/
#define stdin 0 /* what DeSmet does in stdio.h */
#define XSIZE 320 /* number of horizontal pixels */
#define YSIZE 200 /* number of vertical pixels */
#define ASPECT 1.15 /* aspect ratio (adjust for your monitor) */
#define MAXCOUNT 255 /* max number of iterations */
#define SPLITCOUNT 127 /* count above which everything is a solid color */
int count; /* limiting number of iterations */
int zr, zi, br, bi, br2; /* real/imaginary parts of z/c */
double brd= -.24, bid; /* low/high limits of c */
int x_point, y_point, color;/* inputs to put_point() */
char buf[256], fname[256], fbuf[16000];
int xstep, ystep; /* increment for each pixel (1/16k units) */
int fn; /* holds number of open file for r/w display */
main(){
int i, j, k;
printf("Mandelbrot set for the IBM with color graphics\n");
printf("Version 0.0, written by John N. White\n\n");
printf("After displaying something, press <ENTER> to continue.\n");
xstep=180, ystep=180*ASPECT;
loop:
printf("0 exit\n1 display file\n2 set display window\n");
printf("3 save last display\n4 calculate Mandelbrot set\ncommand: ");
fgets(buf, 256, stdin);
i= -1;
sscanf(buf, "%d", &i);
switch(i){
case 0:
printf("ok, bye");
exit(0);
case 1: /* display file */
printf("name of file to display: ");
*fname=0;
fgets(buf, 256, stdin);
sscanf(buf, "%s", fname);
if(*fname==0) break;
if((fn=open(fname, 0)) <= 0){
printf("can't read file %s\n", fname);
break;
}
read(fn, fbuf, 16000);
close(fn);
setup();
for(j=8, k= -1, y_point=0; y_point<YSIZE; y_point++){
for(x_point=0; x_point<XSIZE; x_point++){
if(j>6) j=0, k++;
color= (fbuf[k]>>j) & 3;
j+=2;
put_point();
}
}
fgets(buf, 256, stdin);
restore();
break;
case 2: /* set display window */
printf("real step size (number of 1/16k units per step) [%d]: ", xstep);
fgets(buf, 256, stdin);
i= -1;
sscanf(buf, "%d", &i);
if(i!= -1){
xstep=i;
ystep=xstep*ASPECT;
}
printf("imaginary step size (number of 1/16k units per step) [%d]: ",
ystep);
fgets(buf, 256, stdin);
i= -1;
sscanf(buf, "%d", &i);
if(i!= -1) ystep=i;
printf("enter real and imaginary coordinants of center [%f %f]\n: ",
brd, bid);
fgets(buf, 256, stdin);
if(*buf<' ') break;
sscanf(buf, "%lf%lf", &brd, &bid);
printf("using center %f %f\n", brd, bid);
break;
case 3: /* save last display in a file */
gfnloop:
printf("enter name of save file: ");
fgets(buf, 256, stdin);
if(*buf<' ') break;
sscanf(buf, "%s", fname);
if(*fname){
if(exists()){
printf("file %s already exists\n", fname);
goto gfnloop;
}
if((fn=creat(fname)) <= 0){
printf("can't write file %s\n", fname);
goto gfnloop;
}
}
write(fn, fbuf, 16000);
close(fn);
break;
case 4: /* calculate Mandelbrot set */
br2= round(brd*16384.-(double)xstep*(XSIZE/2));
bi= round(bid*16384.+(double)ystep*(YSIZE/2));
count= MAXCOUNT;
setup();
for(j=8, k= -1, y_point=0; y_point<YSIZE; y_point++, bi-=ystep)
for(x_point=0, br=br2; x_point<XSIZE; x_point++, br+=xstep){
i=do_point();
/* set color to 0-3 */
color= (i<SPLITCOUNT) ? (i&1)+1 : (i==MAXCOUNT ? 0 : 3);
if(j>6) fbuf[++k]=j=0;
fbuf[k]|= (color<<j);
j+=2;
put_point();
}
fgets(buf, 256, stdin); /* wait for <enter> */
restore();
break;
default:
printf("try again\n");
}
goto loop;
}
/* returns num of itterations before z = z*z + b exceeded |z|=2 */
dummy(){
#asm
do_point_:
push bp
xor ax,ax
mov word zr_,ax ;set zr,zi to 0
mov word zi_,ax
mov bp,word count_
cloop_:
mov ax,word zr_
imul ax
xchg ax,si ;set si,di to be low,high part of zr*zr
mov di,dx
mov ax,word zi_
imul ax
mov bx,ax ;set bx,cx to be low,high part of zi*zi
mov cx,dx
add ax,si
adc dx,di
test dh,0C0h
jnz do_point_exit_ ;if |z|>=2, exit
dec bp
jz do_point_exit_ ;if count limit exeeded
sub si,bx
sbb di,cx
shl si,1
rcl di,1
shl si,1
rcl di,1
jo do_point_ov1_
shl si,1
adc di,word br_ ;di= new zr = zr*zr-zi*zi+br
jo do_point_exit_ ;if new zr>=2
do_point_cont1_:
mov ax,word zr_
imul word zi_
shl ax,1
rcl dx,1
shl ax,1
rcl dx,1
shl ax,1
rcl dx,1
jo do_point_ov2_
shl ax,1
adc dx,word bi_ ;dx= new zi = 2*zr*zi+bi
jo do_point_exit_ ;if new zi>=2
do_point_cont2_:
mov word zr_,di
mov word zi_,dx
jmp cloop_ ;loop while <count iterations
do_point_exit_:
mov ax,word count_
sub ax,bp
pop bp
ret ;value in AX is value returned by routine
do_point_ov1_:
shl si,1
adc di,word br_ ;di= new zr = zr*zr-zi*zi+br
jo do_point_cont1_ ;if new zr<2
jmp do_point_exit_
do_point_ov2_:
shl ax,1
adc dx,word bi_ ;dx= new zi = 2*zr*zi+bi
jo do_point_cont2_ ;if new zi<2
jmp do_point_exit_
#
}
/* setup display */
setup(){
#asm
mov ax,4 ;put display in graphics mode 4
int 10h
mov ah,1 ;turn off cursor
mov cx,4000h
int 10h
mov ah,11 ;set up color palette
mov bx,0 ;background color is black
int 10h
mov ah,11
mov bx,101h ;1=cyan, 2=magenta, 3=white
int 10h
mov ax,600h ;clear screen
mov bh,0
mov cx,0
mov dx,1827h
int 10h
#
}
/* restore */
restore(){
#asm
mov ax,2 ;reset display to mode 2
int 10h
#
}
/* ipoint sets pixel x_point,y_point to color */
dummy2(){
#asm
put_point_:
mov ah,12 ;write dot
mov al,byte color_
mov cx,word x_point_
mov dx,word y_point_
int 10h
ret
#
}
/* convert d to int with rounding */
round(d)
double d;{
if(d>=0) d+=.5;
else d-=.5;
return((int)d);
}
/* check to see if filename in bi exists */
exists(){
static i;
#asm
mov dx, offset fname_ ;pointer to filename
mov exists_i_,0 ;init to not exist
mov ax,4300h ;get file mode (to see if exists)
int 21h ;function call
jc no_exists ;if file doesn't exist
mov exists_i_,1 ;if file does exist
no_exists:
#
return(i);
}