lib/cdhc: correct C translation of Cdhc_alnfac() from original Fortra…

…n source (#3873)

The C translation of the Fortran alnfac() function was flawed,
originated from confusion of Fortran 1-based array index and
C's 0-based. An incorrect constant (A0) is also updated. This
corrects the C code to mirror the original Fortran code
(https://lib.stat.cmu.edu/apstat/177).

Source: Royston, J. P. (1982). Algorithm AS 177: Expected Normal
Order Statistics (Exact and Approximate). Journal of the Royal
Statistical Society. Series C (Applied Statistics), 31(2), 161–165.
https://doi.org/10.2307/2347982

lib/cdhc/as177.c

@@ -82,20 +82,20 @@ void init(double work[])
  */
 static double Cdhc_alnfac(int j)
 {
-    static double r[7] = {0.0,           0.0,           0.69314718056,
-                          1.79175946923, 3.17805383035, 4.78749174278,
-                          6.57925121101};
+    static const double r[7] = {0.0,           0.0,           0.69314718056,
+                                1.79175946923, 3.17805383035, 4.78749174278,
+                                6.57925121101};
     double w, z;
 
-    if (j == 1)
-        return (double)1.0;
-    else if (j <= 7)
+    if (j < 0)
+        return 1.0;
+    if (j < 7)
         return r[j];
 
     w = (double)j + 1;
     z = 1.0 / (w * w);
 
-    return (w - 0.5) * log(w) - w + 0.918938522305 +
+    return (w - 0.5) * log(w) - w + 0.918938533205 +
            (((4.0 - 3.0 * z) * z - 14.0) * z + 420.0) / (5040.0 * w);
 }