``````{
SEAN PALMER

I was just toying around With a B-Spline curve routine I got out of an
old issue of Byte, and thought it was pretty neat. I changed it to use
fixed point fractions instead of Reals, and optimized it some...

by Sean Palmer
public domain
}

Var
color : Byte;
Procedure plot(x, y : Word);
begin
mem[\$A000 : y * 320 + x] := color;
end;

Type
coord = Record
x, y : Word;
end;

CurveDataRec = Array [0..65521 div sizeof(coord)] of coord;

Function fracMul(f, f2 : Word) : Word;
Inline(
\$58/                   {pop ax}
\$5B/                   {pop bx}
\$F7/\$E3/               {mul bx}
\$89/\$D0);              {mov ax,dx}

Function mul(f, f2 : Word) : LongInt;
Inline(
\$58/                   {pop ax}
\$5B/                   {pop bx}
\$F7/\$E3);              {mul bx}

Const
nSteps = 1 shl 8;  {about 8 For smoothness (dots), 4 For speed (lines)}

Procedure drawBSpline(Var d0 : coord; nPoints : Word);
Const
nsa  = \$10000 div 6;
nsb  = \$20000 div 3;
step = \$10000 div nSteps;
Var
i, xx, yy,
t1, t2, t3,
c1, c2, c3, c4 : Word;

d : curveDataRec Absolute d0;

begin
t1 := 0;
color := 32 + 2;

For i := 0 to nPoints - 4 do
begin

{algorithm converted from Steve Enns' original Basic subroutine}

Repeat
t2 := fracMul(t1, t1);
t3 := fracMul(t2, t1);
c1 := (Integer(t2 - t1) div 2) + nsa - fracmul(nsa, t3);
c2 := (t3 shr 1) + nsb - t2;
c3 := ((t2 + t1 - t3) shr 1) + nsa;
c4 := fracmul(nsa, t3);
xx := (mul(c1, d[i].x) + mul(c2, d[i + 1].x) +
mul(c3, d[i + 2].x) + mul(c4, d[i + 3].x)) shr 16;
yy := (mul(c1, d[i].y) + mul(c2, d[i + 1].y) +
mul(c3, d[i + 2].y) + mul(c4, d[i + 3].y)) shr 16;
plot(xx, yy);
inc(t1, step);
Until t1 = 0;  {this is why nSteps must be even power of 2}
inc(color);
end;
end;

Const
pts = 24; {number of points} {chose this because of colors}

Var
c : Array [-1..2 + pts] of coord;
i : Integer;
begin
Asm
mov ax, \$13
int \$10
end;  {init vga/mcga Graphics}
randomize;
For i := 1 to pts do
With c[i] do
begin
{x:=i*(319 div pts);}    {for precision demo}
x := random(320);               {for fun demo}
y := random(200);
end;
{for i:=1 to pts div 2 do c[i*2+1].y:=c[i*2].y;}    {fit closer}
For i := 1 to pts do
With c[i] do
begin
color := i + 32;
plot(x, y);
end;
{replicate end points so curves fit to input}
c[-1] := c;
c  := c;
c[pts + 1] := c[pts];
c[pts + 2] := c[pts];
drawBSpline(c[-1], pts + 4);