Lowering and runtime support for F08 transformational intrinsics: BESSEL_JN and BESSEL_YN

[tarunprabhu]

Summary

This patch implements lowering and adds runtime support for F08 transformational variants of 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).

The special case for x == 0 in BESSEL_JN(N1, N2, x) and BESSEL_YN(N1, N2, x) are handled in the same runtime function although the check for x == 0 is carried out in FIR. Similarly, the anchor points for the recursion are also calculated in FIR and passed explicitly to the runtime function.

Tests have been added for the new runtime functions and the lowering.

The existing tests for lowering were modified so the tests for the elemental and transformational variants are stylistically more consistent.

Possible alternatives:

• Pass the elemental Bessel function that is used to calculate the anchor points of the recursion to the runtime instead of calculating it in FIR.

• Split the runtime implementation for the x == 0 and x != 0 cases.

[vzakhari]

flang/lib/Lower/IntrinsicCall.cpp

2695

I think unconditional computation of the recursion bases might be suboptimal at least for the n1 == n2 case. Can you please rewrite the selects with IfOp?

flang/lib/Optimizer/Builder/Runtime/Transformational.cpp

31

For portability, I think it is better to use 32 here and int32_t in the runtime implementation.

flang/runtime/transformational.cpp

168

Should we also verify that n1 and n2 are non-negative?

193

Early return is more preferable.

[tarunprabhu]

flang/runtime/transformational.cpp

168

The Fortran standard specifies that n1 and n2 should be non-negative.

However, gfortran and ifort do not seem to check for this case. The output of this patch as currently implemented is comparable to the behavior of those two compilers when x > 0.0. The x == 0.0 case is a bit messier.

A similar concern arose with one of the other intrinsics (https://reviews.llvm.org/D129316#3638582. See also @klausler's comment in the discourse post linked there). As before, I am happy to adhere more closely to the standard, or be more compatible with other compilers. For the latter case, I could do as @jeanPerier suggested in that patch and document our implementation more clearly.

However, we do diverge from other compilers in the behavior of the intrinsics when compile-time constant negative values are passed. Both gfortran and ifort raise a compile-error in that case, while we do not. The runtime implementation is consistent with the current behavior when the parameters are compile-time constants.

[tarunprabhu]

Significant changes to the patch:

• Use fir::IfOp so anchor points for the recursion are only calculated if necessary.

• Created runtime functions specialized for the x == 0 case. This has cleaned up the signature of the runtime functions so only those arguments that are necessary are passed along.

• Use int32_t in the signature for the runtime functions instead of relying on the width of the system's int datatype.

• Added tests for the genBessel* lowering functions in addition to tests of the runtime functions themselves.

[vzakhari]

Thank you for the explanantion and the changes! We can assert that n1 and n2 are non-negative later, if we decide to be more verbose in runtime.

flang/unittests/Runtime/Transformational.cpp

343

Can you please add a unit test for BesselYnX0_16? We've had issues before that std::numeric_limits was not fully implemented for CppTypeFor<TypeCategory::Real, 16> (e.g. D134496).

[tarunprabhu]

• Added comments in the runtime functions indicating that the non-negativity constraint on n1 and n2 that is required by the standard is not enforced.

• Modified the runtime tests so all kinds (4, 8, 10 and 16) are now tested (with appropriate guards for the KIND=10 and KIND=16 cases).

• Fixed includes in transformational.h and fixed the signature of some of the functions which were incorrect.

[vzakhari]

Thank you, Tarun!

[tarunprabhu]

bef2bb34bfadd61cf399a15f096c84980121279f

[vzakhari]

flang/lib/Lower/IntrinsicCall.cpp

2682

Hello @tarunprabhu, is unordered predicate here intentional?

[tarunprabhu]

flang/lib/Lower/IntrinsicCall.cpp

2682

I think I intended this to replicate the behavior of gfortran in the case when x is NaN. I'll have to double-check that though.

[vzakhari]

flang/lib/Lower/IntrinsicCall.cpp

2682

Yes, please double-check it, because treating NaN argument as zero may not be the right thing to do. We'd better propagate the NaN argument to the result.

flang/include/flang/Optimizer/Builder/Runtime/Transformational.h

 
 namespace fir {
 class ExtendedValue;
 class FirOpBuilder;
 } // namespace fir
 
 namespace fir::runtime {
 
+void genBesselJn(fir::FirOpBuilder &builder, mlir::Location loc,
+                 mlir::Value resultBox, mlir::Value n1, mlir::Value n2,
+                 mlir::Value x, mlir::Value bn2, mlir::Value bn2_1);
+
+void genBesselJnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                   mlir::Type xTy, mlir::Value resultBox, mlir::Value n1,
+                   mlir::Value n2);
+
+void genBesselYn(fir::FirOpBuilder &builder, mlir::Location loc,
+                 mlir::Value resultBox, mlir::Value n1, mlir::Value n2,
+                 mlir::Value x, mlir::Value bn1, mlir::Value bn1_1);
+
+void genBesselYnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                   mlir::Type xTy, mlir::Value resultBox, mlir::Value n1,
+                   mlir::Value n2);
+
 void genCshift(fir::FirOpBuilder &builder, mlir::Location loc,
                mlir::Value resultBox, mlir::Value arrayBox,
                mlir::Value shiftBox, mlir::Value dimBox);
 
 void genCshiftVector(fir::FirOpBuilder &builder, mlir::Location loc,
                      mlir::Value resultBox, mlir::Value arrayBox,
                      mlir::Value shiftBox);
 

