fold-implementation.h

//===-- lib/Evaluate/fold-implementation.h --------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
 
#ifndef FORTRAN_EVALUATE_FOLD_IMPLEMENTATION_H_
#define FORTRAN_EVALUATE_FOLD_IMPLEMENTATION_H_
 
#include "character.h"
#include "host.h"
#include "int-power.h"
#include "flang/Common/indirection.h"
#include "flang/Common/template.h"
#include "flang/Common/unwrap.h"
#include "flang/Evaluate/characteristics.h"
#include "flang/Evaluate/check-expression.h"
#include "flang/Evaluate/common.h"
#include "flang/Evaluate/constant.h"
#include "flang/Evaluate/expression.h"
#include "flang/Evaluate/fold.h"
#include "flang/Evaluate/intrinsics-library.h"
#include "flang/Evaluate/intrinsics.h"
#include "flang/Evaluate/shape.h"
#include "flang/Evaluate/tools.h"
#include "flang/Evaluate/traverse.h"
#include "flang/Evaluate/type.h"
#include "flang/Parser/message.h"
#include "flang/Semantics/scope.h"
#include "flang/Semantics/symbol.h"
#include "flang/Semantics/tools.h"
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <optional>
#include <type_traits>
#include <variant>
 
// Some environments, viz. glibc 2.17 and *BSD, allow the macro HUGE
// to leak out of <math.h>.
#undef HUGE
 
namespace Fortran::evaluate {
 
// Don't use Kahan extended precision summation any more when folding
// transformational intrinsic functions other than SUM, since it is
// not used in the runtime implementations of those functions and we
// want results to match.
static constexpr bool useKahanSummation{false};
 
// Utilities
template <typename T> class Folder {
public:
  explicit Folder(int kind, FoldingContext &c, bool forOptionalArgument = false)
      : kind_{kind}, context_{c}, forOptionalArgument_{forOptionalArgument} {}
  template <typename U = T,
      typename =
          std::enable_if_t<std::is_same_v<U, Type<TypeCategory::Derived>>>>
  explicit Folder(FoldingContext &c, bool forOptionalArgument = false)
      : Folder(/*kind=*/0, c, forOptionalArgument) {}
 
  std::optional<Constant<T>> GetNamedConstant(const Symbol &);
  std::optional<Constant<T>> ApplySubscripts(const Constant<T> &array,
      const std::vector<Constant<SubscriptInteger>> &subscripts);
  std::optional<Constant<T>> ApplyComponent(
      Constant<Type<TypeCategory::Derived>> &&, const Symbol &component,
      const std::vector<Constant<SubscriptInteger>> * = nullptr);
  std::optional<Constant<T>> GetConstantComponent(
      Component &, const std::vector<Constant<SubscriptInteger>> * = nullptr);
  std::optional<Constant<T>> Folding(ArrayRef &);
  std::optional<Constant<T>> Folding(DataRef &);
  Expr<T> Folding(Designator<T> &&);
  Constant<T> *Folding(std::optional<ActualArgument> &);
 
  Expr<T> CSHIFT(FunctionRef<T> &&);
  Expr<T> EOSHIFT(FunctionRef<T> &&);
  Expr<T> MERGE(FunctionRef<T> &&);
  Expr<T> PACK(FunctionRef<T> &&);
  Expr<T> RESHAPE(FunctionRef<T> &&);
  Expr<T> SPREAD(FunctionRef<T> &&);
  Expr<T> TRANSPOSE(FunctionRef<T> &&);
  Expr<T> UNPACK(FunctionRef<T> &&);
 
  Expr<T> TRANSFER(FunctionRef<T> &&);
 
private:
  int kind_;
  FoldingContext &context_;
  bool forOptionalArgument_{false};
};
 
std::optional<Constant<SubscriptInteger>> GetConstantSubscript(
    FoldingContext &, Subscript &, const NamedEntity &, int dim);
 
template <typename TR, typename... TA, std::size_t... I>
std::optional<std::function<Scalar<TR>(FoldingContext &, Scalar<TA>...)>>
GetHostRuntimeWrapperHelper(int resultKind,
    std::array<int, sizeof...(TA)> argKinds, const std::string &name,
    std::index_sequence<I...>) {
  std::vector<DynamicType> argTypes{DynamicType{TA::category, argKinds[I]}...};
  if (auto hostWrapper{GetHostRuntimeWrapper(
          name, DynamicType{TR::category, resultKind}, argTypes)}) {
    return [hostWrapper, resultKind, argKinds](
               FoldingContext &context, Scalar<TA>... args) -> Scalar<TR> {
      std::vector<Expr<SomeType>> genericArgs{
          AsGenericExpr(Constant<TA>{argKinds[I], args})...};
      return GetScalarConstantValue<TR>(
          resultKind, (*hostWrapper)(context, std::move(genericArgs)))
          .value();
    };
  }
  return std::nullopt;
}
 
// Helper to use host runtime on scalars for folding.
template <typename TR, typename... TA>
static std::optional<std::function<Scalar<TR>(FoldingContext &, Scalar<TA>...)>>
GetHostRuntimeWrapper(int resultKind, std::array<int, sizeof...(TA)> argKinds,
    const std::string &name) {
  return GetHostRuntimeWrapperHelper<TR, TA...>(
      resultKind, argKinds, name, std::index_sequence_for<TA...>{});
}
 
// FoldOperation() rewrites expression tree nodes.
// If there is any possibility that the rewritten node will
// not have the same representation type, the result of
// FoldOperation() will be packaged in an Expr<> of the same
// specific type.
 
// no-op base case
template <typename A>
common::IfNoLvalue<Expr<ResultType<A>>, A> FoldOperation(
    FoldingContext &, A &&x) {
  static_assert(!std::is_same_v<A, Expr<ResultType<A>>>,
      "call Fold() instead for Expr<>");
  return Expr<ResultType<A>>{std::move(x)};
}
 
// Forward declarations needed to ensure overload lookup considers all possible
// implementations.
ComplexPart FoldOperation(FoldingContext &, ComplexPart &&);
Expr<Type<TypeCategory::Real>> FoldOperation(
    FoldingContext &, ComplexComponent &&);
Expr<Type<TypeCategory::Complex>> FoldOperation(
    FoldingContext &, ComplexConstructor &&);
Component FoldOperation(FoldingContext &, Component &&);
NamedEntity FoldOperation(FoldingContext &, NamedEntity &&);
Triplet FoldOperation(FoldingContext &, Triplet &&);
Subscript FoldOperation(FoldingContext &, Subscript &&);
ArrayRef FoldOperation(FoldingContext &, ArrayRef &&);
CoarrayRef FoldOperation(FoldingContext &, CoarrayRef &&);
DataRef FoldOperation(FoldingContext &, DataRef &&);
Substring FoldOperation(FoldingContext &, Substring &&);
template <typename T>
Expr<T> FoldOperation(FoldingContext &, FunctionRef<T> &&);
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Designator<T> &&designator) {
  const int kind{designator.kind()};
  return Folder<T>{kind, context}.Folding(std::move(designator));
}
Expr<TypeParamInquiry::Result> FoldOperation(
    FoldingContext &, TypeParamInquiry &&);
Expr<RankOneBoundElement::Result> FoldOperation(
    FoldingContext &, RankOneBoundElement &&);
Expr<ImpliedDoIndex::Result> FoldOperation(
    FoldingContext &context, ImpliedDoIndex &&);
template <typename T>
Expr<T> FoldOperation(FoldingContext &, ArrayConstructor<T> &&);
Expr<Type<TypeCategory::Derived>> FoldOperation(
    FoldingContext &, StructureConstructor &&);
template <typename T>
Expr<T> FoldOperation(FoldingContext &, ConditionalExpr<T> &&);
 
template <typename T>
std::optional<Constant<T>> Folder<T>::GetNamedConstant(const Symbol &symbol0) {
  const Symbol &symbol{ResolveAssociations(symbol0)};
  if (IsNamedConstant(symbol)) {
    if (const auto *object{
            symbol.detailsIf<semantics::ObjectEntityDetails>()}) {
      if (const auto *constant{UnwrapConstantValue<T>(kind_, object->init())}) {
        return *constant;
      }
    }
  }
  return std::nullopt;
}
 
template <typename T>
std::optional<Constant<T>> Folder<T>::Folding(ArrayRef &aRef) {
  std::vector<Constant<SubscriptInteger>> subscripts;
  int dim{0};
  for (Subscript &ss : aRef.subscript()) {
    if (auto constant{GetConstantSubscript(context_, ss, aRef.base(), dim++)}) {
      subscripts.emplace_back(std::move(*constant));
    } else {
      return std::nullopt;
    }
  }
  if (Component * component{aRef.base().UnwrapComponent()}) {
    return GetConstantComponent(*component, &subscripts);
  } else if (std::optional<Constant<T>> array{
                 GetNamedConstant(aRef.base().GetLastSymbol())}) {
    return ApplySubscripts(*array, subscripts);
  } else {
    return std::nullopt;
  }
}
 
template <typename T>
std::optional<Constant<T>> Folder<T>::Folding(DataRef &ref) {
  return common::visit(
      common::visitors{
          [this](SymbolRef &sym) { return GetNamedConstant(*sym); },
          [this](Component &comp) {
            comp = FoldOperation(context_, std::move(comp));
            return GetConstantComponent(comp);
          },
          [this](ArrayRef &aRef) {
            aRef = FoldOperation(context_, std::move(aRef));
            return Folding(aRef);
          },
          [](CoarrayRef &) { return std::optional<Constant<T>>{}; },
      },
      ref.u);
}
 
// TODO: This would be more natural as a member function of Constant<T>.
template <typename T>
std::optional<Constant<T>> Folder<T>::ApplySubscripts(const Constant<T> &array,
    const std::vector<Constant<SubscriptInteger>> &subscripts) {
  const int kind{array.kind()};
  const auto &shape{array.shape()};
  const auto &lbounds{array.lbounds()};
  int rank{GetRank(shape)};
  CHECK(rank == static_cast<int>(subscripts.size()));
  std::size_t elements{1};
  ConstantSubscripts resultShape;
  ConstantSubscripts ssLB;
  for (const auto &ss : subscripts) {
    if (ss.Rank() == 1) {
      resultShape.push_back(static_cast<ConstantSubscript>(ss.size()));
      elements *= ss.size();
      ssLB.push_back(ss.lbounds().front());
    } else if (ss.Rank() > 1) {
      return std::nullopt; // error recovery
    }
  }
  ConstantSubscripts ssAt(rank, 0), at(rank, 0), tmp(1, 0);
  std::vector<Scalar<T>> values;
  while (elements-- > 0) {
    bool increment{true};
    int k{0};
    for (int j{0}; j < rank; ++j) {
      if (subscripts[j].Rank() == 0) {
        at[j] = subscripts[j].GetScalarValue().value().ToInt64();
      } else {
        CHECK(k < GetRank(resultShape));
        tmp[0] = ssLB.at(k) + ssAt.at(k);
        at[j] = subscripts[j].At(tmp).ToInt64();
        if (increment) {
          if (++ssAt[k] == resultShape[k]) {
            ssAt[k] = 0;
          } else {
            increment = false;
          }
        }
        ++k;
      }
      if (at[j] < lbounds[j] || at[j] >= lbounds[j] + shape[j]) {
        context_.messages().Say(
            "Subscript value (%jd) is out of range on dimension %d in reference to a constant array value"_err_en_US,
            at[j], j + 1);
        return std::nullopt;
      }
    }
    values.emplace_back(array.At(at));
    CHECK(!increment || elements == 0);
    CHECK(k == GetRank(resultShape));
  }
  if constexpr (T::category == TypeCategory::Character) {
    return Constant<T>{
        kind, array.LEN(), std::move(values), std::move(resultShape)};
  } else if constexpr (std::is_same_v<T, Type<TypeCategory::Derived>>) {
    CHECK(kind == 0);
    return Constant<T>{array.result().derivedTypeSpec(), std::move(values),
        std::move(resultShape)};
  } else {
    return Constant<T>{kind, std::move(values), std::move(resultShape)};
  }
}
 
