🚀 go-pugleaf

On Wed, 18 Nov 1998 16:08:27 -0500, Ken Laberteaux <klaberte@trc.tellabs.com> wrote:
>I am reading the 2nd Edition of _Numerical Recipes in C_.  I am
>specifically looking at laguer.c.
>
>On p. 373 the line of code which implements Eq 9.5.7 (p. 372) is 
>
>               h=Csub(g2,RCmul(2.0,Cdiv(f,b)));
>
>However Eq 9.5.7 (p. 372) makes no reference of the 2, indicating that
>the correct version should be
>
>               h=Csub(g2,Cdiv(f,b));
>
>Can anyone shed some light into this inconsistency?

There is no inconsistency.  The recurrence relation used to calculate
f in the loop is:

  /* a[i] defines the i-th coefficient of an m degree polynomial */
  /* let p(m) = a[m] * x^m + ... + a[0] */
  f[m] = 0;
  d[m] = 0;
  b[m] = a[m];
  for (j = m-1; j >= 0; j--) {
    f[j] = x * f[j+1] + d[j+1];
    d[j] = x * d[j+1] + b[j+1];
    b[j] = x * b[j+1] + a[j];
  }

You can prove the following loop invariant by induction:

For k = 1 ... m;

  d[m-k](x) = b[m-k](x)'
  f[m-k](x) = (1/2) * b[m-k](x)''

Proof:

  For k = 1;
    b[m-1](x) = x * b[m](x) + a[m-1]
              = a[m] * x + a[m-1];
    
    but,
    d[m-1](x) = x * d[m](x) + b[m](x)
              = a[m]
	      = b[m](x)'
    
    and,
    f[m-1](x) = x * f[m](x) + d[m]
              = 0
	      = (1/2) * b[m-1](x)''

  Therefore, there exists some n > 1 such that hypothesis is true for
  all numbers greater than or equal to 1 but less than n. We will no
  show the hypothesis is true for k = n,

  For k = n;
    b[m-n](x) = x * b[m-(n-1)](x) + a[m-n]

    Thus,
    b[m-n](x)' = x * b[m-(n-1)](x)' + b[m-(n-1)](x)
    
    and
    b[m-n](x)'' = x * b[m-(n-1)](x)'' + 2 * b[m-(n-1)](x)'
    
    But,
    d[m-n](x) = x * d[m-(n-1)](x) + b[m-(n-1)](x)
              = x * b[m-(n-1)](x)' + b[m-(n-1)](x)
   
    which follows from the induction hypothesis.  Thus,
    
    d[m-n](x) = b[m-n](x)'
    
    Also,
    f[m-n](x) = x * f[m-(n-1)] + d[m-(n-1)](x)
              = x * (1/2) * b[m-(n-1)](x)'' + b[m-(n-1)](x)'
	      
    which follows from the induction hypothesis.  Thus,
    
    f[m-n](x) = (1/2) * b[m-n](x)''

  Thus, the hypothesis is true for all n >= 1.

-- 
James C. Hu <jxh@cs.wustl.edu>       Computer Science Doctoral Candidate
http://www.cs.wustl.edu/~jxh/        Washington University in Saint Louis
>>>>>>>>>>>>> I use *SpamBeGone* <URL:http://www.internz.com/SpamBeGone/>