486 lines · cpp
1//===-- lib/Evaluate/fold-real.cpp ----------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9#include "fold-implementation.h"10#include "fold-matmul.h"11#include "fold-reduction.h"12 13namespace Fortran::evaluate {14 15template <typename T>16static Expr<T> FoldTransformationalBessel(17 FunctionRef<T> &&funcRef, FoldingContext &context) {18 CHECK(funcRef.arguments().size() == 3);19 /// Bessel runtime functions use `int` integer arguments. Convert integer20 /// arguments to Int4, any overflow error will be reported during the21 /// conversion folding.22 using Int4 = Type<TypeCategory::Integer, 4>;23 if (auto args{GetConstantArguments<Int4, Int4, T>(24 context, funcRef.arguments(), /*hasOptionalArgument=*/false)}) {25 const std::string &name{std::get<SpecificIntrinsic>(funcRef.proc().u).name};26 if (auto elementalBessel{GetHostRuntimeWrapper<T, Int4, T>(name)}) {27 std::vector<Scalar<T>> results;28 int n1{static_cast<int>(29 std::get<0>(*args)->GetScalarValue().value().ToInt64())};30 int n2{static_cast<int>(31 std::get<1>(*args)->GetScalarValue().value().ToInt64())};32 Scalar<T> x{std::get<2>(*args)->GetScalarValue().value()};33 for (int i{n1}; i <= n2; ++i) {34 results.emplace_back((*elementalBessel)(context, Scalar<Int4>{i}, x));35 }36 return Expr<T>{Constant<T>{37 std::move(results), ConstantSubscripts{std::max(n2 - n1 + 1, 0)}}};38 } else {39 context.Warn(common::UsageWarning::FoldingFailure,40 "%s(integer(kind=4), real(kind=%d)) cannot be folded on host"_warn_en_US,41 name, T::kind);42 }43 }44 return Expr<T>{std::move(funcRef)};45}46 47// NORM248template <int KIND> class Norm2Accumulator {49 using T = Type<TypeCategory::Real, KIND>;50 51public:52 Norm2Accumulator(53 const Constant<T> &array, const Constant<T> &maxAbs, Rounding rounding)54 : array_{array}, maxAbs_{maxAbs}, rounding_{rounding} {};55 void operator()(56 Scalar<T> &element, const ConstantSubscripts &at, bool /*first*/) {57 // Summation of scaled elements:58 // Naively,59 // NORM2(A(:)) = SQRT(SUM(A(:)**2))60 // For any T > 0, we have mathematically61 // SQRT(SUM(A(:)**2))62 // = SQRT(T**2 * (SUM(A(:)**2) / T**2))63 // = SQRT(T**2 * SUM(A(:)**2 / T**2))64 // = SQRT(T**2 * SUM((A(:)/T)**2))65 // = SQRT(T**2) * SQRT(SUM((A(:)/T)**2))66 // = T * SQRT(SUM((A(:)/T)**2))67 // By letting T = MAXVAL(ABS(A)), we ensure that68 // ALL(ABS(A(:)/T) <= 1), so ALL((A(:)/T)**2 <= 1), and the SUM will69 // not overflow unless absolutely necessary.70 auto scale{maxAbs_.At(maxAbsAt_)};71 if (scale.IsZero()) {72 // Maximum value is zero, and so will the result be.73 // Avoid division by zero below.74 element = scale;75 } else {76 auto item{array_.At(at)};77 auto scaled{item.Divide(scale).value};78 auto square{scaled.Multiply(scaled).value};79 if constexpr (useKahanSummation) {80 auto next{square.Subtract(correction_, rounding_)};81 overflow_ |= next.flags.test(RealFlag::Overflow);82 auto sum{element.Add(next.value, rounding_)};83 overflow_ |= sum.flags.test(RealFlag::Overflow);84 correction_ = sum.value.Subtract(element, rounding_)85 .value.Subtract(next.value, rounding_)86 .value;87 element = sum.value;88 } else {89 auto sum{element.Add(square, rounding_)};90 overflow_ |= sum.flags.test(RealFlag::Overflow);91 element = sum.value;92 }93 }94 }95 bool overflow() const { return overflow_; }96 void Done(Scalar<T> &result) {97 // incoming result = SUM((data(:)/maxAbs)**2)98 // outgoing result = maxAbs * SQRT(result)99 auto root{result.SQRT().value};100 auto product{root.Multiply(maxAbs_.At(maxAbsAt_))};101 maxAbs_.IncrementSubscripts(maxAbsAt_);102 overflow_ |= product.flags.test(RealFlag::Overflow);103 result = product.value;104 }105 106private:107 const Constant<T> &array_;108 const Constant<T> &maxAbs_;109 const Rounding rounding_;110 bool overflow_{false};111 Scalar<T> correction_{};112 ConstantSubscripts maxAbsAt_{maxAbs_.lbounds()};113};114 115template <int KIND>116static Expr<Type<TypeCategory::Real, KIND>> FoldNorm2(FoldingContext &context,117 FunctionRef<Type<TypeCategory::Real, KIND>> &&funcRef) {118 using T = Type<TypeCategory::Real, KIND>;119 using Element = typename Constant<T>::Element;120 std::optional<int> dim;121 if (std::optional<ArrayAndMask<T>> arrayAndMask{122 ProcessReductionArgs<T>(context, funcRef.arguments(), dim,123 /*X=*/0, /*DIM=*/1)}) {124 MaxvalMinvalAccumulator<T, /*ABS=*/true> maxAbsAccumulator{125 RelationalOperator::GT, context, arrayAndMask->array};126 const Element identity{};127 Constant<T> maxAbs{DoReduction<T>(arrayAndMask->array, arrayAndMask->mask,128 dim, identity, maxAbsAccumulator)};129 Norm2Accumulator norm2Accumulator{arrayAndMask->array, maxAbs,130 context.targetCharacteristics().roundingMode()};131 Constant<T> result{DoReduction<T>(arrayAndMask->array, arrayAndMask->mask,132 dim, identity, norm2Accumulator)};133 if (norm2Accumulator.overflow()) {134 context.Warn(common::UsageWarning::FoldingException,135 "NORM2() of REAL(%d) data overflowed"_warn_en_US, KIND);136 }137 return Expr<T>{std::move(result)};138 }139 return Expr<T>{std::move(funcRef)};140}141 142template <int KIND>143Expr<Type<TypeCategory::Real, KIND>> FoldIntrinsicFunction(144 FoldingContext &context,145 FunctionRef<Type<TypeCategory::Real, KIND>> &&funcRef) {146 using T = Type<TypeCategory::Real, KIND>;147 using ComplexT = Type<TypeCategory::Complex, KIND>;148 using Int4 = Type<TypeCategory::Integer, 4>;149 ActualArguments &args{funcRef.arguments()};150 auto *intrinsic{std::get_if<SpecificIntrinsic>(&funcRef.proc().u)};151 CHECK(intrinsic);152 std::string name{intrinsic->name};153 if (name == "acos" || name == "acosh" || name == "asin" || name == "asinh" ||154 (name == "atan" && args.size() == 1) || name == "atanh" ||155 name == "bessel_j0" || name == "bessel_j1" || name == "bessel_y0" ||156 name == "bessel_y1" || name == "cos" || name == "cosh" || name == "erf" ||157 name == "erfc" || name == "erfc_scaled" || name == "exp" ||158 name == "gamma" || name == "log" || name == "log10" ||159 name == "log_gamma" || name == "sin" || name == "sinh" || name == "tan" ||160 name == "tanh") {161 CHECK(args.size() == 1);162 if (auto callable{GetHostRuntimeWrapper<T, T>(name)}) {163 return FoldElementalIntrinsic<T, T>(164 context, std::move(funcRef), *callable);165 } else {166 context.Warn(common::UsageWarning::FoldingFailure,167 "%s(real(kind=%d)) cannot be folded on host"_warn_en_US, name, KIND);168 }169 } else if (name == "amax0" || name == "amin0" || name == "amin1" ||170 name == "amax1" || name == "dmin1" || name == "dmax1") {171 return RewriteSpecificMINorMAX(context, std::move(funcRef));172 } else if (name == "atan" || name == "atan2") {173 std::string localName{name == "atan" ? "atan2" : name};174 CHECK(args.size() == 2);175 if (auto callable{GetHostRuntimeWrapper<T, T, T>(localName)}) {176 return FoldElementalIntrinsic<T, T, T>(177 context, std::move(funcRef), *callable);178 } else {179 context.Warn(common::UsageWarning::FoldingFailure,180 "%s(real(kind=%d), real(kind%d)) cannot be folded on host"_warn_en_US,181 name, KIND, KIND);182 }183 } else if (name == "bessel_jn" || name == "bessel_yn") {184 if (args.size() == 2) { // elemental185 // runtime functions use int arg186 if (auto callable{GetHostRuntimeWrapper<T, Int4, T>(name)}) {187 return FoldElementalIntrinsic<T, Int4, T>(188 context, std::move(funcRef), *callable);189 } else {190 context.Warn(common::UsageWarning::FoldingFailure,191 "%s(integer(kind=4), real(kind=%d)) cannot be folded on host"_warn_en_US,192 name, KIND);193 }194 } else {195 return FoldTransformationalBessel<T>(std::move(funcRef), context);196 }197 } else if (name == "abs") { // incl. zabs & cdabs198 // Argument can be complex or real199 if (UnwrapExpr<Expr<SomeReal>>(args[0])) {200 return FoldElementalIntrinsic<T, T>(201 context, std::move(funcRef), &Scalar<T>::ABS);202 } else if (UnwrapExpr<Expr<SomeComplex>>(args[0])) {203 return FoldElementalIntrinsic<T, ComplexT>(context, std::move(funcRef),204 ScalarFunc<T, ComplexT>([&name, &context](205 const Scalar<ComplexT> &z) -> Scalar<T> {206 ValueWithRealFlags<Scalar<T>> y{z.ABS()};207 if (y.flags.test(RealFlag::Overflow)) {208 context.Warn(common::UsageWarning::FoldingException,209 "complex ABS intrinsic folding overflow"_warn_en_US, name);210 }211 return y.value;212 }));213 } else {214 common::die(" unexpected argument type inside abs");215 }216 } else if (name == "aimag") {217 if (auto *zExpr{UnwrapExpr<Expr<ComplexT>>(args[0])}) {218 return Fold(context, Expr<T>{ComplexComponent{true, std::move(*zExpr)}});219 }220 } else if (name == "aint" || name == "anint") {221 // ANINT rounds ties away from zero, not to even222 common::RoundingMode mode{name == "aint"223 ? common::RoundingMode::ToZero224 : common::RoundingMode::TiesAwayFromZero};225 return FoldElementalIntrinsic<T, T>(context, std::move(funcRef),226 ScalarFunc<T, T>(227 [&name, &context, mode](const Scalar<T> &x) -> Scalar<T> {228 ValueWithRealFlags<Scalar<T>> y{x.ToWholeNumber(mode)};229 if (y.flags.test(RealFlag::Overflow)) {230 context.Warn(common::UsageWarning::FoldingException,231 "%s intrinsic folding overflow"_warn_en_US, name);232 }233 return y.value;234 }));235 } else if (name == "dim") {236 return FoldElementalIntrinsic<T, T, T>(context, std::move(funcRef),237 ScalarFunc<T, T, T>([&context](const Scalar<T> &x,238 const Scalar<T> &y) -> Scalar<T> {239 ValueWithRealFlags<Scalar<T>> result{x.DIM(y)};240 if (result.flags.test(RealFlag::Overflow)) {241 context.Warn(common::UsageWarning::FoldingException,242 "DIM intrinsic folding overflow"_warn_en_US);243 }244 return result.value;245 }));246 } else if (name == "dot_product") {247 return FoldDotProduct<T>(context, std::move(funcRef));248 } else if (name == "dprod") {249 // Rewrite DPROD(x,y) -> DBLE(x)*DBLE(y)250 if (args.at(0) && args.at(1)) {251 const auto *xExpr{args[0]->UnwrapExpr()};252 const auto *yExpr{args[1]->UnwrapExpr()};253 if (xExpr && yExpr) {254 return Fold(context,255 ToReal<T::kind>(context, common::Clone(*xExpr)) *256 ToReal<T::kind>(context, common::Clone(*yExpr)));257 }258 }259 } else if (name == "epsilon") {260 return Expr<T>{Scalar<T>::EPSILON()};261 } else if (name == "fraction") {262 return FoldElementalIntrinsic<T, T>(context, std::move(funcRef),263 ScalarFunc<T, T>(264 [](const Scalar<T> &x) -> Scalar<T> { return x.FRACTION(); }));265 } else if (name == "huge") {266 return Expr<T>{Scalar<T>::HUGE()};267 } else if (name == "hypot") {268 CHECK(args.size() == 2);269 return FoldElementalIntrinsic<T, T, T>(context, std::move(funcRef),270 ScalarFunc<T, T, T>(271 [&](const Scalar<T> &x, const Scalar<T> &y) -> Scalar<T> {272 ValueWithRealFlags<Scalar<T>> result{x.HYPOT(y)};273 if (result.flags.test(RealFlag::Overflow)) {274 context.Warn(common::UsageWarning::FoldingException,275 "HYPOT intrinsic folding overflow"_warn_en_US);276 }277 return result.value;278 }));279 } else if (name == "matmul") {280 return FoldMatmul(context, std::move(funcRef));281 } else if (name == "max") {282 return FoldMINorMAX(context, std::move(funcRef), Ordering::Greater);283 } else if (name == "maxval") {284 return FoldMaxvalMinval<T>(context, std::move(funcRef),285 RelationalOperator::GT, T::Scalar::HUGE().Negate());286 } else if (name == "min") {287 return FoldMINorMAX(context, std::move(funcRef), Ordering::Less);288 } else if (name == "minval") {289 return FoldMaxvalMinval<T>(290 context, std::move(funcRef), RelationalOperator::LT, T::Scalar::HUGE());291 } else if (name == "mod") {292 CHECK(args.size() == 2);293 bool badPConst{false};294 if (auto *pExpr{UnwrapExpr<Expr<T>>(args[1])}) {295 *pExpr = Fold(context, std::move(*pExpr));296 if (auto pConst{GetScalarConstantValue<T>(*pExpr)};297 pConst && pConst->IsZero()) {298 context.Warn(common::UsageWarning::FoldingAvoidsRuntimeCrash,299 "MOD: P argument is zero"_warn_en_US);300 badPConst = true;301 }302 }303 return FoldElementalIntrinsic<T, T, T>(context, std::move(funcRef),304 ScalarFunc<T, T, T>([&context, badPConst](const Scalar<T> &x,305 const Scalar<T> &y) -> Scalar<T> {306 auto result{x.MOD(y)};307 if (!badPConst && result.flags.test(RealFlag::DivideByZero)) {308 context.Warn(common::UsageWarning::FoldingAvoidsRuntimeCrash,309 "second argument to MOD must not be zero"_warn_en_US);310 }311 return result.value;312 }));313 } else if (name == "modulo") {314 CHECK(args.size() == 2);315 bool badPConst{false};316 if (auto *pExpr{UnwrapExpr<Expr<T>>(args[1])}) {317 *pExpr = Fold(context, std::move(*pExpr));318 if (auto pConst{GetScalarConstantValue<T>(*pExpr)};319 pConst && pConst->IsZero()) {320 context.Warn(common::UsageWarning::FoldingAvoidsRuntimeCrash,321 "MODULO: P argument is zero"_warn_en_US);322 badPConst = true;323 }324 }325 return FoldElementalIntrinsic<T, T, T>(context, std::move(funcRef),326 ScalarFunc<T, T, T>([&context, badPConst](const Scalar<T> &x,327 const Scalar<T> &y) -> Scalar<T> {328 auto result{x.MODULO(y)};329 if (!badPConst && result.flags.test(RealFlag::DivideByZero)) {330 context.Warn(common::UsageWarning::FoldingAvoidsRuntimeCrash,331 "second argument to MODULO must not be zero"_warn_en_US);332 }333 return result.value;334 }));335 } else if (name == "nearest") {336 if (auto *sExpr{UnwrapExpr<Expr<SomeReal>>(args[1])}) {337 *sExpr = Fold(context, std::move(*sExpr));338 return common::visit(339 [&](const auto &sVal) {340 using TS = ResultType<decltype(sVal)>;341 bool badSConst{false};342 if (auto sConst{GetScalarConstantValue<TS>(sVal)};343 sConst && (sConst->IsZero() || sConst->IsNotANumber())) {344 context.Warn(common::UsageWarning::FoldingValueChecks,345 "NEAREST: S argument is %s"_warn_en_US,346 sConst->IsZero() ? "zero" : "NaN");347 badSConst = true;348 }349 return FoldElementalIntrinsic<T, T, TS>(context, std::move(funcRef),350 ScalarFunc<T, T, TS>([&](const Scalar<T> &x,351 const Scalar<TS> &s) -> Scalar<T> {352 if (!badSConst && (s.IsZero() || s.IsNotANumber())) {353 context.Warn(common::UsageWarning::FoldingValueChecks,354 "NEAREST: S argument is %s"_warn_en_US,355 s.IsZero() ? "zero" : "NaN");356 }357 auto result{x.NEAREST(!s.IsNegative())};358 if (result.flags.test(RealFlag::InvalidArgument)) {359 context.Warn(common::UsageWarning::FoldingException,360 "NEAREST intrinsic folding: bad argument"_warn_en_US);361 }362 return result.value;363 }));364 },365 sExpr->u);366 }367 } else if (name == "norm2") {368 return FoldNorm2<T::kind>(context, std::move(funcRef));369 } else if (name == "product") {370 auto one{Scalar<T>::FromInteger(value::Integer<8>{1}).value};371 return FoldProduct<T>(context, std::move(funcRef), one);372 } else if (name == "real" || name == "dble") {373 if (auto *expr{args[0].value().UnwrapExpr()}) {374 return ToReal<KIND>(context, std::move(*expr));375 }376 } else if (name == "rrspacing") {377 return FoldElementalIntrinsic<T, T>(context, std::move(funcRef),378 ScalarFunc<T, T>(379 [](const Scalar<T> &x) -> Scalar<T> { return x.RRSPACING(); }));380 } else if (name == "scale") {381 if (const auto *byExpr{UnwrapExpr<Expr<SomeInteger>>(args[1])}) {382 return common::visit(383 [&](const auto &byVal) {384 using TBY = ResultType<decltype(byVal)>;385 return FoldElementalIntrinsic<T, T, TBY>(context,386 std::move(funcRef),387 ScalarFunc<T, T, TBY>(388 [&](const Scalar<T> &x, const Scalar<TBY> &y) -> Scalar<T> {389 ValueWithRealFlags<Scalar<T>> result{390 x.391// MSVC chokes on the keyword "template" here in a call to a392// member function template.393#ifndef _MSC_VER394 template395#endif396 SCALE<Scalar<TBY>>(y)};397 if (result.flags.test(RealFlag::Overflow)) {398 context.Warn(common::UsageWarning::FoldingException,399 "SCALE/IEEE_SCALB intrinsic folding overflow"_warn_en_US);400 }401 return result.value;402 }));403 },404 byExpr->u);405 }406 } else if (name == "set_exponent") {407 if (const auto *iExpr{UnwrapExpr<Expr<SomeInteger>>(args[1])}) {408 return common::visit(409 [&](const auto &iVal) {410 using TY = ResultType<decltype(iVal)>;411 return FoldElementalIntrinsic<T, T, TY>(context, std::move(funcRef),412 ScalarFunc<T, T, TY>(413 [&](const Scalar<T> &x, const Scalar<TY> &i) -> Scalar<T> {414 return x.SET_EXPONENT(i.ToInt64());415 }));416 },417 iExpr->u);418 }419 } else if (name == "sign") {420 return FoldElementalIntrinsic<T, T, T>(421 context, std::move(funcRef), &Scalar<T>::SIGN);422 } else if (name == "spacing") {423 return FoldElementalIntrinsic<T, T>(context, std::move(funcRef),424 ScalarFunc<T, T>(425 [](const Scalar<T> &x) -> Scalar<T> { return x.SPACING(); }));426 } else if (name == "sqrt") {427 return FoldElementalIntrinsic<T, T>(context, std::move(funcRef),428 ScalarFunc<T, T>(429 [](const Scalar<T> &x) -> Scalar<T> { return x.SQRT().value; }));430 } else if (name == "sum") {431 return FoldSum<T>(context, std::move(funcRef));432 } else if (name == "tiny") {433 return Expr<T>{Scalar<T>::TINY()};434 } else if (name == "__builtin_fma") {435 CHECK(args.size() == 3);436 } else if (name == "__builtin_ieee_next_after") {437 if (const auto *yExpr{UnwrapExpr<Expr<SomeReal>>(args[1])}) {438 return common::visit(439 [&](const auto &yVal) {440 using TY = ResultType<decltype(yVal)>;441 return FoldElementalIntrinsic<T, T, TY>(context, std::move(funcRef),442 ScalarFunc<T, T, TY>([&](const Scalar<T> &x,443 const Scalar<TY> &y) -> Scalar<T> {444 auto xBig{Scalar<LargestReal>::Convert(x).value};445 auto yBig{Scalar<LargestReal>::Convert(y).value};446 switch (xBig.Compare(yBig)) {447 case Relation::Unordered:448 context.Warn(common::UsageWarning::FoldingValueChecks,449 "IEEE_NEXT_AFTER intrinsic folding: arguments are unordered"_warn_en_US);450 return x.NotANumber();451 case Relation::Equal:452 break;453 case Relation::Less:454 return x.NEAREST(true).value;455 case Relation::Greater:456 return x.NEAREST(false).value;457 }458 return x; // dodge bogus "missing return" GCC warning459 }));460 },461 yExpr->u);462 }463 } else if (name == "__builtin_ieee_next_up" ||464 name == "__builtin_ieee_next_down") {465 bool upward{name == "__builtin_ieee_next_up"};466 const char *iName{upward ? "IEEE_NEXT_UP" : "IEEE_NEXT_DOWN"};467 return FoldElementalIntrinsic<T, T>(context, std::move(funcRef),468 ScalarFunc<T, T>([&](const Scalar<T> &x) -> Scalar<T> {469 auto result{x.NEAREST(upward)};470 if (result.flags.test(RealFlag::InvalidArgument)) {471 context.Warn(common::UsageWarning::FoldingException,472 "%s intrinsic folding: argument is NaN"_warn_en_US, iName);473 }474 return result.value;475 }));476 }477 return Expr<T>{std::move(funcRef)};478}479 480#ifdef _MSC_VER // disable bogus warning about missing definitions481#pragma warning(disable : 4661)482#endif483FOR_EACH_REAL_KIND(template class ExpressionBase, )484template class ExpressionBase<SomeReal>;485} // namespace Fortran::evaluate486