template <typename T>
std::optional<Constant<T>> Folder<T>::ApplyComponent(
    Constant<Type<TypeCategory::Derived>> &&structures, const Symbol &component,
    const std::vector<Constant<SubscriptInteger>> *subscripts) {
  if (auto scalar{structures.GetScalarValue()}) {
    if (std::optional<Expr<SomeType>> expr{scalar->Find(component)}) {
      if (const Constant<T> *value{UnwrapConstantValue<T>(kind_, *expr)}) {
        if (subscripts) {
          return ApplySubscripts(*value, *subscripts);
        } else {
          return *value;
        }
      }
    }
  } else {
    // A(:)%scalar_component & A(:)%array_component(subscripts)
    std::unique_ptr<ArrayConstructor<T>> array;
    if (structures.empty()) {
      return std::nullopt;
    }
    ConstantSubscripts at{structures.lbounds()};
    do {
      StructureConstructor scalar{structures.At(at)};
      if (std::optional<Expr<SomeType>> expr{scalar.Find(component)}) {
        if (const Constant<T> *value{
                UnwrapConstantValue<T>(kind_, expr.value())}) {
          if (!array.get()) {
            // This technique ensures that character length or derived type
            // information is propagated to the array constructor.
            auto *typedExpr{UnwrapExpr<Expr<T>>(kind_, expr.value())};
            CHECK(typedExpr);
            array = std::make_unique<ArrayConstructor<T>>(kind_, *typedExpr);
            if constexpr (T::category == TypeCategory::Character) {
              array->set_LEN(MakeSubscriptIntExpr(value->LEN()));
            }
          }
          if (subscripts) {
            if (auto element{ApplySubscripts(*value, *subscripts)}) {
              CHECK(element->Rank() == 0);
              array->Push(Expr<T>{std::move(*element)});
            } else {
              return std::nullopt;
            }
          } else {
            CHECK(value->Rank() == 0);
            array->Push(Expr<T>{*value});
          }
        } else {
          return std::nullopt;
        }
      }
    } while (structures.IncrementSubscripts(at));
    // Fold the ArrayConstructor<> into a Constant<>.
    CHECK(array);
    Expr<T> result{Fold(context_, Expr<T>{std::move(*array)})};
    if (auto *constant{UnwrapConstantValue<T>(kind_, result)}) {
      return constant->Reshape(common::Clone(structures.shape()));
    }
  }
  return std::nullopt;
}
 
template <typename T>
std::optional<Constant<T>> Folder<T>::GetConstantComponent(Component &component,
    const std::vector<Constant<SubscriptInteger>> *subscripts) {
  if (std::optional<Constant<Type<TypeCategory::Derived>>> structures{
          common::visit(
              common::visitors{
                  [&](const Symbol &symbol) {
                    return Folder<Type<TypeCategory::Derived>>{context_}
                        .GetNamedConstant(symbol);
                  },
                  [&](ArrayRef &aRef) {
                    return Folder<Type<TypeCategory::Derived>>{context_}
                        .Folding(aRef);
                  },
                  [&](Component &base) {
                    return Folder<Type<TypeCategory::Derived>>{context_}
                        .GetConstantComponent(base);
                  },
                  [&](CoarrayRef &) {
                    return std::optional<
                        Constant<Type<TypeCategory::Derived>>>{};
                  },
              },
              component.base().u)}) {
    return ApplyComponent(
        std::move(*structures), component.GetLastSymbol(), subscripts);
  } else {
    return std::nullopt;
  }
}
 
template <typename T> Expr<T> Folder<T>::Folding(Designator<T> &&designator) {
  const int kind{designator.kind()};
  if constexpr (T::category == TypeCategory::Character) {
    if (auto *substring{common::Unwrap<Substring>(designator.u)}) {
      if (std::optional<Expr<Type<TypeCategory::Character>>> folded{
              substring->Fold(context_)};
          folded && folded->kind() == kind) {
        return std::move(*folded);
      }
      // We used to fold zero-length substrings into zero-length
      // constants here, but that led to problems in variable
      // definition contexts.
    }
  } else if constexpr (T::category == TypeCategory::Real) {
    if (auto *zPart{std::get_if<ComplexPart>(&designator.u)}) {
      *zPart = FoldOperation(context_, std::move(*zPart));
      using ComplexT = Type<TypeCategory::Complex>;
      if (auto zConst{
              Folder<ComplexT>{kind, context_}.Folding(zPart->complex())}) {
        return Fold(context_,
            Expr<T>{ComplexComponent{zPart->part() == ComplexPart::Part::IM,
                Expr<ComplexT>{std::move(*zConst)}}});
      } else {
        return Expr<T>{Designator<T>{kind, std::move(*zPart)}};
      }
    }
  }
  return common::visit(
      common::visitors{
          [&](SymbolRef &&symbol) {
            if (auto constant{GetNamedConstant(*symbol)}) {
              return Expr<T>{std::move(*constant)};
            }
            return Expr<T>{std::move(designator)};
          },
          [&, kind](ArrayRef &&aRef) {
            aRef = FoldOperation(context_, std::move(aRef));
            if (auto c{Folding(aRef)}) {
              return Expr<T>{std::move(*c)};
            } else {
              return Expr<T>{Designator<T>{kind, std::move(aRef)}};
            }
          },
          [&, kind](Component &&component) {
            component = FoldOperation(context_, std::move(component));
            if (auto c{GetConstantComponent(component)}) {
              return Expr<T>{std::move(*c)};
            } else {
              return Expr<T>{Designator<T>{kind, std::move(component)}};
            }
          },
          [&, kind](auto &&x) {
            return Expr<T>{
                Designator<T>{kind, FoldOperation(context_, std::move(x))}};
          },
      },
      std::move(designator.u));
}
 
// Apply type conversion and re-folding if necessary.
// This is where BOZ arguments are converted.
template <typename T>
Constant<T> *Folder<T>::Folding(std::optional<ActualArgument> &arg) {
  if (auto *expr{UnwrapExpr<Expr<SomeType>>(arg)}) {
    *expr = Fold(context_, std::move(*expr));
    if constexpr (T::category != TypeCategory::Derived) {
      if (!UnwrapExpr<Expr<T>>(kind_, *expr)) {
        if (const Symbol *
                var{forOptionalArgument_
                        ? UnwrapWholeSymbolOrComponentDataRef(*expr)
                        : nullptr};
            var && (IsOptional(*var) || IsAllocatableOrObjectPointer(var))) {
          // can't safely convert item that may not be present
        } else if (auto converted{ConvertToType(
                       DynamicType{T::category, kind_}, std::move(*expr))}) {
          *expr = Fold(context_, std::move(*converted));
        }
      }
    }
    return UnwrapConstantValue<T>(kind_, *expr);
  }
  return nullptr;
}
 
template <typename... A, std::size_t... I>
std::optional<std::tuple<const Constant<A> *...>> GetConstantArgumentsHelper(
    const std::array<int, sizeof...(A)> &kinds, FoldingContext &context,
    ActualArguments &arguments, bool hasOptionalArgument,
    std::index_sequence<I...>) {
  static_assert(sizeof...(A) > 0);
  std::tuple<const Constant<A> *...> args{
      Folder<A>{kinds[I], context, hasOptionalArgument}.Folding(
          arguments.at(I))...};
  if ((... && (std::get<I>(args)))) {
    return args;
  } else {
    return std::nullopt;
  }
}
 
template <typename... A>
std::optional<std::tuple<const Constant<A> *...>> GetConstantArguments(
    const std::array<int, sizeof...(A)> &kinds, FoldingContext &context,
    ActualArguments &args, bool hasOptionalArgument) {
  return GetConstantArgumentsHelper<A...>(kinds, context, args,
      hasOptionalArgument, std::index_sequence_for<A...>{});
}
 
template <typename... A, std::size_t... I>
std::optional<std::tuple<Scalar<A>...>> GetScalarConstantArgumentsHelper(
    const std::array<int, sizeof...(A)> &kinds, FoldingContext &context,
    ActualArguments &args, bool hasOptionalArgument,
    std::index_sequence<I...>) {
  if (auto constArgs{GetConstantArguments<A...>(
          kinds, context, args, hasOptionalArgument)}) {
    return std::tuple<Scalar<A>...>{
        std::get<I>(*constArgs)->GetScalarValue().value()...};
  } else {
    return std::nullopt;
  }
}
 
template <typename... A>
std::optional<std::tuple<Scalar<A>...>> GetScalarConstantArguments(
    const std::array<int, sizeof...(A)> &kinds, FoldingContext &context,
    ActualArguments &args, bool hasOptionalArgument) {
  return GetScalarConstantArgumentsHelper<A...>(kinds, context, args,
      hasOptionalArgument, std::index_sequence_for<A...>{});
}
 
// helpers to fold intrinsic function references
// Define callable types used in a common utility that
// takes care of array and cast/conversion aspects for elemental intrinsics
 
template <typename TR, typename... TArgs>
using ScalarFunc = std::function<Scalar<TR>(const Scalar<TArgs> &...)>;
template <typename TR, typename... TArgs>
using ScalarFuncWithContext =
    std::function<Scalar<TR>(FoldingContext &, const Scalar<TArgs> &...)>;
 
template <template <typename, typename...> typename WrapperType, typename TR,
    typename... TA, std::size_t... I>
Expr<TR> FoldElementalIntrinsicHelper(int resultKind,
    const std::array<int, sizeof...(TA)> &argKinds, FoldingContext &context,
    FunctionRef<TR> &&funcRef, WrapperType<TR, TA...> func,
    bool hasOptionalArgument, std::index_sequence<I...>) {
  CHECK(funcRef.kind() == resultKind);
  if (std::optional<std::tuple<const Constant<TA> *...>> args{
          GetConstantArguments<TA...>({argKinds[I]...}, context,
              funcRef.arguments(), hasOptionalArgument)}) {
    // Compute the shape of the result based on shapes of arguments
    ConstantSubscripts shape;
    int rank{0};
    const ConstantSubscripts *shapes[]{&std::get<I>(*args)->shape()...};
    const int ranks[]{std::get<I>(*args)->Rank()...};
    for (unsigned int i{0}; i < sizeof...(TA); ++i) {
      if (ranks[i] > 0) {
        if (rank == 0) {
          rank = ranks[i];
          shape = *shapes[i];
        } else {
          if (shape != *shapes[i]) {
            // TODO: Rank compatibility was already checked but it seems to be
            // the first place where the actual shapes are checked to be the
            // same. Shouldn't this be checked elsewhere so that this is also
            // checked for non constexpr call to elemental intrinsics function?
            context.messages().Say(
                "Arguments in elemental intrinsic function are not conformable"_err_en_US);
            return Expr<TR>{std::move(funcRef)};
          }
        }
      }
    }
    CHECK(rank == GetRank(shape));
    // Compute all the scalar values of the results
    std::vector<Scalar<TR>> results;
    std::optional<uint64_t> n{TotalElementCount(shape)};
    if (!n) {
      context.messages().Say(
          "Too many elements in elemental intrinsic function result"_err_en_US);
      return Expr<TR>{std::move(funcRef)};
    }
    if (*n > 0) {
      ConstantBounds bounds{shape};
      ConstantSubscripts resultIndex(rank, 1);
      ConstantSubscripts argIndex[]{std::get<I>(*args)->lbounds()...};
      do {
        if constexpr (std::is_same_v<WrapperType<TR, TA...>,
                          ScalarFuncWithContext<TR, TA...>>) {
          results.emplace_back(
              func(context, std::get<I>(*args)->At(argIndex[I])...));
        } else if constexpr (std::is_same_v<WrapperType<TR, TA...>,
                                 ScalarFunc<TR, TA...>>) {
          results.emplace_back(func(std::get<I>(*args)->At(argIndex[I])...));
        }
        (std::get<I>(*args)->IncrementSubscripts(argIndex[I]), ...);
      } while (bounds.IncrementSubscripts(resultIndex));
    }
    // Build and return constant result
    if constexpr (TR::category == TypeCategory::Character) {
      auto len{static_cast<ConstantSubscript>(
          results.empty() ? 0 : results[0].length())};
      return Expr<TR>{
          Constant<TR>{resultKind, len, std::move(results), std::move(shape)}};
    } else if constexpr (TR::category == TypeCategory::Derived) {
      if (!results.empty()) {
        return Expr<TR>{rank == 0
                ? Constant<TR>{results.front()}
                : Constant<TR>{results.front().derivedTypeSpec(),
                      std::move(results), std::move(shape)}};
      }
    } else {
      return Expr<TR>{
          Constant<TR>{resultKind, std::move(results), std::move(shape)}};
    }
  }
  return Expr<TR>{std::move(funcRef)};
}
 
template <typename TR, typename... TA>
Expr<TR> FoldElementalIntrinsic(int resultKind,
    const std::array<int, sizeof...(TA)> &argKinds, FoldingContext &context,
    FunctionRef<TR> &&funcRef, ScalarFunc<TR, TA...> func,
    bool hasOptionalArgument = false) {
  return FoldElementalIntrinsicHelper<ScalarFunc, TR, TA...>(resultKind,
      argKinds, context, std::move(funcRef), func, hasOptionalArgument,
      std::index_sequence_for<TA...>{});
}
template <typename TR, typename... TA>
Expr<TR> FoldElementalIntrinsic(int resultKind,
    const std::array<int, sizeof...(TA)> &argKinds, FoldingContext &context,
    FunctionRef<TR> &&funcRef, ScalarFuncWithContext<TR, TA...> func,
    bool hasOptionalArgument = false) {
  return FoldElementalIntrinsicHelper<ScalarFuncWithContext, TR, TA...>(
      resultKind, argKinds, context, std::move(funcRef), func,
      hasOptionalArgument, std::index_sequence_for<TA...>{});
}
 
std::optional<std::int64_t> GetInt64ArgOr(
    const std::optional<ActualArgument> &, std::int64_t defaultValue);
 
