Sierpinski (<1K) Recurrance Relations
[Needs typographic polish]
    Let us suppose that we have an object travelling along the X-axis with velocity v starting with x=x0 at time t=0.
At time t the object will have position x = x0 + vt.
Now suppose the the object is not travelling with constant velocity but is accelerating with acceleration a from an initial velocity u so that the its velocity at time t is given by v = u + at.
Its position is then given by x = (u+v)t/2 = ut + at2/2, a quadratic in t.

    If we split the time interval (0,T) into N segments of length l =T/N and iterate
     xnew = xold + lvold ;   vnew = vold + la.

or more formally
     xn+1 = xn + lvn ;  vn+1 = vn + la
N times what will we get? This is an example of a  recurrance relation. Recurrance relations describe a sequence of terms by specifying a given term as an expression involving previous terms in the sequence. The Fibonacci numbers, for example,   1,1,2,3,5,8,13,... are uniquely specified by the recurrance relation xn+1 = xn + xn-1 together with "initial condition" x1 = x0 = 1.

    The advantage of recurrance relations for the programmer is that the position of an object at time t=nl, though theoretically given by a quadratic in nl, can be derived without multiplication from xn-1 (its position last game cycle) and the quantity lv which (if l is sensibly chosen, often l = 1) is easily calculated.
The movement iteration above in fact gives
      xn = x0 + nlu + n(n-1)l2a/2      ;      vn = v0 + nla
since vn is readily solved and then xn follows. We see from this that xn is not the correct position of the body at time t=nl. The correct value is x0 + nlu + n2l2a/2 but this differs by only nl2a/2 which is usually negligible. If the position is required to be exactly correct then there are two ways to achieve this:
(i) Iterate as above but start with v0 = u + la/2.
(ii) Use the following slightly modified iteration:
 vtemp = vold + la/2
 xnew = xold + lvtemp
 vnew = vtemp + la/2
These iterations extend readily to two or three dimensional motion, the position, velocity, and acceleartion terms becoming vectors.

Solving Recurrance Relations
To solve a recurrance relation of the form  g(xn+k,...,xn+1,xn) = f(n) for some linear function g and arbitary function f we find the general solution to g(xn+k,...,xn+1,xn) = 0 and then add to it any specific solution to g(xn+k,...,xn+1,xn) = f(n), which is often a function of similar form to f(n).
The general solution to axn+1 + bxn = 0 (a¹0) is xn= A (-b/a)n.
The form of the general solution to axn+2 + bxn+1 +cxn = 0 (a¹0) depends on the nature of the roots of the  auxiliary equation am2+bm+c=0.
If the roots are repeated (b2=4ac) then the general solution is xn = (A+Bn)mn where m is the repeated root -b/2a.
If the roots are real and distinct (b2>4ac) then the general solution is xn=Am1n + Bm2n where m1 and m2 are the two roots (-b ±  Ö (b2-4ac))/2a.
If the roots are imaginary (b2<4ac) then the general solution is xn = (c/a)n/2(Acos(nq)+Bsin(nq)) or equivalently xn = (c/a)n/2Acos(nq+e) where tanq = - Ö(4ac-b2)/b.

Example - Damped Motion
It is often desirable to damp the motion of an object each game cycle so that objects come to rest if not explicitly accelerated. We will consider motion in N dimensions and write r for the position vector of the point under consideration at time t, v for its velocity vector (the first derivative of r with respect to time) at time t, and a for its constant applied acceleration vector.
We might wish to simulate motion satisfying dv/dt = -sv+a where s is the damping term.

    This  differential equation has solution
r = r0 + (1-e-st)(v0-a/s)/s  + ta/s
v = (v0-a/s)e-st + a/s
for constant a and we see that in a unit time interval v is decreased by a damping factor of e-s. As t®¥ these approximate
r = ta/s + (r0 + (v0-a/s))
v = a/s
We would typically simulate this with the iteration
  rn+1 = rn + lvn
