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

source code of flang/lib/Evaluate/fold-real.cpp