template <typename A, typename B>
std::optional<std::vector<A>> GetIntegerVector(const B &x) {
  static_assert(common::is_integral_v<A>);
  if (const auto *someInteger{
          UnwrapExpr<Expr<Type<TypeCategory::Integer>>>(x)}) {
    if (const auto *constant{
            UnwrapConstantValue<Type<TypeCategory::Integer>>(*someInteger)}) {
      if (constant->Rank() == 1) {
        std::vector<A> result;
        for (const auto &value : constant->values()) {
          result.push_back(static_cast<A>(value.ToInt64()));
        }
        return result;
      }
    }
  }
  return std::nullopt;
}
 
// Transform an intrinsic function reference that contains user errors
// into an intrinsic with the same characteristic but the "invalid" name.
// This to prevent generating warnings over and over if the expression
// gets re-folded.
template <typename T> Expr<T> MakeInvalidIntrinsic(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  SpecificIntrinsic invalid{std::get<SpecificIntrinsic>(funcRef.proc().u)};
  invalid.name = IntrinsicProcTable::InvalidName;
  return Expr<T>{FunctionRef<T>{kind, ProcedureDesignator{std::move(invalid)},
      ActualArguments{std::move(funcRef.arguments())}}};
}
 
template <typename T> Expr<T> Folder<T>::CSHIFT(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  auto args{funcRef.arguments()};
  CHECK(args.size() == 3);
  const auto *array{UnwrapConstantValue<T>(kind, args[0])};
  const auto *shiftExpr{UnwrapExpr<Expr<Type<TypeCategory::Integer>>>(args[1])};
  auto dim{GetInt64ArgOr(args[2], 1)};
  if (!array || !shiftExpr || !dim) {
    return Expr<T>{std::move(funcRef)};
  }
  auto convertedShift{Fold(context_,
      ConvertToType<SubscriptInteger>(SubscriptIntegerKind,
          Expr<Type<TypeCategory::Integer>>{*shiftExpr}))};
  const auto *shift{UnwrapConstantValue<SubscriptInteger>(
      SubscriptIntegerKind, convertedShift)};
  if (!shift) {
    return Expr<T>{std::move(funcRef)};
  }
  // Arguments are constant
  if (*dim < 1 || *dim > array->Rank()) {
    context_.messages().Say("Invalid 'dim=' argument (%jd) in CSHIFT"_err_en_US,
        static_cast<std::intmax_t>(*dim));
  } else if (shift->Rank() > 0 && shift->Rank() != array->Rank() - 1) {
    // message already emitted from intrinsic look-up
  } else {
    int rank{array->Rank()};
    int zbDim{static_cast<int>(*dim) - 1};
    bool ok{true};
    if (shift->Rank() > 0) {
      int k{0};
      for (int j{0}; j < rank; ++j) {
        if (j != zbDim) {
          if (array->shape()[j] != shift->shape()[k]) {
            context_.messages().Say(
                "Invalid 'shift=' argument in CSHIFT: extent on dimension %d is %jd but must be %jd"_err_en_US,
                k + 1, static_cast<std::intmax_t>(shift->shape()[k]),
                static_cast<std::intmax_t>(array->shape()[j]));
            ok = false;
          }
          ++k;
        }
      }
    }
    if (ok) {
      std::vector<Scalar<T>> resultElements;
      ConstantSubscripts arrayLB{array->lbounds()};
      ConstantSubscripts arrayAt{arrayLB};
      ConstantSubscript &dimIndex{arrayAt[zbDim]};
      ConstantSubscript dimLB{dimIndex}; // initial value
      ConstantSubscript dimExtent{array->shape()[zbDim]};
      ConstantSubscripts shiftLB{shift->lbounds()};
      for (auto n{GetSize(array->shape())}; n > 0; --n) {
        ConstantSubscript origDimIndex{dimIndex};
        ConstantSubscripts shiftAt;
        if (shift->Rank() > 0) {
          int k{0};
          for (int j{0}; j < rank; ++j) {
            if (j != zbDim) {
              shiftAt.emplace_back(shiftLB[k++] + arrayAt[j] - arrayLB[j]);
            }
          }
        }
        ConstantSubscript shiftCount{shift->At(shiftAt).ToInt64()};
        dimIndex = dimLB + ((dimIndex - dimLB + shiftCount) % dimExtent);
        if (dimIndex < dimLB) {
          dimIndex += dimExtent;
        } else if (dimIndex >= dimLB + dimExtent) {
          dimIndex -= dimExtent;
        }
        resultElements.push_back(array->At(arrayAt));
        dimIndex = origDimIndex;
        array->IncrementSubscripts(arrayAt);
      }
      return Expr<T>{PackageConstant<T>(
          kind, std::move(resultElements), *array, array->shape())};
    }
  }
  // Invalid, prevent re-folding
  return MakeInvalidIntrinsic(std::move(funcRef));
}
 
template <typename T> Expr<T> Folder<T>::EOSHIFT(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  auto args{funcRef.arguments()};
  CHECK(args.size() == 4);
  const auto *array{UnwrapConstantValue<T>(kind, args[0])};
  const auto *shiftExpr{UnwrapExpr<Expr<Type<TypeCategory::Integer>>>(args[1])};
  auto dim{GetInt64ArgOr(args[3], 1)};
  if (!array || !shiftExpr || !dim) {
    return Expr<T>{std::move(funcRef)};
  }
  // Apply type conversions to the shift= and boundary= arguments.
  auto convertedShift{Fold(context_,
      ConvertToType<SubscriptInteger>(SubscriptIntegerKind,
          Expr<Type<TypeCategory::Integer>>{*shiftExpr}))};
  const auto *shift{UnwrapConstantValue<SubscriptInteger>(
      SubscriptIntegerKind, convertedShift)};
  if (!shift) {
    return Expr<T>{std::move(funcRef)};
  }
  const Constant<T> *boundary{nullptr};
  std::optional<Expr<SomeType>> convertedBoundary;
  if (const auto *boundaryExpr{UnwrapExpr<Expr<SomeType>>(args[2])}) {
    convertedBoundary = Fold(context_,
        ConvertToType(array->GetType(), Expr<SomeType>{*boundaryExpr}));
    boundary = UnwrapExpr<Constant<T>>(kind, convertedBoundary);
    if (!boundary) {
      return Expr<T>{std::move(funcRef)};
    }
  }
  // Arguments are constant
  if (*dim < 1 || *dim > array->Rank()) {
    context_.messages().Say(
        "Invalid 'dim=' argument (%jd) in EOSHIFT"_err_en_US,
        static_cast<std::intmax_t>(*dim));
  } else if (shift->Rank() > 0 && shift->Rank() != array->Rank() - 1) {
    // message already emitted from intrinsic look-up
  } else if (boundary && boundary->Rank() > 0 &&
      boundary->Rank() != array->Rank() - 1) {
    // ditto
  } else {
    int rank{array->Rank()};
    int zbDim{static_cast<int>(*dim) - 1};
    bool ok{true};
    if (shift->Rank() > 0) {
      int k{0};
      for (int j{0}; j < rank; ++j) {
        if (j != zbDim) {
          if (array->shape()[j] != shift->shape()[k]) {
            context_.messages().Say(
                "Invalid 'shift=' argument in EOSHIFT: extent on dimension %d is %jd but must be %jd"_err_en_US,
                k + 1, static_cast<std::intmax_t>(shift->shape()[k]),
                static_cast<std::intmax_t>(array->shape()[j]));
            ok = false;
          }
          ++k;
        }
      }
    }
    if (boundary && boundary->Rank() > 0) {
      int k{0};
      for (int j{0}; j < rank; ++j) {
        if (j != zbDim) {
          if (array->shape()[j] != boundary->shape()[k]) {
            context_.messages().Say(
                "Invalid 'boundary=' argument in EOSHIFT: extent on dimension %d is %jd but must be %jd"_err_en_US,
                k + 1, static_cast<std::intmax_t>(boundary->shape()[k]),
                static_cast<std::intmax_t>(array->shape()[j]));
            ok = false;
          }
          ++k;
        }
      }
    }
    if (ok) {
      std::vector<Scalar<T>> resultElements;
      ConstantSubscripts arrayLB{array->lbounds()};
      ConstantSubscripts arrayAt{arrayLB};
      ConstantSubscript &dimIndex{arrayAt[zbDim]};
      ConstantSubscript dimLB{dimIndex}; // initial value
      ConstantSubscript dimExtent{array->shape()[zbDim]};
      ConstantSubscripts shiftLB{shift->lbounds()};
      ConstantSubscripts boundaryLB;
      if (boundary) {
        boundaryLB = boundary->lbounds();
      }
      for (auto n{GetSize(array->shape())}; n > 0; --n) {
        ConstantSubscript origDimIndex{dimIndex};
        ConstantSubscripts shiftAt;
        if (shift->Rank() > 0) {
          int k{0};
          for (int j{0}; j < rank; ++j) {
            if (j != zbDim) {
              shiftAt.emplace_back(shiftLB[k++] + arrayAt[j] - arrayLB[j]);
            }
          }
        }
        ConstantSubscript shiftCount{shift->At(shiftAt).ToInt64()};
        dimIndex += shiftCount;
        if (dimIndex >= dimLB && dimIndex < dimLB + dimExtent) {
          resultElements.push_back(array->At(arrayAt));
        } else if (boundary) {
          ConstantSubscripts boundaryAt;
          if (boundary->Rank() > 0) {
            for (int j{0}; j < rank; ++j) {
              int k{0};
              if (j != zbDim) {
                boundaryAt.emplace_back(
                    boundaryLB[k++] + arrayAt[j] - arrayLB[j]);
              }
            }
          }
          resultElements.push_back(boundary->At(boundaryAt));
        } else if constexpr (T::category == TypeCategory::Integer ||
            T::category == TypeCategory::Unsigned ||
            T::category == TypeCategory::Real ||
            T::category == TypeCategory::Complex ||
            T::category == TypeCategory::Logical) {
          resultElements.emplace_back(Scalar<T>::Zero(kind));
        } else if constexpr (T::category == TypeCategory::Character) {
          auto len{static_cast<std::size_t>(array->LEN())};
          resultElements.emplace_back(kind, len, ' ');
        } else {
          DIE("no derived type boundary");
        }
        dimIndex = origDimIndex;
        array->IncrementSubscripts(arrayAt);
      }
      return Expr<T>{PackageConstant<T>(
          kind, std::move(resultElements), *array, array->shape())};
    }
  }
  // Invalid, prevent re-folding
  return MakeInvalidIntrinsic(std::move(funcRef));
}
 
template <typename T> Expr<T> Folder<T>::MERGE(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  return FoldElementalIntrinsic<T, T, T, LogicalResult>(kind,
      {kind, kind, LogicalResultKind}, context_, std::move(funcRef),
      ScalarFunc<T, T, T, LogicalResult>(
          [](const Scalar<T> &ifTrue, const Scalar<T> &ifFalse,
              const Scalar<LogicalResult> &predicate) -> Scalar<T> {
            return predicate.IsTrue() ? ifTrue : ifFalse;
          }));
}
 
template <typename T> Expr<T> Folder<T>::PACK(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  auto args{funcRef.arguments()};
  CHECK(args.size() == 3);
  const auto *array{UnwrapConstantValue<T>(kind, args[0])};
  const auto *vector{UnwrapConstantValue<T>(kind, args[2])};
  auto convertedMask{Fold(context_,
      ConvertToType<LogicalResult>(LogicalResultKind,
          Expr<Type<TypeCategory::Logical>>{
              DEREF(UnwrapExpr<Expr<Type<TypeCategory::Logical>>>(args[1]))}))};
  const auto *mask{
      UnwrapConstantValue<LogicalResult>(LogicalResultKind, convertedMask)};
  if (!array || !mask || (args[2] && !vector)) {
    return Expr<T>{std::move(funcRef)};
  }
  // Arguments are constant.
  ConstantSubscript arrayElements{GetSize(array->shape())};
  ConstantSubscript truths{0};
  ConstantSubscripts maskAt{mask->lbounds()};
  if (mask->Rank() == 0) {
    if (mask->At(maskAt).IsTrue()) {
      truths = arrayElements;
    }
  } else if (array->shape() != mask->shape()) {
    // Error already emitted from intrinsic processing
    return MakeInvalidIntrinsic(std::move(funcRef));
  } else {
    for (ConstantSubscript j{0}; j < arrayElements;
         ++j, mask->IncrementSubscripts(maskAt)) {
      if (mask->At(maskAt).IsTrue()) {
        ++truths;
      }
    }
  }
  std::vector<Scalar<T>> resultElements;
  ConstantSubscripts arrayAt{array->lbounds()};
  ConstantSubscript resultSize{truths};
  if (vector) {
    resultSize = vector->shape().at(0);
    if (resultSize < truths) {
      context_.messages().Say(
          "Invalid 'vector=' argument in PACK: the 'mask=' argument has %jd true elements, but the vector has only %jd elements"_err_en_US,
          static_cast<std::intmax_t>(truths),
          static_cast<std::intmax_t>(resultSize));
      return MakeInvalidIntrinsic(std::move(funcRef));
    }
  }
  for (ConstantSubscript j{0}; j < truths;) {
    if (mask->At(maskAt).IsTrue()) {
      resultElements.push_back(array->At(arrayAt));
      ++j;
    }
    array->IncrementSubscripts(arrayAt);
    mask->IncrementSubscripts(maskAt);
  }
  if (vector) {
    ConstantSubscripts vectorAt{vector->lbounds()};
    vectorAt.at(0) += truths;
    for (ConstantSubscript j{truths}; j < resultSize; ++j) {
      resultElements.push_back(vector->At(vectorAt));
      ++vectorAt[0];
    }
  }
  return Expr<T>{PackageConstant<T>(kind, std::move(resultElements), *array,
      ConstantSubscripts{static_cast<ConstantSubscript>(resultSize)})};
}
 
