diff --git a/flang/runtime/Float128Math/CMakeLists.txt b/flang/runtime/Float128Math/CMakeLists.txt
index f11678cd70b769..60d44c78be0faf 100644
--- a/flang/runtime/Float128Math/CMakeLists.txt
+++ b/flang/runtime/Float128Math/CMakeLists.txt
@@ -59,7 +59,9 @@ set(sources
   erf.cpp
   erfc.cpp
   exp.cpp
+  exponent.cpp
   floor.cpp
+  fraction.cpp
   hypot.cpp
   j0.cpp
   j1.cpp
@@ -69,11 +71,18 @@ set(sources
   log.cpp
   log10.cpp
   lround.cpp
+  mod-real.cpp
+  modulo-real.cpp
+  nearest.cpp
   norm2.cpp
   pow.cpp
   round.cpp
+  rrspacing.cpp
+  scale.cpp
+  set-exponent.cpp
   sin.cpp
   sinh.cpp
+  spacing.cpp
   sqrt.cpp
   tan.cpp
   tanh.cpp
diff --git a/flang/runtime/Float128Math/acos.cpp b/flang/runtime/Float128Math/acos.cpp
index 531c79c7444bd3..14ff6944856844 100644
--- a/flang/runtime/Float128Math/acos.cpp
+++ b/flang/runtime/Float128Math/acos.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(AcosF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Acos<RTNAME(AcosF128)>::invoke(x);
+  return Acos<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/acosh.cpp b/flang/runtime/Float128Math/acosh.cpp
index 1495120edd1a07..9d70804e44a470 100644
--- a/flang/runtime/Float128Math/acosh.cpp
+++ b/flang/runtime/Float128Math/acosh.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(AcoshF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Acosh<RTNAME(AcoshF128)>::invoke(x);
+  return Acosh<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/asin.cpp b/flang/runtime/Float128Math/asin.cpp
index 2fb8c6c5e97d71..6781b23f0363db 100644
--- a/flang/runtime/Float128Math/asin.cpp
+++ b/flang/runtime/Float128Math/asin.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(AsinF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Asin<RTNAME(AsinF128)>::invoke(x);
+  return Asin<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/asinh.cpp b/flang/runtime/Float128Math/asinh.cpp
index 3630a77be42b2c..1310bc61c1de0f 100644
--- a/flang/runtime/Float128Math/asinh.cpp
+++ b/flang/runtime/Float128Math/asinh.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(AsinhF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Asinh<RTNAME(AsinhF128)>::invoke(x);
+  return Asinh<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/atan.cpp b/flang/runtime/Float128Math/atan.cpp
index 4609343e9d1273..f01382df90c0ee 100644
--- a/flang/runtime/Float128Math/atan.cpp
+++ b/flang/runtime/Float128Math/atan.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(AtanF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Atan<RTNAME(AtanF128)>::invoke(x);
+  return Atan<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/atan2.cpp b/flang/runtime/Float128Math/atan2.cpp
index c0175e67ec71bd..dd646b0452b115 100644
--- a/flang/runtime/Float128Math/atan2.cpp
+++ b/flang/runtime/Float128Math/atan2.cpp
@@ -15,7 +15,7 @@ extern "C" {
 CppTypeFor<TypeCategory::Real, 16> RTDEF(Atan2F128)(
     CppTypeFor<TypeCategory::Real, 16> x,
     CppTypeFor<TypeCategory::Real, 16> y) {
-  return Atan2<RTNAME(Atan2F128)>::invoke(x, y);
+  return Atan2<true>::invoke(x, y);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/atanh.cpp b/flang/runtime/Float128Math/atanh.cpp
index bfacb967117d70..5fc5ba5debc81a 100644
--- a/flang/runtime/Float128Math/atanh.cpp
+++ b/flang/runtime/Float128Math/atanh.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(AtanhF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Atanh<RTNAME(AtanhF128)>::invoke(x);
+  return Atanh<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/cabs.cpp b/flang/runtime/Float128Math/cabs.cpp
index 827b197a6a81ae..3b8c9d17003c6e 100644
--- a/flang/runtime/Float128Math/cabs.cpp
+++ b/flang/runtime/Float128Math/cabs.cpp
@@ -16,7 +16,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 // NOTE: Flang calls the runtime APIs using C _Complex ABI
 CppTypeFor<TypeCategory::Real, 16> RTDEF(CAbsF128)(CFloat128ComplexType x) {
-  return CAbs<RTNAME(CAbsF128)>::invoke(x);
+  return CAbs<true>::invoke(x);
 }
 #endif
 #endif
diff --git a/flang/runtime/Float128Math/ceil.cpp b/flang/runtime/Float128Math/ceil.cpp
index a53a2c27c616b5..ed4d164a62bedc 100644
--- a/flang/runtime/Float128Math/ceil.cpp
+++ b/flang/runtime/Float128Math/ceil.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(CeilF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Ceil<RTNAME(CeilF128)>::invoke(x);
+  return Ceil<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/cos.cpp b/flang/runtime/Float128Math/cos.cpp
index 845c970bd8e639..b93c92f275f791 100644
--- a/flang/runtime/Float128Math/cos.cpp
+++ b/flang/runtime/Float128Math/cos.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(CosF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Cos<RTNAME(CosF128)>::invoke(x);
+  return Cos<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/cosh.cpp b/flang/runtime/Float128Math/cosh.cpp
index acf6ff4130ee3c..a3662a826dcb1c 100644
--- a/flang/runtime/Float128Math/cosh.cpp
+++ b/flang/runtime/Float128Math/cosh.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(CoshF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Cosh<RTNAME(CoshF128)>::invoke(x);
+  return Cosh<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/erf.cpp b/flang/runtime/Float128Math/erf.cpp
index 862f3b97411873..631f71c76effe7 100644
--- a/flang/runtime/Float128Math/erf.cpp
+++ b/flang/runtime/Float128Math/erf.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(ErfF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Erf<RTNAME(ErfF128)>::invoke(x);
+  return Erf<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/erfc.cpp b/flang/runtime/Float128Math/erfc.cpp
index 0ac0b945563747..ea3cd646d8c4ba 100644
--- a/flang/runtime/Float128Math/erfc.cpp
+++ b/flang/runtime/Float128Math/erfc.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(ErfcF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Erfc<RTNAME(ErfcF128)>::invoke(x);
+  return Erfc<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/exp.cpp b/flang/runtime/Float128Math/exp.cpp
index 50386fdbfb6449..b1161b0f29294c 100644
--- a/flang/runtime/Float128Math/exp.cpp
+++ b/flang/runtime/Float128Math/exp.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(ExpF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Exp<RTNAME(ExpF128)>::invoke(x);
+  return Exp<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/exponent.cpp b/flang/runtime/Float128Math/exponent.cpp
new file mode 100644
index 00000000000000..1be1dd0d0ac8b8
--- /dev/null
+++ b/flang/runtime/Float128Math/exponent.cpp
@@ -0,0 +1,26 @@
+//===-- runtime/Float128Math/exponent.cpp ---------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// EXPONENT (16.9.75)
+CppTypeFor<TypeCategory::Integer, 4> RTDEF(Exponent16_4)(F128Type x) {
+  return Exponent<CppTypeFor<TypeCategory::Integer, 4>>(x);
+}
+CppTypeFor<TypeCategory::Integer, 8> RTDEF(Exponent16_8)(F128Type x) {
+  return Exponent<CppTypeFor<TypeCategory::Integer, 8>>(x);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/floor.cpp b/flang/runtime/Float128Math/floor.cpp
index 48cf4e01448070..78a94984cac8a3 100644
--- a/flang/runtime/Float128Math/floor.cpp
+++ b/flang/runtime/Float128Math/floor.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(FloorF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Floor<RTNAME(FloorF128)>::invoke(x);
+  return Floor<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/fraction.cpp b/flang/runtime/Float128Math/fraction.cpp
new file mode 100644
index 00000000000000..8c9889b7f6871e
--- /dev/null
+++ b/flang/runtime/Float128Math/fraction.cpp
@@ -0,0 +1,21 @@
+//===-- runtime/Float128Math/fraction.cpp ---------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// FRACTION (16.9.80)
+F128Type RTDEF(Fraction16)(F128Type x) { return Fraction(x); }
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/hypot.cpp b/flang/runtime/Float128Math/hypot.cpp
index 33c83a1654993e..b4fa1d66bcfa6a 100644
--- a/flang/runtime/Float128Math/hypot.cpp
+++ b/flang/runtime/Float128Math/hypot.cpp
@@ -15,7 +15,7 @@ extern "C" {
 CppTypeFor<TypeCategory::Real, 16> RTDEF(HypotF128)(
     CppTypeFor<TypeCategory::Real, 16> x,
     CppTypeFor<TypeCategory::Real, 16> y) {
-  return Hypot<RTNAME(HypotF128)>::invoke(x, y);
+  return Hypot<true>::invoke(x, y);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/j0.cpp b/flang/runtime/Float128Math/j0.cpp
index f8f3fe71d8a616..9390a7eeb3c605 100644
--- a/flang/runtime/Float128Math/j0.cpp
+++ b/flang/runtime/Float128Math/j0.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(J0F128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return J0<RTNAME(J0F128)>::invoke(x);
+  return J0<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/j1.cpp b/flang/runtime/Float128Math/j1.cpp
index 9a51b973e1cf88..c54927123388c6 100644
--- a/flang/runtime/Float128Math/j1.cpp
+++ b/flang/runtime/Float128Math/j1.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(J1F128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return J1<RTNAME(J1F128)>::invoke(x);
+  return J1<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/jn.cpp b/flang/runtime/Float128Math/jn.cpp
index 644a66863c0d23..15afd83400c320 100644
--- a/flang/runtime/Float128Math/jn.cpp
+++ b/flang/runtime/Float128Math/jn.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(JnF128)(
     int n, CppTypeFor<TypeCategory::Real, 16> x) {
-  return Jn<RTNAME(JnF128)>::invoke(n, x);
+  return Jn<true>::invoke(n, x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/lgamma.cpp b/flang/runtime/Float128Math/lgamma.cpp
index fff7dfcb9c15db..ac31c89a912b32 100644
--- a/flang/runtime/Float128Math/lgamma.cpp
+++ b/flang/runtime/Float128Math/lgamma.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(LgammaF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Lgamma<RTNAME(LgammaF128)>::invoke(x);
+  return Lgamma<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/llround.cpp b/flang/runtime/Float128Math/llround.cpp
index 00c62818af19db..b77281c507fe7c 100644
--- a/flang/runtime/Float128Math/llround.cpp
+++ b/flang/runtime/Float128Math/llround.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Integer, 8> RTDEF(LlroundF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Llround<RTNAME(LlroundF128)>::invoke(x);
+  return Llround<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/log.cpp b/flang/runtime/Float128Math/log.cpp
index 0cfe329c6f7f59..38e6b581fd849c 100644
--- a/flang/runtime/Float128Math/log.cpp
+++ b/flang/runtime/Float128Math/log.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(LogF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Log<RTNAME(LogF128)>::invoke(x);
+  return Log<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/log10.cpp b/flang/runtime/Float128Math/log10.cpp
index cd8bf27fcb121b..3c89c0e707774f 100644
--- a/flang/runtime/Float128Math/log10.cpp
+++ b/flang/runtime/Float128Math/log10.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(Log10F128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Log10<RTNAME(Log10F128)>::invoke(x);
+  return Log10<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/lround.cpp b/flang/runtime/Float128Math/lround.cpp
index 6ced66a1b2d3af..ce7a228038a1d3 100644
--- a/flang/runtime/Float128Math/lround.cpp
+++ b/flang/runtime/Float128Math/lround.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Integer, 4> RTDEF(LroundF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Lround<RTNAME(LroundF128)>::invoke(x);
+  return Lround<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/math-entries.h b/flang/runtime/Float128Math/math-entries.h
index a0d81d0cbb5407..ad3f6aa18aa9a1 100644
--- a/flang/runtime/Float128Math/math-entries.h
+++ b/flang/runtime/Float128Math/math-entries.h
@@ -13,36 +13,40 @@
 #include "flang/Common/float128.h"
 #include "flang/Runtime/entry-names.h"
 #include <cfloat>
+#include <cmath>
 #include <type_traits>
 
+namespace {
+using namespace Fortran::runtime;
+using F128RetType = CppTypeFor<TypeCategory::Real, 16>;
+using I32RetType = CppTypeFor<TypeCategory::Integer, 4>;
+using I64RetType = CppTypeFor<TypeCategory::Integer, 8>;
+} // namespace
+
 namespace Fortran::runtime {
 
 // Define a class template to gracefully fail, when
 // there is no specialized template that implements
 // the required function via using the third-party
 // implementation.
-#define DEFINE_FALLBACK(caller) \
-  template <auto F> struct caller { \
-    template <typename... ATs> \
-    [[noreturn]] static std::invoke_result_t<decltype(F), ATs...> invoke( \
-        ATs... args) { \
+#define DEFINE_FALLBACK(caller, ret_type) \
+  template <bool = false, typename RT = ret_type> struct caller { \
+    template <typename... ATs> [[noreturn]] static RT invoke(ATs... args) { \
       Terminator terminator{__FILE__, __LINE__}; \
       terminator.Crash("Float128 variant of '%s' is unsupported", #caller); \
     } \
   };
 
 // Define template specialization that is calling the third-party
-// implementation. The template is specialized by a function pointer
-// that is the FortranFloat128Math entry point. The signatures
-// of the caller and the callee must match.
+// implementation.
 //
 // Defining the specialization for any target library requires
 // adding the generic template via DEFINE_FALLBACK, so that
 // a build with another target library that does not define
 // the same alias can gracefully fail in runtime.
 #define DEFINE_SIMPLE_ALIAS(caller, callee) \
-  template <typename RT, typename... ATs, RT (*p)(ATs...)> struct caller<p> { \
-    static RT invoke(ATs... args) { \
+  template <typename RT> struct caller<true, RT> { \
+    template <typename... ATs> static RT invoke(ATs... args) { \
       static_assert(std::is_invocable_r_v<RT, \
           decltype(callee(std::declval<ATs>()...))(ATs...), ATs...>); \
       if constexpr (std::is_same_v<RT, void>) { \
@@ -54,48 +58,58 @@ namespace Fortran::runtime {
   };
 
 // Define fallback callers.
-DEFINE_FALLBACK(Abs)
-DEFINE_FALLBACK(Acos)
-DEFINE_FALLBACK(Acosh)
-DEFINE_FALLBACK(Asin)
-DEFINE_FALLBACK(Asinh)
-DEFINE_FALLBACK(Atan)
-DEFINE_FALLBACK(Atan2)
-DEFINE_FALLBACK(Atanh)
-DEFINE_FALLBACK(Ceil)
-DEFINE_FALLBACK(Cos)
-DEFINE_FALLBACK(Cosh)
-DEFINE_FALLBACK(Erf)
-DEFINE_FALLBACK(Erfc)
-DEFINE_FALLBACK(Exp)
-DEFINE_FALLBACK(Floor)
-DEFINE_FALLBACK(Hypot)
-DEFINE_FALLBACK(J0)
-DEFINE_FALLBACK(J1)
-DEFINE_FALLBACK(Jn)
-DEFINE_FALLBACK(Lgamma)
-DEFINE_FALLBACK(Llround)
-DEFINE_FALLBACK(Lround)
-DEFINE_FALLBACK(Log)
-DEFINE_FALLBACK(Log10)
-DEFINE_FALLBACK(Pow)
-DEFINE_FALLBACK(Round)
-DEFINE_FALLBACK(Sin)
-DEFINE_FALLBACK(Sinh)
-DEFINE_FALLBACK(Sqrt)
-DEFINE_FALLBACK(Tan)
-DEFINE_FALLBACK(Tanh)
-DEFINE_FALLBACK(Tgamma)
-DEFINE_FALLBACK(Trunc)
-DEFINE_FALLBACK(Y0)
-DEFINE_FALLBACK(Y1)
-DEFINE_FALLBACK(Yn)
+#define DEFINE_FALLBACK_F128(caller) DEFINE_FALLBACK(caller, ::F128RetType)
+#define DEFINE_FALLBACK_I32(caller) DEFINE_FALLBACK(caller, ::I32RetType)
+#define DEFINE_FALLBACK_I64(caller) DEFINE_FALLBACK(caller, ::I64RetType)
+
+DEFINE_FALLBACK_F128(Abs)
+DEFINE_FALLBACK_F128(Acos)
+DEFINE_FALLBACK_F128(Acosh)
+DEFINE_FALLBACK_F128(Asin)
+DEFINE_FALLBACK_F128(Asinh)
+DEFINE_FALLBACK_F128(Atan)
+DEFINE_FALLBACK_F128(Atan2)
+DEFINE_FALLBACK_F128(Atanh)
+DEFINE_FALLBACK_F128(Ceil)
+DEFINE_FALLBACK_F128(Cos)
+DEFINE_FALLBACK_F128(Cosh)
+DEFINE_FALLBACK_F128(Erf)
+DEFINE_FALLBACK_F128(Erfc)
+DEFINE_FALLBACK_F128(Exp)
+DEFINE_FALLBACK_F128(Floor)
+DEFINE_FALLBACK_F128(Frexp)
+DEFINE_FALLBACK_F128(Hypot)
+DEFINE_FALLBACK_I32(Ilogb)
+DEFINE_FALLBACK_I32(Isinf)
+DEFINE_FALLBACK_I32(Isnan)
+DEFINE_FALLBACK_F128(J0)
+DEFINE_FALLBACK_F128(J1)
+DEFINE_FALLBACK_F128(Jn)
+DEFINE_FALLBACK_F128(Ldexp)
+DEFINE_FALLBACK_F128(Lgamma)
+DEFINE_FALLBACK_I64(Llround)
+DEFINE_FALLBACK_F128(Log)
+DEFINE_FALLBACK_F128(Log10)
+DEFINE_FALLBACK_I32(Lround)
+DEFINE_FALLBACK_F128(Nextafter)
+DEFINE_FALLBACK_F128(Pow)
+DEFINE_FALLBACK_F128(Qnan)
+DEFINE_FALLBACK_F128(Round)
+DEFINE_FALLBACK_F128(Sin)
+DEFINE_FALLBACK_F128(Sinh)
+DEFINE_FALLBACK_F128(Sqrt)
+DEFINE_FALLBACK_F128(Tan)
+DEFINE_FALLBACK_F128(Tanh)
+DEFINE_FALLBACK_F128(Tgamma)
+DEFINE_FALLBACK_F128(Trunc)
+DEFINE_FALLBACK_F128(Y0)
+DEFINE_FALLBACK_F128(Y1)
+DEFINE_FALLBACK_F128(Yn)
 
 #if HAS_LIBM
-// Define wrapper callers for libm.
-#include <ccomplex>
-#include <cmath>
+#include <limits>
 
+// Define wrapper callers for libm.
 #if LDBL_MANT_DIG == 113
 // Use STD math functions. They provide IEEE-754 128-bit float
 // support either via 'long double' or __float128.
@@ -118,15 +132,21 @@ DEFINE_SIMPLE_ALIAS(Erf, std::erf)
 DEFINE_SIMPLE_ALIAS(Erfc, std::erfc)
 DEFINE_SIMPLE_ALIAS(Exp, std::exp)
 DEFINE_SIMPLE_ALIAS(Floor, std::floor)
+DEFINE_SIMPLE_ALIAS(Frexp, std::frexp)
 DEFINE_SIMPLE_ALIAS(Hypot, std::hypot)
+DEFINE_SIMPLE_ALIAS(Ilogb, std::ilogb)
+DEFINE_SIMPLE_ALIAS(Isinf, std::isinf)
+DEFINE_SIMPLE_ALIAS(Isnan, std::isnan)
 DEFINE_SIMPLE_ALIAS(J0, j0l)
 DEFINE_SIMPLE_ALIAS(J1, j1l)
 DEFINE_SIMPLE_ALIAS(Jn, jnl)
+DEFINE_SIMPLE_ALIAS(Ldexp, std::ldexp)
 DEFINE_SIMPLE_ALIAS(Lgamma, std::lgamma)
 DEFINE_SIMPLE_ALIAS(Llround, std::llround)
-DEFINE_SIMPLE_ALIAS(Lround, std::lround)
 DEFINE_SIMPLE_ALIAS(Log, std::log)
 DEFINE_SIMPLE_ALIAS(Log10, std::log10)
+DEFINE_SIMPLE_ALIAS(Lround, std::lround)
+DEFINE_SIMPLE_ALIAS(Nextafter, std::nextafter)
 DEFINE_SIMPLE_ALIAS(Pow, std::pow)
 DEFINE_SIMPLE_ALIAS(Round, std::round)
 DEFINE_SIMPLE_ALIAS(Sin, std::sin)
@@ -139,6 +159,12 @@ DEFINE_SIMPLE_ALIAS(Trunc, std::trunc)
 DEFINE_SIMPLE_ALIAS(Y0, y0l)
 DEFINE_SIMPLE_ALIAS(Y1, y1l)
 DEFINE_SIMPLE_ALIAS(Yn, ynl)
+
+// Use numeric_limits to produce infinity of the right type.
+#define F128_RT_INFINITY \
+  (std::numeric_limits<CppTypeFor<TypeCategory::Real, 16>>::infinity())
+#define F128_RT_QNAN \
+  (std::numeric_limits<CppTypeFor<TypeCategory::Real, 16>>::quiet_NaN())
 #else // LDBL_MANT_DIG != 113
 #if !HAS_LIBMF128
 // glibc >=2.26 seems to have complete support for __float128
@@ -172,15 +198,21 @@ DEFINE_SIMPLE_ALIAS(Erf, erfq)
 DEFINE_SIMPLE_ALIAS(Erfc, erfcq)
 DEFINE_SIMPLE_ALIAS(Exp, expq)
 DEFINE_SIMPLE_ALIAS(Floor, floorq)
+DEFINE_SIMPLE_ALIAS(Frexp, frexpq)
 DEFINE_SIMPLE_ALIAS(Hypot, hypotq)
+DEFINE_SIMPLE_ALIAS(Ilogb, ilogbq)
+DEFINE_SIMPLE_ALIAS(Isinf, isinfq)
+DEFINE_SIMPLE_ALIAS(Isnan, isnanq)
 DEFINE_SIMPLE_ALIAS(J0, j0q)
 DEFINE_SIMPLE_ALIAS(J1, j1q)
 DEFINE_SIMPLE_ALIAS(Jn, jnq)
+DEFINE_SIMPLE_ALIAS(Ldexp, ldexpq)
 DEFINE_SIMPLE_ALIAS(Lgamma, lgammaq)
 DEFINE_SIMPLE_ALIAS(Llround, llroundq)
-DEFINE_SIMPLE_ALIAS(Lround, lroundq)
 DEFINE_SIMPLE_ALIAS(Log, logq)
 DEFINE_SIMPLE_ALIAS(Log10, log10q)
+DEFINE_SIMPLE_ALIAS(Lround, lroundq)
+DEFINE_SIMPLE_ALIAS(Nextafter, nextafterq)
 DEFINE_SIMPLE_ALIAS(Pow, powq)
 DEFINE_SIMPLE_ALIAS(Round, roundq)
 DEFINE_SIMPLE_ALIAS(Sin, sinq)
@@ -193,19 +225,11 @@ DEFINE_SIMPLE_ALIAS(Trunc, truncq)
 DEFINE_SIMPLE_ALIAS(Y0, y0q)
 DEFINE_SIMPLE_ALIAS(Y1, y1q)
 DEFINE_SIMPLE_ALIAS(Yn, ynq)
-#endif
 
-extern "C" {
-// Declarations of the entry points that might be referenced
-// within the Float128Math library itself.
-// Note that not all of these entry points are actually
-// defined in this library. Some of them are used just
-// as template parameters to call the corresponding callee directly.
-CppTypeFor<TypeCategory::Real, 16> RTDECL(AbsF128)(
-    CppTypeFor<TypeCategory::Real, 16> x);
-CppTypeFor<TypeCategory::Real, 16> RTDECL(SqrtF128)(
-    CppTypeFor<TypeCategory::Real, 16> x);
-} // extern "C"
+// Use cmath INFINITY/NAN definition. Rely on C implicit conversions.
+#define F128_RT_INFINITY (INFINITY)
+#define F128_RT_QNAN (NAN)
+#endif
 
 } // namespace Fortran::runtime
 
diff --git a/flang/runtime/Float128Math/mod-real.cpp b/flang/runtime/Float128Math/mod-real.cpp
new file mode 100644
index 00000000000000..42e6ce76e2fa1b
--- /dev/null
+++ b/flang/runtime/Float128Math/mod-real.cpp
@@ -0,0 +1,24 @@
+//===-- runtime/Float128Math/mod-real.cpp ---------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// MOD (16.9.135)
+F128Type RTDEF(ModReal16)(
+    F128Type x, F128Type p, const char *sourceFile, int sourceLine) {
+  return RealMod<false>(x, p, sourceFile, sourceLine);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/modulo-real.cpp b/flang/runtime/Float128Math/modulo-real.cpp
new file mode 100644
index 00000000000000..13000aba8c8323
--- /dev/null
+++ b/flang/runtime/Float128Math/modulo-real.cpp
@@ -0,0 +1,24 @@
+//===-- runtime/Float128Math/modulo-real.cpp ------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// MODULO (16.9.136)
+F128Type RTDEF(ModuloReal16)(
+    F128Type x, F128Type p, const char *sourceFile, int sourceLine) {
+  return RealMod<true>(x, p, sourceFile, sourceLine);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/nearest.cpp b/flang/runtime/Float128Math/nearest.cpp
new file mode 100644
index 00000000000000..148ac4ef839160
--- /dev/null
+++ b/flang/runtime/Float128Math/nearest.cpp
@@ -0,0 +1,23 @@
+//===-- runtime/Float128Math/nearest.cpp ----------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+CppTypeFor<TypeCategory::Real, 16> RTDEF(Nearest16)(
+    CppTypeFor<TypeCategory::Real, 16> x, bool positive) {
+  return Nextafter<true>::invoke(
+      x, positive ? F128_RT_INFINITY : -F128_RT_INFINITY);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/norm2.cpp b/flang/runtime/Float128Math/norm2.cpp
index 17453bd2d6cbd7..15c482f7f007ce 100644
--- a/flang/runtime/Float128Math/norm2.cpp
+++ b/flang/runtime/Float128Math/norm2.cpp
@@ -7,39 +7,17 @@
 //===----------------------------------------------------------------------===//
 
 #include "math-entries.h"
+#include "numeric-template-specs.h"
 #include "reduction-templates.h"
-#include <cmath>
-
-#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
-
-namespace {
-using namespace Fortran::runtime;
-
-using AccumType = Norm2AccumType<16>;
-
-struct ABSTy {
-  static AccumType compute(AccumType x) {
-    return Sqrt<RTNAME(AbsF128)>::invoke(x);
-  }
-};
-
-struct SQRTTy {
-  static AccumType compute(AccumType x) {
-    return Sqrt<RTNAME(SqrtF128)>::invoke(x);
-  }
-};
-
-using Float128Norm2Accumulator = Norm2Accumulator<16, ABSTy, SQRTTy>;
-} // namespace
 
 namespace Fortran::runtime {
 extern "C" {
 
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(Norm2_16)(
     const Descriptor &x, const char *source, int line, int dim) {
-  auto accumulator{::Float128Norm2Accumulator(x)};
   return GetTotalReduction<TypeCategory::Real, 16>(
-      x, source, line, dim, nullptr, accumulator, "NORM2");
+      x, source, line, dim, nullptr, Norm2Accumulator<16>{x}, "NORM2");
 }
 
 void RTDEF(Norm2DimReal16)(Descriptor &result, const Descriptor &x, int dim,
@@ -49,11 +27,9 @@ void RTDEF(Norm2DimReal16)(Descriptor &result, const Descriptor &x, int dim,
   RUNTIME_CHECK(terminator, type);
   RUNTIME_CHECK(
       terminator, type->first == TypeCategory::Real && type->second == 16);
-  DoMaxMinNorm2<TypeCategory::Real, 16, ::Float128Norm2Accumulator>(
-      result, x, dim, nullptr, "NORM2", terminator);
+  Norm2Helper<16>{}(result, x, dim, nullptr, terminator);
 }
+#endif
 
 } // extern "C"
 } // namespace Fortran::runtime
-
-#endif
diff --git a/flang/runtime/Float128Math/numeric-template-specs.h b/flang/runtime/Float128Math/numeric-template-specs.h
new file mode 100644
index 00000000000000..a0a77230c3e9eb
--- /dev/null
+++ b/flang/runtime/Float128Math/numeric-template-specs.h
@@ -0,0 +1,55 @@
+//===-- runtime/Float128Math/numeric-template-specs.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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef FORTRAN_RUNTIME_FLOAT128MATH_NUMERIC_TEMPLATE_SPECS_H_
+#define FORTRAN_RUNTIME_FLOAT128MATH_NUMERIC_TEMPLATE_SPECS_H_
+
+#include "math-entries.h"
+#include "numeric-templates.h"
+
+namespace Fortran::runtime {
+using F128Type = CppTypeFor<TypeCategory::Real, 16>;
+
+template <> struct ABSTy<F128Type> {
+  static F128Type compute(F128Type x) { return Abs<true>::invoke(x); }
+};
+
+template <> struct FREXPTy<F128Type> {
+  static F128Type compute(F128Type x, int *e) {
+    return Frexp<true>::invoke(x, e);
+  }
+};
+
+template <> struct ILOGBTy<F128Type> {
+  static int compute(F128Type x) { return Ilogb<true>::invoke(x); }
+};
+
+template <> struct ISINFTy<F128Type> {
+  static bool compute(F128Type x) { return Isinf<true>::invoke(x); }
+};
+
+template <> struct ISNANTy<F128Type> {
+  static bool compute(F128Type x) { return Isnan<true>::invoke(x); }
+};
+
+template <> struct LDEXPTy<F128Type> {
+  template <typename ET> static F128Type compute(F128Type x, ET p) {
+    return Ldexp<true>::invoke(x, p);
+  }
+};
+
+template <> struct QNANTy<F128Type> {
+  static F128Type compute() { return F128_RT_QNAN; }
+};
+
+template <> struct SQRTTy<F128Type> {
+  static F128Type compute(F128Type x) { return Sqrt<true>::invoke(x); }
+};
+
+} // namespace Fortran::runtime
+#endif // FORTRAN_RUNTIME_FLOAT128MATH_NUMERIC_TEMPLATE_SPECS_H_
diff --git a/flang/runtime/Float128Math/pow.cpp b/flang/runtime/Float128Math/pow.cpp
index 02958a890e5221..7a48828ee3e765 100644
--- a/flang/runtime/Float128Math/pow.cpp
+++ b/flang/runtime/Float128Math/pow.cpp
@@ -15,7 +15,7 @@ extern "C" {
 CppTypeFor<TypeCategory::Real, 16> RTDEF(PowF128)(
     CppTypeFor<TypeCategory::Real, 16> x,
     CppTypeFor<TypeCategory::Real, 16> y) {
-  return Pow<RTNAME(PowF128)>::invoke(x, y);
+  return Pow<true>::invoke(x, y);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/round.cpp b/flang/runtime/Float128Math/round.cpp
index 43ab57768cb77a..6420c1bc9cd25d 100644
--- a/flang/runtime/Float128Math/round.cpp
+++ b/flang/runtime/Float128Math/round.cpp
@@ -18,7 +18,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(RoundF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Round<RTNAME(RoundF128)>::invoke(x);
+  return Round<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/rrspacing.cpp b/flang/runtime/Float128Math/rrspacing.cpp
new file mode 100644
index 00000000000000..feddac418eec39
--- /dev/null
+++ b/flang/runtime/Float128Math/rrspacing.cpp
@@ -0,0 +1,21 @@
+//===-- runtime/Float128Math/rrspacing.cpp --------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// FRACTION (16.9.80)
+F128Type RTDEF(RRSpacing16)(F128Type x) { return RRSpacing<113>(x); }
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/scale.cpp b/flang/runtime/Float128Math/scale.cpp
new file mode 100644
index 00000000000000..0be958bd9f2a72
--- /dev/null
+++ b/flang/runtime/Float128Math/scale.cpp
@@ -0,0 +1,28 @@
+//===-- runtime/Float128Math/scale.cpp ------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+#include <limits>
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+F128Type RTDEF(Scale16)(F128Type x, std::int64_t p) {
+  auto ip{static_cast<int>(p)};
+  if (ip != p) {
+    ip = p < 0 ? std::numeric_limits<int>::min()
+               : std::numeric_limits<int>::max();
+  }
+  return LDEXPTy<F128Type>::compute(x, ip);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/set-exponent.cpp b/flang/runtime/Float128Math/set-exponent.cpp
new file mode 100644
index 00000000000000..99c34af7962b9a
--- /dev/null
+++ b/flang/runtime/Float128Math/set-exponent.cpp
@@ -0,0 +1,23 @@
+//===-- runtime/Float128Math/set-exponent.cpp -----------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// SET_EXPONENT (16.9.171)
+F128Type RTDEF(SetExponent16)(F128Type x, std::int64_t p) {
+  return SetExponent(x, p);
+}
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/sin.cpp b/flang/runtime/Float128Math/sin.cpp
index 013eb9d119a6a3..8ebc3f9881586e 100644
--- a/flang/runtime/Float128Math/sin.cpp
+++ b/flang/runtime/Float128Math/sin.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(SinF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Sin<RTNAME(SinF128)>::invoke(x);
+  return Sin<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/sinh.cpp b/flang/runtime/Float128Math/sinh.cpp
index 9c907041fd7eb4..aa716a3e51ef5a 100644
--- a/flang/runtime/Float128Math/sinh.cpp
+++ b/flang/runtime/Float128Math/sinh.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(SinhF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Sinh<RTNAME(SinhF128)>::invoke(x);
+  return Sinh<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/spacing.cpp b/flang/runtime/Float128Math/spacing.cpp
new file mode 100644
index 00000000000000..a86c0b30e567ab
--- /dev/null
+++ b/flang/runtime/Float128Math/spacing.cpp
@@ -0,0 +1,21 @@
+//===-- runtime/Float128Math/spacing.cpp ----------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "math-entries.h"
+#include "numeric-template-specs.h"
+
+namespace Fortran::runtime {
+extern "C" {
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+// SPACING (16.9.180)
+F128Type RTDEF(Spacing16)(F128Type x) { return Spacing<113>(x); }
+#endif
+
+} // extern "C"
+} // namespace Fortran::runtime
diff --git a/flang/runtime/Float128Math/sqrt.cpp b/flang/runtime/Float128Math/sqrt.cpp
index aafbd850ca973a..83165a4c623191 100644
--- a/flang/runtime/Float128Math/sqrt.cpp
+++ b/flang/runtime/Float128Math/sqrt.cpp
@@ -7,15 +7,13 @@
 //===----------------------------------------------------------------------===//
 
 #include "math-entries.h"
+#include "numeric-template-specs.h"
 
 namespace Fortran::runtime {
 extern "C" {
 
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
-CppTypeFor<TypeCategory::Real, 16> RTDEF(SqrtF128)(
-    CppTypeFor<TypeCategory::Real, 16> x) {
-  return Sqrt<RTNAME(SqrtF128)>::invoke(x);
-}
+F128Type RTDEF(SqrtF128)(F128Type x) { return SQRTTy<F128Type>::compute(x); }
 #endif
 
 } // extern "C"
diff --git a/flang/runtime/Float128Math/tan.cpp b/flang/runtime/Float128Math/tan.cpp
index 01d3c7bdd2e85d..8f4b723ca977bd 100644
--- a/flang/runtime/Float128Math/tan.cpp
+++ b/flang/runtime/Float128Math/tan.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(TanF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Tan<RTNAME(TanF128)>::invoke(x);
+  return Tan<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/tanh.cpp b/flang/runtime/Float128Math/tanh.cpp
index fedc1a4120caf5..b43a89520b6797 100644
--- a/flang/runtime/Float128Math/tanh.cpp
+++ b/flang/runtime/Float128Math/tanh.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(TanhF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Tanh<RTNAME(TanhF128)>::invoke(x);
+  return Tanh<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/tgamma.cpp b/flang/runtime/Float128Math/tgamma.cpp
index 329defff38cf91..93f97800bdc966 100644
--- a/flang/runtime/Float128Math/tgamma.cpp
+++ b/flang/runtime/Float128Math/tgamma.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(TgammaF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Tgamma<RTNAME(TgammaF128)>::invoke(x);
+  return Tgamma<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/trunc.cpp b/flang/runtime/Float128Math/trunc.cpp
index 3cab219ce31c2d..ca15a739c030e8 100644
--- a/flang/runtime/Float128Math/trunc.cpp
+++ b/flang/runtime/Float128Math/trunc.cpp
@@ -18,7 +18,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(TruncF128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Trunc<RTNAME(TruncF128)>::invoke(x);
+  return Trunc<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/y0.cpp b/flang/runtime/Float128Math/y0.cpp
index f3e2ee454aeab5..d6f39aac1053a8 100644
--- a/flang/runtime/Float128Math/y0.cpp
+++ b/flang/runtime/Float128Math/y0.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(Y0F128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Y0<RTNAME(Y0F128)>::invoke(x);
+  return Y0<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/y1.cpp b/flang/runtime/Float128Math/y1.cpp
index c117bbcb2b5a86..477d36a9ea3c66 100644
--- a/flang/runtime/Float128Math/y1.cpp
+++ b/flang/runtime/Float128Math/y1.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(Y1F128)(
     CppTypeFor<TypeCategory::Real, 16> x) {
-  return Y1<RTNAME(Y1F128)>::invoke(x);
+  return Y1<true>::invoke(x);
 }
 #endif
 
diff --git a/flang/runtime/Float128Math/yn.cpp b/flang/runtime/Float128Math/yn.cpp
index 237bc2866a0d5b..3a040cc8858970 100644
--- a/flang/runtime/Float128Math/yn.cpp
+++ b/flang/runtime/Float128Math/yn.cpp
@@ -14,7 +14,7 @@ extern "C" {
 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 CppTypeFor<TypeCategory::Real, 16> RTDEF(YnF128)(
     int n, CppTypeFor<TypeCategory::Real, 16> x) {
-  return Yn<RTNAME(YnF128)>::invoke(n, x);
+  return Yn<true>::invoke(n, x);
 }
 #endif
 
diff --git a/flang/runtime/extrema.cpp b/flang/runtime/extrema.cpp
index fc2b4e165cb269..61afb0458430db 100644
--- a/flang/runtime/extrema.cpp
+++ b/flang/runtime/extrema.cpp
@@ -424,62 +424,6 @@ RT_EXT_API_GROUP_END
 
 // MAXVAL and MINVAL
 
-template <TypeCategory CAT, int KIND, bool IS_MAXVAL, typename Enable = void>
-struct MaxOrMinIdentity {
-  using Type = CppTypeFor<CAT, KIND>;
-  static constexpr RT_API_ATTRS Type Value() {
-    return IS_MAXVAL ? std::numeric_limits<Type>::lowest()
-                     : std::numeric_limits<Type>::max();
-  }
-};
-
-// std::numeric_limits<> may not know int128_t
-template <bool IS_MAXVAL>
-struct MaxOrMinIdentity<TypeCategory::Integer, 16, IS_MAXVAL> {
-  using Type = CppTypeFor<TypeCategory::Integer, 16>;
-  static constexpr RT_API_ATTRS Type Value() {
-    return IS_MAXVAL ? Type{1} << 127 : ~Type{0} >> 1;
-  }
-};
-
-#if HAS_FLOAT128
-// std::numeric_limits<> may not support __float128.
-//
-// Usage of GCC quadmath.h's FLT128_MAX is complicated by the fact that
-// even GCC complains about 'Q' literal suffix under -Wpedantic.
-// We just recreate FLT128_MAX ourselves.
-//
-// This specialization must engage only when
-// CppTypeFor<TypeCategory::Real, 16> is __float128.
-template <bool IS_MAXVAL>
-struct MaxOrMinIdentity<TypeCategory::Real, 16, IS_MAXVAL,
-    typename std::enable_if_t<
-        std::is_same_v<CppTypeFor<TypeCategory::Real, 16>, __float128>>> {
-  using Type = __float128;
-  static RT_API_ATTRS Type Value() {
-    // Create a buffer to store binary representation of __float128 constant.
-    constexpr std::size_t alignment =
-        std::max(alignof(Type), alignof(std::uint64_t));
-    alignas(alignment) char data[sizeof(Type)];
-
-    // First, verify that our interpretation of __float128 format is correct,
-    // e.g. by checking at least one known constant.
-    *reinterpret_cast<Type *>(data) = Type(1.0);
-    if (*reinterpret_cast<std::uint64_t *>(data) != 0 ||
-        *(reinterpret_cast<std::uint64_t *>(data) + 1) != 0x3FFF000000000000) {
-      Terminator terminator{__FILE__, __LINE__};
-      terminator.Crash("not yet implemented: no full support for __float128");
-    }
-
-    // Recreate FLT128_MAX.
-    *reinterpret_cast<std::uint64_t *>(data) = 0xFFFFFFFFFFFFFFFF;
-    *(reinterpret_cast<std::uint64_t *>(data) + 1) = 0x7FFEFFFFFFFFFFFF;
-    Type max = *reinterpret_cast<Type *>(data);
-    return IS_MAXVAL ? -max : max;
-  }
-};
-#endif // HAS_FLOAT128
-
 template <TypeCategory CAT, int KIND, bool IS_MAXVAL>
 class NumericExtremumAccumulator {
 public:
@@ -773,42 +717,25 @@ RT_EXT_API_GROUP_END
 
 // NORM2
 
-template <int KIND> struct Norm2Helper {
-  RT_API_ATTRS void operator()(Descriptor &result, const Descriptor &x, int dim,
-      const Descriptor *mask, Terminator &terminator) const {
-    DoMaxMinNorm2<TypeCategory::Real, KIND,
-        typename Norm2AccumulatorGetter<KIND>::Type>(
-        result, x, dim, mask, "NORM2", terminator);
-  }
-};
-
 extern "C" {
 RT_EXT_API_GROUP_BEGIN
 
 // TODO: REAL(2 & 3)
 CppTypeFor<TypeCategory::Real, 4> RTDEF(Norm2_4)(
     const Descriptor &x, const char *source, int line, int dim) {
-  return GetTotalReduction<TypeCategory::Real, 4>(x, source, line, dim, nullptr,
-      Norm2AccumulatorGetter<4>::create(x), "NORM2");
+  return GetTotalReduction<TypeCategory::Real, 4>(
+      x, source, line, dim, nullptr, Norm2Accumulator<4>{x}, "NORM2");
 }
 CppTypeFor<TypeCategory::Real, 8> RTDEF(Norm2_8)(
     const Descriptor &x, const char *source, int line, int dim) {
-  return GetTotalReduction<TypeCategory::Real, 8>(x, source, line, dim, nullptr,
-      Norm2AccumulatorGetter<8>::create(x), "NORM2");
+  return GetTotalReduction<TypeCategory::Real, 8>(
+      x, source, line, dim, nullptr, Norm2Accumulator<8>{x}, "NORM2");
 }
 #if LDBL_MANT_DIG == 64
 CppTypeFor<TypeCategory::Real, 10> RTDEF(Norm2_10)(
     const Descriptor &x, const char *source, int line, int dim) {
-  return GetTotalReduction<TypeCategory::Real, 10>(x, source, line, dim,
-      nullptr, Norm2AccumulatorGetter<10>::create(x), "NORM2");
-}
-#endif
-#if LDBL_MANT_DIG == 113
-// The __float128 implementation resides in FortranFloat128Math library.
-CppTypeFor<TypeCategory::Real, 16> RTDEF(Norm2_16)(
-    const Descriptor &x, const char *source, int line, int dim) {
-  return GetTotalReduction<TypeCategory::Real, 16>(x, source, line, dim,
-      nullptr, Norm2AccumulatorGetter<16>::create(x), "NORM2");
+  return GetTotalReduction<TypeCategory::Real, 10>(
+      x, source, line, dim, nullptr, Norm2Accumulator<10>{x}, "NORM2");
 }
 #endif
 
diff --git a/flang/runtime/numeric-templates.h b/flang/runtime/numeric-templates.h
new file mode 100644
index 00000000000000..b16440dbc2241a
--- /dev/null
+++ b/flang/runtime/numeric-templates.h
@@ -0,0 +1,339 @@
+//===-- runtime/numeric-templates.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
+//
+//===----------------------------------------------------------------------===//
+
+// Generic class and function templates used for implementing
+// various numeric intrinsics (EXPONENT, FRACTION, etc.).
+//
+// This header file also defines generic templates for "basic"
+// math operations like abs, isnan, etc. The Float128Math
+// library provides specializations for these templates
+// for the data type corresponding to CppTypeFor<TypeCategory::Real, 16>
+// on the target.
+
+#ifndef FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_
+#define FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_
+
+#include "terminator.h"
+#include "tools.h"
+#include "flang/Common/float128.h"
+#include <cstdint>
+#include <limits>
+
+namespace Fortran::runtime {
+
+// MAX/MIN/LOWEST values for different data types.
+
+// MaxOrMinIdentity returns MAX or LOWEST value of the given type.
+template <TypeCategory CAT, int KIND, bool IS_MAXVAL, typename Enable = void>
+struct MaxOrMinIdentity {
+  using Type = CppTypeFor<CAT, KIND>;
+  static constexpr RT_API_ATTRS Type Value() {
+    return IS_MAXVAL ? std::numeric_limits<Type>::lowest()
+                     : std::numeric_limits<Type>::max();
+  }
+};
+
+// std::numeric_limits<> may not know int128_t
+template <bool IS_MAXVAL>
+struct MaxOrMinIdentity<TypeCategory::Integer, 16, IS_MAXVAL> {
+  using Type = CppTypeFor<TypeCategory::Integer, 16>;
+  static constexpr RT_API_ATTRS Type Value() {
+    return IS_MAXVAL ? Type{1} << 127 : ~Type{0} >> 1;
+  }
+};
+
+#if HAS_FLOAT128
+// std::numeric_limits<> may not support __float128.
+//
+// Usage of GCC quadmath.h's FLT128_MAX is complicated by the fact that
+// even GCC complains about 'Q' literal suffix under -Wpedantic.
+// We just recreate FLT128_MAX ourselves.
+//
+// This specialization must engage only when
+// CppTypeFor<TypeCategory::Real, 16> is __float128.
+template <bool IS_MAXVAL>
+struct MaxOrMinIdentity<TypeCategory::Real, 16, IS_MAXVAL,
+    typename std::enable_if_t<
+        std::is_same_v<CppTypeFor<TypeCategory::Real, 16>, __float128>>> {
+  using Type = __float128;
+  static RT_API_ATTRS Type Value() {
+    // Create a buffer to store binary representation of __float128 constant.
+    constexpr std::size_t alignment =
+        std::max(alignof(Type), alignof(std::uint64_t));
+    alignas(alignment) char data[sizeof(Type)];
+
+    // First, verify that our interpretation of __float128 format is correct,
+    // e.g. by checking at least one known constant.
+    *reinterpret_cast<Type *>(data) = Type(1.0);
+    if (*reinterpret_cast<std::uint64_t *>(data) != 0 ||
+        *(reinterpret_cast<std::uint64_t *>(data) + 1) != 0x3FFF000000000000) {
+      Terminator terminator{__FILE__, __LINE__};
+      terminator.Crash("not yet implemented: no full support for __float128");
+    }
+
+    // Recreate FLT128_MAX.
+    *reinterpret_cast<std::uint64_t *>(data) = 0xFFFFFFFFFFFFFFFF;
+    *(reinterpret_cast<std::uint64_t *>(data) + 1) = 0x7FFEFFFFFFFFFFFF;
+    Type max = *reinterpret_cast<Type *>(data);
+    return IS_MAXVAL ? -max : max;
+  }
+};
+#endif // HAS_FLOAT128
+
+// Minimum finite representable value.
+// For floating-point types, returns minimum positive normalized value.
+template <typename T> struct MinValue {
+  static RT_API_ATTRS T get() { return std::numeric_limits<T>::min(); }
+};
+
+#if HAS_FLOAT128
+template <> struct MinValue<CppTypeFor<TypeCategory::Real, 16>> {
+  using Type = CppTypeFor<TypeCategory::Real, 16>;
+  static RT_API_ATTRS Type get() {
+    // Create a buffer to store binary representation of __float128 constant.
+    constexpr std::size_t alignment =
+        std::max(alignof(Type), alignof(std::uint64_t));
+    alignas(alignment) char data[sizeof(Type)];
+
+    // First, verify that our interpretation of __float128 format is correct,
+    // e.g. by checking at least one known constant.
+    *reinterpret_cast<Type *>(data) = Type(1.0);
+    if (*reinterpret_cast<std::uint64_t *>(data) != 0 ||
+        *(reinterpret_cast<std::uint64_t *>(data) + 1) != 0x3FFF000000000000) {
+      Terminator terminator{__FILE__, __LINE__};
+      terminator.Crash("not yet implemented: no full support for __float128");
+    }
+
+    // Recreate FLT128_MIN.
+    *reinterpret_cast<std::uint64_t *>(data) = 0;
+    *(reinterpret_cast<std::uint64_t *>(data) + 1) = 0x1000000000000;
+    return *reinterpret_cast<Type *>(data);
+  }
+};
+#endif // HAS_FLOAT128
+
+template <typename T> struct ABSTy {
+  static constexpr RT_API_ATTRS T compute(T x) { return std::abs(x); }
+};
+
+template <typename T> struct FREXPTy {
+  static constexpr RT_API_ATTRS T compute(T x, int *e) {
+    return std::frexp(x, e);
+  }
+};
+
+template <typename T> struct ILOGBTy {
+  static constexpr RT_API_ATTRS int compute(T x) { return std::ilogb(x); }
+};
+
+template <typename T> struct ISINFTy {
+  static constexpr RT_API_ATTRS bool compute(T x) { return std::isinf(x); }
+};
+
+template <typename T> struct ISNANTy {
+  static constexpr RT_API_ATTRS bool compute(T x) { return std::isnan(x); }
+};
+
+template <typename T> struct LDEXPTy {
+  template <typename ET> static constexpr RT_API_ATTRS T compute(T x, ET e) {
+    return std::ldexp(x, e);
+  }
+};
+
+template <typename T> struct MAXTy {
+  static constexpr RT_API_ATTRS T compute() {
+    return std::numeric_limits<T>::max();
+  }
+};
+
+#if LDBL_MANT_DIG == 113 || HAS_FLOAT128
+template <> struct MAXTy<CppTypeFor<TypeCategory::Real, 16>> {
+  static CppTypeFor<TypeCategory::Real, 16> compute() {
+    return MaxOrMinIdentity<TypeCategory::Real, 16, true>::Value();
+  }
+};
+#endif
+
+template <typename T> struct MINTy {
+  static constexpr RT_API_ATTRS T compute() { return MinValue<T>::get(); }
+};
+
+template <typename T> struct QNANTy {
+  static constexpr RT_API_ATTRS T compute() {
+    return std::numeric_limits<T>::quiet_NaN();
+  }
+};
+
+template <typename T> struct SQRTTy {
+  static constexpr RT_API_ATTRS T compute(T x) { return std::sqrt(x); }
+};
+
+// EXPONENT (16.9.75)
+template <typename RESULT, typename ARG>
+inline RT_API_ATTRS RESULT Exponent(ARG x) {
+  if (ISINFTy<ARG>::compute(x) || ISNANTy<ARG>::compute(x)) {
+    return MAXTy<RESULT>::compute(); // +/-Inf, NaN -> HUGE(0)
+  } else if (x == 0) {
+    return 0; // 0 -> 0
+  } else {
+    return ILOGBTy<ARG>::compute(x) + 1;
+  }
+}
+
+// Suppress the warnings about calling __host__-only std::frexp,
+// defined in C++ STD header files, from __device__ code.
+RT_DIAG_PUSH
+RT_DIAG_DISABLE_CALL_HOST_FROM_DEVICE_WARN
+
+// FRACTION (16.9.80)
+template <typename T> inline RT_API_ATTRS T Fraction(T x) {
+  if (ISNANTy<T>::compute(x)) {
+    return x; // NaN -> same NaN
+  } else if (ISINFTy<T>::compute(x)) {
+    return QNANTy<T>::compute(); // +/-Inf -> NaN
+  } else if (x == 0) {
+    return x; // 0 -> same 0
+  } else {
+    int ignoredExp;
+    return FREXPTy<T>::compute(x, &ignoredExp);
+  }
+}
+
+RT_DIAG_POP
+
+// SET_EXPONENT (16.9.171)
+template <typename T> inline RT_API_ATTRS T SetExponent(T x, std::int64_t p) {
+  if (ISNANTy<T>::compute(x)) {
+    return x; // NaN -> same NaN
+  } else if (ISINFTy<T>::compute(x)) {
+    return QNANTy<T>::compute(); // +/-Inf -> NaN
+  } else if (x == 0) {
+    return x; // return negative zero if x is negative zero
+  } else {
+    int expo{ILOGBTy<T>::compute(x) + 1};
+    auto ip{static_cast<int>(p - expo)};
+    if (ip != p - expo) {
+      ip = p < 0 ? std::numeric_limits<int>::min()
+                 : std::numeric_limits<int>::max();
+    }
+    return LDEXPTy<T>::compute(x, ip); // x*2**(p-e)
+  }
+}
+
+// MOD & MODULO (16.9.135, .136)
+template <bool IS_MODULO, typename T>
+inline RT_API_ATTRS T RealMod(
+    T a, T p, const char *sourceFile, int sourceLine) {
+  if (p == 0) {
+    Terminator{sourceFile, sourceLine}.Crash(
+        IS_MODULO ? "MODULO with P==0" : "MOD with P==0");
+  }
+  if (ISNANTy<T>::compute(a) || ISNANTy<T>::compute(p) ||
+      ISINFTy<T>::compute(a)) {
+    return QNANTy<T>::compute();
+  } else if (ISINFTy<T>::compute(p)) {
+    return a;
+  }
+  T aAbs{ABSTy<T>::compute(a)};
+  T pAbs{ABSTy<T>::compute(p)};
+  if (aAbs <= static_cast<T>(std::numeric_limits<std::int64_t>::max()) &&
+      pAbs <= static_cast<T>(std::numeric_limits<std::int64_t>::max())) {
+    if (auto aInt{static_cast<std::int64_t>(a)}; a == aInt) {
+      if (auto pInt{static_cast<std::int64_t>(p)}; p == pInt) {
+        // Fast exact case for integer operands
+        auto mod{aInt - (aInt / pInt) * pInt};
+        if (IS_MODULO && (aInt > 0) != (pInt > 0)) {
+          mod += pInt;
+        }
+        return static_cast<T>(mod);
+      }
+    }
+  }
+  if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double> ||
+      std::is_same_v<T, long double>) {
+    // std::fmod() semantics on signed operands seems to match
+    // the requirements of MOD().  MODULO() needs adjustment.
+    T result{std::fmod(a, p)};
+    if constexpr (IS_MODULO) {
+      if ((a < 0) != (p < 0)) {
+        if (result == 0.) {
+          result = -result;
+        } else {
+          result += p;
+        }
+      }
+    }
+    return result;
+  } else {
+    // The standard defines MOD(a,p)=a-AINT(a/p)*p and
+    // MODULO(a,p)=a-FLOOR(a/p)*p, but those definitions lose
+    // precision badly due to cancellation when ABS(a) is
+    // much larger than ABS(p).
+    // Insights:
+    //  - MOD(a,p)=MOD(a-n*p,p) when a>0, p>0, integer n>0, and a>=n*p
+    //  - when n is a power of two, n*p is exact
+    //  - as a>=n*p, a-n*p does not round.
+    // So repeatedly reduce a by all n*p in decreasing order of n;
+    // what's left is the desired remainder.  This is basically
+    // the same algorithm as arbitrary precision binary long division,
+    // discarding the quotient.
+    T tmp{aAbs};
+    for (T adj{SetExponent(pAbs, Exponent<int>(aAbs))}; tmp >= pAbs; adj /= 2) {
+      if (tmp >= adj) {
+        tmp -= adj;
+        if (tmp == 0) {
+          break;
+        }
+      }
+    }
+    if (a < 0) {
+      tmp = -tmp;
+    }
+    if constexpr (IS_MODULO) {
+      if ((a < 0) != (p < 0)) {
+        tmp += p;
+      }
+    }
+    return tmp;
+  }
+}
+
+// RRSPACING (16.9.164)
+template <int PREC, typename T> inline RT_API_ATTRS T RRSpacing(T x) {
+  if (ISNANTy<T>::compute(x)) {
+    return x; // NaN -> same NaN
+  } else if (ISINFTy<T>::compute(x)) {
+    return QNANTy<T>::compute(); // +/-Inf -> NaN
+  } else if (x == 0) {
+    return 0; // 0 -> 0
+  } else {
+    return LDEXPTy<T>::compute(
+        ABSTy<T>::compute(x), PREC - (ILOGBTy<T>::compute(x) + 1));
+  }
+}
+
+// SPACING (16.9.180)
+template <int PREC, typename T> inline RT_API_ATTRS T Spacing(T x) {
+  if (ISNANTy<T>::compute(x)) {
+    return x; // NaN -> same NaN
+  } else if (ISINFTy<T>::compute(x)) {
+    return QNANTy<T>::compute(); // +/-Inf -> NaN
+  } else if (x == 0) {
+    // The standard-mandated behavior seems broken, since TINY() can't be
+    // subnormal.
+    return MINTy<T>::compute(); // 0 -> TINY(x)
+  } else {
+    T result{LDEXPTy<T>::compute(
+        static_cast<T>(1.0), ILOGBTy<T>::compute(x) + 1 - PREC)}; // 2**(e-p)
+    return result == 0 ? /*TINY(x)*/ MINTy<T>::compute() : result;
+  }
+}
+
+} // namespace Fortran::runtime
+
+#endif // FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_
diff --git a/flang/runtime/numeric.cpp b/flang/runtime/numeric.cpp
index ede00d69f20e25..abd3e500029fe4 100644
--- a/flang/runtime/numeric.cpp
+++ b/flang/runtime/numeric.cpp
@@ -7,6 +7,7 @@
 //===----------------------------------------------------------------------===//
 
 #include "flang/Runtime/numeric.h"
+#include "numeric-templates.h"
 #include "terminator.h"
 #include "flang/Common/float128.h"
 #include <cfloat>
@@ -68,58 +69,6 @@ inline RT_API_ATTRS RESULT Floor(ARG x) {
   return std::floor(x);
 }
 
-// EXPONENT (16.9.75)
-template <typename RESULT, typename ARG>
-inline RT_API_ATTRS RESULT Exponent(ARG x) {
-  if (std::isinf(x) || std::isnan(x)) {
-    return std::numeric_limits<RESULT>::max(); // +/-Inf, NaN -> HUGE(0)
-  } else if (x == 0) {
-    return 0; // 0 -> 0
-  } else {
-    return std::ilogb(x) + 1;
-  }
-}
-
-// Suppress the warnings about calling __host__-only std::frexp,
-// defined in C++ STD header files, from __device__ code.
-RT_DIAG_PUSH
-RT_DIAG_DISABLE_CALL_HOST_FROM_DEVICE_WARN
-
-// FRACTION (16.9.80)
-template <typename T> inline RT_API_ATTRS T Fraction(T x) {
-  if (std::isnan(x)) {
-    return x; // NaN -> same NaN
-  } else if (std::isinf(x)) {
-    return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN
-  } else if (x == 0) {
-    return x; // 0 -> same 0
-  } else {
-    int ignoredExp;
-    return std::frexp(x, &ignoredExp);
-  }
-}
-
-RT_DIAG_POP
-
-// SET_EXPONENT (16.9.171)
-template <typename T> inline RT_API_ATTRS T SetExponent(T x, std::int64_t p) {
-  if (std::isnan(x)) {
-    return x; // NaN -> same NaN
-  } else if (std::isinf(x)) {
-    return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN
-  } else if (x == 0) {
-    return x; // return negative zero if x is negative zero
-  } else {
-    int expo{std::ilogb(x) + 1};
-    auto ip{static_cast<int>(p - expo)};
-    if (ip != p - expo) {
-      ip = p < 0 ? std::numeric_limits<int>::min()
-                 : std::numeric_limits<int>::max();
-    }
-    return std::ldexp(x, ip); // x*2**(p-e)
-  }
-}
-
 // MOD & MODULO (16.9.135, .136)
 template <bool IS_MODULO, typename T>
 inline RT_API_ATTRS T IntMod(T x, T p, const char *sourceFile, int sourceLine) {
@@ -133,94 +82,6 @@ inline RT_API_ATTRS T IntMod(T x, T p, const char *sourceFile, int sourceLine) {
   }
   return mod;
 }
-template <bool IS_MODULO, typename T>
-inline RT_API_ATTRS T RealMod(
-    T a, T p, const char *sourceFile, int sourceLine) {
-  if (p == 0) {
-    Terminator{sourceFile, sourceLine}.Crash(
-        IS_MODULO ? "MODULO with P==0" : "MOD with P==0");
-  }
-  if (std::isnan(a) || std::isnan(p) || std::isinf(a)) {
-    return std::numeric_limits<T>::quiet_NaN();
-  } else if (std::isinf(p)) {
-    return a;
-  }
-  T aAbs{std::abs(a)};
-  T pAbs{std::abs(p)};
-  if (aAbs <= static_cast<T>(std::numeric_limits<std::int64_t>::max()) &&
-      pAbs <= static_cast<T>(std::numeric_limits<std::int64_t>::max())) {
-    if (auto aInt{static_cast<std::int64_t>(a)}; a == aInt) {
-      if (auto pInt{static_cast<std::int64_t>(p)}; p == pInt) {
-        // Fast exact case for integer operands
-        auto mod{aInt - (aInt / pInt) * pInt};
-        if (IS_MODULO && (aInt > 0) != (pInt > 0)) {
-          mod += pInt;
-        }
-        return static_cast<T>(mod);
-      }
-    }
-  }
-  if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double> ||
-      std::is_same_v<T, long double>) {
-    // std::fmod() semantics on signed operands seems to match
-    // the requirements of MOD().  MODULO() needs adjustment.
-    T result{std::fmod(a, p)};
-    if constexpr (IS_MODULO) {
-      if ((a < 0) != (p < 0)) {
-        if (result == 0.) {
-          result = -result;
-        } else {
-          result += p;
-        }
-      }
-    }
-    return result;
-  } else {
-    // The standard defines MOD(a,p)=a-AINT(a/p)*p and
-    // MODULO(a,p)=a-FLOOR(a/p)*p, but those definitions lose
-    // precision badly due to cancellation when ABS(a) is
-    // much larger than ABS(p).
-    // Insights:
-    //  - MOD(a,p)=MOD(a-n*p,p) when a>0, p>0, integer n>0, and a>=n*p
-    //  - when n is a power of two, n*p is exact
-    //  - as a>=n*p, a-n*p does not round.
-    // So repeatedly reduce a by all n*p in decreasing order of n;
-    // what's left is the desired remainder.  This is basically
-    // the same algorithm as arbitrary precision binary long division,
-    // discarding the quotient.
-    T tmp{aAbs};
-    for (T adj{SetExponent(pAbs, Exponent<int>(aAbs))}; tmp >= pAbs; adj /= 2) {
-      if (tmp >= adj) {
-        tmp -= adj;
-        if (tmp == 0) {
-          break;
-        }
-      }
-    }
-    if (a < 0) {
-      tmp = -tmp;
-    }
-    if constexpr (IS_MODULO) {
-      if ((a < 0) != (p < 0)) {
-        tmp += p;
-      }
-    }
-    return tmp;
-  }
-}
-
-// RRSPACING (16.9.164)
-template <int PREC, typename T> inline RT_API_ATTRS T RRSpacing(T x) {
-  if (std::isnan(x)) {
-    return x; // NaN -> same NaN
-  } else if (std::isinf(x)) {
-    return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN
-  } else if (x == 0) {
-    return 0; // 0 -> 0
-  } else {
-    return std::ldexp(std::abs(x), PREC - (std::ilogb(x) + 1));
-  }
-}
 
 // SCALE (16.9.166)
 template <typename T> inline RT_API_ATTRS T Scale(T x, std::int64_t p) {
@@ -229,7 +90,7 @@ template <typename T> inline RT_API_ATTRS T Scale(T x, std::int64_t p) {
     ip = p < 0 ? std::numeric_limits<int>::min()
                : std::numeric_limits<int>::max();
   }
-  return std::ldexp(x, p); // x*2**p
+  return std::ldexp(x, ip); // x*2**p
 }
 
 // SELECTED_INT_KIND (16.9.169)
@@ -300,23 +161,6 @@ inline RT_API_ATTRS CppTypeFor<TypeCategory::Integer, 4> SelectedRealKind(
   return error ? error : kind;
 }
 
-// SPACING (16.9.180)
-template <int PREC, typename T> inline RT_API_ATTRS T Spacing(T x) {
-  if (std::isnan(x)) {
-    return x; // NaN -> same NaN
-  } else if (std::isinf(x)) {
-    return std::numeric_limits<T>::quiet_NaN(); // +/-Inf -> NaN
-  } else if (x == 0) {
-    // The standard-mandated behavior seems broken, since TINY() can't be
-    // subnormal.
-    return std::numeric_limits<T>::min(); // 0 -> TINY(x)
-  } else {
-    T result{
-        std::ldexp(static_cast<T>(1.0), std::ilogb(x) + 1 - PREC)}; // 2**(e-p)
-    return result == 0 ? /*TINY(x)*/ std::numeric_limits<T>::min() : result;
-  }
-}
-
 // NEAREST (16.9.139)
 template <int PREC, typename T>
 inline RT_API_ATTRS T Nearest(T x, bool positive) {
@@ -480,15 +324,6 @@ CppTypeFor<TypeCategory::Integer, 8> RTDEF(Exponent10_8)(
     CppTypeFor<TypeCategory::Real, 10> x) {
   return Exponent<CppTypeFor<TypeCategory::Integer, 8>>(x);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Integer, 4> RTDEF(Exponent16_4)(
-    CppTypeFor<TypeCategory::Real, 16> x) {
-  return Exponent<CppTypeFor<TypeCategory::Integer, 4>>(x);
-}
-CppTypeFor<TypeCategory::Integer, 8> RTDEF(Exponent16_8)(
-    CppTypeFor<TypeCategory::Real, 16> x) {
-  return Exponent<CppTypeFor<TypeCategory::Integer, 8>>(x);
-}
 #endif
 
 CppTypeFor<TypeCategory::Integer, 1> RTDEF(Floor4_1)(
@@ -596,11 +431,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(Fraction10)(
     CppTypeFor<TypeCategory::Real, 10> x) {
   return Fraction(x);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(Fraction16)(
-    CppTypeFor<TypeCategory::Real, 16> x) {
-  return Fraction(x);
-}
 #endif
 
 bool RTDEF(IsFinite4)(CppTypeFor<TypeCategory::Real, 4> x) {
@@ -683,12 +513,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(ModReal10)(
     const char *sourceFile, int sourceLine) {
   return RealMod<false>(x, p, sourceFile, sourceLine);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(ModReal16)(
-    CppTypeFor<TypeCategory::Real, 16> x, CppTypeFor<TypeCategory::Real, 16> p,
-    const char *sourceFile, int sourceLine) {
-  return RealMod<false>(x, p, sourceFile, sourceLine);
-}
 #endif
 
 CppTypeFor<TypeCategory::Integer, 1> RTDEF(ModuloInteger1)(
@@ -739,12 +563,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(ModuloReal10)(
     const char *sourceFile, int sourceLine) {
   return RealMod<true>(x, p, sourceFile, sourceLine);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(ModuloReal16)(
-    CppTypeFor<TypeCategory::Real, 16> x, CppTypeFor<TypeCategory::Real, 16> p,
-    const char *sourceFile, int sourceLine) {
-  return RealMod<true>(x, p, sourceFile, sourceLine);
-}
 #endif
 
 CppTypeFor<TypeCategory::Real, 4> RTDEF(Nearest4)(
@@ -760,11 +578,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(Nearest10)(
     CppTypeFor<TypeCategory::Real, 10> x, bool positive) {
   return Nearest<64>(x, positive);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(Nearest16)(
-    CppTypeFor<TypeCategory::Real, 16> x, bool positive) {
-  return Nearest<113>(x, positive);
-}
 #endif
 
 CppTypeFor<TypeCategory::Integer, 1> RTDEF(Nint4_1)(
@@ -872,11 +685,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(RRSpacing10)(
     CppTypeFor<TypeCategory::Real, 10> x) {
   return RRSpacing<64>(x);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(RRSpacing16)(
-    CppTypeFor<TypeCategory::Real, 16> x) {
-  return RRSpacing<113>(x);
-}
 #endif
 
 CppTypeFor<TypeCategory::Real, 4> RTDEF(SetExponent4)(
@@ -892,11 +700,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(SetExponent10)(
     CppTypeFor<TypeCategory::Real, 10> x, std::int64_t p) {
   return SetExponent(x, p);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(SetExponent16)(
-    CppTypeFor<TypeCategory::Real, 16> x, std::int64_t p) {
-  return SetExponent(x, p);
-}
 #endif
 
 CppTypeFor<TypeCategory::Real, 4> RTDEF(Scale4)(
@@ -912,11 +715,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(Scale10)(
     CppTypeFor<TypeCategory::Real, 10> x, std::int64_t p) {
   return Scale(x, p);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(Scale16)(
-    CppTypeFor<TypeCategory::Real, 16> x, std::int64_t p) {
-  return Scale(x, p);
-}
 #endif
 
 // SELECTED_INT_KIND
@@ -971,11 +769,6 @@ CppTypeFor<TypeCategory::Real, 10> RTDEF(Spacing10)(
     CppTypeFor<TypeCategory::Real, 10> x) {
   return Spacing<64>(x);
 }
-#elif LDBL_MANT_DIG == 113
-CppTypeFor<TypeCategory::Real, 16> RTDEF(Spacing16)(
-    CppTypeFor<TypeCategory::Real, 16> x) {
-  return Spacing<113>(x);
-}
 #endif
 
 CppTypeFor<TypeCategory::Real, 4> RTDEF(FPow4i)(
diff --git a/flang/runtime/reduction-templates.h b/flang/runtime/reduction-templates.h
index 0891bc021ff753..5b793deb2a123d 100644
--- a/flang/runtime/reduction-templates.h
+++ b/flang/runtime/reduction-templates.h
@@ -1,4 +1,4 @@
-//===-- runtime/reduction-templates.h -------------------------------------===//
+//===-- runtime/reduction-templates.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.
@@ -21,6 +21,7 @@
 #ifndef FORTRAN_RUNTIME_REDUCTION_TEMPLATES_H_
 #define FORTRAN_RUNTIME_REDUCTION_TEMPLATES_H_
 
+#include "numeric-templates.h"
 #include "terminator.h"
 #include "tools.h"
 #include "flang/Runtime/cpp-type.h"
@@ -385,7 +386,7 @@ template <int KIND>
 using Norm2AccumType =
     CppTypeFor<TypeCategory::Real, std::clamp(KIND, 8, Norm2LargestLDKind)>;
 
-template <int KIND, typename ABS, typename SQRT> class Norm2Accumulator {
+template <int KIND> class Norm2Accumulator {
 public:
   using Type = CppTypeFor<TypeCategory::Real, KIND>;
   using AccumType = Norm2AccumType<KIND>;
@@ -395,10 +396,10 @@ template <int KIND, typename ABS, typename SQRT> class Norm2Accumulator {
   template <typename A>
   RT_API_ATTRS void GetResult(A *p, int /*zeroBasedDim*/ = -1) const {
     // m * sqrt(1 + sum((others(:)/m)**2))
-    *p = static_cast<Type>(max_ * SQRT::compute(1 + sum_));
+    *p = static_cast<Type>(max_ * SQRTTy<AccumType>::compute(1 + sum_));
   }
   RT_API_ATTRS bool Accumulate(Type x) {
-    auto absX{ABS::compute(static_cast<AccumType>(x))};
+    auto absX{ABSTy<AccumType>::compute(static_cast<AccumType>(x))};
     if (!max_) {
       max_ = absX;
     } else if (absX > max_) {
@@ -424,27 +425,12 @@ template <int KIND, typename ABS, typename SQRT> class Norm2Accumulator {
   AccumType sum_{0}; // sum((others(:)/m)**2)
 };
 
-// Helper class for creating Norm2Accumulator instance
-// based on the given KIND. This helper returns and instance
-// that uses std::abs and std::sqrt for the computations.
-template <int KIND> class Norm2AccumulatorGetter {
-  using AccumType = Norm2AccumType<KIND>;
-
-public:
-  struct ABSTy {
-    static constexpr RT_API_ATTRS AccumType compute(AccumType &&x) {
-      return std::abs(std::forward<AccumType>(x));
-    }
-  };
-  struct SQRTTy {
-    static constexpr RT_API_ATTRS AccumType compute(AccumType &&x) {
-      return std::sqrt(std::forward<AccumType>(x));
-    }
-  };
-
-  using Type = Norm2Accumulator<KIND, ABSTy, SQRTTy>;
-
-  static RT_API_ATTRS Type create(const Descriptor &x) { return Type(x); }
+template <int KIND> struct Norm2Helper {
+  RT_API_ATTRS void operator()(Descriptor &result, const Descriptor &x, int dim,
+      const Descriptor *mask, Terminator &terminator) const {
+    DoMaxMinNorm2<TypeCategory::Real, KIND, Norm2Accumulator<KIND>>(
+        result, x, dim, mask, "NORM2", terminator);
+  }
 };
 
 } // namespace Fortran::runtime