Sign, Sign, Everywhere a Sign
Reference [1] describes a numerical method for computing the derivative of a 1D function, using Neville's algorithm to extrapolate from a sequence of simple polynomial approximations based on interpolating points within specific bounds of a given point. The original code in the published paper in 1980 was written in ALGOL 60. It seems as if a Fortran 77 translation was produced by David Kahaner at NIST circa 1989. I found the file on this NIST server where it has been sitting quietly since 1992 (originally this was an ftp server when I first found it). It's an interesting algorithm, and I included a modernized version in my NumDiff numerical differentiation library. I also note that a version of the code is embedded within the Dataplot package.
But, there is a problem!
Here is a line of code in the original faccuracy function:
IF 16*ABS(H1)>ABS(H0) THEN H1:=SIGN(H1)*ABS(H0)/16;
Now, I know nothing about ALGOL 60, except that apparently John Backus (the creator of Fortran) was one of the committee members that developed it. Notice the SIGN function in this code.
This line was translated into Fortran 77 (where the function is renamed as FACCUR) as:
IF(16.*ABS(H1) .GT. ABS(H0)) H1 = SIGN(H1,1.)*ABS(H0)/16.
They are very similar, but note the SIGN function again.
In Fortran, the SIGN function (added in Fortran 77) is a little different since it has two arguments. While more flexible, it is a never-ending source of confusion, and I always have to look up the meaning of the two arguments. The following table summarizes the difference between the two:
| Language | Syntax | Description | Reference |
|---|---|---|---|
| ALGOL 60 | sign(E) |
the sign of the value of E (+1 for E>0, 0 for E=0, -1 for E<0) | [2] |
| Fortran | sign(A, B) |
the value of A with the sign of B | [3] |
So, do you see the bug? SIGN(H1,1.) doesn't return the same thing that SIGN(H1) did in the original code. SIGN(H1,1.) returns the value of H1 with the sign of 1.0, which is always +H1, not the original intent at all! It really should have been translated to SIGN(1.0,H1), to give the +1/-1 of the original code. So, this bug is over 30 years old. Technically, the 0 case is still not handled the same, but that's a degenerate case that would never happen or give meaningful results in this context. But to be totally faithful we would need to use a function like this:
elemental real(wp) function sign_algol(x)
real(wp),intent(in) :: x
sign_algol = merge(0.0_wp, sign(1.0_wp,x), x==0.0_wp)
end function sign_algol
Other Languages
Python doesn't have a built-in sign function, but the numpy one behaves exactly like the ALGOL one. The advantage of the Fortran one is that it does support signed zeros, which the Numpy one does not appear to support (for that you need to use copysign). C++ also has a copysign function that does the same thing.
References
- J. Oliver, "Algorithm 017: An algorithm for numerical differentiation of a function of one real variable", Journal of Computational and Applied Mathematics 6 (2) (1980) 145-160. Fortran 77 code from NIST.
- J. W. Backus, F. L. Bauer, J. Green, C. Katz, J. McCarthy, A. J. Perlis, H. Rutishauser, K. Samelson, B. Vauquois, J. H. Wegstein, A. van Wijngaarden, M. Woodger Editor: P. Naur, "Revised Report on the Algorithmic Language ALGOL 60", Communications of the ACM, Volume 6, Issue 1 Pages 1 - 17, 01 January 1963.
- SIGN — Sign copying function [gcc.gnu.org]
- Numerical Differentiation, degenerateconic.com, Dec 04, 2016