template <typename T> Expr<T> Folder<T>::RESHAPE(FunctionRef<T> &&funcRef) {
  auto args{funcRef.arguments()};
  CHECK(args.size() == 4);
  const auto *source{UnwrapConstantValue<T>(funcRef.kind(), args[0])};
  const auto *pad{UnwrapConstantValue<T>(funcRef.kind(), args[2])};
  std::optional<std::vector<ConstantSubscript>> shape{
      GetIntegerVector<ConstantSubscript>(args[1])};
  std::optional<std::vector<int>> order{GetIntegerVector<int>(args[3])};
  std::optional<uint64_t> optResultElement;
  std::optional<std::vector<int>> dimOrder;
  bool ok{true};
  if (shape) {
    if (shape->size() > common::maxRank) {
      context_.messages().Say(
          "Size of 'shape=' argument (%zd) must not be greater than %d"_err_en_US,
          shape->size(), common::maxRank);
      ok = false;
    } else if (HasNegativeExtent(*shape)) {
      context_.messages().Say(
          "'shape=' argument (%s) must not have a negative extent"_err_en_US,
          DEREF(args[1]->UnwrapExpr()).AsFortran());
      ok = false;
    } else {
      optResultElement = TotalElementCount(*shape);
      if (!optResultElement) {
        context_.messages().Say(
            "'shape=' argument (%s) specifies an array with too many elements"_err_en_US,
            DEREF(args[1]->UnwrapExpr()).AsFortran());
        ok = false;
      }
    }
    if (order) {
      dimOrder = ValidateDimensionOrder(GetRank(*shape), *order);
      if (!dimOrder) {
        context_.messages().Say(
            "Invalid 'order=' argument (%s) in RESHAPE"_err_en_US,
            DEREF(args[3]->UnwrapExpr()).AsFortran());
        ok = false;
      }
    }
  }
  if (!ok) {
    // convert into an invalid intrinsic procedure call below
  } else if (!source || !shape || (args[2] && !pad) || (args[3] && !order)) {
    return Expr<T>{std::move(funcRef)}; // Non-constant arguments
  } else {
    uint64_t resultElements{*optResultElement};
    std::vector<int> *dimOrderPtr{dimOrder ? &dimOrder.value() : nullptr};
    if (resultElements > source->size() && (!pad || pad->empty())) {
      context_.messages().Say(
          "Too few elements in 'source=' argument and 'pad=' "
          "argument is not present or has null size"_err_en_US);
      ok = false;
    } else {
      Constant<T> result{!source->empty() || !pad
              ? source->Reshape(std::move(shape.value()))
              : pad->Reshape(std::move(shape.value()))};
      ConstantSubscripts subscripts{result.lbounds()};
      auto copied{result.CopyFrom(*source,
          std::min(static_cast<uint64_t>(source->size()), resultElements),
          subscripts, dimOrderPtr)};
      if (copied < resultElements) {
        CHECK(pad);
        copied += result.CopyFrom(
            *pad, resultElements - copied, subscripts, dimOrderPtr);
      }
      CHECK(copied == resultElements);
      return Expr<T>{std::move(result)};
    }
  }
  // Invalid, prevent re-folding
  return MakeInvalidIntrinsic(std::move(funcRef));
}
 
template <typename T> Expr<T> Folder<T>::SPREAD(FunctionRef<T> &&funcRef) {
  auto args{funcRef.arguments()};
  CHECK(args.size() == 3);
  const Constant<T> *source{UnwrapConstantValue<T>(funcRef.kind(), args[0])};
  auto dim{ToInt64(args[1])};
  auto ncopies{ToInt64(args[2])};
  if (!source || !dim) {
    return Expr<T>{std::move(funcRef)};
  }
  int sourceRank{source->Rank()};
  if (sourceRank >= common::maxRank) {
    context_.messages().Say(
        "SOURCE= argument to SPREAD has rank %d but must have rank less than %d"_err_en_US,
        sourceRank, common::maxRank);
  } else if (*dim < 1 || *dim > sourceRank + 1) {
    context_.messages().Say(
        "DIM=%d argument to SPREAD must be between 1 and %d"_err_en_US, *dim,
        sourceRank + 1);
  } else if (!ncopies) {
    return Expr<T>{std::move(funcRef)};
  } else {
    if (*ncopies < 0) {
      ncopies = 0;
    }
    // TODO: Consider moving this implementation (after the user error
    // checks), along with other transformational intrinsics, into
    // constant.h (or a new header) so that the transformationals
    // are available for all Constant<>s without needing to be packaged
    // as references to intrinsic functions for folding.
    ConstantSubscripts shape{source->shape()};
    shape.insert(shape.begin() + *dim - 1, *ncopies);
    Constant<T> spread{source->Reshape(std::move(shape))};
    std::optional<uint64_t> n{TotalElementCount(spread.shape())};
    if (!n) {
      context_.messages().Say("Too many elements in SPREAD result"_err_en_US);
    } else {
      std::vector<int> dimOrder;
      for (int j{0}; j < sourceRank; ++j) {
        dimOrder.push_back(j < *dim - 1 ? j : j + 1);
      }
      dimOrder.push_back(*dim - 1);
      ConstantSubscripts at{spread.lbounds()}; // all 1
      spread.CopyFrom(*source, *n, at, &dimOrder);
      return Expr<T>{std::move(spread)};
    }
  }
  // Invalid, prevent re-folding
  return MakeInvalidIntrinsic(std::move(funcRef));
}
 
template <typename T> Expr<T> Folder<T>::TRANSPOSE(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  auto args{funcRef.arguments()};
  CHECK(args.size() == 1);
  const auto *matrix{UnwrapConstantValue<T>(kind, args[0])};
  if (!matrix) {
    return Expr<T>{std::move(funcRef)};
  }
  // Argument is constant.  Traverse its elements in transposed order.
  std::vector<Scalar<T>> resultElements;
  ConstantSubscripts at(2);
  for (ConstantSubscript j{0}; j < matrix->shape()[0]; ++j) {
    at[0] = matrix->lbounds()[0] + j;
    for (ConstantSubscript k{0}; k < matrix->shape()[1]; ++k) {
      at[1] = matrix->lbounds()[1] + k;
      resultElements.push_back(matrix->At(at));
    }
  }
  at = matrix->shape();
  std::swap(at[0], at[1]);
  return Expr<T>{
      PackageConstant<T>(kind, std::move(resultElements), *matrix, at)};
}
 
template <typename T> Expr<T> Folder<T>::UNPACK(FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  auto args{funcRef.arguments()};
  CHECK(args.size() == 3);
  const auto *vector{UnwrapConstantValue<T>(kind, args[0])};
  auto convertedMask{Fold(context_,
      ConvertToType<LogicalResult>(LogicalResultKind,
          Expr<Type<TypeCategory::Logical>>{
              DEREF(UnwrapExpr<Expr<Type<TypeCategory::Logical>>>(args[1]))}))};
  const auto *mask{
      UnwrapConstantValue<LogicalResult>(LogicalResultKind, convertedMask)};
  const auto *field{UnwrapConstantValue<T>(kind, args[2])};
  if (!vector || !mask || !field) {
    return Expr<T>{std::move(funcRef)};
  }
  // Arguments are constant.
  if (field->Rank() > 0 && field->shape() != mask->shape()) {
    // Error already emitted from intrinsic processing
    return MakeInvalidIntrinsic(std::move(funcRef));
  }
  ConstantSubscript maskElements{GetSize(mask->shape())};
  ConstantSubscript truths{0};
  ConstantSubscripts maskAt{mask->lbounds()};
  for (ConstantSubscript j{0}; j < maskElements;
       ++j, mask->IncrementSubscripts(maskAt)) {
    if (mask->At(maskAt).IsTrue()) {
      ++truths;
    }
  }
  if (truths > GetSize(vector->shape())) {
    context_.messages().Say(
        "Invalid 'vector=' argument in UNPACK: the 'mask=' argument has %jd true elements, but the vector has only %jd elements"_err_en_US,
        static_cast<std::intmax_t>(truths),
        static_cast<std::intmax_t>(GetSize(vector->shape())));
    return MakeInvalidIntrinsic(std::move(funcRef));
  }
  std::vector<Scalar<T>> resultElements;
  ConstantSubscripts vectorAt{vector->lbounds()};
  ConstantSubscripts fieldAt{field->lbounds()};
  for (ConstantSubscript j{0}; j < maskElements; ++j) {
    if (mask->At(maskAt).IsTrue()) {
      resultElements.push_back(vector->At(vectorAt));
      vector->IncrementSubscripts(vectorAt);
    } else {
      resultElements.push_back(field->At(fieldAt));
    }
    mask->IncrementSubscripts(maskAt);
    field->IncrementSubscripts(fieldAt);
  }
  return Expr<T>{PackageConstant<T>(
      kind, std::move(resultElements), *vector, mask->shape())};
}
 
std::optional<Expr<SomeType>> FoldTransfer(
    FoldingContext &, const ActualArguments &);
 
template <typename T> Expr<T> Folder<T>::TRANSFER(FunctionRef<T> &&funcRef) {
  if (auto folded{FoldTransfer(context_, funcRef.arguments())}) {
    return DEREF(UnwrapExpr<Expr<T>>(funcRef.kind(), *folded));
  } else {
    return Expr<T>{std::move(funcRef)};
  }
}
 
// TODO: Once the backend supports character extremums we could support
// min/max with non-optional arguments to trees of extremum operations.
template <typename T>
Expr<T> FoldMINorMAX(
    FoldingContext &context, FunctionRef<T> &&funcRef, Ordering order) {
  const int kind{funcRef.kind()};
  static_assert(T::category == TypeCategory::Integer ||
      T::category == TypeCategory::Unsigned ||
      T::category == TypeCategory::Real ||
      T::category == TypeCategory::Character);
 
  // Lots of constraints:
  // - We want Extremum<T> generated by semantics to compare equal to
  //   Extremum<T> written out to module files as max or min calls.
  // - Users can also write min/max calls that must also compare equal
  //   to min/max calls that wind up being written to module files.
  // - Extremeum<T> is binary and can't currently handle processing
  //   optional arguments that may show up in 3rd + argument.
  // - The code below only accepts more than 2 arguments if all the
  //   arguments are constant (and hence known to be present).
  // - ConvertExprToHLFIR can't currently handle Extremum<Character>
  // - Semantics doesn't currently generate Extremum<Character>
  // The original code did the folding of arguments and the overall extremum
  // operation in a single pass. This was shorter code-wise, but took me
  // a while to tease out all the logic and was doing redundant work.
  // So I split it into two passes:
  // 1) fold the arguments and check if they are constant,
  // 2) Decide if we:
  //    - can constant-fold the min/max operation, or
  //    - need to generate an extremum anyway,
  //    and do it if so.
  //    Otherwise, return the original call.
  auto &args{funcRef.arguments()};
  std::size_t nargs{args.size()};
  bool allArgsConstant{true};
  bool extremumAnyway{nargs == 2 && T::category != TypeCategory::Character};
  // 1a)Fold the first two arguments.
  {
    Folder<T> folder{kind, context, /*forOptionalArgument=*/false};
    if (!folder.Folding(args[0])) {
      allArgsConstant = false;
    }
    if (!folder.Folding(args[1])) {
      allArgsConstant = false;
    }
  }
  // 1b) Fold any optional arguments.
  if (nargs > 2) {
    Folder<T> folder{kind, context, /*forOptionalArgument=*/true};
    for (std::size_t i{2}; i < nargs; ++i) {
      if (args[i]) {
        if (!folder.Folding(args[i])) {
          allArgsConstant = false;
        }
      }
    }
  }
  // 2) If we can fold the result or the call to min/max may compare equal to
  // an extremum generated by semantics go ahead and convert to an extremum,
  // and try to fold the result.
  if (allArgsConstant || extremumAnyway) {
    // Folding updates the argument expressions in place, no need to call
    // Fold() on each argument again.
    if (const auto *resultp{UnwrapExpr<Expr<T>>(kind, args[0])}) {
      Expr<T> result{*resultp};
      for (std::size_t i{1}; i < nargs; ++i) {
        if (const auto *tExpr{UnwrapExpr<Expr<T>>(kind, args[i])}) {
          result = FoldOperation(
              context, Extremum<T>{order, std::move(result), *tExpr});
        } else {
          // This should never happen, but here is a value to return.
          return Expr<T>{std::move(funcRef)};
        }
      }
      return result;
    }
  }
  // If we decided to not generate an extremum just return the original call,
  // with the arguments folded.
  return Expr<T>{std::move(funcRef)};
}
 
