[flang] Lowering and runtime support for F08 transformational intrinsics: BESSEL_JN and BESSEL_YN

The runtime implementation uses the recurrence relations

`J(n-1, x) = (2.0 / x) * n * J(n, x) - J(n+1, x)`
`Y(n+1, x) = (2.0 / x) * n * Y(n, x) - Y(n-1, x)`

(see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).

Although the standard requires that `N1` and `N2` in `BESSEL_JN(N1, N2, x)`
and `BESSEL_YN(N1, N2, x)` be non-negative, this is not checked in the
runtime functions. This is in keeping with some other compilers which also
return some results when `N1` and/or `N2` are negative.

The special case for `x == 0` is  handled in different runtime functions
for each of `BESSEL_JN` and `BESSEL_YN`. The lowering code checks for this
case and inserts the checks and the appropriate runtime calls in FIR.

The existing tests for the two intrinsics was modified to keep the style
consistent with the additional lowering tests that were added.

2 files changed, 291 insertions, 0 deletions

diff --git a/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp b/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
index 37e13ccb06b5..d5884ae3febb 100644
--- a/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
+++ b/flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp
@@ -10,6 +10,92 @@
 #include "RuntimeCallTestBase.h"
 #include "gtest/gtest.h"
 
+void testGenBesselJn(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value x = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn1 = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn2 = builder.create<fir::UndefOp>(loc, realTy);
+  fir::runtime::genBesselJn(builder, loc, result, n1, n2, x, bn1, bn2);
+  checkCallOpFromResultBox(result, fctName, 6);
+}
+
+TEST_F(RuntimeCallTest, genBesselJnTest) {
+  testGenBesselJn(*firBuilder, f32Ty, "_FortranABesselJn_4");
+  testGenBesselJn(*firBuilder, f64Ty, "_FortranABesselJn_8");
+  testGenBesselJn(*firBuilder, f80Ty, "_FortranABesselJn_10");
+  testGenBesselJn(*firBuilder, f128Ty, "_FortranABesselJn_16");
+}
+
+void testGenBesselJnX0(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  fir::runtime::genBesselJnX0(builder, loc, realTy, result, n1, n2);
+  checkCallOpFromResultBox(result, fctName, 3);
+}
+
+TEST_F(RuntimeCallTest, genBesselJnX0Test) {
+  testGenBesselJnX0(*firBuilder, f32Ty, "_FortranABesselJnX0_4");
+  testGenBesselJnX0(*firBuilder, f64Ty, "_FortranABesselJnX0_8");
+  testGenBesselJnX0(*firBuilder, f80Ty, "_FortranABesselJnX0_10");
+  testGenBesselJnX0(*firBuilder, f128Ty, "_FortranABesselJnX0_16");
+}
+
+void testGenBesselYn(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value x = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn1 = builder.create<fir::UndefOp>(loc, realTy);
+  mlir::Value bn2 = builder.create<fir::UndefOp>(loc, realTy);
+  fir::runtime::genBesselYn(builder, loc, result, n1, n2, x, bn1, bn2);
+  checkCallOpFromResultBox(result, fctName, 6);
+}
+
+TEST_F(RuntimeCallTest, genBesselYnTest) {
+  testGenBesselYn(*firBuilder, f32Ty, "_FortranABesselYn_4");
+  testGenBesselYn(*firBuilder, f64Ty, "_FortranABesselYn_8");
+  testGenBesselYn(*firBuilder, f80Ty, "_FortranABesselYn_10");
+  testGenBesselYn(*firBuilder, f128Ty, "_FortranABesselYn_16");
+}
+
+void testGenBesselYnX0(
+    fir::FirOpBuilder &builder, mlir::Type realTy, llvm::StringRef fctName) {
+  mlir::Location loc = builder.getUnknownLoc();
+  mlir::Type i32Ty = builder.getIntegerType(32);
+  mlir::Type seqTy =
+      fir::SequenceType::get(fir::SequenceType::Shape(1, 10), realTy);
+  mlir::Value result = builder.create<fir::UndefOp>(loc, seqTy);
+  mlir::Value n1 = builder.create<fir::UndefOp>(loc, i32Ty);
+  mlir::Value n2 = builder.create<fir::UndefOp>(loc, i32Ty);
+  fir::runtime::genBesselYnX0(builder, loc, realTy, result, n1, n2);
+  checkCallOpFromResultBox(result, fctName, 3);
+}
+
+TEST_F(RuntimeCallTest, genBesselYnX0Test) {
+  testGenBesselYnX0(*firBuilder, f32Ty, "_FortranABesselYnX0_4");
+  testGenBesselYnX0(*firBuilder, f64Ty, "_FortranABesselYnX0_8");
+  testGenBesselYnX0(*firBuilder, f80Ty, "_FortranABesselYnX0_10");
+  testGenBesselYnX0(*firBuilder, f128Ty, "_FortranABesselYnX0_16");
+}
+
 TEST_F(RuntimeCallTest, genCshiftTest) {
   auto loc = firBuilder->getUnknownLoc();
   mlir::Type seqTy =
diff --git a/flang/unittests/Runtime/Transformational.cpp b/flang/unittests/Runtime/Transformational.cpp
index a4672369f782..70ab42416564 100644
--- a/flang/unittests/Runtime/Transformational.cpp
+++ b/flang/unittests/Runtime/Transformational.cpp
@@ -10,10 +10,215 @@
 #include "gtest/gtest.h"
 #include "tools.h"
 #include "flang/Runtime/type-code.h"
+#include <vector>
 
 using namespace Fortran::runtime;
 using Fortran::common::TypeCategory;
 
+template <int KIND>
+using BesselFuncType = std::function<void(Descriptor &, int32_t, int32_t,
+    CppTypeFor<TypeCategory::Real, KIND>, CppTypeFor<TypeCategory::Real, KIND>,
+    CppTypeFor<TypeCategory::Real, KIND>, const char *, int)>;
+
+template <int KIND>
+using BesselX0FuncType =
+    std::function<void(Descriptor &, int32_t, int32_t, const char *, int)>;
+
+template <int KIND>
+constexpr CppTypeFor<TypeCategory::Real, KIND>
+    besselEpsilon = CppTypeFor<TypeCategory::Real, KIND>(1e-4);
+
+template <int KIND>
+static void testBesselJn(BesselFuncType<KIND> rtFunc, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, KIND> x,
+    const std::vector<CppTypeFor<TypeCategory::Real, KIND>> &expected) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+  unsigned len = expected.size();
+
+  CppTypeFor<TypeCategory::Real, KIND> anc0 = len > 0 ? expected[len - 1] : 0.0;
+  CppTypeFor<TypeCategory::Real, KIND> anc1 = len > 1 ? expected[len - 2] : 0.0;
+
+  rtFunc(result, n1, n2, x, anc0, anc1, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(
+      result.ElementBytes(), sizeof(CppTypeFor<TypeCategory::Real, KIND>));
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), len);
+
+  for (size_t j{0}; j < len; ++j) {
+    EXPECT_NEAR(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        expected[j], besselEpsilon<KIND>);
+  }
+}
+
+template <int KIND>
+static void testBesselJnX0(
+    BesselX0FuncType<KIND> rtFunc, int32_t n1, int32_t n2) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+
+  rtFunc(result, n1, n2, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), n2 >= n1 ? n2 - n1 + 1 : 0);
+
+  if (n2 < n1) {
+    return;
+  }
+
+  EXPECT_NEAR(
+      (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+          0)),
+      (n1 == 0) ? 1.0 : 0.0, 1e-5);
+
+  for (int j{1}; j < (n2 - n1 + 1); ++j) {
+    EXPECT_NEAR(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        0.0, besselEpsilon<KIND>);
+  }
+}
+
+template <int KIND> static void testBesselJn(BesselFuncType<KIND> rtFunc) {
+  testBesselJn<KIND>(rtFunc, 1, 0, 1.0, {});
+  testBesselJn<KIND>(rtFunc, 0, 0, 1.0, {0.765197694});
+  testBesselJn<KIND>(rtFunc, 0, 1, 1.0, {0.765197694, 0.440050572});
+  testBesselJn<KIND>(
+      rtFunc, 0, 2, 1.0, {0.765197694, 0.440050572, 0.114903487});
+  testBesselJn<KIND>(rtFunc, 1, 5, 5.0,
+      {-0.327579111, 0.046565145, 0.364831239, 0.391232371, 0.261140555});
+}
+
+template <int KIND> static void testBesselJnX0(BesselX0FuncType<KIND> rtFunc) {
+  testBesselJnX0<KIND>(rtFunc, 1, 0);
+  testBesselJnX0<KIND>(rtFunc, 0, 0);
+  testBesselJnX0<KIND>(rtFunc, 1, 1);
+  testBesselJnX0<KIND>(rtFunc, 0, 3);
+  testBesselJnX0<KIND>(rtFunc, 1, 4);
+}
+
+static void testBesselJn() {
+  testBesselJn<4>(RTNAME(BesselJn_4));
+  testBesselJn<8>(RTNAME(BesselJn_8));
+#if LDBL_MANT_DIG == 64
+  testBesselJn<10>(RTNAME(BesselJn_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselJn<16>(RTNAME(BesselJn_16));
+#endif
+
+  testBesselJnX0<4>(RTNAME(BesselJnX0_4));
+  testBesselJnX0<8>(RTNAME(BesselJnX0_8));
+#if LDBL_MANT_DIG == 64
+  testBesselJnX0<10>(RTNAME(BesselJnX0_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselJnX0<16>(RTNAME(BesselJnX0_16));
+#endif
+}
+
+TEST(Transformational, BesselJn) { testBesselJn(); }
+
+template <int KIND>
+static void testBesselYn(BesselFuncType<KIND> rtFunc, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, KIND> x,
+    const std::vector<CppTypeFor<TypeCategory::Real, KIND>> &expected) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+  unsigned len = expected.size();
+
+  CppTypeFor<TypeCategory::Real, KIND> anc0 = len > 0 ? expected[0] : 0.0;
+  CppTypeFor<TypeCategory::Real, KIND> anc1 = len > 1 ? expected[1] : 0.0;
+
+  rtFunc(result, n1, n2, x, anc0, anc1, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(
+      result.ElementBytes(), sizeof(CppTypeFor<TypeCategory::Real, KIND>));
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), len);
+
+  for (size_t j{0}; j < len; ++j) {
+    EXPECT_NEAR(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        expected[j], besselEpsilon<KIND>);
+  }
+}
+
+template <int KIND>
+static void testBesselYnX0(
+    BesselX0FuncType<KIND> rtFunc, int32_t n1, int32_t n2) {
+  StaticDescriptor<1> desc;
+  Descriptor &result{desc.descriptor()};
+
+  rtFunc(result, n1, n2, __FILE__, __LINE__);
+
+  EXPECT_EQ(result.type(), (TypeCode{TypeCategory::Real, KIND}));
+  EXPECT_EQ(result.rank(), 1);
+  EXPECT_EQ(result.GetDimension(0).LowerBound(), 1);
+  EXPECT_EQ(result.GetDimension(0).Extent(), n2 >= n1 ? n2 - n1 + 1 : 0);
+
+  if (n2 < n1) {
+    return;
+  }
+
+  for (int j{0}; j < (n2 - n1 + 1); ++j) {
+    EXPECT_EQ(
+        (*result.ZeroBasedIndexedElement<CppTypeFor<TypeCategory::Real, KIND>>(
+            j)),
+        (-std::numeric_limits<
+            CppTypeFor<TypeCategory::Real, KIND>>::infinity()));
+  }
+}
+
+template <int KIND> static void testBesselYn(BesselFuncType<KIND> rtFunc) {
+  testBesselYn<KIND>(rtFunc, 1, 0, 1.0, {});
+  testBesselYn<KIND>(rtFunc, 0, 0, 1.0, {0.08825695});
+  testBesselYn<KIND>(rtFunc, 0, 1, 1.0, {0.08825695, -0.7812128});
+  testBesselYn<KIND>(rtFunc, 0, 2, 1.0, {0.0882569555, -0.7812128, -1.6506826});
+  testBesselYn<KIND>(rtFunc, 1, 5, 1.0,
+      {-0.7812128, -1.6506826, -5.8215175, -33.278423, -260.40588});
+}
+
+template <int KIND> static void testBesselYnX0(BesselX0FuncType<KIND> rtFunc) {
+  testBesselYnX0<KIND>(rtFunc, 1, 0);
+  testBesselYnX0<KIND>(rtFunc, 0, 0);
+  testBesselYnX0<KIND>(rtFunc, 1, 1);
+  testBesselYnX0<KIND>(rtFunc, 0, 3);
+  testBesselYnX0<KIND>(rtFunc, 1, 4);
+}
+
+static void testBesselYn() {
+  testBesselYn<4>(RTNAME(BesselYn_4));
+  testBesselYn<8>(RTNAME(BesselYn_8));
+#if LDBL_MANT_DIG == 64
+  testBesselYn<10>(RTNAME(BesselYn_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselYn<16>(RTNAME(BesselYn_16));
+#endif
+
+  testBesselYnX0<4>(RTNAME(BesselYnX0_4));
+  testBesselYnX0<8>(RTNAME(BesselYnX0_8));
+#if LDBL_MANT_DIG == 64
+  testBesselYnX0<10>(RTNAME(BesselYnX0_10));
+#endif
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+  testBesselYnX0<16>(RTNAME(BesselYnX0_16));
+#endif
+}
+
+TEST(Transformational, BesselYn) { testBesselYn(); }
+
 TEST(Transformational, Shifts) {
   // ARRAY  1 3 5
   //        2 4 6