Index: flang/include/flang/Optimizer/Builder/Runtime/Transformational.h =================================================================== --- flang/include/flang/Optimizer/Builder/Runtime/Transformational.h +++ flang/include/flang/Optimizer/Builder/Runtime/Transformational.h @@ -19,6 +19,22 @@ 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); Index: flang/include/flang/Runtime/transformational.h =================================================================== --- flang/include/flang/Runtime/transformational.h +++ flang/include/flang/Runtime/transformational.h @@ -17,7 +17,9 @@ #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 namespace Fortran::runtime { @@ -31,6 +33,98 @@ 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); Index: flang/lib/Lower/IntrinsicCall.cpp =================================================================== --- flang/lib/Lower/IntrinsicCall.cpp +++ flang/lib/Lower/IntrinsicCall.cpp @@ -465,7 +465,10 @@ genCommandArgumentCount(mlir::Type, llvm::ArrayRef); fir::ExtendedValue genAssociated(mlir::Type, llvm::ArrayRef); - + fir::ExtendedValue genBesselJn(mlir::Type, + llvm::ArrayRef); + fir::ExtendedValue genBesselYn(mlir::Type, + llvm::ArrayRef); /// Lower a bitwise comparison intrinsic using the given comparator. template mlir::Value genBitwiseCompare(mlir::Type resultType, @@ -735,6 +738,14 @@ &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}, {"bgt", &I::genBitwiseCompare}, {"ble", &I::genBitwiseCompare}, @@ -2640,6 +2651,196 @@ return Fortran::lower::genAssociated(builder, loc, pointerBox, targetBox); } +// BESSEL_JN +fir::ExtendedValue +IntrinsicLibrary::genBesselJn(mlir::Type resultType, + llvm::ArrayRef 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( + loc, mlir::arith::CmpFPredicate::UEQ, x, zero); + mlir::Value cmpN1LtN2 = builder.create( + loc, mlir::arith::CmpIPredicate::slt, n1, n2); + mlir::Value cmpN1EqN2 = builder.create( + 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(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 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( + loc, mlir::arith::CmpFPredicate::UEQ, x, zero); + mlir::Value cmpN1LtN2 = builder.create( + loc, mlir::arith::CmpIPredicate::slt, n1, n2); + mlir::Value cmpN1EqN2 = builder.create( + 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(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::Value Index: flang/lib/Optimizer/Builder/Runtime/Transformational.cpp =================================================================== --- flang/lib/Optimizer/Builder/Runtime/Transformational.cpp +++ flang/lib/Optimizer/Builder/Runtime/Transformational.cpp @@ -19,6 +19,256 @@ 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()(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()(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()(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()(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()(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()(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()(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()(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(loc, builder); + else if (xTy.isF64()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF80()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF128()) + func = fir::runtime::getRuntimeFunc(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(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(loc, builder); + else if (xTy.isF64()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF80()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF128()) + func = fir::runtime::getRuntimeFunc(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(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(loc, builder); + else if (xTy.isF64()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF80()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF128()) + func = fir::runtime::getRuntimeFunc(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(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(loc, builder); + else if (xTy.isF64()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF80()) + func = fir::runtime::getRuntimeFunc(loc, builder); + else if (xTy.isF128()) + func = fir::runtime::getRuntimeFunc(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(loc, func, args); +} + /// Generate call to Cshift intrinsic void fir::runtime::genCshift(fir::FirOpBuilder &builder, mlir::Location loc, mlir::Value resultBox, mlir::Value arrayBox, Index: flang/runtime/transformational.cpp =================================================================== --- flang/runtime/transformational.cpp +++ flang/runtime/transformational.cpp @@ -133,8 +133,319 @@ return elementLen; } +template +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 +static inline void DoBesselJn(Descriptor &result, int32_t n1, int32_t n2, + CppTypeFor x, CppTypeFor bn2, + CppTypeFor bn2_1, const char *sourceFile, int line) { + Terminator terminator{sourceFile, line}; + AllocateBesselResult(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>(at) = bn2; + + if (n2 == n1) { + return; + } + + at[0] = n2 - n1; + *result.Element>(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 bn_2 = bn2; + CppTypeFor bn_1 = bn2_1; + CppTypeFor 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>(at) = bn; + + bn_2 = bn_1; + bn_1 = bn; + } +} + +template +static inline void DoBesselJnX0(Descriptor &result, int32_t n1, int32_t n2, + const char *sourceFile, int line) { + Terminator terminator{sourceFile, line}; + AllocateBesselResult(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>(at) = (n1 == 0) ? 1.0 : 0.0; + for (int j{2}; j <= n2 - n1 + 1; ++j) { + at[0] = j; + *result.Element>(at) = 0.0; + } +} + +template +static inline void DoBesselYn(Descriptor &result, int32_t n1, int32_t n2, + CppTypeFor x, CppTypeFor bn1, + CppTypeFor bn1_1, const char *sourceFile, int line) { + Terminator terminator{sourceFile, line}; + AllocateBesselResult(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>(at) = bn1; + + if (n2 == n1) { + return; + } + + at[0] = 2; + *result.Element>(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 bn_2 = bn1; + CppTypeFor bn_1 = bn1_1; + CppTypeFor 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>(at) = bn; + + bn_2 = bn_1; + bn_1 = bn; + } +} + +template +static inline void DoBesselYnX0(Descriptor &result, int32_t n1, int32_t n2, + const char *sourceFile, int line) { + Terminator terminator{sourceFile, line}; + AllocateBesselResult(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>(at) = + -std::numeric_limits>::infinity(); + } +} + extern "C" { +// BESSEL_JN +// TODO: REAL(2 & 3) +void RTNAME(BesselJn_4)(Descriptor &result, int32_t n1, int32_t n2, + CppTypeFor x, CppTypeFor bn2, + CppTypeFor bn2_1, const char *sourceFile, int line) { + DoBesselJn( + result, n1, n2, x, bn2, bn2_1, sourceFile, line); +} + +void RTNAME(BesselJn_8)(Descriptor &result, int32_t n1, int32_t n2, + CppTypeFor x, CppTypeFor bn2, + CppTypeFor bn2_1, const char *sourceFile, int line) { + DoBesselJn( + 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 x, + CppTypeFor bn2, + CppTypeFor bn2_1, const char *sourceFile, + int line) { + DoBesselJn( + 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 x, + CppTypeFor bn2, + CppTypeFor bn2_1, const char *sourceFile, + int line) { + DoBesselJn( + 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(result, n1, n2, sourceFile, line); +} + +void RTNAME(BesselJnX0_8)(Descriptor &result, int32_t n1, int32_t n2, + const char *sourceFile, int line) { + DoBesselJnX0(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(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(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 x, CppTypeFor bn1, + CppTypeFor bn1_1, const char *sourceFile, int line) { + DoBesselYn( + result, n1, n2, x, bn1, bn1_1, sourceFile, line); +} + +void RTNAME(BesselYn_8)(Descriptor &result, int32_t n1, int32_t n2, + CppTypeFor x, CppTypeFor bn1, + CppTypeFor bn1_1, const char *sourceFile, int line) { + DoBesselYn( + 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 x, + CppTypeFor bn1, + CppTypeFor bn1_1, const char *sourceFile, + int line) { + DoBesselYn( + 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 x, + CppTypeFor bn1, + CppTypeFor bn1_1, const char *sourceFile, + int line) { + DoBesselYn( + 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(result, n1, n2, sourceFile, line); +} + +void RTNAME(BesselYnX0_8)(Descriptor &result, int32_t n1, int32_t n2, + const char *sourceFile, int line) { + DoBesselYnX0(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(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(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) { Index: flang/test/Lower/Intrinsics/bessel_jn.f90 =================================================================== --- flang/test/Lower/Intrinsics/bessel_jn.f90 +++ flang/test/Lower/Intrinsics/bessel_jn.f90 @@ -5,20 +5,112 @@ ! 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{{.*}}, %[[argn:.*]]: !fir.ref{{.*}}) -> f32 function test_real4(x, n) real :: x, test_real4 integer :: n + + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref + ! 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{{.*}}, %[[argn:.*]]: !fir.ref{{.*}}) -> f64 function test_real8(x, n) real(8) :: x, test_real8 integer :: n + + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref + ! 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{{.*}}, %[[argn1:.*]]: !fir.ref{{.*}}, %[[argn2:.*]]: !fir.ref{{.*}} +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>> + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref + ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref + ! 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>, i32, i32, !fir.ref, 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>, i32, i32, f32, f32, f32, !fir.ref, 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>, i32, i32, f32, f32, f32, !fir.ref, i32) -> none + ! ALL-NEXT: } else { + ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}} + ! ALL: fir.call @_FortranABesselJn_4(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref>, i32, i32, f32, f32, f32, !fir.ref, i32) -> none + ! ALL-NEXT: } + ! ALL-NEXT: } + ! ALL-NEXT: } + r = bessel_jn(n1, n2, x) + ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref>>> + ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box>>) -> !fir.heap> + ! ALL: fir.freemem %[[addr]] : !fir.heap> +end subroutine test_transformational_real4 + +! ALL-LABEL: @_QPtest_transformational_real8 +! ALL-SAME: %[[argx:.*]]: !fir.ref{{.*}}, %[[argn1:.*]]: !fir.ref{{.*}}, %[[argn2:.*]]: !fir.ref{{.*}} +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>> + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref + ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref + ! 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>, i32, i32, !fir.ref, 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>, i32, i32, f64, f64, f64, !fir.ref, 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>, i32, i32, f64, f64, f64, !fir.ref, i32) -> none + ! ALL-NEXT: } else { + ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}} + ! ALL: fir.call @_FortranABesselJn_8(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref>, i32, i32, f64, f64, f64, !fir.ref, i32) -> none + ! ALL-NEXT: } + ! ALL-NEXT: } + ! ALL-NEXT: } + r = bessel_jn(n1, n2, x) + ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref>>> + ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box>>) -> !fir.heap> + ! ALL: fir.freemem %[[addr]] : !fir.heap> +end subroutine test_transformational_real8 Index: flang/test/Lower/Intrinsics/bessel_yn.f90 =================================================================== --- flang/test/Lower/Intrinsics/bessel_yn.f90 +++ flang/test/Lower/Intrinsics/bessel_yn.f90 @@ -5,20 +5,112 @@ ! 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{{.*}}, %[[argn:.*]]: !fir.ref{{.*}}) -> f32 function test_real4(x, n) real :: x, test_real4 integer :: n + + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref + ! 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{{.*}}, %[[argn:.*]]: !fir.ref{{.*}}) -> f64 function test_real8(x, n) real(8) :: x, test_real8 integer :: n + + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n:.*]] = fir.load %[[argn]] : !fir.ref + ! 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{{.*}}, %[[argn1:.*]]: !fir.ref{{.*}}, %[[argn2:.*]]: !fir.ref{{.*}} +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>> + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref + ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref + ! 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>, i32, i32, !fir.ref, 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>, i32, i32, f32, f32, f32, !fir.ref, 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>, i32, i32, f32, f32, f32, !fir.ref, i32) -> none + ! ALL-NEXT: } else { + ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}} + ! ALL: fir.call @_FortranABesselYn_4(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref>, i32, i32, f32, f32, f32, !fir.ref, i32) -> none + ! ALL-NEXT: } + ! ALL-NEXT: } + ! ALL-NEXT: } + r = bessel_yn(n1, n2, x) + ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref>>> + ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box>>) -> !fir.heap> + ! ALL: fir.freemem %[[addr]] : !fir.heap> +end subroutine test_transformational_real4 + +! ALL-LABEL: @_QPtest_transformational_real8 +! ALL-SAME: %[[argx:.*]]: !fir.ref{{.*}}, %[[argn1:.*]]: !fir.ref{{.*}}, %[[argn2:.*]]: !fir.ref{{.*}} +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>> + ! ALL-DAG: %[[x:.*]] = fir.load %[[argx]] : !fir.ref + ! ALL-DAG: %[[n1:.*]] = fir.load %[[argn1]] : !fir.ref + ! ALL-DAG: %[[n2:.*]] = fir.load %[[argn2]] : !fir.ref + ! 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>, i32, i32, !fir.ref, 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>, i32, i32, f64, f64, f64, !fir.ref, 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>, i32, i32, f64, f64, f64, !fir.ref, i32) -> none + ! ALL-NEXT: } else { + ! ALL-DAG: %[[resn1gtn2:.*]] = fir.convert %[[r]] {{.*}} + ! ALL: fir.call @_FortranABesselYn_8(%[[resn1gtn2]], %[[n1]], %[[n2]], %[[x]], %[[zero]], %[[zero]], {{.*}}, {{.*}}) {{.*}} : (!fir.ref>, i32, i32, f64, f64, f64, !fir.ref, i32) -> none + ! ALL-NEXT: } + ! ALL-NEXT: } + ! ALL-NEXT: } + r = bessel_yn(n1, n2, x) + ! ALL: %[[box:.*]] = fir.load %[[r]] : !fir.ref>>> + ! ALL: %[[addr:.*]] = fir.box_addr %[[box]] : (!fir.box>>) -> !fir.heap> + ! ALL: fir.freemem %[[addr]] : !fir.heap> +end subroutine test_transformational_real8 Index: flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp =================================================================== --- flang/unittests/Optimizer/Builder/Runtime/TransformationalTest.cpp +++ 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(loc, seqTy); + mlir::Value n1 = builder.create(loc, i32Ty); + mlir::Value n2 = builder.create(loc, i32Ty); + mlir::Value x = builder.create(loc, realTy); + mlir::Value bn1 = builder.create(loc, realTy); + mlir::Value bn2 = builder.create(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(loc, seqTy); + mlir::Value n1 = builder.create(loc, i32Ty); + mlir::Value n2 = builder.create(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(loc, seqTy); + mlir::Value n1 = builder.create(loc, i32Ty); + mlir::Value n2 = builder.create(loc, i32Ty); + mlir::Value x = builder.create(loc, realTy); + mlir::Value bn1 = builder.create(loc, realTy); + mlir::Value bn2 = builder.create(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(loc, seqTy); + mlir::Value n1 = builder.create(loc, i32Ty); + mlir::Value n2 = builder.create(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 = Index: flang/unittests/Runtime/Transformational.cpp =================================================================== --- flang/unittests/Runtime/Transformational.cpp +++ flang/unittests/Runtime/Transformational.cpp @@ -10,10 +10,215 @@ #include "gtest/gtest.h" #include "tools.h" #include "flang/Runtime/type-code.h" +#include using namespace Fortran::runtime; using Fortran::common::TypeCategory; +template +using BesselFuncType = std::function, CppTypeFor, + CppTypeFor, const char *, int)>; + +template +using BesselX0FuncType = + std::function; + +template +constexpr CppTypeFor + besselEpsilon = CppTypeFor(1e-4); + +template +static void testBesselJn(BesselFuncType rtFunc, int32_t n1, int32_t n2, + CppTypeFor x, + const std::vector> &expected) { + StaticDescriptor<1> desc; + Descriptor &result{desc.descriptor()}; + unsigned len = expected.size(); + + CppTypeFor anc0 = len > 0 ? expected[len - 1] : 0.0; + CppTypeFor 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)); + 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>( + j)), + expected[j], besselEpsilon); + } +} + +template +static void testBesselJnX0( + BesselX0FuncType 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>( + 0)), + (n1 == 0) ? 1.0 : 0.0, 1e-5); + + for (int j{1}; j < (n2 - n1 + 1); ++j) { + EXPECT_NEAR( + (*result.ZeroBasedIndexedElement>( + j)), + 0.0, besselEpsilon); + } +} + +template static void testBesselJn(BesselFuncType rtFunc) { + testBesselJn(rtFunc, 1, 0, 1.0, {}); + testBesselJn(rtFunc, 0, 0, 1.0, {0.765197694}); + testBesselJn(rtFunc, 0, 1, 1.0, {0.765197694, 0.440050572}); + testBesselJn( + rtFunc, 0, 2, 1.0, {0.765197694, 0.440050572, 0.114903487}); + testBesselJn(rtFunc, 1, 5, 5.0, + {-0.327579111, 0.046565145, 0.364831239, 0.391232371, 0.261140555}); +} + +template static void testBesselJnX0(BesselX0FuncType rtFunc) { + testBesselJnX0(rtFunc, 1, 0); + testBesselJnX0(rtFunc, 0, 0); + testBesselJnX0(rtFunc, 1, 1); + testBesselJnX0(rtFunc, 0, 3); + testBesselJnX0(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 +static void testBesselYn(BesselFuncType rtFunc, int32_t n1, int32_t n2, + CppTypeFor x, + const std::vector> &expected) { + StaticDescriptor<1> desc; + Descriptor &result{desc.descriptor()}; + unsigned len = expected.size(); + + CppTypeFor anc0 = len > 0 ? expected[0] : 0.0; + CppTypeFor 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)); + 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>( + j)), + expected[j], besselEpsilon); + } +} + +template +static void testBesselYnX0( + BesselX0FuncType 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>( + j)), + (-std::numeric_limits< + CppTypeFor>::infinity())); + } +} + +template static void testBesselYn(BesselFuncType rtFunc) { + testBesselYn(rtFunc, 1, 0, 1.0, {}); + testBesselYn(rtFunc, 0, 0, 1.0, {0.08825695}); + testBesselYn(rtFunc, 0, 1, 1.0, {0.08825695, -0.7812128}); + testBesselYn(rtFunc, 0, 2, 1.0, {0.0882569555, -0.7812128, -1.6506826}); + testBesselYn(rtFunc, 1, 5, 1.0, + {-0.7812128, -1.6506826, -5.8215175, -33.278423, -260.40588}); +} + +template static void testBesselYnX0(BesselX0FuncType rtFunc) { + testBesselYnX0(rtFunc, 1, 0); + testBesselYnX0(rtFunc, 0, 0); + testBesselYnX0(rtFunc, 1, 1); + testBesselYnX0(rtFunc, 0, 3); + testBesselYnX0(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