[flang-commits] [flang] [flang] Fold ERFC_SCALED accurately for all kinds, including REAL(16) (PR #222756)
Eugene Epshteyn via flang-commits
flang-commits at lists.llvm.org
Thu Sep 10 14:01:19 PDT 2026
https://github.com/eugeneepshteyn updated https://github.com/llvm/llvm-project/pull/222756
>From f9d9eb2d1d2e41cc23528c151fac893e96df82e6 Mon Sep 17 00:00:00 2001
From: Eugene Epshteyn <eepshteyn at nvidia.com>
Date: Thu, 10 Sep 2026 12:53:20 -0700
Subject: [PATCH 1/3] [flang] Fold ERFC_SCALED accurately for all kinds,
including REAL(16)
ERFC_SCALED could not be folded at REAL(16) on x86-64: the __float128
host-runtime table had no erfc_scaled entry, so a constant argument drew
a -Wfolding-failure warning and a REAL(16) named constant initialized
with erfc_scaled() was a hard error. Folding at the other kinds went
through the netlib-derived approximation in
flang/Common/erfc-scaled.h, which returns HUGE(x) for arguments in
(-26.6287357, -26.628) at kind 8 where the true value is representable
(up to ~3e14 ulps off), and is limited to about double precision.
This extracts the accurate algorithm contributed for the runtime in
llvm/llvm-project#219697 into policy-parameterized template functions in
a new header, flang/include/flang/Common/erfc-scaled-accurate.h, usable
by both the compiler and flang-rt, and switches host folding to it for
all kinds:
* the generic host-runtime table (float, double, long double) now folds
erfc_scaled through ErfcScaledAccurate with the <cmath> policy, and
* the __float128 table gains an erfc_scaled entry backed by a
libquadmath policy (expq/erfcq/fmaq/HUGE_VALQ).
The header generalizes the original algorithm per type and applies
corrections from the review of #219697: a per-type direct/series
threshold (16 in general, 9 for float, where expf overflows and erfcf
underflows well below the general threshold), a negative-branch cutoff
at 107 rather than 106 (true binary128 values remain representable
through x = -106.567; the cutoff itself stays load-bearing because
letting x*x overflow corrupts the fma compensation into an infinity of
the wrong sign), and a series division order that avoids flushing to
zero near the top of the format.
User-visible changes beyond the REAL(16) fix: large-negative arguments
whose true value exceeds the format's range now fold to +infinity
rather than HUGE(x) at kinds 4 and 8, and folded values elsewhere move
by up to 2 ulps toward the correctly rounded result (measured against a
300-bit mpmath 1.4.1 oracle over a 160-point sweep per kind: worst
ordinary deviation 3 ulps before, 2 ulps after; kind-8's categorical
HUGE window eliminated).
The runtime keeps its current implementations until #219697 lands;
after that the same header can serve flang-rt, keeping compile-time and
run-time results consistent.
Co-authored-by: Dmitry Mikushin <dmitry at kernelgen.org>
---
.../flang/Common/erfc-scaled-accurate.h | 148 ++++++++++++++++++
flang/lib/Evaluate/intrinsics-library.cpp | 28 +++-
flang/test/CMakeLists.txt | 9 ++
flang/test/Evaluate/fold-erfc-scaled.f90 | 50 +++++-
flang/test/Evaluate/fold-erfc-scaled16.f90 | 43 +++++
flang/test/lit.cfg.py | 4 +
flang/test/lit.site.cfg.py.in | 1 +
flang/unittests/Evaluate/folding.cpp | 83 ++++++++++
8 files changed, 358 insertions(+), 8 deletions(-)
create mode 100644 flang/include/flang/Common/erfc-scaled-accurate.h
create mode 100644 flang/test/Evaluate/fold-erfc-scaled16.f90
diff --git a/flang/include/flang/Common/erfc-scaled-accurate.h b/flang/include/flang/Common/erfc-scaled-accurate.h
new file mode 100644
index 0000000000000..4b222f26ae1ff
--- /dev/null
+++ b/flang/include/flang/Common/erfc-scaled-accurate.h
@@ -0,0 +1,148 @@
+//===-- include/flang/Common/erfc-scaled-accurate.h -------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+// Implements ERFC_SCALED (erfcx) to nearly full precision of the argument
+// type, templated over both the type and the math functions it needs, so the
+// same algorithm can serve compile-time constant folding (see
+// flang/lib/Evaluate/intrinsics-library.cpp) and the runtime (flang-rt), and
+// types with no <cmath> overloads (__float128 via libquadmath) alongside
+// float, double, and long double.
+//
+// The algorithm was contributed in
+// https://github.com/llvm/llvm-project/pull/219697 for binary128; this header
+// generalizes it per type:
+//
+// x >= threshold: the asymptotic expansion
+// erfcx(x) ~ 1/(x*sqrt(pi)) * (1 - 1/(2x^2) + 3/(4x^4) -
+// ...) truncated after kTerms terms. Host-libm-free, so
+// results are reproducible across hosts with the same IEEE
+// format, rounding mode, and no FP contraction.
+// 0 <= x < threshold: exp(x*x)*erfc(x) directly, with the squaring error
+// compensated: u = x*x, e = fma(x,x,-u) (the exact low
+// part), and exp(u+e) expanded as exp(u)*(1+e). Without
+// the compensation, exp() amplifies the squaring rounding
+// error by a factor of x*x (hundreds of ulps near the
+// threshold).
+// x < 0: the reflection erfc(-x) = 2 - erfc(x), so
+// erfcx(x) = 2*exp(x*x) - erfcx(-x), which grows and
+// overflows to +infinity where 2*exp(x*x) does.
+//
+// Two constants deserve their shared story told once:
+//
+// * The direct/series threshold is per-type: 16 in general, 9 for float.
+// The series needs x large enough for kTerms truncation to be accurate
+// (with 22 terms: ~3.1e5 binary128 eps at x = 12, <= 0.59 eps at x = 16, and
+// below one eps above; at float's 24-bit precision, x = 9 is already exact
+// to the last ulp). The direct form needs exp(x*x) finite and erfc(x)
+// normal: for binary128 that holds through the threshold with a wide
+// margin, but expf overflows at x*x > 88.72 (x ~ 9.42) and erfcf is
+// subnormal there and exactly zero by x ~ 10 - a threshold of 16 would
+// return Inf then NaN over [9.42, 16) at float. Threshold 9 keeps float on
+// the series exactly where its direct form breaks down.
+//
+// * The negative-branch cutoff 107 is one universal constant. It is sized
+// for binary128, where the true value first exceeds the format's range at
+// x*x > 16384*ln2 (|x| ~ 106.567): returning +Inf from |x| >= 107 on is
+// correct there, and the cutoff must not be larger than ~1.09e2466, where
+// x*x itself overflows and the compensation term e = fma(x,x,-Inf) becomes
+// -Inf, turning exp(u)*(1+e) into Inf*(-Inf) = -Inf - an infinity of the
+// wrong sign. For narrower types the same constant is safe without
+// adjustment: their 2*exp(x*x) overflows to +Inf - the correctly rounded
+// answer - far below 107 (x ~ 9.4 for float, ~26.63 for double), the
+// per-type threshold routes the subtrahend erfcx(-x) to the never-
+// overflowing series, and x*x stays finite below 107 in every supported
+// format, so the fma hazard is unreachable.
+//
+// Do not reach for std::numeric_limits<T> here beyond the default policy:
+// libstdc++ leaves numeric_limits<__float128> unspecialized (all members
+// zero) while libc++ and MSVC differ again, so any use of it in the shared
+// code would compile everywhere and be wrong somewhere. Everything
+// type-specific comes in through the policy.
+
+#ifndef FORTRAN_COMMON_ERFC_SCALED_ACCURATE_H_
+#define FORTRAN_COMMON_ERFC_SCALED_ACCURATE_H_
+
+#include "flang/Common/api-attrs.h"
+#include <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace Fortran::common {
+
+// Default math policy for types <cmath> covers. Consumers whose type has no
+// std:: overloads (e.g. __float128) supply their own policy with the same
+// four members instead of specializing anything.
+template <typename T> struct ErfcScaledStdHostPolicy {
+ static inline RT_API_ATTRS T Exp(T x) { return std::exp(x); }
+ static inline RT_API_ATTRS T Erfc(T x) { return std::erfc(x); }
+ static inline RT_API_ATTRS T Fma(T x, T y, T z) { return std::fma(x, y, z); }
+ static inline RT_API_ATTRS T Infinity() {
+ return std::numeric_limits<T>::infinity();
+ }
+};
+
+// sqrt(pi) as the sum of three exactly representable doubles: correctly
+// rounded up to and including binary128 with no long double literal capping
+// the precision and no Q suffix (a GNU extension). Keep the grouping; do not
+// fold the sum through a narrower type.
+template <typename T> constexpr RT_API_ATTRS T ErfcScaledSqrtPi() {
+ return static_cast<T>(1.772453850905516) +
+ static_cast<T>(-7.666586499825799e-17) +
+ static_cast<T>(-1.3058334907945429e-33);
+}
+
+// exp(x*x) with the squaring rounding error compensated (see file comment).
+template <typename T, typename P = ErfcScaledStdHostPolicy<T>>
+inline RT_API_ATTRS T ErfcScaledExpOfSquare(T x) {
+ T u{x * x};
+ T e{P::Fma(x, x, -u)};
+ return P::Exp(u) * (T{1} + e);
+}
+
+template <typename T, typename P = ErfcScaledStdHostPolicy<T>>
+inline RT_API_ATTRS T ErfcScaledPositive(T x) {
+ // Per-type direct/series threshold; rationale in the file comment.
+ constexpr T threshold{std::is_same_v<T, float> ? T{9} : T{16}};
+ if (x < threshold) {
+ return ErfcScaledExpOfSquare<T, P>(x) * P::Erfc(x);
+ }
+ // Asymptotic series, truncated after kTerms terms. Truncation error at the
+ // binary128 threshold x = 16 measures <= 0.59 eps and shrinks above it;
+ // do not read the count as giving correctly rounded results.
+ constexpr int kTerms{22};
+ const T inv2x2{T{1} / (T{2} * x * x)};
+ T term{1};
+ T sum{term};
+ for (int k{1}; k < kTerms; ++k) {
+ term *= -static_cast<T>(2 * k - 1) * inv2x2;
+ sum += term;
+ }
+ // Divide by sqrt(pi) before dividing by x: x*sqrt(pi) can overflow near
+ // the top of the format (x >~ 6.7e4931 in binary128) where the true result
+ // is a representable subnormal.
+ return (sum / ErfcScaledSqrtPi<T>()) / x;
+}
+
+template <typename T, typename P = ErfcScaledStdHostPolicy<T>>
+inline RT_API_ATTRS T ErfcScaledAccurate(T x) {
+ if (x < T{0}) {
+ T ax{-x};
+ if (ax >= T{107}) { // universal cutoff; rationale in the file comment
+ return P::Infinity();
+ }
+ // erfc(-x) = 2 - erfc(x), so erfcx(-ax) = 2*exp(ax*ax) - erfcx(ax).
+ // No catastrophic cancellation: the minuend 2*exp(ax*ax) >= 2 while the
+ // subtrahend erfcx(ax) <= 1.
+ return T{2} * ErfcScaledExpOfSquare<T, P>(ax) -
+ ErfcScaledPositive<T, P>(ax);
+ }
+ return ErfcScaledPositive<T, P>(x);
+}
+
+} // namespace Fortran::common
+#endif // FORTRAN_COMMON_ERFC_SCALED_ACCURATE_H_
diff --git a/flang/lib/Evaluate/intrinsics-library.cpp b/flang/lib/Evaluate/intrinsics-library.cpp
index f2c1a7bfaf50b..d76adb7b25883 100644
--- a/flang/lib/Evaluate/intrinsics-library.cpp
+++ b/flang/lib/Evaluate/intrinsics-library.cpp
@@ -14,7 +14,7 @@
#include "flang/Evaluate/intrinsics-library.h"
#include "fold-implementation.h"
#include "host.h"
-#include "flang/Common/erfc-scaled.h"
+#include "flang/Common/erfc-scaled-accurate.h"
#include "flang/Common/idioms.h"
#include "flang/Common/static-multimap-view.h"
#include "flang/Evaluate/expression.h"
@@ -235,7 +235,8 @@ struct HostRuntimeLibrary<HostT, LibraryVersion::Libm> {
FolderFactory<F, F{std::cosh}>::Create("cosh"),
FolderFactory<F, F{std::erf}>::Create("erf"),
FolderFactory<F, F{std::erfc}>::Create("erfc"),
- FolderFactory<F, F{common::ErfcScaled}>::Create("erfc_scaled"),
+ FolderFactory<F, F{common::ErfcScaledAccurate<HostT>}>::Create(
+ "erfc_scaled"),
FolderFactory<F, F{std::exp}>::Create("exp"),
FolderFactory<F, F{std::tgamma}>::Create("gamma"),
FolderFactory<F, F{std::log}>::Create("log"),
@@ -247,7 +248,8 @@ struct HostRuntimeLibrary<HostT, LibraryVersion::Libm> {
FolderFactory<F, F{std::tan}>::Create("tan"),
FolderFactory<F, F{std::tanh}>::Create("tanh"),
};
- // Note: cmath does not have modulo and erfc_scaled equivalent
+ // Note: cmath does not have a modulo equivalent; erfc_scaled is folded
+ // through flang/Common/erfc-scaled-accurate.h.
// Note regarding lack of bessel function support:
// C++17 defined standard Bessel math functions std::cyl_bessel_j
@@ -435,7 +437,8 @@ struct HostRuntimeLibrary<std::complex<double>, LibraryVersion::Libm> {
#endif // _AIX
// Note regarding cmath:
-// - cmath does not have modulo and erfc_scaled equivalent
+// - cmath does not have a modulo equivalent (erfc_scaled, which it also
+// lacks, is folded through flang/Common/erfc-scaled-accurate.h)
// - C++17 defined standard Bessel math functions std::cyl_bessel_j
// and std::cyl_neumann that can be used for Fortran j and y
// bessel functions. However, they are not yet implemented in
@@ -516,6 +519,22 @@ template <> struct HostRuntimeLibrary<double, LibraryVersion::LibmExtensions> {
#endif // _WIN32
#if HAS_QUADMATHLIB
+// Math policy routing flang/Common/erfc-scaled-accurate.h to libquadmath.
+// Do not use std::numeric_limits<__float128> here: libstdc++ leaves it
+// unspecialized (all members zero).
+struct ErfcScaledQuadPolicy {
+ static __float128 Exp(__float128 x) { return ::expq(x); }
+ static __float128 Erfc(__float128 x) { return ::erfcq(x); }
+ static __float128 Fma(__float128 x, __float128 y, __float128 z) {
+ return ::fmaq(x, y, z);
+ }
+ // Not HUGE_VALQ: it expands to __builtin_huge_valq(), which clang does not
+ // implement. INFINITY converts exactly (same as flang-rt's math-entries.h).
+ static __float128 Infinity() { return INFINITY; }
+};
+static __float128 ErfcScaledF128(__float128 x) {
+ return common::ErfcScaledAccurate<__float128, ErfcScaledQuadPolicy>(x);
+}
template <> struct HostRuntimeLibrary<__float128, LibraryVersion::Libm> {
using F = FuncPointer<__float128, __float128>;
using F2 = FuncPointer<__float128, __float128, __float128>;
@@ -538,6 +557,7 @@ template <> struct HostRuntimeLibrary<__float128, LibraryVersion::Libm> {
FolderFactory<F, F{::coshq}>::Create("cosh"),
FolderFactory<F, F{::erfq}>::Create("erf"),
FolderFactory<F, F{::erfcq}>::Create("erfc"),
+ FolderFactory<F, F{&ErfcScaledF128}>::Create("erfc_scaled"),
FolderFactory<F, F{::expq}>::Create("exp"),
FolderFactory<F, F{::tgammaq}>::Create("gamma"),
FolderFactory<F, F{::logq}>::Create("log"),
diff --git a/flang/test/CMakeLists.txt b/flang/test/CMakeLists.txt
index 50980a521b241..a9e290e7e95b0 100644
--- a/flang/test/CMakeLists.txt
+++ b/flang/test/CMakeLists.txt
@@ -61,6 +61,15 @@ check_cxx_source_compiles(
"
HAVE_LDBL_MANT_DIG_113)
+# Check whether the flang compiler itself can fold REAL(16) intrinsic calls
+# through libquadmath. This mirrors the HAS_QUADMATHLIB condition in
+# flang/lib/Evaluate/CMakeLists.txt (FOUND_QUADMATH_LIB is computed and
+# cached there, which runs before this directory).
+set(FLANG_FOLD_REAL16_QUADMATH "")
+if (FLANG_INCLUDE_QUADMATH_H AND FOUND_QUADMATH_LIB)
+ set(FLANG_FOLD_REAL16_QUADMATH "TRUE")
+endif ()
+
# FIXME In out-of-tree builds, "SHLIBDIR" is undefined and passing it to
# `configure_lit_site_cfg` leads to a configuration error. This is currently
# only required by plugins/examples, which are not supported in out-of-tree
diff --git a/flang/test/Evaluate/fold-erfc-scaled.f90 b/flang/test/Evaluate/fold-erfc-scaled.f90
index b38cd0157d0ba..75140c28ad7b4 100644
--- a/flang/test/Evaluate/fold-erfc-scaled.f90
+++ b/flang/test/Evaluate/fold-erfc-scaled.f90
@@ -1,7 +1,49 @@
! RUN: %python %S/test_folding.py %s %flang_fc1
+! Folding of ERFC_SCALED for kinds 4 and 8.
+!
+! On the asymptotic-series branch (x >= 16 in general, x >= 9 for kind 4) the
+! computation involves no math-library calls, so the folded values are
+! host-independent and compared exactly; the expected literals are the exact
+! decimal expansions of the algorithm's results, which measure within 2 ulps
+! of a 300-bit mpmath oracle. Direct-branch and negative results pass through
+! host exp/erfc, whose last-ulp rounding may vary across hosts, and are
+! compared to a small relative tolerance instead.
module m
- real(4), parameter :: x20_4 = erfc_scaled(20._4)
- logical, parameter :: t20_4 = x20_4 == 0.02817435003817081451416015625_4
- real(8), parameter :: x20_8 = erfc_scaled(20._8)
- logical, parameter :: t20_8 = x20_8 == 0.0281743487410513193669459042212110944092273712158203125_8
+ ! Series branch: exact.
+ logical, parameter :: test_p9_4 = &
+ erfc_scaled(9._4) == 0.06230773031711578369140625_4
+ logical, parameter :: test_p95_4 = &
+ erfc_scaled(9.5_4) == 0.0590646751224994659423828125_4
+ logical, parameter :: test_p12_4 = &
+ erfc_scaled(12._4) == 0.0468542166054248809814453125_4
+ logical, parameter :: test_p20_4 = &
+ erfc_scaled(20._4) == 0.02817435003817081451416015625_4
+ logical, parameter :: test_p20_8 = erfc_scaled(20._8) == &
+ 0.028174348741051312428052000313982716761529445648193359375_8
+
+ ! Direct branch: tolerance (about four ulps).
+ real(4), parameter :: tol_4 = 5.0e-7_4
+ real(8), parameter :: tol_8 = 1.0e-15_8
+ real(4), parameter :: ref_p85_4 = 0.065925121307373046875_4
+ logical, parameter :: test_p85_4 = &
+ abs(erfc_scaled(8.5_4) - ref_p85_4) <= tol_4 * ref_p85_4
+
+ ! Negative arguments with representable results.
+ real(4), parameter :: ref_n9_4 = 3.01219472345922556007081388223234048e35_4
+ logical, parameter :: test_n9_4 = &
+ abs(erfc_scaled(-9._4) - ref_n9_4) <= tol_4 * ref_n9_4
+ ! At kind 8 the last finite results lie just below |x| = 26.6287357; the
+ ! previous implementation returned HUGE(x) from x < -26.628 on.
+ real(8), parameter :: ref_sliver_8 = 1.7378490118139657e308_8
+ logical, parameter :: test_sliver_8 = &
+ abs(erfc_scaled(-26.6281_8) - ref_sliver_8) <= tol_8 * ref_sliver_8
+
+ ! Negative arguments whose true value exceeds the format's range fold to
+ ! +infinity.
+ !WARN: warning: overflow on evaluation of intrinsic function or operation [-Wfolding-exception]
+ logical, parameter :: test_n12_4 = erfc_scaled(-12._4) > huge(0._4)
+ !WARN: warning: overflow on evaluation of intrinsic function or operation [-Wfolding-exception]
+ logical, parameter :: test_n27_4 = erfc_scaled(-27._4) > huge(0._4)
+ !WARN: warning: overflow on evaluation of intrinsic function or operation [-Wfolding-exception]
+ logical, parameter :: test_n27_8 = erfc_scaled(-27._8) > huge(0._8)
end
diff --git a/flang/test/Evaluate/fold-erfc-scaled16.f90 b/flang/test/Evaluate/fold-erfc-scaled16.f90
new file mode 100644
index 0000000000000..bd1e479d4c3cf
--- /dev/null
+++ b/flang/test/Evaluate/fold-erfc-scaled16.f90
@@ -0,0 +1,43 @@
+! REQUIRES: flang-fold-real16-quadmath
+! RUN: %python %S/test_folding.py %s %flang_fc1
+! Folding of ERFC_SCALED for REAL(16) on hosts where the compiler folds
+! REAL(16) through libquadmath.
+!
+! The asymptotic-series results (x >= 16) involve no math-library calls and
+! are compared exactly; the expected literals are the exact decimal
+! expansions of the algorithm's results, which sit within one ulp of a
+! 300-bit mpmath oracle. The remaining points go through expq/erfcq, whose
+! last-ulp rounding may vary across libquadmath versions, so they are
+! compared to a few-ulp relative tolerance (binary128 epsilon is 1.93e-34);
+! their reference literals are the correctly rounded true values.
+module m
+ integer, parameter :: qp = selected_real_kind(33)
+ real(qp), parameter :: tol = 5.0e-33_qp
+
+ ! Series branch: exact.
+ logical, parameter :: test_p20 = erfc_scaled(20._qp) == &
+ 2.817434874105131931864915453447075566e-2_qp
+ logical, parameter :: test_p1e6 = erfc_scaled(1000000._qp) == &
+ 5.641895835474741921563059965594862664e-7_qp
+
+ ! Direct branch.
+ real(qp), parameter :: ref_half = 0.615690344192925874870793422683741924_qp
+ logical, parameter :: test_phalf = &
+ abs(erfc_scaled(0.5_qp) - ref_half) <= tol * ref_half
+
+ ! Negative arguments with representable binary128 results. The previous
+ ! implementation returned HUGE(x) for x < -26.628; the true values stay
+ ! representable down to x ~ -106.567.
+ real(qp), parameter :: ref_n27 = &
+ 7.97457052408519312709372209466870032e316_qp
+ logical, parameter :: test_n27 = &
+ abs(erfc_scaled(-27._qp) - ref_n27) <= tol * ref_n27
+ real(qp), parameter :: ref_n1063 = &
+ 4.94361012472232021713994535196543716e4907_qp
+ logical, parameter :: test_n1063 = &
+ abs(erfc_scaled(-106.3_qp) - ref_n1063) <= tol * ref_n1063
+
+ ! Beyond the last representable result, +infinity.
+ !WARN: warning: overflow on evaluation of intrinsic function or operation [-Wfolding-exception]
+ logical, parameter :: test_n107 = erfc_scaled(-107._qp) > huge(0._qp)
+end
diff --git a/flang/test/lit.cfg.py b/flang/test/lit.cfg.py
index e7aa50e001b92..c273ed9c92c88 100644
--- a/flang/test/lit.cfg.py
+++ b/flang/test/lit.cfg.py
@@ -280,6 +280,10 @@ def get_resource_module_intrinsic_dir(modfile):
else:
config.substitutions.append(("%f128-lib", "NONE"))
+# The compiler itself folds REAL(16) intrinsic calls through libquadmath.
+if config.flang_fold_real16_quadmath:
+ config.available_features.add("flang-fold-real16-quadmath")
+
# Set OBJECT_MODE=64 as tools on AIX default to 32-bit.
if "system-aix" in config.available_features:
config.environment["OBJECT_MODE"] = "64"
diff --git a/flang/test/lit.site.cfg.py.in b/flang/test/lit.site.cfg.py.in
index ca94ef4153390..5996d27171184 100644
--- a/flang/test/lit.site.cfg.py.in
+++ b/flang/test/lit.site.cfg.py.in
@@ -32,6 +32,7 @@ config.default_sysroot = "@DEFAULT_SYSROOT@"
config.clang_default_unwindlib = "@CLANG_DEFAULT_UNWINDLIB@"
config.flang_runtime_f128_math_lib = "@FLANG_RUNTIME_F128_MATH_LIB@"
config.have_ldbl_mant_dig_113 = "@HAVE_LDBL_MANT_DIG_113@"
+config.flang_fold_real16_quadmath = "@FLANG_FOLD_REAL16_QUADMATH@"
config.perf_helper_dir = "@LLVM_MAIN_SRC_DIR@/../clang/utils/perf-training"
config.flang_bolt_mode = "@FLANG_BOLT@"
diff --git a/flang/unittests/Evaluate/folding.cpp b/flang/unittests/Evaluate/folding.cpp
index 832e55d44316d..61dabc0f13d6a 100644
--- a/flang/unittests/Evaluate/folding.cpp
+++ b/flang/unittests/Evaluate/folding.cpp
@@ -79,8 +79,91 @@ void TestHostRuntimeSubnormalFlushing() {
}
}
+// Host folding of ERFC_SCALED at REAL(16), on hosts that support it (either
+// __float128 through libquadmath or a binary128 long double). Reference
+// values were computed with mpmath 1.4.1 at 300-bit precision and rounded to
+// binary128 (round-to-nearest-even); inputs and expected results are given
+// as raw words so that no host-dependent decimal conversion contaminates the
+// comparison.
+void TestErfcScaledFoldingReal16() {
+ using R16 = Type<TypeCategory::Real, 16>;
+ using Word = typename Scalar<R16>::Word;
+ const auto real16{[](std::uint64_t hi, std::uint64_t lo) {
+ return Scalar<R16>{Word{hi}.SHIFTL(64).IOR(Word{lo})};
+ }};
+ // True iff a and b, both finite and of the same sign, are at most maxUlps
+ // representable values apart.
+ const auto withinUlps{
+ [](const Scalar<R16> &a, const Scalar<R16> &b, std::uint64_t maxUlps) {
+ Word wa{a.RawBits()};
+ Word wb{b.RawBits()};
+ if (wa.CompareUnsigned(wb) == Fortran::evaluate::Ordering::Less) {
+ std::swap(wa, wb);
+ }
+ Word diff{wa.SubtractSigned(wb).value};
+ return diff.CompareUnsigned(Word{maxUlps}) !=
+ Fortran::evaluate::Ordering::Greater;
+ }};
+ Fortran::parser::CharBlock src;
+ Fortran::parser::ContextualMessages messages{src, nullptr};
+ Fortran::common::IntrinsicTypeDefaultKinds defaults;
+ auto intrinsics{Fortran::evaluate::IntrinsicProcTable::Configure(defaults)};
+ TargetCharacteristics targetCharacteristics;
+ Fortran::common::LanguageFeatureControl languageFeatures;
+ std::set<std::string> tempNames;
+ FoldingContext context{messages, defaults, intrinsics, targetCharacteristics,
+ languageFeatures, tempNames};
+ DynamicType r16{R16{}.GetType()};
+ auto callable{GetHostRuntimeWrapper("erfc_scaled", r16, {r16})};
+ if (!callable) {
+ return; // This host cannot fold REAL(16); folding is legitimately absent.
+ }
+ const auto fold{[&](const Scalar<R16> &x) {
+ return CallHostRt<R16>(*callable, context, x);
+ }};
+ // The results of the asymptotic-series branch (x >= 16) involve no
+ // math-library calls and must reproduce exactly on every host. The
+ // expected words are the algorithm's own results; each sits one ulp from
+ // the correctly rounded value (series truncation): ...1fcf for 20.0 and
+ // ...be12 for 1.0e6.
+ TEST(fold(real16(0x4003400000000000, 0)) == // erfc_scaled(20.0)
+ real16(0x3ff9cd9bc89b7354, 0x7fc6fb33fcba1fce));
+ TEST(fold(real16(0x4012e84800000000, 0)) == // erfc_scaled(1.0e6)
+ real16(0x3fea2ee5a03c7620, 0xf1a7ebc5b5aebe13));
+ // Results that pass through exp/erfc are compared to the correctly
+ // rounded true values with a small allowance, since math-library rounding
+ // differs across implementations (libquadmath versions, long double
+ // libm). Measured deviation on x86-64 with GCC 9.3 libquadmath: 0 ulps at
+ // all three points.
+ constexpr std::uint64_t maxUlps{16};
+ TEST(withinUlps(fold(real16(0x3ffe000000000000, 0)), // erfc_scaled(0.5)
+ real16(0x3ffe3b3bc3c98b0f, 0x2caaf529dbcefe16), maxUlps));
+ // The true values below stay representable in binary128 well past the
+ // former implementation's cutoff at -26.628, which returned HUGE from
+ // there on.
+ TEST(withinUlps(fold(real16(0xc003b00000000000, 0)), // erfc_scaled(-27.0)
+ real16(0x441ba70cd52262b7, 0x9eae660a818fc57b), maxUlps));
+ TEST(withinUlps( // erfc_scaled(-106.3)
+ fold(real16(0xc005a93333333333, 0x3333333333333333)),
+ real16(0x7fae0132469c65b8, 0xbc8f583f6c46a408), maxUlps));
+ // Past the last representable result (x ~ -106.567), +infinity.
+ Scalar<R16> atNeg107{fold(real16(0xc005ac0000000000, 0))};
+ TEST(atNeg107.IsInfinite() && !atNeg107.IsNegative());
+ // Specials.
+ TEST(fold(real16(0x7fff800000000000, 0)).IsNotANumber()); // NaN
+ Scalar<R16> atPosInf{fold(real16(0x7fff000000000000, 0))};
+ TEST(atPosInf.IsZero() && !atPosInf.IsNegative()); // erfc_scaled(+Inf)->+0
+ Scalar<R16> atNegInf{fold(real16(0xffff000000000000, 0))};
+ TEST(atNegInf.IsInfinite() && !atNegInf.IsNegative());
+ TEST(fold(real16(0x0000000000000000, 0)) == // erfc_scaled(+0) == 1
+ real16(0x3fff000000000000, 0));
+ TEST(fold(real16(0x8000000000000000, 0)) == // erfc_scaled(-0) == 1
+ real16(0x3fff000000000000, 0));
+}
+
int main() {
RunOnTypes<TestGetScalarConstantValue, AllIntrinsicTypes>::Run();
TestHostRuntimeSubnormalFlushing();
+ TestErfcScaledFoldingReal16();
return testing::Complete();
}
>From 09ac590a5562835f19e9bd0f3e02f3da6121d192 Mon Sep 17 00:00:00 2001
From: Eugene Epshteyn <eepshteyn at nvidia.com>
Date: Thu, 10 Sep 2026 13:50:42 -0700
Subject: [PATCH 2/3] [flang][test] Print compiler diagnostics when
test_folding.py's compile fails
When the flang -fc1 invocation exits non-zero, test_folding.py raised
CalledProcessError, which swallows the compiler's stderr - the CI log
then shows only a Python traceback with no clue what the compiler
diagnosed. Print the diagnostics and exit instead.
---
flang/test/Evaluate/test_folding.py | 7 ++++++-
1 file changed, 6 insertions(+), 1 deletion(-)
diff --git a/flang/test/Evaluate/test_folding.py b/flang/test/Evaluate/test_folding.py
index dcd1541997c2a..f346218e1180e 100755
--- a/flang/test/Evaluate/test_folding.py
+++ b/flang/test/Evaluate/test_folding.py
@@ -95,10 +95,15 @@ def set_executable(exe):
cmd,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
- check=True,
universal_newlines=True,
cwd=tmpdir,
)
+ if proc.returncode != 0:
+ # Show the compiler's diagnostics rather than swallowing them in a
+ # CalledProcessError traceback.
+ print(f"{cmd} failed with exit status {proc.returncode}:")
+ sys.stdout.write(proc.stderr)
+ sys.exit(1)
src1 = proc.stdout
messages = proc.stderr
>From 2e1b474fd03cbafd784187149f9233275c70e995 Mon Sep 17 00:00:00 2001
From: Eugene Epshteyn <eepshteyn at nvidia.com>
Date: Thu, 10 Sep 2026 14:01:07 -0700
Subject: [PATCH 3/3] [flang][test] Gate REAL(16) folding test on frontend
REAL(16) support
The flang-fold-real16-quadmath lit feature only checked that the
compiler was built against libquadmath, but per
flang/include/flang/Tools/TargetSetup.h the frontend disables REAL(16)
altogether unless FLANG_RUNTIME_F128_MATH_LIB was configured or long
double is binary128. On such hosts (e.g. the Linux premerge bot)
selected_real_kind(33) is -1 and fold-erfc-scaled16.f90 failed with
'REAL(KIND=-1) is not a supported type' instead of being skipped.
Folding coverage on those hosts remains via the folding.cpp unit test,
which calls GetHostRuntimeWrapper below the frontend's type gating.
---
flang/test/CMakeLists.txt | 13 +++++++++----
1 file changed, 9 insertions(+), 4 deletions(-)
diff --git a/flang/test/CMakeLists.txt b/flang/test/CMakeLists.txt
index a9e290e7e95b0..220f4852b065e 100644
--- a/flang/test/CMakeLists.txt
+++ b/flang/test/CMakeLists.txt
@@ -62,11 +62,16 @@ check_cxx_source_compiles(
HAVE_LDBL_MANT_DIG_113)
# Check whether the flang compiler itself can fold REAL(16) intrinsic calls
-# through libquadmath. This mirrors the HAS_QUADMATHLIB condition in
-# flang/lib/Evaluate/CMakeLists.txt (FOUND_QUADMATH_LIB is computed and
-# cached there, which runs before this directory).
+# through libquadmath. Two conditions must hold: the compiler was built
+# against libquadmath (mirrors the HAS_QUADMATHLIB condition in
+# flang/lib/Evaluate/CMakeLists.txt; FOUND_QUADMATH_LIB is computed and
+# cached there, which runs before this directory), and the frontend accepts
+# REAL(16) at all - per flang/include/flang/Tools/TargetSetup.h it disables
+# the type unless FLANG_RUNTIME_F128_MATH_LIB was set or long double is
+# binary128.
set(FLANG_FOLD_REAL16_QUADMATH "")
-if (FLANG_INCLUDE_QUADMATH_H AND FOUND_QUADMATH_LIB)
+if (FLANG_INCLUDE_QUADMATH_H AND FOUND_QUADMATH_LIB
+ AND (FLANG_RUNTIME_F128_MATH_LIB OR HAVE_LDBL_MANT_DIG_113))
set(FLANG_FOLD_REAL16_QUADMATH "TRUE")
endif ()
More information about the flang-commits
mailing list