vn+1 = dvn + la
where l is time step simulated by one iteration cycle and d = e-sl is the damping factor for one iteration period.
We will consider here the more general iteration
rn+1 = rn + l1vn
vn+1 = dvn + l2a.
(by taking l1=l, l2=dl, this covers the alternative iteration rn+1 = rn + lvn  ;   vn+1 = d(vn + la) ).

    To find the explicit forms for this iteration we first find vn as the solution of vn+1 - dvn = l2a. The general solution is thus Adn where A is an arbitary vector. The form of f(n) is a constant vector so we look for a constant specific solution. That is, we try to solve  x-dx = l2a and obtain x=l2a/(1-d). Our full solution is therefore vn = Adn + l2a/(1-d) and it remains only to  satisfy the initial conditions, ie. to find a value for A that gives the desired value for v0 .
Once vn is solved we use it to find rn using results A1.1 and A1.9.
The results are
l1(1-dn) l1l2 æ (1-dn) ö
rn=r0+¾¾¾¾v0+¾¾ ç n-¾¾ ÷ a
(1-d) (1-d) è (1-d) ø


which for large n approximates
l1l2 æ l2 l1l2 ö
rn=n¾¾a+ ç r0+¾¾v0-¾¾a ÷
(1-d) è (1-d) (1-d) ø


so to approximate the behaviour of the continuous solution for large n we need l1=l2 and (1-d)/l2 = s .

Successive Evaluations of Polynomials Via Iteration
Evaluating the position of a body moving with constant acceleration using iteration is an example of the evaluation of polynomials via recurrence relations - sometimes known as  the method of differences. Since such motion is a quadratic equation (a polynomial of degree two) in t, we keep track of two counts (xn and vn).
In general if we iterate M+1 counts x0n, x1n,...xMn thus:
x0n = x00  " n ³ 0
x1n = x1(n-1) + lx00 " n ³ 1
x2n = x2(n-1) + lx1(n-1) " n ³ 1
x3n = x3(n-1) + lx2(n-1) " n ³ 1
xMn = xM(n-1) + lx(M-1)(n-1) " n ³ 1
then we obtain
x0n = x00
x1n = x10 + lnx00
x2n = x20 + lnx10 + l2n(n-1)x00/2
x3n = x30 + ....... + (ln)3x00/6
xMn = xM0 + ....... + (ln)Mx00/M!.

xmn is a polynomial of degree m in (ln) whose coefficients are determined by the initial values x00,x10,...xm0.
Though we could determine these coefficients in terms of the initial values, it is more common to want to know what initial values must be used in order to get specific coefficients. If, for example, we want x2n = a(nl)2+b(nl)+c then we set x20=2a, x10=b+la, x00=c.  These values were obtained by inspection of the formula for x2n above but we can obtain the initial values for polynomials of greater degree by  back calculation from the the first N desired values.

Example - Iteration of Bezier Curves
We are interseted in the cubic a(nl)3+b(nl)2+c(nl)+d for nl ranging from zero to one. We evaluate it for n=0,1,2,3 to give x30,x31,x32,x33 and work backwards to get x30,x20,x10, shown in the table below.
0dal2+ bl+c6al+2b6a
1al3+ bl2+ cl+d7al2+3bl +c12al +2b

In the above table, the left column contains the desired values for x3n. The other entries are successively back calculated to derive the initial count values in the top row. x22, for example, was calculated as (x33-x32)/l.

    Each stage of the iterative process requires M additions of the form x=y+lz. If l is a "friendly" value such as 2-p this is okay but if l is simply 1/N for some general N we would rather not incorporate it into the iteration if at all possible. Suppose we wish to evaluate the X­coordinate of a bezier curve cubic x(t)=at3+bt2+ct+d for t=n/N " n Î ZN.
To do this we iterate
x0n' = x00   " n ³ 0
x1n' = x1(n-1)' + x00' " n ³ 1
x2n' = x2(n-1)' + x1(n-1)' " n ³ 1
x3n' = x3(n-1)' + x2(n-1)' " n ³ 1

which will work if x3n' = x3n   ; x2n' = lx2n ; x1n' = l2x1n ; x0n' = l3x0n  " n ³ 1 so we start from the initial values
  x00' = l3x00, x10' = l2x10, x20' = lx20, x30' = x30 where l =1/N .

    The price of loosing l from the iteration is that the initial count values become very small for small l but any inaccuracies in them will ripple up into increasing inaccuracies in x3n'.

    This method extends readily to polynomials of higher degree and can be used to fit polynomials to given data.