// For AMAX0, AMIN0, AMAX1, AMIN1, DMAX1, DMIN1, MAX0, MIN0, MAX1, and MIN1
// a special care has to be taken to insert the conversion on the result
// of the MIN/MAX. This is made slightly more complex by the extension
// supported by f18 that arguments may have different kinds. This implies
// that the created MIN/MAX result type cannot be deduced from the standard but
// has to be deduced from the arguments.
// e.g. AMAX0(int8, int4) is rewritten to REAL(MAX(int8, INT(int4, 8)))).
template <typename T>
Expr<T> RewriteSpecificMINorMAX(
    FoldingContext &context, FunctionRef<T> &&funcRef) {
  int kind{funcRef.kind()}; // not const: GCC 7 ICE in generic lambdas
  ActualArguments &args{funcRef.arguments()};
  auto &intrinsic{DEREF(std::get_if<SpecificIntrinsic>(&funcRef.proc().u))};
  // Rewrite MAX1(args) to INT(MAX(args)) and fold. Same logic for MIN1.
  // Find result type for max/min based on the arguments.
  std::optional<DynamicType> resultType;
  ActualArgument *resultTypeArg{nullptr};
  for (auto j{args.size()}; j-- > 0;) {
    if (args[j]) {
      DynamicType type{args[j]->GetType().value()};
      // Handle mixed real/integer arguments: all the previous arguments were
      // integers and this one is real. The type of the MAX/MIN result will
      // be the one of the real argument.
      if (!resultType ||
          (type.category() == resultType->category() &&
              type.kind() > resultType->kind()) ||
          resultType->category() == TypeCategory::Integer) {
        resultType = type;
        resultTypeArg = &*args[j];
      }
    }
  }
  if (!resultType) { // error recovery
    return Expr<T>{std::move(funcRef)};
  }
  intrinsic.name =
      intrinsic.name.find("max") != std::string::npos ? "max"s : "min"s;
  intrinsic.characteristics.value().functionResult.value().SetType(*resultType);
  auto insertConversion{[&](const auto &x) -> Expr<T> {
    using TR = ResultType<decltype(x)>;
    const int rKind{x.kind()};
    FunctionRef<TR> maxRef{
        rKind, ProcedureDesignator{funcRef.proc()}, ActualArguments{args}};
    return Fold(
        context, ConvertToType<T>(kind, AsCategoryExpr(std::move(maxRef))));
  }};
  if (auto *sx{UnwrapExpr<Expr<Type<TypeCategory::Real>>>(*resultTypeArg)}) {
    return insertConversion(*sx);
  } else if (auto *sx{UnwrapExpr<Expr<Type<TypeCategory::Integer>>>(
                 *resultTypeArg)}) {
    return insertConversion(*sx);
  } else {
    return Expr<T>{std::move(funcRef)}; // error recovery
  }
}
 
// FoldIntrinsicFunction()
Expr<Type<TypeCategory::Integer>> FoldIntrinsicFunction(
    FoldingContext &context, FunctionRef<Type<TypeCategory::Integer>> &&);
Expr<Type<TypeCategory::Unsigned>> FoldIntrinsicFunction(
    FoldingContext &context, FunctionRef<Type<TypeCategory::Unsigned>> &&);
Expr<Type<TypeCategory::Real>> FoldIntrinsicFunction(
    FoldingContext &context, FunctionRef<Type<TypeCategory::Real>> &&);
Expr<Type<TypeCategory::Complex>> FoldIntrinsicFunction(
    FoldingContext &context, FunctionRef<Type<TypeCategory::Complex>> &&);
Expr<Type<TypeCategory::Logical>> FoldIntrinsicFunction(
    FoldingContext &context, FunctionRef<Type<TypeCategory::Logical>> &&);
 
// Fold NEXT(enum) / PREVIOUS(enum) to a constant enumerator (scalar or array)
// when the argument is a constant enumeration value.  Returns the folded
// constant, or the original reference left unfolded (e.g. at a boundary,
// where error termination is deferred to run time).
Expr<Type<TypeCategory::Derived>> FoldEnumerationNextOrPrevious(
    FoldingContext &context, FunctionRef<Type<TypeCategory::Derived>> &&funcRef,
    bool isNext);
 
// Returns true when the j-th dummy argument of an intrinsic contributes only
// its declared type, kind type parameters, and rank to the result, such as X
// in PRECISION(X). The intrinsic table marks these dummies with
// ArgFlag::onlyConstantInquiry, and IsConstantExpr() uses the same attribute.
// STORAGE_SIZE and LEN aren't marked, because their results also depend on
// the dynamic type or on length type parameters.
bool IsTypeOnlyInquiryArg(const SpecificIntrinsic &intrinsic, std::size_t j);
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, FunctionRef<T> &&funcRef) {
  const int kind{funcRef.kind()};
  ActualArguments &args{funcRef.arguments()};
  const auto *intrinsic{std::get_if<SpecificIntrinsic>(&funcRef.proc().u)};
  if (!intrinsic || intrinsic->name != "kind") {
    // Don't fold the argument to KIND(); it might be a TypeParamInquiry
    // with a forced result type that doesn't match the parameter.
    for (std::optional<ActualArgument> &arg : args) {
      if (arg && arg->GetConditionalArg()) {
        FoldConditionalArg(context, arg);
      } else if (auto *expr{UnwrapExpr<Expr<SomeType>>(arg)}) {
        if (!intrinsic && IsNamedConstantDesignator(*expr)) {
          // A designator of a named constant retained as the actual argument
          // of a nonintrinsic function reference stays in designator form so
          // that lowering can associate the dummy argument with the named
          // constant's storage (in particular for element sequence
          // association); folding it here would replace it with its scalar
          // value.  Intrinsic references still fold, since intrinsic folding
          // inspects constant values.
          continue;
        }
        *expr = Fold(context, std::move(*expr));
      }
    }
  }
  if (intrinsic) {
    // A conditional arg that survives the argument folding above has a
    // condition that isn't a compile-time constant, or belongs to KIND, whose
    // arguments aren't folded. Either way, skip intrinsic folding unless every
    // such argument corresponds to a type-only inquiry dummy. F2023 C1538
    // (declared type and kind type parameters) and C1539 (rank) make every
    // consequent-arg agree on what such an inquiry examines, so the first
    // consequent-arg stands in for the whole argument in a copy of the
    // reference.
    //
    // TODO:
    // For elemental/pure intrinsics, distribute the call over each
    // consequent of the conditional arg and fold each branch independently:
    //   abs((c1 ? a : c2 ? b : c))
    //     → (c1 ? abs(a) : c2 ? abs(b) : abs(c))
    // Use ForEachConsequent to walk the chain, clone the call per
    // consequent, fold each clone, and reassemble into a new ConditionalArg.
    // When multiple arguments are conditional args, distribute one at a
    // time to avoid a combinatorial cross-product expansion.
    // This is not valid for non-elemental intrinsics such as RESHAPE or
    // TRANSFER, whose results depend on seeing all arguments together.
    std::optional<FunctionRef<T>> viaRepresentative;
    for (std::size_t j{0}; j < args.size(); ++j) {
      if (args[j] && args[j]->isConditionalArg()) {
        const ActualArgument::ConditionalArg &condArg{
            DEREF(args[j]->GetConditionalArg())};
        // C1540 rules out .NIL. here in conforming code: no type-only inquiry
        // dummy is OPTIONAL in the standard (the IEEE_SUPPORT_* inquiries are
        // generic with and without X). Bail out defensively.
        if (condArg.HasNilConsequent() ||
            !IsTypeOnlyInquiryArg(*intrinsic, j)) {
          return Expr<T>{std::move(funcRef)};
        }
        if (!viaRepresentative) {
          viaRepresentative.emplace(funcRef);
        }
        // Assigning the Expr keeps the argument's keyword and dummy intent.
        *viaRepresentative->arguments()[j] =
            common::Clone(DEREF(condArg.FirstNonNilConsequent()));
      }
    }
    if (viaRepresentative) {
      Expr<T> folded{FoldOperation(context, std::move(*viaRepresentative))};
      // Anything but a constant still refers to the representative alone, so
      // keep the conditional reference. RANK of consequent-args that are all
      // assumed-rank, which C1539 allows, folds to a DescriptorInquiry.
      return std::holds_alternative<Constant<T>>(folded.u)
          ? std::move(folded)
          : Expr<T>{std::move(funcRef)};
    }
    const std::string name{intrinsic->name};
    if (name == "cshift") {
      return Folder<T>{kind, context}.CSHIFT(std::move(funcRef));
    } else if (name == "eoshift") {
      return Folder<T>{kind, context}.EOSHIFT(std::move(funcRef));
    } else if (name == "merge") {
      return Folder<T>{kind, context}.MERGE(std::move(funcRef));
    } else if (name == "pack") {
      return Folder<T>{kind, context}.PACK(std::move(funcRef));
    } else if (name == "reshape") {
      return Folder<T>{kind, context}.RESHAPE(std::move(funcRef));
    } else if (name == "spread") {
      return Folder<T>{kind, context}.SPREAD(std::move(funcRef));
    } else if (name == "transfer") {
      return Folder<T>{kind, context}.TRANSFER(std::move(funcRef));
    } else if (name == "transpose") {
      return Folder<T>{kind, context}.TRANSPOSE(std::move(funcRef));
    } else if (name == "unpack") {
      return Folder<T>{kind, context}.UNPACK(std::move(funcRef));
    }
    // TODO: extends_type_of, same_type_as
    if constexpr (std::is_same_v<T, Type<TypeCategory::Derived>>) {
      // Fold enumeration type intrinsics: HUGE(enum), NEXT(enum),
      // PREVIOUS(enum)
      if (name == "huge") {
        // HUGE was eagerly folded — the first arg is the constant result
        if (args.size() >= 1 && args[0]) {
          if (auto *expr{
                  UnwrapExpr<Expr<Type<TypeCategory::Derived>>>(args[0])}) {
            return std::move(*expr);
          }
        }
      } else if (name == "next" || name == "previous") {
        return FoldEnumerationNextOrPrevious(
            context, std::move(funcRef), name == "next");
      }
    } else {
      return FoldIntrinsicFunction(context, std::move(funcRef));
    }
  }
  return Expr<T>{std::move(funcRef)};
}
 
// Array constructor folding
template <typename T> class ArrayConstructorFolder {
public:
  explicit ArrayConstructorFolder(int kind, FoldingContext &c)
      : context_{c}, resultInfo_{kind} {}
 
  Expr<T> FoldArray(ArrayConstructor<T> &&array) {
    const int kind{array.kind()};
    if constexpr (T::category == TypeCategory::Character) {
      if (const auto *len{array.LEN()}) {
        charLength_ = ToInt64(Fold(context_, common::Clone(*len)));
        knownCharLength_ = charLength_.has_value();
      }
    }
    // Calls FoldArray(const ArrayConstructorValues<T> &) below
    if (FoldArray(array)) {
      auto n{static_cast<ConstantSubscript>(elements_.size())};
      if constexpr (std::is_same_v<T, Type<TypeCategory::Derived>>) {
        CHECK(kind == 0);
        return Expr<T>{Constant<T>{array.GetType().GetDerivedTypeSpec(),
            std::move(elements_), ConstantSubscripts{n}}};
      } else if constexpr (T::category == TypeCategory::Character) {
        if (charLength_) {
          return Expr<T>{Constant<T>{
              kind, *charLength_, std::move(elements_), ConstantSubscripts{n}}};
        }
      } else {
        CHECK(resultInfo_.kind() == kind);
        return Expr<T>{Constant<T>{
            kind, std::move(elements_), ConstantSubscripts{n}, resultInfo_}};
      }
    }
    return Expr<T>{std::move(array)};
  }
 
private:
  bool FoldArray(const Expr<T> &expr) {
    const int kind{expr.kind()};
    Expr<T> folded{Fold(context_, common::Clone(expr))};
    if (const auto *c{UnwrapConstantValue<T>(kind, folded)}) {
      // Copy elements in Fortran array element order
      if (!c->empty()) {
        ConstantSubscripts index{c->lbounds()};
        do {
          elements_.emplace_back(c->At(index));
        } while (c->IncrementSubscripts(index));
      }
      if constexpr (T::category == TypeCategory::Character) {
        if (!knownCharLength_) {
          charLength_ = std::max(c->LEN(), charLength_.value_or(-1));
        }
      } else if constexpr (T::category == TypeCategory::Real ||
          T::category == TypeCategory::Complex) {
        if (c->result().isFromInexactLiteralConversion()) {
          resultInfo_.set_isFromInexactLiteralConversion();
        }
      }
      return true;
    } else {
      return false;
    }
  }
  bool FoldArray(const common::CopyableIndirection<Expr<T>> &expr) {
    return FoldArray(expr.value());
  }
  bool FoldArray(const ImpliedDo<T> &iDo) {
    Expr<SubscriptInteger> lower{
        Fold(context_, Expr<SubscriptInteger>{iDo.lower()})};
    Expr<SubscriptInteger> upper{
        Fold(context_, Expr<SubscriptInteger>{iDo.upper()})};
    Expr<SubscriptInteger> stride{
        Fold(context_, Expr<SubscriptInteger>{iDo.stride()})};
    std::optional<ConstantSubscript> start{ToInt64(lower)}, end{ToInt64(upper)},
        step{ToInt64(stride)};
    if (start && end && step && *step != 0) {
      bool result{true};
      ConstantSubscript &j{context_.StartImpliedDo(iDo.name(), *start)};
      if (*step > 0) {
        for (; j <= *end; j += *step) {
          result &= FoldArray(iDo.values());
        }
      } else {
        for (; j >= *end; j += *step) {
          result &= FoldArray(iDo.values());
        }
      }
      context_.EndImpliedDo(iDo.name());
      return result;
    } else {
      return false;
    }
  }
  bool FoldArray(const ArrayConstructorValue<T> &x) {
    return common::visit([&](const auto &y) { return FoldArray(y); }, x.u);
  }
  bool FoldArray(const ArrayConstructorValues<T> &xs) {
    for (const auto &x : xs) {
      if (!FoldArray(x)) {
        return false;
      }
    }
    return true;
  }
 
  FoldingContext &context_;
  std::vector<Scalar<T>> elements_;
  std::optional<ConstantSubscript> charLength_;
  bool knownCharLength_{false};
  typename Constant<T>::Result resultInfo_;
};
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, ArrayConstructor<T> &&array) {
  const int kind{array.kind()};
  return ArrayConstructorFolder<T>{kind, context}.FoldArray(std::move(array));
}
 
