🚀 go-pugleaf

I  have  the 1989 version  of   Numerical Recipees  in  C, both   in book  and
electronic form. I believe I have found a bug---or at least  a nuisance.  When
a matrix  contains some  extreme values,  it   is  possible  to get  an access
violation.

(Does anyone have a bug list of NRinC?)

Note: The code below is copyrighted, but since it is only a small  fragment, I
presume it is ok to post it to better explain the bug.

----------------

Description: On Chapter 2, page 70, the code looks like

  for (k=n;k>=1;k--) {
    for (its=1;its<=30;its++) {
      flag=1;
      for (l=k;l>=1;l--) {
        nm=l-1;
        if (fabs(rv1[l])+anorm == anorm) {
          flag=0;
          break;
        }
        if (fabs(w[nm])+anorm == anorm) break;
      }

# note that l ends with a value of 0 if the loop exits. Typically,
# rv1[l]=0, so the loop is exited early---no problem then. However,
# if the matrix is poorly setup (and I have had such a case), rv1[1]
# can become NaN (0/0), and the loop completes rather than breaks.
# I believe rv1[1] turned into NaN, since my matrix contained an IEEE
# -Inf.

# I presume this problem can be fixed if the condition reads l>=2,
# and nm=l-1 is moved into both the start and iteration.

      if (flag) {
        c=0.0;
        s=1.0;
        for (i=l;i<=k;i++) {

# note that i starts with a value of l. If l can be 0
# this means that the subsequent indices refer to 0, not 1.
# Since NRinC allocates all matrices/vectors from 1, not 0,
# this results in an access violation error---if you are lucky.

          f=s*rv1[i];
          if (fabs(f)+anorm != anorm) {
            g=w[i];
            h=PYTHAG(f,g);
            w[i]=h;
            h=1.0/h;
            c=g*h;
            s=(-f*h);
            for (j=1;j<=m;j++) {
              y=a[j][nm];
              z=a[j][i];
              a[j][nm]=y*c+z*s;
              a[j][i]=z*c-y*s;
            }
          }
        }
      }