Succesive Evaluations Via The Midpoint Algorithm
Suppose we wish to iterate a function u=U(x) for x=0,1,2,.... and that the gradient u'=du/dx is known to be DU - d(x) over a given subrange of x, where DU is a constant integer and d(x) is known to lie between 0 and 1.
    DU is sometimes refered to as the  "integral gear" of the iteration.
    we know that Int[u(x+1)] will be either Int[u(x)]+DU-1 or Int[u(x)]+DU so we can construct iteratively from x and Int[u(x)] .
    Given two candidate points (x+1,u+DU-1) and (x+1,u+DU). we decide which of these is closest to the true curve { (x,u) : u=U(x) } by rewriting u=U(x) in the form f(x,u)=0 and considering  the sign of a  "decision function"  g(x,u) = mf(x+1,u+DU-1+0.5) where m is   a positive integer of convenience often chosen to make g(x,u) integral.
    g corresponds to f evaluated half way between the two candidate points.
    If g ³ 0 we select (x+1,u+DU) and if g < 0 we select (x+1,u+DU-1). The idea is of course to iterate g rather than evaluate it anew for each x,u pair.

Example - Division-Free 2D Line Draw
u(x)(DY/DX)x + c where 0<DY<DX
f(x,u)=uDX - xDY - cDX
g(x,u)=2f(x+1,u+DU-0.5) where DU=1
 =2uDX - 2xDY + DX -2DY -2cDX

    If g(x,u) <0 then we move to x+1,u. The next g is then given by g(x+1,u) = g(x,u)-2DY.
    If g(x,u) >0 then we instead move to x+1,u+1. The next g is then given by g(x+1,u+1) = g(x,u)+2DX-2DY.

Example - Division-Free Perspective Texture Mapping
[As proposed by John Mears.]
f(x,u)=Cxu + Du - Ax -B

    If g(x,u) <0 then we move to x+1,u+DU-1. The next g is then given by
    g(x+1,u+DU-1) = g(x,u) + k1(x,u)     where k1(u,x) = 2Cx(DU-1) + 2Cu + C(6DU-5) + 2D(DU-1) - 2A is maintained by iteration.