// Array operation elemental application: When all operands to an operation
// are constant arrays, array constructors without any implied DO loops,
// &/or expanded scalars, pull the operation "into" the array result by
// applying it in an elementwise fashion.  For example, [A,1]+[B,2]
// is rewritten into [A+B,1+2] and then partially folded to [A+B,3].
 
// If possible, restructures an array expression into an array constructor
// that comprises a "flat" ArrayConstructorValues with no implied DO loops.
template <typename T>
bool ArrayConstructorIsFlat(const ArrayConstructorValues<T> &values) {
  for (const ArrayConstructorValue<T> &x : values) {
    if (!std::holds_alternative<Expr<T>>(x.u)) {
      return false;
    }
  }
  return true;
}
 
template <typename T>
std::optional<Expr<T>> AsFlatArrayConstructor(const Expr<T> &expr) {
  const int kind{expr.kind()};
  if (const auto *c{UnwrapConstantValue<T>(kind, expr)}) {
    ArrayConstructor<T> result{kind, expr};
    if (!c->empty()) {
      ConstantSubscripts at{c->lbounds()};
      do {
        result.Push(MakeConstantExpr<T>(kind, c->At(at)));
      } while (c->IncrementSubscripts(at));
    }
    return std::make_optional<Expr<T>>(std::move(result));
  } else if (const auto *a{UnwrapExpr<ArrayConstructor<T>>(expr)}) {
    if (ArrayConstructorIsFlat(*a)) {
      return std::make_optional<Expr<T>>(expr);
    }
  } else if (const auto *p{UnwrapExpr<Parentheses<T>>(expr)}) {
    return AsFlatArrayConstructor(Expr<T>{p->left()});
  }
  return std::nullopt;
}
 
// FromArrayConstructor is a subroutine for MapOperation() below.
// Given a flat ArrayConstructor<T> and a shape, it wraps the array
// into an Expr<T>, folds it, and returns the resulting wrapped
// array constructor or constant array value.
template <typename T>
std::optional<Expr<T>> FromArrayConstructor(
    FoldingContext &context, ArrayConstructor<T> &&values, const Shape &shape) {
  const int kind{values.kind()};
  if (auto constShape{AsConstantExtents(context, shape)};
      constShape && !HasNegativeExtent(*constShape)) {
    Expr<T> result{Fold(context, Expr<T>{std::move(values)})};
    if (auto *constant{UnwrapConstantValue<T>(kind, result)}) {
      // Elements and shape are both constant.
      return Expr<T>{constant->Reshape(std::move(*constShape))};
    }
    if (constShape->size() == 1) {
      if (auto elements{GetShape(context, result)}) {
        if (auto constElements{AsConstantExtents(context, *elements)}) {
          if (constElements->size() == 1 &&
              constElements->at(0) == constShape->at(0)) {
            // Elements are not constant, but array constructor has
            // the right known shape and can be simply returned as is.
            return std::move(result);
          }
        }
      }
    }
  }
  return std::nullopt;
}
 
// MapOperation is a utility for various specializations of ApplyElementwise()
// that follow.  Given one or two flat ArrayConstructor<OPERAND> (wrapped in an
// Expr<OPERAND>) for some specific operand type(s), apply a given function f
// to each of their corresponding elements to produce a flat
// ArrayConstructor<RESULT> (wrapped in an Expr<RESULT>).
// Preserves shape.
 
// Unary case
template <typename RESULT, typename OPERAND>
std::optional<Expr<RESULT>> MapOperation(int kind, FoldingContext &context,
    std::function<Expr<RESULT>(Expr<OPERAND> &&)> &&f, const Shape &shape,
    [[maybe_unused]] std::optional<Expr<SubscriptInteger>> &&length,
    Expr<OPERAND> &&values) {
  ArrayConstructor<RESULT> result{kind, values};
  auto &aConst{std::get<ArrayConstructor<OPERAND>>(values.u)};
  for (auto &acValue : aConst) {
    auto &scalar{std::get<Expr<OPERAND>>(acValue.u)};
    result.Push(Fold(context, f(std::move(scalar))));
  }
  if constexpr (RESULT::category == TypeCategory::Character) {
    if (length) {
      result.set_LEN(std::move(*length));
    }
  }
  return FromArrayConstructor(context, std::move(result), shape);
}
 
template <typename RESULT, typename A>
ArrayConstructor<RESULT> ArrayConstructorFromMold(int kind, const A &prototype,
    std::optional<Expr<SubscriptInteger>> &&length) {
  ArrayConstructor<RESULT> result{kind, prototype};
  if constexpr (RESULT::category == TypeCategory::Character) {
    if (length) {
      result.set_LEN(std::move(*length));
    }
  }
  return result;
}
 
template <typename LEFT, typename RIGHT>
bool ShapesMatch(FoldingContext &context,
    const ArrayConstructor<LEFT> &leftArrConst,
    const ArrayConstructor<RIGHT> &rightArrConst) {
  auto rightIter{rightArrConst.begin()};
  for (auto &leftValue : leftArrConst) {
    CHECK(rightIter != rightArrConst.end());
    auto &leftExpr{std::get<Expr<LEFT>>(leftValue.u)};
    auto &rightExpr{std::get<Expr<RIGHT>>(rightIter->u)};
    if (leftExpr.Rank() != rightExpr.Rank()) {
      return false;
    }
    std::optional<Shape> leftShape{GetShape(context, leftExpr)};
    std::optional<Shape> rightShape{GetShape(context, rightExpr)};
    if (!leftShape || !rightShape || *leftShape != *rightShape) {
      return false;
    }
    ++rightIter;
  }
  return true;
}
 
// array * array case
template <typename RESULT, typename LEFT, typename RIGHT>
auto MapOperation(int kind, FoldingContext &context,
    std::function<Expr<RESULT>(Expr<LEFT> &&, Expr<RIGHT> &&)> &&f,
    const Shape &shape, std::optional<Expr<SubscriptInteger>> &&length,
    Expr<LEFT> &&leftValues, Expr<RIGHT> &&rightValues)
    -> std::optional<Expr<RESULT>> {
  auto result{
      ArrayConstructorFromMold<RESULT>(kind, leftValues, std::move(length))};
  auto &leftArrConst{std::get<ArrayConstructor<LEFT>>(leftValues.u)};
  auto &rightArrConst{std::get<ArrayConstructor<RIGHT>>(rightValues.u)};
  if (!ShapesMatch(context, leftArrConst, rightArrConst)) {
    return std::nullopt;
  }
  auto rightIter{rightArrConst.begin()};
  for (auto &leftValue : leftArrConst) {
    CHECK(rightIter != rightArrConst.end());
    auto &leftScalar{std::get<Expr<LEFT>>(leftValue.u)};
    auto &rightScalar{std::get<Expr<RIGHT>>(rightIter->u)};
    result.Push(
        Fold(context, f(std::move(leftScalar), std::move(rightScalar))));
    ++rightIter;
  }
  return FromArrayConstructor(context, std::move(result), shape);
}
 
// array * scalar case
template <typename RESULT, typename LEFT, typename RIGHT>
auto MapOperation(int kind, FoldingContext &context,
    std::function<Expr<RESULT>(Expr<LEFT> &&, Expr<RIGHT> &&)> &&f,
    const Shape &shape, std::optional<Expr<SubscriptInteger>> &&length,
    Expr<LEFT> &&leftValues, const Expr<RIGHT> &rightScalar)
    -> std::optional<Expr<RESULT>> {
  auto result{
      ArrayConstructorFromMold<RESULT>(kind, leftValues, std::move(length))};
  auto &leftArrConst{std::get<ArrayConstructor<LEFT>>(leftValues.u)};
  for (auto &leftValue : leftArrConst) {
    auto &leftScalar{std::get<Expr<LEFT>>(leftValue.u)};
    result.Push(
        Fold(context, f(std::move(leftScalar), Expr<RIGHT>{rightScalar})));
  }
  return FromArrayConstructor(context, std::move(result), shape);
}
 
// scalar * array case
template <typename RESULT, typename LEFT, typename RIGHT>
auto MapOperation(int kind, FoldingContext &context,
    std::function<Expr<RESULT>(Expr<LEFT> &&, Expr<RIGHT> &&)> &&f,
    const Shape &shape, std::optional<Expr<SubscriptInteger>> &&length,
    const Expr<LEFT> &leftScalar, Expr<RIGHT> &&rightValues)
    -> std::optional<Expr<RESULT>> {
  auto result{
      ArrayConstructorFromMold<RESULT>(kind, leftScalar, std::move(length))};
  auto &rightArrConst{std::get<ArrayConstructor<RIGHT>>(rightValues.u)};
  for (auto &rightValue : rightArrConst) {
    auto &rightScalar{std::get<Expr<RIGHT>>(rightValue.u)};
    result.Push(
        Fold(context, f(Expr<LEFT>{leftScalar}, std::move(rightScalar))));
  }
  return FromArrayConstructor(context, std::move(result), shape);
}
 
template <typename DERIVED, typename RESULT, typename... OPD>
std::optional<Expr<SubscriptInteger>> ComputeResultLength(
    Operation<DERIVED, RESULT, OPD...> &operation) {
  if constexpr (RESULT::category == TypeCategory::Character) {
    return Expr<RESULT>{operation.derived()}.LEN();
  }
  return std::nullopt;
}
 
// ApplyElementwise() recursively folds the operand expression(s) of an
// operation, then attempts to apply the operation to the (corresponding)
// scalar element(s) of those operands.  Returns std::nullopt for scalars
// or unlinearizable operands.
template <typename DERIVED, typename RESULT, typename OPERAND>
auto ApplyElementwise(FoldingContext &context,
    Operation<DERIVED, RESULT, OPERAND> &operation,
    std::function<Expr<RESULT>(Expr<OPERAND> &&)> &&f)
    -> std::optional<Expr<RESULT>> {
  const int kind{operation.kind()};
  auto &expr{operation.left()};
  expr = Fold(context, std::move(expr));
  if (expr.Rank() > 0) {
    if (std::optional<Shape> shape{GetShape(context, expr)}) {
      if (auto values{AsFlatArrayConstructor(expr)}) {
        return MapOperation(kind, context, std::move(f), *shape,
            ComputeResultLength(operation), std::move(*values));
      }
    }
  }
  return std::nullopt;
}
 
template <typename DERIVED, typename RESULT, typename OPERAND>
auto ApplyElementwise(
    FoldingContext &context, Operation<DERIVED, RESULT, OPERAND> &operation)
    -> std::optional<Expr<RESULT>> {
  const int kind{operation.kind()};
  return ApplyElementwise(context, operation,
      std::function<Expr<RESULT>(Expr<OPERAND> &&)>{
          [kind](Expr<OPERAND> &&operand) {
            return Expr<RESULT>{DERIVED{kind, std::move(operand)}};
          }});
}
 
