Successive squares and cubes
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
by Alain Brobecker (aka baah/Arm's Tech)
(abrobecker@yahoo.com or rte de Dardagny; 01630 CHALLEX; FRANCE)
This trick i found in an iterative Bezier drawing routine by Jan Vlietinck
(Jan/BASS), provided without comments (as usual, isn't it Jan? ;).
----[ Successive squares ]------------------------------------------------------
Suppose we want to compute the successive squares: 0^2, 1^2, 2^2, 3^2...
Let's consider the formulae (k+1)^2=k^2+2*k+1. We know k^2 from the previous
iteration, and 2*k+1 can easily be computed by adding 2 at each iteration:
square=0
increment=1
FOR k=0 TO 15
PRINT k;"^2=";square
square+=increment
increment+=2
NEXT
Just fine, but squares of integer quantities are not that usefull, we prefer
something like 0^2, 0.1^2, 0.2^2, 0.3^2... Let "step" be the increment of the
quantity to square (0.1 in our example). To compute squares of k*step we notice
that ((k+1)*step)^2=(k*step)^2+(2*k+1)*step^2. Again (k*step)^2 comes from last
iteration and (2*k+1)*step^2 can be computed incrementally:
step=0.1
square=0
inc1=step^2
inc2=2*inc1
FOR k=0 TO 10
PRINT "(";k*step;")^2=";square
square+=inc1
inc1+=inc2
NEXT
----[ Practical implementation ]------------------------------------------------
Some algorithms (such as iterative Bezier curves) consider a quadratic
polynomial P(x)=a2*x^2+a1*x+a0 and need successive values of P(k/n) where
0<=k<=n. So we have P(0)=a0 and we remark that:
P((k+1)/n)=a2*((k+1)/n)^2+a1*(k+1)/n+a0
=a2*(k/n)^2+a1*k/n+a0+a2*(2*k+1)/n^2+a1/n
=P(k/n)+a2*(2*k+1)/n^2+a1/n
As usual we know P(k/n) and the rest can be computed incrementally, with
the following program:
value=a0
inc1=a2/n^2+a1/n
inc2=2*a2/n^2
FOR k=0 TO n
PRINT "P(";k/n;")^2=";value
value+=inc1
inc1+=inc2
NEXT
Now let's think a bit about a pure assembly version. To achieve high speed
we will of course use fixed point values to implement the 1/n and 1/n^2 stuff.
Moreover if we choose n=2^m we can perform EXACT calculii, because arithmetics
tells us the binary development of k/n^2 is then finite in base 2. So, in the
following program "value" contains P((k/n)^2)<<2*m, where the "<<" stands for
a logical shift left (LSL):
n=2^m
value=a0<<(2*m)
inc1=a2+(a1<