[ Proof :  
=2(C(x+1)(u+2DU-1.5)+D(u+2DU-1.5)-Ax-B) + 2(C(u+2DU-1.5)-A)
=2(C(x+1)(u+DU-0.5)+D(u+DU-0.5)-Ax-B) + 2(C(x+1)(DU-1) +D(DU-1)) + 2(C(u+2DU-1.5)-A)
=g(x,u) + 2(C(x+1)(DU-1) + D(DU-1) + C(u+2DU-1.5)-A)
=g(x,u) + 2Cx(DU-1) + C(6DU-5) + 2D(DU-1) + 2Cu -2A
=g(x,u) + k1(x,u)


    If g(x,u) >0 then we instead move to x+1,u+DU. The next g is then given by
    g(x+1,u+DU) = g(x,u) + k2(x,u)     where k2(x,u) =  2CxDU + 2Cu + C(6DU-1) + 2DDU - 2A is maintained by iteration.
[ Proof :
=2(C(x+1)(u+2DU-0.5)+D(u+2DU-0.5)-Ax-B) + 2(C(u+2DU-0.5)-A)
=2(C(x+1)(u+DU-0.5)+D(u+DU-0.5)-Ax-B) + 2(C(x+1)DU +DDU) + 2(C(u+2DU-0.5)-A)
=g(x,u) + 2(C(x+1)DU + DDU + C(u+2DU-0.5)-A)
=g(x,u) + 2(CxDU + 3CDU + DDU + Cu - C/2 -A)
=g(x,u) + 2CxDU + 2Cu C(6DU-1) + 2DDU - 2A
=g(x,u) + k2(x,u)


    Note that k1(x,u) = k2(x,u)-2Cx-4C-2D

    At a given stage in the iteration we consider g(x,u).
If g(x,u) < 0 we move to x+1,u'=u+DU-1 setting
k1(x+1,u')=k1(x,u) + 4C(DU-1)
k2(x+1,u')=k2(x,u) + 4CDU - 2C

If g(x,u) ³ 0 we move to x+1,u'=u+DU setting
k1(x+1,u')=k1(x,u)+4CDU - 2C

    It only remains to detect when and how the "gear" DU needs to be changed.  
u'(x) = du/dx = (A(Cx+D)-(Ax+B)C)/(Cx+D)2 = (AD-BC)/(Cx+D)2.
u"(x) = d2u/dx2 = -2C(AD-BC)/(Cx+D)3
Thus we see that u' is either always >0, always <0, or always =0 and that it increases while (Cx+D)<0.
 Our two candidate points (x+1,u+DU-1) and (x+1,u+DU) are only the correct ones if g takes opposite signs when evaluated at them (or is zero at one of them).
Thus the gear needs increasing if both g+k1 and g+k2 are < 0, and decreasing if both are > 0.  
Gear changes will be upwards while Cx+D<0.

    Mears asserts that if 4CDU - 2C > 0 then all subsequent gear changes will be upwards. This is incorrect.

Inverse Functions Via Iterative Approach
In Square Roots the iteration xn+1 = (xn + y/xn)/2 was shown to approach Öy extremely rapidly.  This is one example of the general iteration xn+1 = xn + (y-f(xn))/f '(xn) which tends to x¥=f -1(y). Infact, if dn=xn-x¥ then dn+1=f"(x)dn2/2f'(x) + O(dn3), and so the convergence is often extremely rapid.
The following table shows iterations converging to integral roots, reciprocals, and arctangents.
f(x)xn+1x ¥ dn+1
xm ((m-1)x+ y/xnm-1)/m  mÖy  (m-1)dn2/2x¥
p/x 2xn - yxn2/p p/y -ydn2/p
tan(x)  xn + cos2xn(y-tanxn)  tan-1(y) -sinx¥dn2
tan(x) ((2-sin2xn)xn + (1+cos2xn)y)/2  tan-1(y) -sinx¥dn2

General Functions Via Iterative Approach
By inspecting the  Taylor Series   f(a+h) =åi=0¥ f(i)(a)hi/(i!) for a function f(x) it is often possible to deduce an iteration converging to f(a+h) when |h|<1. Such iterations usually have xn =åi=0n f(i)(a)hi/(i!) which can be achieved crudely by iterating xn+1 = xn + Cn+1hn+1 where Cn = f(n)(a)/n! . More efficient iterations are often possible, sometimes of the form yn+1 = g(yn,n) ; xn+1 = xn + yn+1. These iterations give dn = xn-f(a+h) = Cn+1hn+1/Cn = f(n+1)(a)hn+1/(n+1)! + O(hn+2)
so that dn+1 »  (hCn+1/Cn)dn = (hf(n+1)(a)/nf(n)(a))dn. This method will therefore usually converge far less repidly than the Inverse Function Iteration described above; typically only giving one extra bit of accuracy per iteration.
f(x)xn+1dn dn+1
1/(1-x) 1 + xxn -xn+1 xdn
sinx -x2xn/2n(2n-1) (-1)nx2n+1/(2n+1)! -x2dn/2n(2n-1)

Circle Drawing Using Simultaneous Recurrance Relations
It is possible to draw circles by conducting two iterations in parallel. Suppose we wish to draw a circle having centre (0,0) and passing through the point (x0,y0).
Consider the iteration:
xn+1  = xn + ayn
yn+1  = yn + bxn+1.
To solve this we "expand out" xn+1 thus:
xn+1 = xn + a(yn-1+bxn)  = xn + abn + a(yn-2+bxn-1)
      = xn + ab(xn+xn-1+...+x1) + ay0

    Writing Xn for åi=1n xi we have Xn+1-Xn = Xn-Xn-1 + abXn + ay0, ie. Xn+1­(2+ab)Xn+Xn-1 = ay0 which is a second order recurrance relation having specific solution Xn = -y0/b and an auxilliary equation m2-(2+ab)m+1=0 having roots ((2+abÖ (ab(4+ab))/2.

    If ab(4+ab)<0 then we have Xn = Acos(nq)+Bsin(nq)-y0/b where q= tan-1( Öab(4+ab))/(2+ab)).
We then have
xn  = Xn-Xn-1 = A(cos(nq)-cos((n-1)q))+B(sin(nq)­sin((n-1)q))
= -2Asin(q/2)sin(nq)+2Bcos(q/2)cos(nq)
= Ccos(nq) + Dsin(nq).
From this general solution for xn we can obtain the specific solution in terms of x0 and y0. We can derive yn similarly.
xn = x0cos(nq) + [(x0+ay0-x0cosq)/sinq]sin(nq)
yn = y0cos(nq) + [(bx0+(1-ab-cos+q)y0)/sinq]sin(nq).

If a=-b and a is small then  q= tan-1( Ö (a2(4-a2))/(2-a2)) »  tan-1(a) so that sinq » a and cosq » 1-a2/2 then these become
xn = x0cos(nq) + y0sin(nq) + x0asin(nq)/2
yn = y0cos(nq) - x0sin(nq) - y0asin(nq)/2.

Typically we might have x0=R, y0=0 giving
xn = Rcos(nq) + Rasin(nq)/2
yn = -Rsin(nq).

This is our desired circle but with a small "error term" (x0asin(nq)/2,-y0asin(nq)/2). If the units used are screen pixels, then for this error term be less than half a pixel and so have no effect we need a <Min{1/x0,1/y0}.
Thus taking a < 1/R ensures that xn and yn never change by more than one.
To map out the full circle we need 2p/q » 2p/a iterations.
If the reader is wondering why the general result is asymmetric in x and y it is because the iteration
xn+1  = xn + ayn
yn+1  = yn + bxn+1
is itself asymmetric. The symmetric iteration
xn+1  = xn + ayn
yn+1  = yn + bxn
is far less useful, giving a spiral of increasing radius (if ab<0):
xn = (1-ab)n/2(x0cos(nq) + ((x0(1-cosq)+ay0)/sinq)sin(nq))
yn = (1-ab)n/2(y0cos(nq) + ((y0(1-cosq)+bx0)/sinq)sin(nq))
where q  =  tan-1(Ö(ab)).

Ellipse Drawing Using Circle Iterations
It is possible to generate a general (ie. non axi-aligned) ellipse without using any multiplications other than those implicit in two circle iterations of the type described above.
An axi-aligned ellipse having centre (0,0) has parametric form:
 X = acosq
 Y = bsinq
and we see immediately that by running two parallel circle iterations, one from (a,0) and one from (b,0) we can pick the X term from the X­coord of the first circle and the Y term from the Y­coord of the second.
This same ellipse rotated by an angle a about the origin has parametric form:
 X = a cosacosq ­ b sinasinq
 Y = a sinacosq + b cosasinq
which we can rewrite using A?.? as
 X = Acos(q+a) + Bcos(q-a)
 Y = Asin(q+a) - Bsin(q-a)
where A=(a+b)/2 ; B=(a-b)/2. So if we iterate one circle of radius A and another of radius B such that the second circle is 2a behind the first, then adding the X-coords of the two circles and subtracting the Y-coord of the second from that of the first will give, respectively, the X and Y­coords of the desired ellipes. One such pair of circles is obtained by starting the first from (A,0) and the second from (Bcos(2a),-Bsin(2a)); which corresponds to starting with q=-a.

    If we wish to start from the point on the ellipse having maximal Y­coordinate then we need to begin with q=tan-1((bcota)/a). This corresponds to starting the first circle at
((a+b)(a­b)sin2a/4C , (a+b)(bcos2a + asin2a)/2C)
and the second at
((a+b)(a­b)sin2a/4C , (a-b)(bcos2a-asin2a)/2C)
where C= Ö (b2cos2a + a2sin2a).

    Note that q is  not the angle subtended with the major axis. That angle is given by tan-1( b *  tanq /a) .

Next : Bezier Curves

Glossary   Contents   Author
Copyright (c) Ian C G Bell 198
Web Source: or
18 Nov 2006.