BUG: amos/amos.h: fix translation bugs in binu and rati

include/xsf/amos/amos.h

@@ -2666,7 +2666,7 @@ namespace amos {
         dfnu = fnu + n - 1;
         if ((az <= 2.) || (az * az * 0.25 <= (dfnu + 1.0))) {
             /* GOTO 10 */
-            nw = seri(z, fnu, kode, n, cy, tol, elim, alim);
+            nw = seri(z, fnu, kode, nn, cy, tol, elim, alim);
             inw = std::abs(nw);
             nz += inw;
             nn -= inw;
@@ -2689,7 +2689,7 @@ namespace amos {
                 //
                 // MILLER ALGORITHM NORMALIZED BY THE SERIES
                 //
-                nw = mlri(z, fnu, kode, n, cy, tol);
+                nw = mlri(z, fnu, kode, nn, cy, tol);
                 if (nw < 0) {
                     nz = -1;
                     if (nw == -2) {
@@ -2706,7 +2706,7 @@ namespace amos {
                 //
                 // ASYMPTOTIC EXPANSION FOR LARGE Z
                 //
-                nw = asyi(z, fnu, kode, n, cy, rl, tol, elim, alim);
+                nw = asyi(z, fnu, kode, nn, cy, rl, tol, elim, alim);
                 if (nw < 0) {
                     nz = -1;
                     if (nw == -2) {
@@ -2760,7 +2760,7 @@ namespace amos {
         /* 60 */
         if (az <= rl) {
             /* 70 */
-            nw = mlri(z, fnu, kode, n, cy, tol);
+            nw = mlri(z, fnu, kode, nn, cy, tol);
             if (nw < 0) {
                 nz = -1;
                 if (nw == -2) {
@@ -4296,7 +4296,7 @@ namespace amos {
         }
 
         if (p1 == 0.) {
-            p1 = std::complex<double>(0, 0);
+            p1 = std::complex<double>(tol, tol);
         }
 
         cy[n - 1] = p2 / p1;
@@ -4309,7 +4309,7 @@ namespace amos {
         cdfnu = fnu * rz;
 
         for (i = 2; i < (n + 1); i++) {
-            pt = cdfnu + t1 * rz * cy[k];
+            pt = cdfnu + t1 * rz + cy[k];
             if (pt == 0.0) {
                 pt = std::complex<double>(tol, tol);
             }

tests/xsf_tests/test_amos.cpp

@@ -0,0 +1,51 @@
+#include "../testing_utils.h"
+#include <xsf/amos/amos.h>
+
+TEST_CASE("amos besj", "[amos][xsf_tests]") {
+    using std::complex;
+
+    const complex<double> z{70.0, -0.7};
+    const double fnu = 0.0;
+    const int kode = 1;
+    const int n = 15;
+
+    std::array<complex<double>, n> cy{};
+    int ierr = 0;
+
+    const int nz = xsf::amos::besj(z, fnu, kode, n, cy.data(), &ierr);
+
+    REQUIRE(ierr == 0);
+    REQUIRE(nz == 0);
+
+    // generated using:
+    // from mpmath import mp
+    // mp.dps = 50
+    // n = 15
+    // z = 70.0 - 0.7j
+    // for n in range(n):
+    //     result = mp.besselj(n, z)
+    //     print(f"complex<double>{{{result.real:.20f}, {result.imag:.20f}}},")
+    const std::array<complex<double>, n> ref = {
+        complex<double>{0.11908800062250928549, 0.00765768917663659615},
+        complex<double>{0.01289516192962112526, -0.07187585288729215033},
+        complex<double>{-0.11869907035957860579, -0.00970739567078706507},
+        complex<double>{-0.01967174120900008651, 0.07125337877045338558},
+        complex<double>{0.11695202149331912143, 0.01579735764816333185},
+        complex<double>{0.03301829774384231522, -0.06931450090245772537},
+        complex<double>{-0.11213658264393736068, -0.02565127481716347501},
+        complex<double>{-0.05219582054512202598, 0.06472536427721690399},
+        complex<double>{0.10156902456216669875, 0.03849067209178914920},
+        complex<double>{0.07532130636025452583, -0.05569624151237091633},
+        complex<double>{-0.08205942069459905224, -0.05261746672772313695},
+        complex<double>{-0.09861419022253167551, 0.04042975075154153429},
+        complex<double>{0.05094243624931449501, 0.06501278997661191646},
+        complex<double>{0.11585554323564695001, -0.01796723780085772449},
+        complex<double>{-0.00784795132702038573, -0.07125539059627898735},
+    };
+
+    for (int i = 0; i < n; ++i) {
+        const auto rel_error = xsf::extended_relative_error(cy[i], ref[i]);
+        CAPTURE(i, cy[i], ref[i], rel_error);
+        REQUIRE(rel_error <= 1e-14);
+    }
+}