This file lists known or suspected bugs that were reported or discovered after the deadline for the current release of Numerical Recipes. Not all the reports listed here are fully validated, so this listing should not be relied on as definitive. All the entries here will be further investigated before the next release. However, users of the current release who encounter bugs may wish to see if their bugs are already in this listing and, if so, whether they have additional information that may be useful for a fix. (If so, we encourage email reports to "bugs@nr.com".)
Note that this file does not include bugs already fixed in the current release. If you want information on those, look at these upgrade patch files.
Version 2.10
No new bugs yet reported!
Version 2.08: These bugs have all been fixed in Version 2.10.
1. In Fortran 77 caldat and julday change all the real constants to double by appending "d0". In the Fortran 90 version, change all _sp's to _dp's. Fixes a roundoff problem present on some machines. No dates plus or minus 300 years of the present are affected. No change in C version is required.
2. In Fortran 77 and 90 julday change the obvious line to read
julday=365*jy+int(0.25*jy+2000.)+int(30.6001*jm)+id+1718995
No change in C version is required. Bug affects only some B.C. dates.
3. For icrc() in C, documentation should state clearly that the string to be checksummed is in locations bufptr[1..len], not bufptr[0..len-1].
4. In the main books, equations 18.3.14 and 18.3.16 are missing a minus sign before the ln. The programs and figure are correct.
5. In the C version of dftcor the obvious line should be
if (a >= b || th < 0.0e0 || th > 3.1416e0) ... etc.
6. medfit encounters a divide by zero if the input
vectors are a perfect straight line. The fix is to put a line
if (sigb>0.0) {
just after the line f1=rofunc(b1), after sigb is computed,
and then to put the corresponding line
}
just before the line *a=aa;
7. In mpinv.f90 (Fortran 90 version only), some
Windows distribution media contain the incorrect line
allocate(rr(max(n+m)+n+1),s(n))
instead of the intended and correct
allocate(rr(max(n,m)+n+1),s(n))
8. In the vector version of bessjy.f90
(Fortran 90 version only) the line
where (h < FPMIN) h=FPMIN
in bessjy_xltxmin should be deleted. (It's already appeared,
correctly, in the calling routine bessjy_v.) The bug causes
many, but not all, values to be incorrect for x < 2.
9. shoot and shootf in the Fortran 90 version (ONLY)
can give a runtime error because of the statement
if (associated(xp)) deallocate(xp,yp)
The associated function is called even if xp is not defined,
which is an error. The fix is to move the statement
nullify(xp,yp)
in odeint.f90 to just before the first occurrence of the statement
if (save_steps) then
10. In rj.c, the lines
if (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(x+y,x+z),FMIN(y+z,fabs(p))) < TINY
|| FMAX(FMAX(x,y),FMAX(z,fabs(p))) > BIG)
nrerror("invalid arguments in rj");
should be
if (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(FMIN(x+y,x+z),y+z),fabs(p)) < TINY
|| FMAX(FMAX(FMAX(x,y),z),fabs(p)) > BIG)
nrerror("invalid arguments in rj");
and similarly the lines
} while (FMAX(FMAX(fabs(delx),fabs(dely)),
FMAX(fabs(delz),fabs(delp))) > ERRTOL);
should be
} while (FMAX(FMAX(FMAX(fabs(delx),fabs(dely)),
fabs(delz)),fabs(delp)) > ERRTOL);
11. In indexx_i4b (Fortran 90 only) change do i=j-1,1,-1 to do i=j-1,l,-1 (i.e., change the limit from "one" to "ell"). There is no corresponding change needed in indexx.for or indexx.c. However, the routines iindexx.for and iindex.c (integer versions of indexx) need several changes to bring them into conformity with indexx.for and indexx.c. These are obvious from the respective codes.
12. Professor Christian Reinsch (reinsch@mathematik.tu-muenchen.de) has distributed the following error notice, which affects the hqr routine in Numerical Recipes (all versions, all languages):
The original Handbook routine HQR for the computation of all eigenvalues
of a real, non-symmetric matrix by Francis' method of double-QR-steps
contains an error. Since this routine has also been included in the
Eispack-software and in the "Numerical Recipes", you are kindly asked to
help inform potential users.
The error occurs only in rare situations: If a subdiagonal element is
negligible in the sense of routine HQR, viz.
As a correction, it is therefore suggested to explicitly reset to zero
a subdiagonal entry which satisfies (*). This makes the breakpoint permanent
and the splitting feasible.
Another remedy would be to do the column modifications on the righthand block
[i:n] * [1:n], (as it is done in the Handbook routine HQR2).
The error was detected when a C-version of routine HQR was applied to a large
series of projectors A. Projectors have A*A = A,
so that their eigenvalues
must be from {0, 1} and in linear elementary divisors (1*1 Jordan blocks).
For example, the 5*5 lower Hessenberg matrix
C. Reinsch, Munich, Feb 16, 2000
|H[i,i-1]| <= epsmach * (|H[i-1,i-1]| + |H[i,i]|) (*)
then the lower, righthand block [i:n]*[i:n] alone
undergoes the next
transformation. This saves computational labor. However, the new
H[i,i]
could be smaller, so that at the beginning of a subsequent iteration, the
criterion (*) need no longer be satisfied. (The other two entries,
H[i-1,i-1] and H[i,i-1],
are not changed.) In such a case, a larger block
(or even the entire matrix) will be treated, and matrix entries come back
into play which are invalid since they missed previous QR steps.
7 0 0 0 0 1 7 49 0 0
-5 8 4 0 0 1 7 -135 0 2
A = 10 -2 -1 0 0 X = 1 7 95 Y = 0 -4
-10 -2 3 8 8 1 -33 -1 1 1
5 3 1 -1 -1 1 12 -8 -1 1
has A*X = 7*X and A*Y = 0.
The routine HQR in IEEE double precision
would deliver +-O( 3e-8) instead of two times zero or O(1e-15).
dist(span(X),span(Y)) = 0.368...
shows that this matrix is well-conditioned.