template <typename DERIVED, typename RESULT, typename LEFT, typename RIGHT>
auto ApplyElementwise(FoldingContext &context,
    Operation<DERIVED, RESULT, LEFT, RIGHT> &operation,
    std::function<Expr<RESULT>(Expr<LEFT> &&, Expr<RIGHT> &&)> &&f)
    -> std::optional<Expr<RESULT>> {
  const int kind{operation.kind()};
  auto resultLength{ComputeResultLength(operation)};
  auto &leftExpr{operation.left()};
  auto &rightExpr{operation.right()};
  if (leftExpr.Rank() != rightExpr.Rank() && leftExpr.Rank() != 0 &&
      rightExpr.Rank() != 0) {
    return std::nullopt; // error recovery
  }
  leftExpr = Fold(context, std::move(leftExpr));
  rightExpr = Fold(context, std::move(rightExpr));
  if (leftExpr.Rank() > 0) {
    if (std::optional<Shape> leftShape{GetShape(context, leftExpr)}) {
      if (auto left{AsFlatArrayConstructor(leftExpr)}) {
        if (rightExpr.Rank() > 0) {
          if (std::optional<Shape> rightShape{GetShape(context, rightExpr)}) {
            if (auto right{AsFlatArrayConstructor(rightExpr)}) {
              if (CheckConformance(context.messages(), *leftShape, *rightShape,
                      CheckConformanceFlags::EitherScalarExpandable)
                      .value_or(false /*fail if not known now to conform*/)) {
                return MapOperation(kind, context, std::move(f), *leftShape,
                    std::move(resultLength), std::move(*left),
                    std::move(*right));
              } else {
                return std::nullopt;
              }
              return MapOperation(kind, context, std::move(f), *leftShape,
                  std::move(resultLength), std::move(*left), std::move(*right));
            }
          }
        } else if (IsExpandableScalar(rightExpr, context, *leftShape)) {
          return MapOperation(kind, context, std::move(f), *leftShape,
              std::move(resultLength), std::move(*left), rightExpr);
        }
      }
    }
  } else if (rightExpr.Rank() > 0) {
    if (std::optional<Shape> rightShape{GetShape(context, rightExpr)}) {
      if (IsExpandableScalar(leftExpr, context, *rightShape)) {
        if (auto right{AsFlatArrayConstructor(rightExpr)}) {
          return MapOperation(kind, context, std::move(f), *rightShape,
              std::move(resultLength), leftExpr, std::move(*right));
        }
      }
    }
  }
  return std::nullopt;
}
 
template <typename DERIVED, typename RESULT, typename LEFT, typename RIGHT>
auto ApplyElementwise(
    FoldingContext &context, Operation<DERIVED, RESULT, LEFT, RIGHT> &operation)
    -> std::optional<Expr<RESULT>> {
  const int kind{operation.kind()};
  return ApplyElementwise(context, operation,
      std::function<Expr<RESULT>(Expr<LEFT> &&, Expr<RIGHT> &&)>{
          [kind](Expr<LEFT> &&left, Expr<RIGHT> &&right) {
            return Expr<RESULT>{
                DERIVED{kind, std::move(left), std::move(right)}};
          }});
}
 
// Unary operations
 
template <typename TO, typename FROM>
common::IfNoLvalue<std::optional<TO>, FROM> ConvertString(
    int toKind, FROM &&s) {
  const int fromKind{s.kind()};
  if (std::is_same_v<TO, FROM> && toKind == fromKind) {
    return std::make_optional<TO>(std::move(s));
  } else {
    return s.ToAscii(toKind);
  }
}
 
template <typename TO, TypeCategory FROMCAT>
Expr<TO> FoldOperation(
    FoldingContext &context, Convert<TO, FROMCAT> &&convert) {
  const int toKind{convert.kind()};
  if (auto array{ApplyElementwise(context, convert)}) {
    return *array;
  }
  using Operand = Type<FROMCAT>;
  auto &operand{convert.left()};
  const int operandKind{operand.kind()};
  // This variable is a workaround for msvc which emits an error when
  // using the FROMCAT template parameter below.
  TypeCategory constexpr FromCat{FROMCAT};
  if (auto value{GetScalarConstantValue<Operand>(operand)}) {
    FoldingContext &ctx{context};
    const int fromKind{value->kind()};
    if constexpr (TO::category == TypeCategory::Integer) {
      if constexpr (FromCat == TypeCategory::Integer) {
        auto converted{Scalar<TO>::ConvertSigned(toKind, *value)};
        if (converted.overflow) {
          ctx.Warn(common::UsageWarning::FoldingException,
              "conversion of %s_%d to INTEGER(%d) overflowed; result is %s"_warn_en_US,
              value->SignedDecimal(), fromKind, toKind,
              converted.value.SignedDecimal());
        }
        return MakeConstantExpr<TO>(toKind, std::move(converted.value));
      } else if constexpr (FromCat == TypeCategory::Unsigned) {
        auto converted{Scalar<TO>::ConvertUnsigned(toKind, *value)};
        if ((converted.overflow || converted.value.IsNegative())) {
          ctx.Warn(common::UsageWarning::FoldingException,
              "conversion of %s_U%d to INTEGER(%d) overflowed; result is %s"_warn_en_US,
              value->UnsignedDecimal(), fromKind, toKind,
              converted.value.SignedDecimal());
        }
        return MakeConstantExpr<TO>(toKind, std::move(converted.value));
      } else if constexpr (FromCat == TypeCategory::Real) {
        auto converted{value->ToInteger(
            common::RoundingMode::ToZero, Scalar<TO>::bits(toKind))};
        if (converted.flags.test(RealFlag::InvalidArgument)) {
          ctx.Warn(common::UsageWarning::FoldingException,
              "REAL(%d) to INTEGER(%d) conversion: invalid argument"_warn_en_US,
              fromKind, toKind);
        } else if (converted.flags.test(RealFlag::Overflow)) {
          ctx.Warn(common::UsageWarning::FoldingException,
              "REAL(%d) to INTEGER(%d) conversion overflowed"_warn_en_US,
              fromKind, toKind);
        }
        return MakeConstantExpr<TO>(toKind, std::move(converted.value));
      }
    } else if constexpr (TO::category == TypeCategory::Unsigned) {
      if constexpr (FromCat == TypeCategory::Integer ||
          FromCat == TypeCategory::Unsigned) {
        return MakeConstantExpr<TO>(
            toKind, Scalar<TO>::ConvertUnsigned(toKind, *value).value);
      } else if constexpr (FromCat == TypeCategory::Real) {
        return MakeConstantExpr<TO>(toKind,
            value
                ->ToInteger(
                    common::RoundingMode::ToZero, Scalar<TO>::bits(toKind))
                .value);
      }
    } else if constexpr (TO::category == TypeCategory::Real) {
      if constexpr (FromCat == TypeCategory::Integer ||
          FromCat == TypeCategory::Unsigned) {
        auto converted{Scalar<TO>::FromInteger(
            toKind, *value, FromCat == TypeCategory::Unsigned)};
        if (!converted.flags.empty()) {
          char buffer[64];
          std::snprintf(buffer, sizeof buffer,
              "INTEGER(%d) to REAL(%d) conversion", fromKind, toKind);
          ctx.RealFlagWarnings(converted.flags, buffer);
        }
        return MakeConstantExpr<TO>(toKind, std::move(converted.value));
      } else if constexpr (FromCat == TypeCategory::Real) {
        auto converted{Scalar<TO>::Convert(toKind, *value)};
        char buffer[64];
        if (!converted.flags.empty()) {
          std::snprintf(buffer, sizeof buffer,
              "REAL(%d) to REAL(%d) conversion", fromKind, toKind);
          ctx.RealFlagWarnings(converted.flags, buffer);
        }
        if (ctx.targetCharacteristics().areSubnormalsFlushedToZero()) {
          converted.value = converted.value.FlushSubnormalToZero();
        }
        return MakeConstantExpr<TO>(toKind, std::move(converted.value));
      }
    } else if constexpr (TO::category == TypeCategory::Complex) {
      if constexpr (FromCat == TypeCategory::Complex) {
        auto real{value->REAL()};
        auto imag{value->AIMAG()};
        return FoldOperation(ctx,
            ComplexConstructor{toKind,
                Fold(ctx,
                    AsExpr(Convert<typename TO::Part>(toKind,
                        AsCategoryExpr(Constant<typename Operand::Part>{
                            real.kind(), real})))),
                Fold(ctx,
                    AsExpr(Convert<typename TO::Part>(toKind,
                        AsCategoryExpr(Constant<typename Operand::Part>{
                            imag.kind(), imag}))))});
      }
    } else if constexpr (TO::category == TypeCategory::Character &&
        FromCat == TypeCategory::Character) {
      if (auto converted{
              ConvertString<Scalar<TO>>(toKind, std::move(*value))}) {
        return ScalarConstantToExpr(std::move(*converted));
      }
    } else if constexpr (TO::category == TypeCategory::Logical &&
        FromCat == TypeCategory::Logical) {
      return MakeConstantExpr<TO>(toKind, value->IsTrue());
    }
  } else if constexpr (TO::category == FromCat &&
      FromCat != TypeCategory::Character) {
    // Conversion of non-constant in same type category
    if (Operand::category == TO::category && toKind == operandKind) {
      return std::move(operand); // remove needless conversion
    } else if constexpr (TO::category == TypeCategory::Logical ||
        TO::category == TypeCategory::Integer) {
      if (auto *innerConv{
              std::get_if<Convert<Operand, TO::category>>(&operand.u)}) {
        // Conversion of conversion of same category & kind
        if (auto *x{&innerConv->left()}; x->kind() == toKind) {
          if (TO::category == TypeCategory::Logical || toKind <= operandKind) {
            return std::move(*x); // no-op Logical or Integer
                                  // widening/narrowing conversion pair
          } else if constexpr (std::is_same_v<TO, DescriptorInquiry::Result>) {
            if (toKind == DescriptorInquiry::kind()) {
              if (std::holds_alternative<DescriptorInquiry>(x->u) ||
                  std::holds_alternative<TypeParamInquiry>(x->u)) {
                // int(int(size(...),kind=k),kind=8) -> size(...)
                return std::move(*x);
              }
            }
          }
        }
      }
    }
  }
  return Expr<TO>{std::move(convert)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Parentheses<T> &&x) {
  const int kind{x.kind()};
  auto &operand{x.left()};
  operand = Fold(context, std::move(operand));
  if (auto value{GetScalarConstantValue<T>(kind, operand)}) {
    // Preserve parentheses, even around constants.
    return Expr<T>{Parentheses<T>{kind, MakeConstantExpr<T>(kind, *value)}};
  } else if (std::holds_alternative<Parentheses<T>>(operand.u)) {
    // ((x)) -> (x)
    return std::move(operand);
  } else {
    return Expr<T>{Parentheses<T>{std::move(operand)}};
  }
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Negate<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  auto &operand{x.left()};
  if (auto *nn{std::get_if<Negate<T>>(&x.left().u)}; nn && nn->kind() == kind) {
    // -(-x) -> (x)
    if (IsVariable(nn->left())) {
      return FoldOperation(context, Parentheses<T>{std::move(nn->left())});
    } else {
      return std::move(nn->left());
    }
  } else if (auto value{GetScalarConstantValue<T>(kind, operand)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      auto negated{value->Negate()};
      if (negated.overflow) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) negation overflowed"_warn_en_US, kind);
      }
      return Expr<T>{Constant<T>{kind, std::move(negated.value)}};
    } else if constexpr (T::category == TypeCategory::Unsigned) {
      return Expr<T>{Constant<T>{kind, std::move(value->Negate().value)}};
    } else {
      // REAL & COMPLEX negation: no exceptions possible
      return Expr<T>{Constant<T>{kind, value->Negate()}};
    }
  }
  return Expr<T>{std::move(x)};
}
 
// Binary (dyadic) operations
 
template <typename LEFT, typename RIGHT>
std::optional<std::pair<Scalar<LEFT>, Scalar<RIGHT>>> OperandsAreConstants(
    const Expr<LEFT> &x, const Expr<RIGHT> &y) {
  if (auto xvalue{GetScalarConstantValue<LEFT>(x)}) {
    if (auto yvalue{GetScalarConstantValue<RIGHT>(y)}) {
      return {std::make_pair(*xvalue, *yvalue)};
    }
  }
  return std::nullopt;
}
 