flang/include/flang/Runtime/transformational.h

 // These are naive allocating implementations; optimized forms that manipulate
 // pointer descriptors or that supply functional views of arrays remain to
 // be defined and may instead be part of lowering (see docs/ArrayComposition.md)
 // for details).
 
 #ifndef FORTRAN_RUNTIME_TRANSFORMATIONAL_H_
 #define FORTRAN_RUNTIME_TRANSFORMATIONAL_H_
 
+#include "flang/Runtime/cpp-type.h"
 #include "flang/Runtime/entry-names.h"
+#include "flang/Runtime/float128.h"
 #include <cinttypes>
 
 namespace Fortran::runtime {
 
 class Descriptor;
 
 extern "C" {
 
 void RTNAME(Reshape)(Descriptor &result, const Descriptor &source,
     const Descriptor &shape, const Descriptor *pad = nullptr,
     const Descriptor *order = nullptr, const char *sourceFile = nullptr,
     int line = 0);
 
+void RTNAME(BesselJn_2)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn2, float bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJn_3)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn2, float bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJn_4)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn2, float bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJn_8)(Descriptor &result, int32_t n1, int32_t n2, double x,
+    double bn2, double bn2_1, const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    long double x, long double bn2, long double bn2_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppFloat128Type x, CppFloat128Type bn2, CppFloat128Type bn2_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+void RTNAME(BesselJnX0_2)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJnX0_3)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselJnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+void RTNAME(BesselYn_2)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn1, float bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYn_3)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn1, float bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYn_4)(Descriptor &result, int32_t n1, int32_t n2, float x,
+    float bn1, float bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYn_8)(Descriptor &result, int32_t n1, int32_t n2, double x,
+    double bn1, double bn1_1, const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    long double x, long double bn1, long double bn1_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppFloat128Type x, CppFloat128Type bn1, CppFloat128Type bn1_1,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+void RTNAME(BesselYnX0_2)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYnX0_3)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+void RTNAME(BesselYnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile = nullptr, int line = 0);
+#endif
+
 void RTNAME(Cshift)(Descriptor &result, const Descriptor &source,
     const Descriptor &shift, int dim = 1, const char *sourceFile = nullptr,
     int line = 0);
 void RTNAME(CshiftVector)(Descriptor &result, const Descriptor &source,
     std::int64_t shift, const char *sourceFile = nullptr, int line = 0);
 
 void RTNAME(Eoshift)(Descriptor &result, const Descriptor &source,
     const Descriptor &shift, const Descriptor *boundary = nullptr, int dim = 1,

flang/lib/Lower/IntrinsicCall.cpp

struct IntrinsicLibrary {
   fir::ExtendedValue genAllocated(mlir::Type,
                                   llvm::ArrayRef<fir::ExtendedValue>);
   mlir::Value genAnint(mlir::Type, llvm::ArrayRef<mlir::Value>);
   fir::ExtendedValue genAny(mlir::Type, llvm::ArrayRef<fir::ExtendedValue>);
   fir::ExtendedValue
       genCommandArgumentCount(mlir::Type, llvm::ArrayRef<fir::ExtendedValue>);
   fir::ExtendedValue genAssociated(mlir::Type,
                                    llvm::ArrayRef<fir::ExtendedValue>);
-
+  fir::ExtendedValue genBesselJn(mlir::Type,
+                                 llvm::ArrayRef<fir::ExtendedValue>);
+  fir::ExtendedValue genBesselYn(mlir::Type,
+                                 llvm::ArrayRef<fir::ExtendedValue>);
   /// Lower a bitwise comparison intrinsic using the given comparator.
   template <mlir::arith::CmpIPredicate pred>
   mlir::Value genBitwiseCompare(mlir::Type resultType,
                                 llvm::ArrayRef<mlir::Value> args);
 
   mlir::Value genBtest(mlir::Type, llvm::ArrayRef<mlir::Value>);
   mlir::Value genCeiling(mlir::Type, llvm::ArrayRef<mlir::Value>);
   fir::ExtendedValue genChar(mlir::Type, llvm::ArrayRef<fir::ExtendedValue>);
static constexpr IntrinsicHandler handlers[]{
     {"any",
      &I::genAny,
      {{{"mask", asAddr}, {"dim", asValue}}},
      /*isElemental=*/false},
     {"associated",
      &I::genAssociated,
      {{{"pointer", asInquired}, {"target", asInquired}}},
      /*isElemental=*/false},
+    {"bessel_jn",
+     &I::genBesselJn,
+     {{{"n1", asValue}, {"n2", asValue}, {"x", asValue}}},
+     /*isElemental=*/false},
+    {"bessel_yn",
+     &I::genBesselYn,
+     {{{"n1", asValue}, {"n2", asValue}, {"x", asValue}}},
+     /*isElemental=*/false},
     {"bge", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::uge>},
     {"bgt", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::ugt>},
     {"ble", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::ule>},
     {"blt", &I::genBitwiseCompare<mlir::arith::CmpIPredicate::ult>},
     {"btest", &I::genBtest},
     {"c_associated_c_funptr",
      &I::genCAssociatedCFunPtr,
      {{{"c_ptr_1", asAddr}, {"c_ptr_2", asAddr, handleDynamicOptional}}},
if (fir::valueHasFirAttribute(fir::getBase(target),
     targetBox = builder.createBox(loc, target);
   }
   mlir::Value pointerBoxRef =
       fir::factory::getMutableIRBox(builder, loc, *pointer);
   auto pointerBox = builder.create<fir::LoadOp>(loc, pointerBoxRef);
   return Fortran::lower::genAssociated(builder, loc, pointerBox, targetBox);
 }
 
+// BESSEL_JN
+fir::ExtendedValue
+IntrinsicLibrary::genBesselJn(mlir::Type resultType,
+                              llvm::ArrayRef<fir::ExtendedValue> args) {
+  assert(args.size() == 2 || args.size() == 3);
+
+  mlir::Value x = fir::getBase(args.back());
+
+  if (args.size() == 2) {
+    mlir::Value n = fir::getBase(args[0]);
+
+    return genRuntimeCall("bessel_jn", resultType, {n, x});
+  } else {
+    mlir::Value n1 = fir::getBase(args[0]);
+    mlir::Value n2 = fir::getBase(args[1]);
+
+    mlir::Type intTy = n1.getType();
+    mlir::Type floatTy = x.getType();
+    mlir::Value zero = builder.createRealZeroConstant(loc, floatTy);
+    mlir::Value one = builder.createIntegerConstant(loc, intTy, 1);
+
+    mlir::Type resultArrayType = builder.getVarLenSeqTy(resultType, 1);
+    fir::MutableBoxValue resultMutableBox =
+        fir::factory::createTempMutableBox(builder, loc, resultArrayType);
+    mlir::Value resultBox =
+        fir::factory::getMutableIRBox(builder, loc, resultMutableBox);
+
+    mlir::Value cmpXEq0 = builder.create<mlir::arith::CmpFOp>(
+        loc, mlir::arith::CmpFPredicate::UEQ, x, zero);
+    mlir::Value cmpN1LtN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::slt, n1, n2);
+    mlir::Value cmpN1EqN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::eq, n1, n2);
+
+    auto genXEq0 = [&]() {
+      fir::runtime::genBesselJnX0(builder, loc, floatTy, resultBox, n1, n2);
+    };
+
+    auto genN1LtN2 = [&]() {
+      // The runtime generates the values in the range using a backward
+      // recursion from n2 to n1. (see https://dlmf.nist.gov/10.74.iv and
+      // https://dlmf.nist.gov/10.6.E1). When n1 < n2, this requires
+      // the values of BESSEL_JN(n2) and BESSEL_JN(n2 - 1) since they
+      // are the anchors of the recursion.
+      mlir::Value n2_1 = builder.create<mlir::arith::SubIOp>(loc, n2, one);
+      mlir::Value bn2 = genRuntimeCall("bessel_jn", resultType, {n2, x});
+      mlir::Value bn2_1 = genRuntimeCall("bessel_jn", resultType, {n2_1, x});
+      fir::runtime::genBesselJn(builder, loc, resultBox, n1, n2, x, bn2, bn2_1);
+    };
+
+    auto genN1EqN2 = [&]() {
+      // When n1 == n2, only BESSEL_JN(n2) is needed.
+      mlir::Value bn2 = genRuntimeCall("bessel_jn", resultType, {n2, x});
+      fir::runtime::genBesselJn(builder, loc, resultBox, n1, n2, x, bn2, zero);
+    };
+
+    auto genN1GtN2 = [&]() {
+      // The standard requires n1 <= n2. However, we still need to allocate
+      // a zero-length array and return it when n1 > n2, so we do need to call
+      // the runtime function.
+      fir::runtime::genBesselJn(builder, loc, resultBox, n1, n2, x, zero, zero);
+    };
+
+    auto genN1GeN2 = [&] {
+      builder.genIfThenElse(loc, cmpN1EqN2)
+          .genThen(genN1EqN2)
+          .genElse(genN1GtN2)
+          .end();
+    };
+
+    auto genXNeq0 = [&]() {
+      builder.genIfThenElse(loc, cmpN1LtN2)
+          .genThen(genN1LtN2)
+          .genElse(genN1GeN2)
+          .end();
+    };
+
+    builder.genIfThenElse(loc, cmpXEq0)
+        .genThen(genXEq0)
+        .genElse(genXNeq0)
+        .end();
+
+    fir::ExtendedValue res =
+        fir::factory::genMutableBoxRead(builder, loc, resultMutableBox);
+    return res.match(
+        [&](const fir::ArrayBoxValue &box) -> fir::ExtendedValue {
+          addCleanUpForTemp(loc, box.getAddr());
+          return box;
+        },
+        [&](const auto &) -> fir::ExtendedValue {
+          fir::emitFatalError(loc, "unexpected result for BESSEL_JN");
+        });
+  }
+}
+
+// BESSEL_YN
+fir::ExtendedValue
+IntrinsicLibrary::genBesselYn(mlir::Type resultType,
+                              llvm::ArrayRef<fir::ExtendedValue> args) {
+  assert(args.size() == 2 || args.size() == 3);
+
+  mlir::Value x = fir::getBase(args.back());
+
+  if (args.size() == 2) {
+    mlir::Value n = fir::getBase(args[0]);
+
+    return genRuntimeCall("bessel_yn", resultType, {n, x});
+  } else {
+    mlir::Value n1 = fir::getBase(args[0]);
+    mlir::Value n2 = fir::getBase(args[1]);
+
+    mlir::Type floatTy = x.getType();
+    mlir::Type intTy = n1.getType();
+    mlir::Value zero = builder.createRealZeroConstant(loc, floatTy);
+    mlir::Value one = builder.createIntegerConstant(loc, intTy, 1);
+
+    mlir::Type resultArrayType = builder.getVarLenSeqTy(resultType, 1);
+    fir::MutableBoxValue resultMutableBox =
+        fir::factory::createTempMutableBox(builder, loc, resultArrayType);
+    mlir::Value resultBox =
+        fir::factory::getMutableIRBox(builder, loc, resultMutableBox);
+
+    mlir::Value cmpXEq0 = builder.create<mlir::arith::CmpFOp>(
+        loc, mlir::arith::CmpFPredicate::UEQ, x, zero);
+    mlir::Value cmpN1LtN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::slt, n1, n2);
+    mlir::Value cmpN1EqN2 = builder.create<mlir::arith::CmpIOp>(
+        loc, mlir::arith::CmpIPredicate::eq, n1, n2);
+
+    auto genXEq0 = [&]() {
+      fir::runtime::genBesselYnX0(builder, loc, floatTy, resultBox, n1, n2);
+    };
+
+    auto genN1LtN2 = [&]() {
+      // The runtime generates the values in the range using a forward
+      // recursion from n1 to n2. (see https://dlmf.nist.gov/10.74.iv and
+      // https://dlmf.nist.gov/10.6.E1). When n1 < n2, this requires
+      // the values of BESSEL_YN(n1) and BESSEL_YN(n1 + 1) since they
+      // are the anchors of the recursion.
+      mlir::Value n1_1 = builder.create<mlir::arith::AddIOp>(loc, n1, one);
+      mlir::Value bn1 = genRuntimeCall("bessel_yn", resultType, {n1, x});
+      mlir::Value bn1_1 = genRuntimeCall("bessel_yn", resultType, {n1_1, x});
+      fir::runtime::genBesselYn(builder, loc, resultBox, n1, n2, x, bn1, bn1_1);
+    };
+
+    auto genN1EqN2 = [&]() {
+      // When n1 == n2, only BESSEL_YN(n1) is needed.
+      mlir::Value bn1 = genRuntimeCall("bessel_yn", resultType, {n1, x});
+      fir::runtime::genBesselYn(builder, loc, resultBox, n1, n2, x, bn1, zero);
+    };
+
+    auto genN1GtN2 = [&]() {
+      // The standard requires n1 <= n2. However, we still need to allocate
+      // a zero-length array and return it when n1 > n2, so we do need to call
+      // the runtime function.
+      fir::runtime::genBesselYn(builder, loc, resultBox, n1, n2, x, zero, zero);
+    };
+
+    auto genN1GeN2 = [&] {
+      builder.genIfThenElse(loc, cmpN1EqN2)
+          .genThen(genN1EqN2)
+          .genElse(genN1GtN2)
+          .end();
+    };
+
+    auto genXNeq0 = [&]() {
+      builder.genIfThenElse(loc, cmpN1LtN2)
+          .genThen(genN1LtN2)
+          .genElse(genN1GeN2)
+          .end();
+    };
+
+    builder.genIfThenElse(loc, cmpXEq0)
+        .genThen(genXEq0)
+        .genElse(genXNeq0)
+        .end();
+
+    fir::ExtendedValue res =
+        fir::factory::genMutableBoxRead(builder, loc, resultMutableBox);
+    return res.match(
+        [&](const fir::ArrayBoxValue &box) -> fir::ExtendedValue {
+          addCleanUpForTemp(loc, box.getAddr());
+          return box;
+        },
+        [&](const auto &) -> fir::ExtendedValue {
+          fir::emitFatalError(loc, "unexpected result for BESSEL_YN");
+        });
+  }
+}
+
 // BGE, BGT, BLE, BLT
 template <mlir::arith::CmpIPredicate pred>
 mlir::Value
 IntrinsicLibrary::genBitwiseCompare(mlir::Type resultType,
                                     llvm::ArrayRef<mlir::Value> args) {
   assert(args.size() == 2);
 
   mlir::Value arg0 = args[0];

flang/lib/Optimizer/Builder/Runtime/Transformational.cpp

 #include "flang/Optimizer/Builder/Runtime/RTBuilder.h"
 #include "flang/Optimizer/Builder/Todo.h"
 #include "flang/Runtime/matmul.h"
 #include "flang/Runtime/transformational.h"
 #include "mlir/Dialect/Func/IR/FuncOps.h"
 
 using namespace Fortran::runtime;
 
+/// Placeholder for real*10 version of BesselJn intrinsic.
+struct ForcedBesselJn_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJn_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF80(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselJn intrinsic.
+struct ForcedBesselJn_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJn_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF128(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*10 version of BesselJn intrinsic when `x == 0.0`.
+struct ForcedBesselJnX0_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJnX0_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselJn intrinsic when `x == 0.0`.
+struct ForcedBesselJnX0_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselJnX0_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*10 version of BesselYn intrinsic.
+struct ForcedBesselYn_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYn_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF80(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselYn intrinsic.
+struct ForcedBesselYn_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYn_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto ty = mlir::FloatType::getF128(ctx);
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(
+          ctx, {boxTy, intTy, intTy, ty, ty, ty, strTy, intTy}, {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*10 version of BesselYn intrinsic when `x == 0.0`.
+struct ForcedBesselYnX0_10 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYnX0_10));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Placeholder for real*16 version of BesselYn intrinsic when `x == 0.0`.
+struct ForcedBesselYnX0_16 {
+  static constexpr const char *name = ExpandAndQuoteKey(RTNAME(BesselYnX0_16));
+  static constexpr fir::runtime::FuncTypeBuilderFunc getTypeModel() {
+    return [](mlir::MLIRContext *ctx) {
+      auto boxTy =
+          fir::runtime::getModel<Fortran::runtime::Descriptor &>()(ctx);
+      auto strTy = fir::ReferenceType::get(mlir::IntegerType::get(ctx, 8));
+      auto intTy = mlir::IntegerType::get(ctx, 32);
+      auto noneTy = mlir::NoneType::get(ctx);
+      return mlir::FunctionType::get(ctx, {boxTy, intTy, intTy, strTy, intTy},
+                                     {noneTy});
+    };
+  }
+};
+
+/// Generate call to `BesselJn` intrinsic.
+void fir::runtime::genBesselJn(fir::FirOpBuilder &builder, mlir::Location loc,
+                               mlir::Value resultBox, mlir::Value n1,
+                               mlir::Value n2, mlir::Value x, mlir::Value bn2,
+                               mlir::Value bn2_1) {
+  mlir::func::FuncOp func;
+  auto xTy = x.getType();
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_JN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJn_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJn_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJn_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJn_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_JN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(7));
+  auto args =
+      fir::runtime::createArguments(builder, loc, fTy, resultBox, n1, n2, x,
+                                    bn2, bn2_1, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
+/// Generate call to `BesselJn` intrinsic. This is used when `x == 0.0`.
+void fir::runtime::genBesselJnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                                 mlir::Type xTy, mlir::Value resultBox,
+                                 mlir::Value n1, mlir::Value n2) {
+  mlir::func::FuncOp func;
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_JN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJnX0_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselJnX0_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJnX0_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselJnX0_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_JN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(4));
+  auto args = fir::runtime::createArguments(builder, loc, fTy, resultBox, n1,
+                                            n2, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
+/// Generate call to `BesselYn` intrinsic.
+void fir::runtime::genBesselYn(fir::FirOpBuilder &builder, mlir::Location loc,
+                               mlir::Value resultBox, mlir::Value n1,
+                               mlir::Value n2, mlir::Value x, mlir::Value bn1,
+                               mlir::Value bn1_1) {
+  mlir::func::FuncOp func;
+  auto xTy = x.getType();
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_YN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYn_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYn_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYn_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYn_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_YN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(7));
+  auto args =
+      fir::runtime::createArguments(builder, loc, fTy, resultBox, n1, n2, x,
+                                    bn1, bn1_1, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
+/// Generate call to `BesselYn` intrinsic. This is used when `x == 0.0`.
+void fir::runtime::genBesselYnX0(fir::FirOpBuilder &builder, mlir::Location loc,
+                                 mlir::Type xTy, mlir::Value resultBox,
+                                 mlir::Value n1, mlir::Value n2) {
+  mlir::func::FuncOp func;
+
+  if (xTy.isF16() || xTy.isBF16())
+    TODO(loc, "half-precision BESSEL_YN");
+  else if (xTy.isF32())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYnX0_4)>(loc, builder);
+  else if (xTy.isF64())
+    func = fir::runtime::getRuntimeFunc<mkRTKey(BesselYnX0_8)>(loc, builder);
+  else if (xTy.isF80())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYnX0_10>(loc, builder);
+  else if (xTy.isF128())
+    func = fir::runtime::getRuntimeFunc<ForcedBesselYnX0_16>(loc, builder);
+  else
+    fir::emitFatalError(loc, "invalid type in BESSEL_YN");
+
+  auto fTy = func.getFunctionType();
+  auto sourceFile = fir::factory::locationToFilename(builder, loc);
+  auto sourceLine =
+      fir::factory::locationToLineNo(builder, loc, fTy.getInput(4));
+  auto args = fir::runtime::createArguments(builder, loc, fTy, resultBox, n1,
+                                            n2, sourceFile, sourceLine);
+  builder.create<fir::CallOp>(loc, func, args);
+}
+
 /// Generate call to Cshift intrinsic
 void fir::runtime::genCshift(fir::FirOpBuilder &builder, mlir::Location loc,
                              mlir::Value resultBox, mlir::Value arrayBox,
                              mlir::Value shiftBox, mlir::Value dimBox) {
   auto cshiftFunc = fir::runtime::getRuntimeFunc<mkRTKey(Cshift)>(loc, builder);
   auto fTy = cshiftFunc.getFunctionType();
   auto sourceFile = fir::factory::locationToFilename(builder, loc);
   auto sourceLine =

flang/runtime/transformational.cpp

static inline std::size_t AllocateResult(Descriptor &result,
   }
   if (int stat{result.Allocate()}) {
     terminator.Crash(
         "%s: Could not allocate memory for result (stat=%d)", function, stat);
   }
   return elementLen;
 }
 
+template <TypeCategory CAT, int KIND>
+static inline std::size_t AllocateBesselResult(Descriptor &result, int32_t n1,
+    int32_t n2, Terminator &terminator, const char *function) {
+  int rank{1};
+  SubscriptValue extent[maxRank];
+  for (int j{0}; j < maxRank; j++) {
+    extent[j] = 0;
+  }
+  if (n1 <= n2) {
+    extent[0] = n2 - n1 + 1;
+  }
+
+  std::size_t elementLen{Descriptor::BytesFor(CAT, KIND)};
+  result.Establish(TypeCode{CAT, KIND}, elementLen, nullptr, rank, extent,
+      CFI_attribute_allocatable, false);
+  for (int j{0}; j < rank; ++j) {
+    result.GetDimension(j).SetBounds(1, extent[j]);
+  }
+  if (int stat{result.Allocate()}) {
+    terminator.Crash(
+        "%s: Could not allocate memory for result (stat=%d)", function, stat);
+  }
+  return elementLen;
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselJn(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<CAT, KIND> x, CppTypeFor<CAT, KIND> bn2,
+    CppTypeFor<CAT, KIND> bn2_1, const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_JN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // if n2 >= n1, there will be at least one element in the result.
+  at[0] = n2 - n1 + 1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn2;
+
+  if (n2 == n1) {
+    return;
+  }
+
+  at[0] = n2 - n1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn2_1;
+
+  // Bessel functions of the first kind are stable for a backward recursion
+  // (see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).
+  //
+  //     J(n-1, x) = (2.0 / x) * n * J(n, x) - J(n+1, x)
+  //
+  // which is equivalent to
+  //
+  //     J(n, x) = (2.0 / x) * (n + 1) * J(n+1, x) - J(n+2, x)
+  //
+  CppTypeFor<CAT, KIND> bn_2 = bn2;
+  CppTypeFor<CAT, KIND> bn_1 = bn2_1;
+  CppTypeFor<CAT, KIND> twoOverX = 2.0 / x;
+  for (int n{n2 - 2}; n >= n1; --n) {
+    auto bn = twoOverX * (n + 1) * bn_1 - bn_2;
+
+    at[0] = n - n1 + 1;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) = bn;
+
+    bn_2 = bn_1;
+    bn_1 = bn;
+  }
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselJnX0(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_JN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // J(0, 0.0) = 1.0, when n == 0.
+  // J(n, 0.0) = 0.0, when n > 0.
+  at[0] = 1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = (n1 == 0) ? 1.0 : 0.0;
+  for (int j{2}; j <= n2 - n1 + 1; ++j) {
+    at[0] = j;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) = 0.0;
+  }
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselYn(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<CAT, KIND> x, CppTypeFor<CAT, KIND> bn1,
+    CppTypeFor<CAT, KIND> bn1_1, const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_YN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // if n2 >= n1, there will be at least one element in the result.
+  at[0] = 1;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn1;
+
+  if (n2 == n1) {
+    return;
+  }
+
+  at[0] = 2;
+  *result.Element<CppTypeFor<CAT, KIND>>(at) = bn1_1;
+
+  // Bessel functions of the second kind are stable for a forward recursion
+  // (see https://dlmf.nist.gov/10.74.iv and https://dlmf.nist.gov/10.6.E1).
+  //
+  //     Y(n+1, x) = (2.0 / x) * n * Y(n, x) - Y(n-1, x)
+  //
+  // which is equivalent to
+  //
+  //     Y(n, x) = (2.0 / x) * (n - 1) * Y(n-1, x) - Y(n-2, x)
+  //
+  CppTypeFor<CAT, KIND> bn_2 = bn1;
+  CppTypeFor<CAT, KIND> bn_1 = bn1_1;
+  CppTypeFor<CAT, KIND> twoOverX = 2.0 / x;
+  for (int n{n1 + 2}; n <= n2; ++n) {
+    auto bn = twoOverX * (n - 1) * bn_1 - bn_2;
+
+    at[0] = n - n1 + 1;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) = bn;
+
+    bn_2 = bn_1;
+    bn_1 = bn;
+  }
+}
+
+template <TypeCategory CAT, int KIND>
+static inline void DoBesselYnX0(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  Terminator terminator{sourceFile, line};
+  AllocateBesselResult<CAT, KIND>(result, n1, n2, terminator, "BESSEL_YN");
+
+  // The standard requires that n1 and n2 be non-negative. However, some other
+  // compilers generate results even when n1 and/or n2 are negative. For now,
+  // we also do not enforce the non-negativity constraint.
+  if (n2 < n1) {
+    return;
+  }
+
+  SubscriptValue at[maxRank];
+  for (int j{0}; j < maxRank; ++j) {
+    at[j] = 0;
+  }
+
+  // Y(n, 0.0) = -Inf, when n >= 0
+  for (int j{1}; j <= n2 - n1 + 1; ++j) {
+    at[0] = j;
+    *result.Element<CppTypeFor<CAT, KIND>>(at) =
+        -std::numeric_limits<CppTypeFor<CAT, KIND>>::infinity();
+  }
+}
+
 extern "C" {
 
+// BESSEL_JN
+// TODO: REAL(2 & 3)
+void RTNAME(BesselJn_4)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 4> x, CppTypeFor<TypeCategory::Real, 4> bn2,
+    CppTypeFor<TypeCategory::Real, 4> bn2_1, const char *sourceFile, int line) {
+  DoBesselJn<TypeCategory::Real, 4>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+
+void RTNAME(BesselJn_8)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 8> x, CppTypeFor<TypeCategory::Real, 8> bn2,
+    CppTypeFor<TypeCategory::Real, 8> bn2_1, const char *sourceFile, int line) {
+  DoBesselJn<TypeCategory::Real, 8>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 10> x,
+    CppTypeFor<TypeCategory::Real, 10> bn2,
+    CppTypeFor<TypeCategory::Real, 10> bn2_1, const char *sourceFile,
+    int line) {
+  DoBesselJn<TypeCategory::Real, 10>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 16> x,
+    CppTypeFor<TypeCategory::Real, 16> bn2,
+    CppTypeFor<TypeCategory::Real, 16> bn2_1, const char *sourceFile,
+    int line) {
+  DoBesselJn<TypeCategory::Real, 16>(
+      result, n1, n2, x, bn2, bn2_1, sourceFile, line);
+}
+#endif
+
+// TODO: REAL(2 & 3)
+void RTNAME(BesselJnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 4>(result, n1, n2, sourceFile, line);
+}
+
+void RTNAME(BesselJnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 8>(result, n1, n2, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselJnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 10>(result, n1, n2, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselJnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselJnX0<TypeCategory::Real, 16>(result, n1, n2, sourceFile, line);
+}
+#endif
+
+// BESSEL_YN
+// TODO: REAL(2 & 3)
+void RTNAME(BesselYn_4)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 4> x, CppTypeFor<TypeCategory::Real, 4> bn1,
+    CppTypeFor<TypeCategory::Real, 4> bn1_1, const char *sourceFile, int line) {
+  DoBesselYn<TypeCategory::Real, 4>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+
+void RTNAME(BesselYn_8)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 8> x, CppTypeFor<TypeCategory::Real, 8> bn1,
+    CppTypeFor<TypeCategory::Real, 8> bn1_1, const char *sourceFile, int line) {
+  DoBesselYn<TypeCategory::Real, 8>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYn_10)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 10> x,
+    CppTypeFor<TypeCategory::Real, 10> bn1,
+    CppTypeFor<TypeCategory::Real, 10> bn1_1, const char *sourceFile,
+    int line) {
+  DoBesselYn<TypeCategory::Real, 10>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYn_16)(Descriptor &result, int32_t n1, int32_t n2,
+    CppTypeFor<TypeCategory::Real, 16> x,
+    CppTypeFor<TypeCategory::Real, 16> bn1,
+    CppTypeFor<TypeCategory::Real, 16> bn1_1, const char *sourceFile,
+    int line) {
+  DoBesselYn<TypeCategory::Real, 16>(
+      result, n1, n2, x, bn1, bn1_1, sourceFile, line);
+}
+#endif
+
+// TODO: REAL(2 & 3)
+void RTNAME(BesselYnX0_4)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 4>(result, n1, n2, sourceFile, line);
+}
+
+void RTNAME(BesselYnX0_8)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 8>(result, n1, n2, sourceFile, line);
+}
+
+#if LDBL_MANT_DIG == 64
+void RTNAME(BesselYnX0_10)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 10>(result, n1, n2, sourceFile, line);
+}
+#endif
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+void RTNAME(BesselYnX0_16)(Descriptor &result, int32_t n1, int32_t n2,
+    const char *sourceFile, int line) {
+  DoBesselYnX0<TypeCategory::Real, 16>(result, n1, n2, sourceFile, line);
+}
+#endif
+
 // CSHIFT where rank of ARRAY argument > 1
 void RTNAME(Cshift)(Descriptor &result, const Descriptor &source,
     const Descriptor &shift, int dim, const char *sourceFile, int line) {
   Terminator terminator{sourceFile, line};
   int rank{source.rank()};
   RUNTIME_CHECK(terminator, rank > 1);
   if (dim < 1 || dim > rank) {
     terminator.Crash(

flang/test/Lower/Intrinsics/bessel_jn.f90

 ! RUN: bbc -emit-fir %s -o - --math-runtime=fast | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=fast %s -o - | FileCheck --check-prefixes=ALL %s
 ! RUN: bbc -emit-fir %s -o - --math-runtime=relaxed | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=relaxed %s -o - | FileCheck --check-prefixes=ALL %s
 ! RUN: bbc -emit-fir %s -o - --math-runtime=precise | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=precise %s -o - | FileCheck --check-prefixes=ALL %s
 
+! ALL-LABEL: @_QPtest_real4
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f32
 function test_real4(x, n)
   real :: x, test_real4
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @jnf(%[[n]], %[[x]]) {{.*}} : (i32, f32) -> f32
   test_real4 = bessel_jn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real4
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @jnf({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f32) -> f32
-
+! ALL-LABEL: @_QPtest_real8
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f64
 function test_real8(x, n)
   real(8) :: x, test_real8
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @jn(%[[n]], %[[x]]) {{.*}} : (i32, f64) -> f64
   test_real8 = bessel_jn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real8
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @jn({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f64) -> f64
+! ALL-LABEL: @_QPtest_transformational_real4
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real4(x, n1, n2, r)
+  real(4) :: x
+  integer :: n1, n2
+  real(4) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f32
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf32>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f32
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJnX0_4(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n2_1:.*]] = arith.subi %[[n2]], %[[one]] : i32
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jnf(%[[n2]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[bn2_1:.*]] = fir.call @jnf(%[[n2_1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_4(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[bn2_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jnf(%[[n2]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_4(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_4(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_jn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf32>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf32>>>) -> !fir.heap<!fir.array<?xf32>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf32>>
+end subroutine test_transformational_real4
+
+! ALL-LABEL: @_QPtest_transformational_real8
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real8(x, n1, n2, r)
+  real(8) :: x
+  integer :: n1, n2
+  real(8) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f64
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf64>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f64
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJnX0_8(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n2_1:.*]] = arith.subi %[[n2]], %[[one]] : i32
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jn(%[[n2]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[bn2_1:.*]] = fir.call @jn(%[[n2_1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_8(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[bn2_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn2:.*]] = fir.call @jn(%[[n2]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_8(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn2]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselJn_8(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_jn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf64>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf64>>>) -> !fir.heap<!fir.array<?xf64>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf64>>
+end subroutine test_transformational_real8

flang/test/Lower/Intrinsics/bessel_yn.f90

 ! RUN: bbc -emit-fir %s -o - --math-runtime=fast | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=fast %s -o - | FileCheck --check-prefixes=ALL %s
 ! RUN: bbc -emit-fir %s -o - --math-runtime=relaxed | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=relaxed %s -o - | FileCheck --check-prefixes=ALL %s
 ! RUN: bbc -emit-fir %s -o - --math-runtime=precise | FileCheck --check-prefixes=ALL %s
 ! RUN: %flang_fc1 -emit-fir -mllvm -math-runtime=precise %s -o - | FileCheck --check-prefixes=ALL %s
 
+! ALL-LABEL: @_QPtest_real4
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f32
 function test_real4(x, n)
   real :: x, test_real4
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @ynf(%[[n]], %[[x]]) {{.*}} : (i32, f32) -> f32
   test_real4 = bessel_yn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real4
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @ynf({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f32) -> f32
-
+! ALL-LABEL: @_QPtest_real8
+! ALL-SAME: (%[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn:.*]]: !fir.ref<i32>{{.*}}) -> f64
 function test_real8(x, n)
   real(8) :: x, test_real8
   integer :: n
+
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref<i32>
+  ! ALL: fir.call @yn(%[[n]], %[[x]]) {{.*}} : (i32, f64) -> f64
   test_real8 = bessel_yn(n, x)
 end function
 
-! ALL-LABEL: @_QPtest_real8
-! ALL: {{%[A-Za-z0-9._]+}} = fir.call @yn({{%[A-Za-z0-9._]+}}, {{%[A-Za-z0-9._]+}}) {{.*}}: (i32, f64) -> f64
+! ALL-LABEL: @_QPtest_transformational_real4
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f32>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real4(x, n1, n2, r)
+  real(4) :: x
+  integer :: n1, n2
+  real(4) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f32
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf32>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f32>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f32
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYnX0_4(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n1_1:.*]] = arith.addi %[[n1]], %[[one]] : i32
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @ynf(%[[n1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[bn1_1:.*]] = fir.call @ynf(%[[n1_1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_4(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[bn1_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @ynf(%[[n1]], %[[x]]) {{.*}} : (i32, f32) -> f32
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_4(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_4(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f32, f32, f32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_yn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf32>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf32>>>) -> !fir.heap<!fir.array<?xf32>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf32>>
+end subroutine test_transformational_real4
+
+! ALL-LABEL: @_QPtest_transformational_real8
+! ALL-SAME: %[[argx:.*]]: !fir.ref<f64>{{.*}}, %[[argn1:.*]]: !fir.ref<i32>{{.*}}, %[[argn2:.*]]: !fir.ref<i32>{{.*}}
+subroutine test_transformational_real8(x, n1, n2, r)
+  real(8) :: x
+  integer :: n1, n2
+  real(8) :: r(:)
+
+  ! ALL-DAG: %[[zero:.*]] = arith.constant 0{{.*}} : f64
+  ! ALL-DAG: %[[one:.*]] = arith.constant 1 : i32
+  ! ALL-DAG: %[[r:.*]] = fir.alloca !fir.box<!fir.heap<!fir.array<?xf64>>>
+  ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref<f64>
+  ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref<i32>
+  ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref<i32>
+  ! ALL-DAG: %[[xeq0:.*]] = arith.cmpf ueq, %[[x]], %[[zero]] : f64
+  ! ALL-DAG: %[[n1ltn2:.*]] = arith.cmpi slt, %[[n1]], %[[n2]] : i32
+  ! ALL-DAG: %[[n1eqn2:.*]] = arith.cmpi eq, %[[n1]], %[[n2]] : i32
+  ! ALL: fir.if %[[xeq0]] {
+  ! ALL: %[[resxeq0:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYnX0_8(%[[resxeq0]], %[[n1]], %[[n2]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1ltn2]] {
+  ! ALL-DAG: %[[n1_1:.*]] = arith.addi %[[n1]], %[[one]] : i32
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @yn(%[[n1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[bn1_1:.*]] = fir.call @yn(%[[n1_1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1ltn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_8(%[[resn1ltn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[bn1_1]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-NEXT: fir.if %[[n1eqn2]] {
+  ! ALL-DAG: %[[bn1:.*]] = fir.call @yn(%[[n1]], %[[x]]) {{.*}} : (i32, f64) -> f64
+  ! ALL-DAG: %[[resn1eqn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_8(%[[resn1eqn2]], %[[n1]], %[[n2]], %[[x]], %[[bn1]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: } else {
+  ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}}
+  ! ALL: fir.call @_FortranABesselYn_8(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref<!fir.box<none>>, i32, i32, f64, f64, f64, !fir.ref<i8>, i32) -> none
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  ! ALL-NEXT: }
+  r = bessel_yn(n1, n2, x)
+  ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref<!fir.box<!fir.heap<!fir.array<?xf64>>>>
+  ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box<!fir.heap<!fir.array<?xf64>>>) -> !fir.heap<!fir.array<?xf64>>
+  ! ALL: fir.freemem %[[addr]] : !fir.heap<!fir.array<?xf64>>
+end subroutine test_transformational_real8

flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp

 //===- TransformationalTest.cpp -- Transformational intrinsic generation --===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "flang/Optimizer/Builder/Runtime/Transformational.h"
 #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 =
       fir::SequenceType::get(fir::SequenceType::Shape(1, 10), i32Ty);
   mlir::Value result = firBuilder->create<fir::UndefOp>(loc, seqTy);
   mlir::Value array = firBuilder->create<fir::UndefOp>(loc, seqTy);
   mlir::Value shift = firBuilder->create<fir::UndefOp>(loc, seqTy);
   mlir::Value dim = firBuilder->create<fir::UndefOp>(loc, seqTy);

flang/unittests/Runtime/Transformational.cpp

 //===-- flang/unittests/Runtime/Transformational.cpp ----------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "flang/Runtime/transformational.h"
 #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
   auto array{MakeArray<TypeCategory::Integer, 4>(
       std::vector<int>{2, 3}, std::vector<std::int32_t>{1, 2, 3, 4, 5, 6})};
   array->GetDimension(0).SetLowerBound(0); // shouldn't matter
   array->GetDimension(1).SetLowerBound(-1);
   StaticDescriptor<2, true> statDesc;
TEST(Transformational, Shifts) {
   }
   result.Destroy();
 
   // VECTOR EOSHIFT
   StaticDescriptor<0> boundaryDescriptor;
   Descriptor vectorBoundary{boundaryDescriptor.descriptor()};
   std::int32_t boundaryValue{343};
   vectorBoundary.Establish(TypeCategory::Integer, 4,
       const_cast<void *>(reinterpret_cast<const void *>(&boundaryValue)), 0);
   RTNAME(EoshiftVector)
   (vectorResult, *vector, 2, &vectorBoundary, __FILE__, __LINE__);
   EXPECT_EQ(vectorResult.type(), array->type());
   EXPECT_EQ(vectorResult.rank(), 1);
   EXPECT_EQ(vectorResult.GetDimension(0).LowerBound(), 1);
   EXPECT_EQ(vectorResult.GetDimension(0).Extent(), 6);
   EXPECT_EQ(vectorResult.type(), (TypeCode{TypeCategory::Integer, 4}));
   static std::int32_t eoshiftVectorExpect[6]{3, 4, 5, 6, 343, 343};