[flang-commits] [flang] [llvm] [flang-rt] Implement ERFC_SCALED for REAL(16) (PR #219697)
Dmitry Mikushin via flang-commits
flang-commits at lists.llvm.org
Sat Aug 29 23:54:58 PDT 2026
https://github.com/dmikushin updated https://github.com/llvm/llvm-project/pull/219697
>From 3a14758a2b7bd231a74387f7d7c3e6c7c8701dc9 Mon Sep 17 00:00:00 2001
From: Dmitry Mikushin <dmitry at kernelgen.org>
Date: Fri, 28 Aug 2026 13:59:36 +0200
Subject: [PATCH 1/4] [flang-rt] Implement ERFC_SCALED at REAL(16) over the
quadmath entry points
ErfcScaled16 was defined in lib/runtime/numeric.cpp under HAS_LDBL128 alone,
while every neighbouring intrinsic there uses HAS_LDBL128 || HAS_FLOAT128. On
x86-64 the symbol was therefore absent and a REAL(16) program using
ERFC_SCALED compiled cleanly and failed at link.
Widening that guard looks like the fix and is not. The definition forwards to
the Cody approximation in flang/Common/erfc-scaled.h, which calls std::exp on
its working type; for __float128 that overload is ambiguous rather than absent,
so the template does not compile there at all. And if it did, it should not be
used: its coefficients carry about seventeen significant digits and its
sqrtpi/rsqrtpi are long double literals, so at binary128 it is wrong by up to
2.8e-17 relative - roughly 1.5e17 epsilons, the accuracy of double precision in
a container that promises thirty-three digits. Measured against mpmath at 50
decimal places; gfortran on the same points is correctly rounded.
So REAL(16) gets its own implementation, over entry points this build already
has, in the directory where the other kind-16 numerics live.
Below x = 16 it is exp(x*x)*erfc(x) with the square taken exactly - x*x = u + e
via FMA, and exp(u)*(1 + e) - because exp amplifies the rounding of a squaring
by a factor of x*x, some 250 epsilons near the threshold. Above x = 16 it is
the asymptotic expansion to 22 terms. The window for that switch is [12, 106.5]
and is dictated, not chosen: the series is still 1.3e8 epsilons out at x = 8
and reaches one epsilon by x = 12, while exp(x*x) overflows binary128 at
x = 106.567 and erfc underflows at 106.536. Sixteen sits with margin on both
sides. Negative arguments use 2*exp(x*x) - ERFC_SCALED(|x|) and saturate where
the format does.
sqrt(pi) is the sum of three doubles rather than one literal: unsuffixed it
would be a double, an L suffix would cap it at the 64-bit mantissa of x86-64
long double and floor every result near 1e-19, and a Q suffix would be wrong
where this type is long double.
Measured over sixteen points including both sides of the threshold, negative
arguments and 1e5: worst 0.51 epsilons, against 1.07 for gfortran on the same
grid. The grid is deliberately not a regular one - a step of 0.5 gives exactly
representable squares, where the x*x amplification vanishes and the defect this
compensation exists for is invisible.
---
flang-rt/lib/quadmath/CMakeLists.txt | 1 +
flang-rt/lib/quadmath/erfc-scaled.cpp | 111 ++++++++++++++++++++++++++
flang-rt/lib/runtime/numeric.cpp | 10 +--
3 files changed, 116 insertions(+), 6 deletions(-)
create mode 100644 flang-rt/lib/quadmath/erfc-scaled.cpp
diff --git a/flang-rt/lib/quadmath/CMakeLists.txt b/flang-rt/lib/quadmath/CMakeLists.txt
index 353c72def24de..ae8c6fba4e005 100644
--- a/flang-rt/lib/quadmath/CMakeLists.txt
+++ b/flang-rt/lib/quadmath/CMakeLists.txt
@@ -31,6 +31,7 @@ set(sources
cosh.cpp
erf.cpp
erfc.cpp
+ erfc-scaled.cpp
exp.cpp
exponent.cpp
floor.cpp
diff --git a/flang-rt/lib/quadmath/erfc-scaled.cpp b/flang-rt/lib/quadmath/erfc-scaled.cpp
new file mode 100644
index 0000000000000..8177f6945fe0a
--- /dev/null
+++ b/flang-rt/lib/quadmath/erfc-scaled.cpp
@@ -0,0 +1,111 @@
+//===-- lib/quadmath/erfc-scaled.cpp -----------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+//
+// ERFC_SCALED at REAL(16).
+//
+// The other kinds share a Cody rational approximation in
+// flang/include/flang/Common/erfc-scaled.h. It is not used here, deliberately.
+// Its coefficient arrays carry about seventeen significant digits and its
+// sqrtpi/rsqrtpi are long double literals, so evaluated at __float128 it
+// returns a plausible answer that is wrong by up to 2.8e-17 relative - some
+// 1.5e17 epsilons of binary128, which is the accuracy of double precision
+// delivered in a container that promises thirty-three digits. That was measured
+// against mpmath at 50 decimal places before this file was written; gfortran on
+// the same points is correctly rounded to within 0.91 eps.
+//
+// This uses the exact entry points the quadmath build already provides.
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if HAS_LDBL128 || HAS_FLOAT128
+
+using F128Ty = CppTypeFor<TypeCategory::Real, 16>;
+
+// Where the two branches meet.
+//
+// Below it, ERFC_SCALED is exp(x*x)*erfc(x) evaluated directly, which is only
+// possible while neither factor leaves the format: exp(x*x) overflows binary128
+// above x = 106.567 (x*x > 16384*ln2) and erfc underflows above x = 106.536.
+//
+// Above it, the asymptotic expansion is used. That reaches one epsilon of
+// binary128 from about x = 12 upwards - 1.4e-6 eps at x = 12, but still 1.3e8
+// eps at x = 8 - so it must not be used below there.
+//
+// The usable window is therefore [12, 106.5], and sixteen sits with margin on
+// both sides: the series is already good to 0.59 eps there, and the direct form
+// is a factor of six short of overflowing.
+static constexpr F128Ty kThreshold{16};
+
+// Terms of the asymptotic series. Twenty-two are needed at the threshold; the
+// series does not begin to diverge until k ~ x*x = 256, so the same count is
+// safe everywhere above it and only grows more accurate.
+static constexpr int kTerms{22};
+
+// sqrt(pi), as the sum of three doubles.
+//
+// Not one literal: an unsuffixed literal is a double, and an L suffix would cap
+// this at the 64-bit mantissa of x86-64 long double and put a 1e-19 floor under
+// every result. A Q suffix would be wrong when this type is long double. Three
+// exactly representable doubles sum to sqrt(pi) within 7.6e-17 eps of binary128
+// and need no suffix at all.
+static constexpr F128Ty kSqrtPi{static_cast<F128Ty>(1.772453850905516) +
+ static_cast<F128Ty>(-7.666586499825799e-17) +
+ static_cast<F128Ty>(-1.3058334907945429e-33)};
+
+// exp() amplifies the rounding error of a squaring by a factor of x*x - some
+// 250 epsilons near the threshold. So the square is taken exactly: x*x = u + e
+// with e the FMA residual, and exp(u + e) = exp(u)*exp(e) is expanded as
+// exp(u)*(1 + e), which suffices because e is below 1e-30 here and the dropped
+// e*e/2 is below 1e-60.
+//
+// Without this the error reaches 123 eps on arguments whose square is not
+// exactly representable - and stays near 0.4 eps on arguments whose square is,
+// which is exactly why an accuracy test must not step by 0.5.
+static F128Ty ScaledExpOfSquare(F128Ty x) {
+ F128Ty u{x * x};
+ F128Ty e{Fma<true>::invoke(x, x, -u)};
+ return Exp<true>::invoke(u) * (1 + e);
+}
+
+static F128Ty ErfcScaledPositive(F128Ty x) {
+ if (x < kThreshold) {
+ return ScaledExpOfSquare(x) * Erfc<true>::invoke(x);
+ }
+ // 1/(x*sqrt(pi)) * (1 - 1/(2x^2) + 3/(4x^4) - 15/(8x^6) + ...)
+ F128Ty inv2x2{1 / (2 * x * x)};
+ F128Ty sum{1}, term{1};
+ for (int k{1}; k <= kTerms; ++k) {
+ term *= -static_cast<F128Ty>(2 * k - 1) * inv2x2;
+ sum += term;
+ }
+ return sum / (x * kSqrtPi);
+}
+
+F128Ty RTDEF(ErfcScaled16)(F128Ty x) {
+ if (x >= 0) {
+ return ErfcScaledPositive(x);
+ }
+ // erfc(-x) = 2 - erfc(x), so ERFC_SCALED(-x) = 2*exp(x*x) - ERFC_SCALED(x).
+ // The asymptotic branch describes the decaying tail only and must never see a
+ // negative argument; this side grows instead, and leaves the format at the
+ // same place the direct branch does.
+ F128Ty ax{-x};
+ if (ax >= 106) {
+ return F128_RT_INFINITY;
+ }
+ return 2 * ScaledExpOfSquare(ax) - ErfcScaledPositive(ax);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang-rt/lib/runtime/numeric.cpp b/flang-rt/lib/runtime/numeric.cpp
index 78f148dbc5d8c..4514e311fff38 100644
--- a/flang-rt/lib/runtime/numeric.cpp
+++ b/flang-rt/lib/runtime/numeric.cpp
@@ -356,12 +356,10 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(ErfcScaled10)(
return ErfcScaled(x);
}
#endif
-#if HAS_LDBL128
-CppTypeFor<TypeCategory::Real, 16> RTDEF(ErfcScaled16)(
- CppTypeFor<TypeCategory::Real, 16> x) {
- return ErfcScaled(x);
-}
-#endif
+// ErfcScaled16 is not defined here. The shared approximation this file uses is
+// accurate to about seventeen digits, which is a thousand times short of what
+// binary128 holds, so REAL(16) has its own implementation over the exact
+// quadmath entry points - see flang-rt/lib/quadmath/erfc-scaled.cpp.
CppTypeFor<TypeCategory::Integer, 4> RTDEF(Exponent4_4)(
CppTypeFor<TypeCategory::Real, 4> x) {
>From 349fc48e1df3d6fe07d21e6540b34e00b43d70b9 Mon Sep 17 00:00:00 2001
From: Dmitry Mikushin <dmitry at kernelgen.org>
Date: Fri, 28 Aug 2026 14:06:37 +0200
Subject: [PATCH 2/4] [flang-rt][test] Cover ERFC_SCALED at REAL(16), and give
flang-rt an f128 gate
flang-rt/test had no REAL(16) tests at all and no way to express that a test
needs them: the feature flang-supports-f128-math existed only in flang/test.
This adds the same feature here, derived from the same two CMake values, so
that the compiler-side and runtime-side suites cannot disagree about whether
this build has quad precision. The dependency hook already existed
(if (TARGET flang_rt.quadmath)).
The test checks accuracy, not merely that the symbol resolves. References come
from mpmath at 50 decimal places rather than from another compiler:
ERFC_SCALED is not a libm entry point, so every implementation is somebody's
approximation, and two of them agreeing would show common ancestry rather than
correctness.
The grid is deliberately irregular and the source says so at length, because
the natural instinct is to tidy it. A step of 0.5 - and the tidy values 0.125,
0.5, 1.0, 2.0, 5.0, 13.75 - all have exactly representable squares, where the
x*x amplification in exp(x*x) vanishes. An uncompensated implementation
measures 0.4 epsilons on such points and 123 off them, so a regular grid would
report the defect this code exists to avoid as absent. Both sides of the
branch threshold are present against the same references, so the switch at
x = 16 cannot hide a discontinuity; negative arguments and 1e3 are covered too.
The tolerance is two epsilons rather than one. Both independently chosen grids
came in under one during development, but erfc at binary128 comes from an
external library whose future versions may move the last place.
Seen red before being trusted: perturbing one reference digit gives 4.29
epsilons, "TOO LARGE", and FileCheck exits 1.
---
flang-rt/test/Runtime/erfc-scaled-real16.f90 | 71 ++++++++++++++++++++
flang-rt/test/lit.cfg.py | 8 +++
flang-rt/test/lit.site.cfg.py.in | 2 +
3 files changed, 81 insertions(+)
create mode 100644 flang-rt/test/Runtime/erfc-scaled-real16.f90
diff --git a/flang-rt/test/Runtime/erfc-scaled-real16.f90 b/flang-rt/test/Runtime/erfc-scaled-real16.f90
new file mode 100644
index 0000000000000..ba05dab8ab04f
--- /dev/null
+++ b/flang-rt/test/Runtime/erfc-scaled-real16.f90
@@ -0,0 +1,71 @@
+! Accuracy of ERFC_SCALED at REAL(16), against values computed independently.
+!
+! REQUIRES: flang-supports-f128-math
+! RUN: %flang %s -o %t && %t | FileCheck %s
+!
+! The reference values below are from mpmath at 50 decimal places, not from
+! another Fortran compiler. ERFC_SCALED is not a libm entry point: every
+! implementation is somebody's approximation, so two of them agreeing would
+! show common ancestry rather than correctness. Only an independent ground
+! truth distinguishes "both right" from "both wrong the same way".
+!
+! THE GRID IS DELIBERATELY NOT REGULAR, AND MUST STAY THAT WAY. A step of 0.5
+! -- or any of the tidy values 0.125, 0.5, 1.0, 2.0, 5.0, 13.75 -- gives an
+! exactly representable square, and exp(x*x) amplifies the rounding of that
+! squaring by a factor of x*x. On such points the error is zero by
+! construction and an uncompensated implementation measures 0.4 epsilons; off
+! them it measures 123. Tidying this grid would blind the test to the defect
+! the implementation exists to avoid.
+!
+! Both sides of the branch threshold are here, against the same references, so
+! that the switch at x = 16 cannot hide a discontinuity.
+!
+! The tolerance is 2 epsilons rather than 1. Both independently chosen grids
+! measured under 1 during development, but erfc at binary128 comes from an
+! external library whose future versions may move the last place.
+
+program erfc_scaled_real16
+ implicit none
+ integer, parameter :: qp = selected_real_kind(33)
+ integer, parameter :: n = 11
+ real(qp), parameter :: delta = 1.0e-10_qp
+ ! Not tidy on purpose: see above.
+ real(qp), parameter :: x(n) = [ &
+ 0.1_qp, 1.7_qp, 15.9_qp, &
+ 16.0_qp - delta, 16.0_qp + delta, &
+ 20.0_qp, 64.0_qp, 1000.0_qp, &
+ -1.7_qp, -13.75_qp, 0.0_qp]
+ ! exp(x*x)*erfc(x), mpmath, 50 dps, rounded to binary128.
+ real(qp), parameter :: want(n) = [ &
+ 0.89645697996912664193188374864404227_qp, &
+ 0.291663297075343466454843377681505751_qp, &
+ 0.0354138554979901787367777955827486539_qp, &
+ 0.0351933778251499452361469135833690837_qp, &
+ 0.0351933778247117298966017592289243255_qp, &
+ 0.0281743487410513193186491545344707584_qp, &
+ 0.00881438653054441357566924341796977317_qp, &
+ 0.000564189301453387654199745028061695727_qp, &
+ 35.6949559060252863357458060477745282_qp, &
+ 2.56939266805496674962846043540050838e+82_qp, &
+ 1.0_qp]
+ real(qp) :: got, rel, worst
+ integer :: i
+
+ worst = 0.0_qp
+ do i = 1, n
+ got = erfc_scaled(x(i))
+ if (want(i) == 0.0_qp) then
+ rel = abs(got)
+ else
+ rel = abs(got - want(i))/abs(want(i))
+ end if
+ worst = max(worst, rel)
+ end do
+
+ ! CHECK: worst deviation in epsilons:
+ ! CHECK-SAME: within tolerance
+ write(*,'(A,ES12.5,A)') 'worst deviation in epsilons: ', &
+ worst/epsilon(1.0_qp), &
+ merge(' within tolerance', ' TOO LARGE ', &
+ worst <= 2.0_qp*epsilon(1.0_qp))
+end program erfc_scaled_real16
diff --git a/flang-rt/test/lit.cfg.py b/flang-rt/test/lit.cfg.py
index 74958d8155ef6..60d43f7af883a 100644
--- a/flang-rt/test/lit.cfg.py
+++ b/flang-rt/test/lit.cfg.py
@@ -111,3 +111,11 @@ def shjoin(args, sep=" "):
# Set OBJECT_MODE=64 as tools on AIX default to 32-bit.
if "system-aix" in config.available_features:
config.environment["OBJECT_MODE"] = "64"
+
+# Tests that need REAL(16) require the runtime to have been built with a
+# quad-precision math library. The condition is the same one flang/test uses
+# (see flang/test/lit.cfg.py), and it must stay the same: the compiler and the
+# runtime derive REAL(16) support from this single CMake value, and a test that
+# disagreed with either would report a configuration difference as a defect.
+if config.flang_runtime_f128_math_lib or config.have_ldbl_mant_dig_113:
+ config.available_features.add("flang-supports-f128-math")
diff --git a/flang-rt/test/lit.site.cfg.py.in b/flang-rt/test/lit.site.cfg.py.in
index 8e3c902d4a3a6..b10ad7b4515c3 100644
--- a/flang-rt/test/lit.site.cfg.py.in
+++ b/flang-rt/test/lit.site.cfg.py.in
@@ -15,6 +15,8 @@ config.cc = "@CMAKE_C_COMPILER@"
config.flang = "@CMAKE_Fortran_COMPILER@"
config.osx_sysroot = path(r"@CMAKE_OSX_SYSROOT@")
config.target_triple = "@LLVM_TARGET_TRIPLE@"
+config.flang_runtime_f128_math_lib = "@FLANG_RUNTIME_F128_MATH_LIB@"
+config.have_ldbl_mant_dig_113 = "@HAVE_LDBL_MANT_DIG_113@"
import lit.llvm
lit.llvm.initialize(lit_config, config)
>From 90e11a358367f0c0ff44643f5a99253a614560ff Mon Sep 17 00:00:00 2001
From: Dmitry Mikushin <dmitry at kernelgen.org>
Date: Sat, 29 Aug 2026 19:43:28 +0200
Subject: [PATCH 3/4] [flang-rt] Fit the file banner in 80 columns
clang-format flagged the header line at 81 characters and wanted to wrap
it, which would have broken the LLVM banner across two lines. One dash
fewer is the fix; the banner stays on one line and the format check is
clean.
---
flang-rt/lib/quadmath/erfc-scaled.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/flang-rt/lib/quadmath/erfc-scaled.cpp b/flang-rt/lib/quadmath/erfc-scaled.cpp
index 8177f6945fe0a..9f51ab07487a4 100644
--- a/flang-rt/lib/quadmath/erfc-scaled.cpp
+++ b/flang-rt/lib/quadmath/erfc-scaled.cpp
@@ -1,4 +1,4 @@
-//===-- lib/quadmath/erfc-scaled.cpp -----------------------------*- C++ -*-===//
+//===-- lib/quadmath/erfc-scaled.cpp ----------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
>From 95cb4c0eb7ec8092197c0ec95d44e40342b8255a Mon Sep 17 00:00:00 2001
From: Dmitry Mikushin <dmitry at kernelgen.org>
Date: Sun, 30 Aug 2026 08:54:29 +0200
Subject: [PATCH 4/4] [flang][docs] Release note for ERFC_SCALED at REAL(16)
The failure mode is worth stating in the note rather than just the fix:
it was not a diagnostic but a run-time abort, so anyone who hit it saw a
crash in the runtime with no indication that the intrinsic was simply
absent at that kind.
---
flang/docs/ReleaseNotes.md | 6 ++++++
1 file changed, 6 insertions(+)
diff --git a/flang/docs/ReleaseNotes.md b/flang/docs/ReleaseNotes.md
index 079a8e5ba0e94..c4292f0d9130f 100644
--- a/flang/docs/ReleaseNotes.md
+++ b/flang/docs/ReleaseNotes.md
@@ -31,6 +31,12 @@ page](https://llvm.org/releases/).
## Bug Fixes
+- `ERFC_SCALED` no longer aborts at run time when called on a `REAL(16)`
+ argument. The entry point was only compiled where `long double` is itself
+ binary128, so on x86-64 a program using it built and linked and then died
+ inside the runtime. It is now implemented over the quad-precision math
+ library, to an accuracy the existing shared approximation could not reach.
+
## Non-comprehensive list of changes in this release
- Added support for the OpenMP implementation-defined extension sentinels
More information about the flang-commits
mailing list