template <typename DERIVED, typename RESULT, typename LEFT, typename RIGHT>
std::optional<std::pair<Scalar<LEFT>, Scalar<RIGHT>>> OperandsAreConstants(
    const Operation<DERIVED, RESULT, LEFT, RIGHT> &operation) {
  return OperandsAreConstants(operation.left(), operation.right());
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Add<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      auto sum{folded->first.AddSigned(folded->second)};
      if (sum.overflow) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) addition overflowed"_warn_en_US, kind);
      }
      return Expr<T>{Constant<T>{kind, sum.value}};
    } else if constexpr (T::category == TypeCategory::Unsigned) {
      return Expr<T>{
          Constant<T>{kind, folded->first.AddUnsigned(folded->second).value}};
    } else {
      auto sum{folded->first.Add(
          folded->second, context.targetCharacteristics().roundingMode())};
      context.RealFlagWarnings(sum.flags, "addition");
      if (context.targetCharacteristics().areSubnormalsFlushedToZero()) {
        sum.value = sum.value.FlushSubnormalToZero();
      }
      return Expr<T>{Constant<T>{kind, sum.value}};
    }
  } else if constexpr (T::category == TypeCategory::Integer ||
      T::category == TypeCategory::Unsigned) {
    if (auto c{GetScalarConstantValue<T>(kind, x.right())}) {
      if (c->IsZero() && x.left().Rank() == 0) {
        if (IsVariable(x.left())) {
          return FoldOperation(context, Parentheses<T>{std::move(x.left())});
        } else {
          return std::move(x.left());
        }
      }
    } else if (auto c{GetScalarConstantValue<T>(kind, x.left())}) {
      if (c->IsZero() && x.right().Rank() == 0) {
        if (IsVariable(x.right())) {
          return FoldOperation(context, Parentheses<T>{std::move(x.right())});
        } else {
          return std::move(x.right());
        }
      }
    }
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Subtract<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      auto difference{folded->first.SubtractSigned(folded->second)};
      if (difference.overflow) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) subtraction overflowed"_warn_en_US, kind);
      }
      return Expr<T>{Constant<T>{kind, difference.value}};
    } else if constexpr (T::category == TypeCategory::Unsigned) {
      return Expr<T>{Constant<T>{
          kind, folded->first.SubtractSigned(folded->second).value}};
    } else {
      auto difference{folded->first.Subtract(
          folded->second, context.targetCharacteristics().roundingMode())};
      context.RealFlagWarnings(difference.flags, "subtraction");
      if (context.targetCharacteristics().areSubnormalsFlushedToZero()) {
        difference.value = difference.value.FlushSubnormalToZero();
      }
      return Expr<T>{Constant<T>{kind, difference.value}};
    }
  } else if constexpr (T::category == TypeCategory::Integer ||
      T::category == TypeCategory::Unsigned) {
    if (auto c{GetScalarConstantValue<T>(kind, x.right())}) {
      if (c->IsZero() && x.left().Rank() == 0) {
        if (IsVariable(x.left())) {
          return FoldOperation(context, Parentheses<T>{std::move(x.left())});
        } else {
          return std::move(x.left());
        }
      }
    }
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Multiply<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      auto product{folded->first.MultiplySigned(folded->second)};
      if (product.SignedMultiplicationOverflowed()) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) multiplication overflowed"_warn_en_US, kind);
      }
      return Expr<T>{Constant<T>{kind, product.lower}};
    } else if constexpr (T::category == TypeCategory::Unsigned) {
      return Expr<T>{Constant<T>{
          kind, folded->first.MultiplyUnsigned(folded->second).lower}};
    } else {
      auto product{folded->first.Multiply(
          folded->second, context.targetCharacteristics().roundingMode())};
      context.RealFlagWarnings(product.flags, "multiplication");
      if (context.targetCharacteristics().areSubnormalsFlushedToZero()) {
        product.value = product.value.FlushSubnormalToZero();
      }
      return Expr<T>{Constant<T>{kind, product.value}};
    }
  } else if constexpr (T::category == TypeCategory::Integer) {
    if (auto c{GetScalarConstantValue<T>(kind, x.right())}) {
      x.right() = std::move(x.left());
      x.left() = MakeConstantExpr<T>(kind, *c);
    }
    if (auto c{GetScalarConstantValue<T>(kind, x.left())}) {
      if (c->IsZero() && x.right().Rank() == 0) {
        return std::move(x.left());
      } else if (c->CompareSigned(value::IntegerValue{kind, 1}) ==
          Ordering::Equal) {
        if (IsVariable(x.right())) {
          return FoldOperation(context, Parentheses<T>{std::move(x.right())});
        } else {
          return std::move(x.right());
        }
      } else if (c->CompareSigned(value::IntegerValue{kind, -1}) ==
          Ordering::Equal) {
        return FoldOperation(context, Negate<T>{std::move(x.right())});
      }
    }
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Divide<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      auto quotAndRem{folded->first.DivideSigned(folded->second)};
      if (quotAndRem.divisionByZero) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) division by zero"_warn_en_US, kind);
        return Expr<T>{std::move(x)};
      }
      if (quotAndRem.overflow) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) division overflowed"_warn_en_US, kind);
      }
      return Expr<T>{Constant<T>{kind, quotAndRem.quotient}};
    } else if constexpr (T::category == TypeCategory::Unsigned) {
      auto quotAndRem{folded->first.DivideUnsigned(folded->second)};
      if (quotAndRem.divisionByZero) {
        context.Warn(common::UsageWarning::FoldingException,
            "UNSIGNED(%d) division by zero"_warn_en_US, kind);
        return Expr<T>{std::move(x)};
      }
      return Expr<T>{Constant<T>{kind, quotAndRem.quotient}};
    } else {
      auto quotient{folded->first.Divide(
          folded->second, context.targetCharacteristics().roundingMode())};
      // Don't warn about -1./0., 0./0., or 1./0. from a module file
      // they are interpreted as canonical Fortran representations of -Inf,
      // NaN, and Inf respectively.
      bool isCanonicalNaNOrInf{false};
      if constexpr (T::category == TypeCategory::Real) {
        if (folded->second.IsZero() && context.moduleFileName().has_value()) {
          using IntType = typename T::Scalar::Word;
          auto intNumerator{folded->first.ToInteger()};
          const int intKind{intNumerator.value.kind()};
          isCanonicalNaNOrInf = intNumerator.flags == RealFlags{} &&
              intNumerator.value >= IntType{intKind, -1} &&
              intNumerator.value <= IntType{intKind, 1};
        }
      }
      if (!isCanonicalNaNOrInf) {
        context.RealFlagWarnings(quotient.flags, "division");
      }
      if (context.targetCharacteristics().areSubnormalsFlushedToZero()) {
        quotient.value = quotient.value.FlushSubnormalToZero();
      }
      return Expr<T>{Constant<T>{kind, quotient.value}};
    }
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Power<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      auto power{folded->first.Power(folded->second)};
      if (power.divisionByZero) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) zero to negative power"_warn_en_US, kind);
      } else if (power.overflow) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) power overflowed"_warn_en_US, kind);
      } else if (power.zeroToZero) {
        context.Warn(common::UsageWarning::FoldingException,
            "INTEGER(%d) 0**0 is not defined"_warn_en_US, kind);
      }
      return Expr<T>{Constant<T>{kind, power.power}};
    } else {
      if (folded->first.IsZero()) {
        if (folded->second.IsZero()) {
          context.Warn(common::UsageWarning::FoldingException,
              "REAL/COMPLEX 0**0 is not defined"_warn_en_US);
        } else {
          return Expr<T>(
              Constant<T>{kind, folded->first}); // 0. ** nonzero -> 0.
        }
      } else if (auto callable{GetHostRuntimeWrapper<T, T, T>(
                     kind, {kind, kind}, "pow")}) {
        return Expr<T>{Constant<T>{
            kind, (*callable)(context, folded->first, folded->second)}};
      } else {
        context.Warn(common::UsageWarning::FoldingFailure,
            "Power for %s cannot be folded on host"_warn_en_US,
            DynamicType{T::category, kind}.AsFortran());
      }
    }
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, RealToIntPower<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x)}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x.left(), x.right())}) {
    ValueWithRealFlags<Scalar<T>> power;
    if constexpr (T::category == TypeCategory::Complex) {
      if (kind == 4) {
        // Match flang-rt's cpowi/cpowk (complex-powi.cpp): accumulate in
        // COMPLEX(8) and round to single once, rather than rounding every
        // intermediate product (which costs a few single-precision ULPs).
        using DoubleComplex = Scalar<Type<TypeCategory::Complex>>;
        using DoublePart = DoubleComplex::Part;
        using SinglePart = typename Scalar<T>::Part;
        // Single -> double is exact, so these conversions raise no flags.
        DoubleComplex wideBase{
            DoublePart::Convert(8, folded->first.REAL()).value,
            DoublePart::Convert(8, folded->first.AIMAG()).value};
        auto wide{evaluate::IntPower(wideBase, folded->second)};
        power.flags = wide.flags;
        power.value = Scalar<T>{SinglePart::Convert(kind, wide.value.REAL())
                                    .AccumulateFlags(power.flags),
            SinglePart::Convert(kind, wide.value.AIMAG())
                .AccumulateFlags(power.flags)};
      } else {
        power = evaluate::IntPower(folded->first, folded->second);
      }
    } else {
      power = evaluate::IntPower(folded->first, folded->second);
    }
    context.RealFlagWarnings(power.flags, "power with INTEGER exponent");
    if (context.targetCharacteristics().areSubnormalsFlushedToZero()) {
      power.value = power.value.FlushSubnormalToZero();
    }
    return Expr<T>{Constant<T>{kind, power.value}};
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, ConditionalExpr<T> &&x) {
  x.condition() = Fold(context, std::move(x.condition()));
  // If the condition is a scalar logical constant, select the branch.
  if (auto cst{GetScalarConstantValue<LogicalResult>(
          LogicalResultKind, x.condition())}) {
    return cst->IsTrue() ? Fold(context, std::move(x.thenValue()))
                         : Fold(context, std::move(x.elseValue()));
  }
  return Expr<T>{std::move(x)};
}
 
template <typename T>
Expr<T> FoldOperation(FoldingContext &context, Extremum<T> &&x) {
  const int kind{x.kind()};
  if (auto array{ApplyElementwise(context, x,
          std::function<Expr<T>(Expr<T> &&, Expr<T> &&)>{
              [ordering{x.ordering}](Expr<T> &&l, Expr<T> &&r) {
                return Expr<T>{
                    Extremum<T>{ordering, std::move(l), std::move(r)}};
              }})}) {
    return *array;
  }
  if (auto folded{OperandsAreConstants(x)}) {
    if constexpr (T::category == TypeCategory::Integer) {
      if (folded->first.CompareSigned(folded->second) == x.ordering) {
        return Expr<T>{Constant<T>{kind, folded->first}};
      }
    } else if constexpr (T::category == TypeCategory::Unsigned) {
      if (folded->first.CompareUnsigned(folded->second) == x.ordering) {
        return Expr<T>{Constant<T>{kind, folded->first}};
      }
    } else if constexpr (T::category == TypeCategory::Real) {
      if (folded->first.IsNotANumber() ||
          (folded->first.Compare(folded->second) == Relation::Less) ==
              (x.ordering == Ordering::Less)) {
        return Expr<T>{Constant<T>{kind, folded->first}};
      }
    } else {
      static_assert(T::category == TypeCategory::Character);
      // Result of MIN and MAX on character has the length of
      // the longest argument.
      auto maxLen{std::max(folded->first.length(), folded->second.length())};
      bool isFirst{x.ordering == Compare(folded->first, folded->second)};
      auto res{isFirst ? std::move(folded->first) : std::move(folded->second)};
      res = res.length() == maxLen ? std::move(res)
                                   : CharacterUtils::Resize(res, maxLen);
      return Expr<T>{Constant<T>{kind, std::move(res)}};
    }
    return Expr<T>{Constant<T>{kind, folded->second}};
  }
  // max(c, max(c, x)) -> max(c, x), and likewise for min, in any operand
  // order, when c is a constant.  An extent max(0,ub) whose ub serializes as
  // ubound(a,dim) is re-folded on read-back into max(0,max(0,ub)); collapsing
  // it makes the module-file round trip a fixed point, so structural
  // comparisons of function results and character lengths across module files
  // still match.  Only a repeated constant is absorbed: Expr::operator== is
  // defined for comparing procedure characteristics, where two dummy arguments
  // at the same position compare equal regardless of identity, so a repeated
  // variable is not reliably the same variable, and a repeated impure function
  // reference or VOLATILE read must not be dropped.
  if constexpr (T::category == TypeCategory::Integer ||
      T::category == TypeCategory::Unsigned) {
    auto absorbs{[&](const Expr<T> &operand, const Expr<T> &nested) {
      const auto *c{std::get_if<Constant<T>>(&operand.u)};
      const auto *inner{std::get_if<Extremum<T>>(&nested.u)};
      if (!c || !inner || inner->ordering != x.ordering) {
        return false;
      }
      auto isSameConstant{[&](const Expr<T> &y) {
        const auto *ic{std::get_if<Constant<T>>(&y.u)};
        return ic && *ic == *c;
      }};
      return isSameConstant(inner->left()) || isSameConstant(inner->right());
    }};
    if (absorbs(x.left(), x.right())) {
      return Expr<T>{std::move(x.right())};
    } else if (absorbs(x.right(), x.left())) {
      return Expr<T>{std::move(x.left())};
    }
  }
  return Expr<T>{std::move(x)};
}
 
Expr<Type<TypeCategory::Real>> ToReal(
    int kind, FoldingContext &context, Expr<SomeType> &&expr);
 
template <typename T>
Expr<T> ExpressionBase<T>::Rewrite(FoldingContext &context, Expr<T> &&expr) {
  return common::visit(
      [&](auto &&x) -> Expr<T> {
        if constexpr (IsSpecificIntrinsicType<T>) {
          return FoldOperation(context, std::move(x));
        } else if constexpr (std::is_same_v<T, Type<TypeCategory::Derived>>) {
          return FoldOperation(context, std::move(x));
        } else if constexpr (common::HasMember<decltype(x),
                                 TypelessExpression>) {
          return std::move(expr);
        } else {
          return Expr<T>{Fold(context, std::move(x))};
        }
      },
      std::move(expr.u));
}
 
FOR_EACH_TYPE_AND_KIND(extern template class ExpressionBase, )
} // namespace Fortran::evaluate
#endif // FORTRAN_EVALUATE_FOLD_IMPLEMENTATION_H_