04

Drawing Circles

I had a friend who had Windows 3.11, and it included some screensavers (or maybe one from After Dark) that I wanted to try to reproduce. One of them drew lots of circles on the screen, but when I tried to code this myself it was insanely slow, and I was determined to solve it.

I was completely oblivious to the optimal algorithms that had been developed decades earlier and had no internet to consult, so I went about it using the math available in my head.

The obvious starting place was the Pythagorean Theorem. For each row, I knew I could a2 + b2 = c2 my way to place the expected point. It worked, but it was insanely slow.

My first clever idea was to calculate just one quadrant and mirror it to the others. Soon after I realized I could get away with half of that, mirroring one octant to the rest of the circle. It was faster, but not fast enough.

I was stuck for a while, until one morning I woke up from a dream where I had solved it with calculus. I was able to replace almost all of the square roots with a differential approach. This got me fast enough for what I wanted, and I was convinced I had perfected it.

I always had this nag in the back of my head wondering if I had invented something novel back at age 17, and this was part of what inspired me to come dig out these files. Of course, while porting it to JavaScript, Claude quickly informed me that Bresenham had an even better solution in 1962.

methoditerationssqrtdividesmax radial error
Costs are for one filled circle of radius 180. The 1995 code ran on a 386 with Turbo Pascal's 6-byte software real; no FPU, so every sqrt and every divide was a subroutine call.

The differential circle

This is the one piece of real mathematics invented here. To walk a circle you normally need trigonometry per point, or Bresenham's integer error term. This method needs neither:

a := a - c/a

One sqrt for the seed, one divide per point, no trig at all. It falls out of differentiating x² + y² = r², and it is stable only over one octant, which is why the octant mirroring had to come first.

The most-reused code here

This routine appears in eight files: TIME, CIRCLE, FUN, FUN1, FUN1.BAK, FUN2, ALL and ALL.BAK - verbatim in all but FUN2, where the radius argument is gone and the loop runs for c:=0 to 1, so the same code draws a three-pixel blob. Nothing else in the archive travels like that.

TIME2 is measuring something else entirely

This one was looking for the fastest way to draw a point, and the version in the archive is testing with 172,800 calls to the library's set_point routine, timed end to end. The file also contains an alternative procedure p(x,y,c:integer); begin Mem[$a000:480*y+x]:=c; end; that I had likely also tested in a version that did not survive.

That is the direct-memory alternative to set_point; write the byte yourself and skip the library. It would not have worked as written: the addressing is linear, which is mode 13h thinking, and in Mode X the four planes are selected through a port. 480*y+x also runs past 64 KB before y reaches 140. But it is the right instinct, two months early - by December the FUN family is writing to $A000 directly and never calls set_point again.

CIRCLE.PAS

1995-07-23 12:30 · 690 bytes · One procedure, circ, drawn once and timed.

uses crt,jmodex,dos;var q:integer;d,e,f,g:word;

procedure circ(x,y,r,l:integer);var b,c:integer;a:real;
begin if r=0 then b:=0 else b:=round(sqrt(r*r-1));a:=b;
for c:=0to round(r/sqrt(2))do begin if a=0 then a:=1;
a:=(a-(c/a));b:=round(a);draw_line(x+b,y+c,x-b,y+c,l);
draw_line(x+b,y-c,x-b,y-c,l);draw_line(x+c,y+b,x-c,y+b,l);
draw_line(x+c,y-b,x-c,y-b,l);end;end;

begin if set_vga_modex(7,360,480,1)=0 then
halt(0);clear_vga_screen(0);for q:=1 to 64 do set_dac_register(q,q,0,0);
gettime(d,d,d,e);
circ(180,240,180,63);
gettime(f,f,f,g);
repeat until keypressed;asm mov ah,0;mov al,3;int 16;end;
writeln(d,'.',e,' to ',f,'.',g);writeln((f-d)*100+g-e);readln;readln;end.
TIME.PAS

1995-07-22 06:09 · 595 bytes · The benchmark.

uses crt,jmodex,dos;var d,e,f,g:word;

procedure circ(x,y,r,l:integer);var b,c:integer;a:real;
begin b:=round(sqrt(r*r-1));a:=b;for c:=0to round(r/sqrt(2))do begin
a:=(a-(c/a));b:=round(a);draw_line(x+b,y+c,x-b,y+c,l);
draw_line(x+b,y-c,x-b,y-c,l);draw_line(x+c,y+b,x-c,y+b,l);
draw_line(x+c,y-b,x-c,y-b,l);end;end;

begin if set_vga_modex(7,360,480,1)=0 then halt(0);clear_vga_screen(0);
gettime(d,d,d,e);
circ(180,240,180,30);
gettime(f,f,f,g);
repeat until keypressed;asm mov ah,0;mov al,3;int 16;end;
writeln(d,'.',e,' to ',f,'.',g);writeln((f-d)*100+g-e);readln;readln;
end.
TIME2.PAS

1995-09-16 16:41 · 433 bytes · A different benchmark: 172,800 set_point calls, timed.

uses crt,jmodex,dos;var a,b,c,d,e,f,g,h:word;X,Y:integer;

procedure p(x,y,c:integer);
begin
Mem[$a000:480*y+x]:=c;
end;


begin if set_vga_modex(7,360,480,1)=0 then halt(0);clear_vga_screen(0);
gettime(a,b,c,d);
for x:=1to 360do for y:=1to 480do set_point(x,y,x);

gettime(e,f,g,h);
repeat until keypressed;asm mov ah,0;mov al,3;int 16;end;
writeln(c,'.',d,' to ',g,'.',h);writeln((g-c)*100+h-d);readln;readln;
end.