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<<m)
            inc2=a2<<1
            FOR k=0 TO n
             PRINT "P(";k/n;")^2=";value/n^2
             value+=inc1
             inc1+=inc2
            NEXT


----[ Higher order ]------------------------------------------------------------

  Higher order polynomials can be computed similarly, though i've never needed
to go beyond the 3rd degree. Here's the program for cubic polynomials: observe
that in (k+1)^3=k^3+3*k^2+3*k+1 the k^3 comes from previous iteration and other
terms can also be computed incrementally since we just gave a method to compute
k^2 incrementally. Then you do the whole work again, and considering the
polynomial P(x)=a3*x^3+a2*x^2+a1*x+a0 the code now is:
            n=2^m
            value=a0<<(3*m)
            inc1=a3+(a2<<m)+(a1<<(2*m))
            inc2=6*a3+(2*a2<<m)
            inc3=6*a3
            FOR k=0 TO n
             PRINT "P(";k/n;")=";value/n^3
             value+=inc1
             inc1+=inc2
             inc2+=inc3
            NEXT


----[ Go forth and multiply ]---------------------------------------------------

  This method allows to compute successive values of a polynomial expression
with adds only, so at a very low processor power cost. Here's a small list of
routines i have designed and using this trick:
            * Circle routine
            * Iterative Bezier curves drawing
            * (bi)quadratic & (bi)cubic interpolation

  The last routines have not been seen frequently and gives much better result
than stupid linear interpolation, at approximately the same expense (except for
the number of registers).