flang/runtime/numeric-templates.h

   1 //===-- runtime/numeric-templates.h -----------------------------*- C++ -*-===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 // Generic class and function templates used for implementing
  10 // various numeric intrinsics (EXPONENT, FRACTION, etc.).
  11 //
  12 // This header file also defines generic templates for "basic"
  13 // math operations like abs, isnan, etc. The Float128Math
  14 // library provides specializations for these templates
  15 // for the data type corresponding to CppTypeFor<TypeCategory::Real, 16>
  16 // on the target.
  17
  18 #ifndef FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_
  19 #define FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_
  20
  21 #include "terminator.h"
  22 #include "tools.h"
  23 #include "flang/Common/api-attrs.h"
  24 #include "flang/Common/float128.h"
  25 #include <cstdint>
  26 #include <limits>
  27
  28 namespace Fortran::runtime {
  29
  30 // MAX/MIN/LOWEST values for different data types.
  31
  32 // MaxOrMinIdentity returns MAX or LOWEST value of the given type.
  33 template <TypeCategory CAT, int KIND, bool IS_MAXVAL, typename Enable = void>
  34 struct MaxOrMinIdentity {
  35   using Type = CppTypeFor<CAT, KIND>;
  36   static constexpr RT_API_ATTRS Type Value() {
  37     return IS_MAXVAL ? std::numeric_limits<Type>::lowest()
  38                      : std::numeric_limits<Type>::max();
  39   }
  40 };
  41
  42 // std::numeric_limits<> may not know int128_t
  43 template <bool IS_MAXVAL>
  44 struct MaxOrMinIdentity<TypeCategory::Integer, 16, IS_MAXVAL> {
  45   using Type = CppTypeFor<TypeCategory::Integer, 16>;
  46   static constexpr RT_API_ATTRS Type Value() {
  47     return IS_MAXVAL ? Type{1} << 127 : ~Type{0} >> 1;
  48   }
  49 };
  50
  51 #if HAS_FLOAT128
  52 // std::numeric_limits<> may not support __float128.
  53 //
  54 // Usage of GCC quadmath.h's FLT128_MAX is complicated by the fact that
  55 // even GCC complains about 'Q' literal suffix under -Wpedantic.
  56 // We just recreate FLT128_MAX ourselves.
  57 //
  58 // This specialization must engage only when
  59 // CppTypeFor<TypeCategory::Real, 16> is __float128.
  60 template <bool IS_MAXVAL>
  61 struct MaxOrMinIdentity<TypeCategory::Real, 16, IS_MAXVAL,
  62     typename std::enable_if_t<
  63         std::is_same_v<CppTypeFor<TypeCategory::Real, 16>, __float128>>> {
  64   using Type = __float128;
  65   static RT_API_ATTRS Type Value() {
  66     // Create a buffer to store binary representation of __float128 constant.
  67     constexpr std::size_t alignment =
  68         std::max(alignof(Type), alignof(std::uint64_t));
  69     alignas(alignment) char data[sizeof(Type)];
  70
  71     // First, verify that our interpretation of __float128 format is correct,
  72     // e.g. by checking at least one known constant.
  73     *reinterpret_cast<Type *>(data) = Type(1.0);
  74     if (*reinterpret_cast<std::uint64_t *>(data) != 0 ||
  75         *(reinterpret_cast<std::uint64_t *>(data) + 1) != 0x3FFF000000000000) {
  76       Terminator terminator{__FILE__, __LINE__};
  77       terminator.Crash("not yet implemented: no full support for __float128");
  78     }
  79
  80     // Recreate FLT128_MAX.
  81     *reinterpret_cast<std::uint64_t *>(data) = 0xFFFFFFFFFFFFFFFF;
  82     *(reinterpret_cast<std::uint64_t *>(data) + 1) = 0x7FFEFFFFFFFFFFFF;
  83     Type max = *reinterpret_cast<Type *>(data);
  84     return IS_MAXVAL ? -max : max;
  85   }
  86 };
  87 #endif // HAS_FLOAT128
  88
  89 // Minimum finite representable value.
  90 // For floating-point types, returns minimum positive normalized value.
  91 template <typename T> struct MinValue {
  92   static RT_API_ATTRS T get() { return std::numeric_limits<T>::min(); }
  93 };
  94
  95 #if HAS_FLOAT128
  96 template <> struct MinValue<CppTypeFor<TypeCategory::Real, 16>> {
  97   using Type = CppTypeFor<TypeCategory::Real, 16>;
  98   static RT_API_ATTRS Type get() {
  99     // Create a buffer to store binary representation of __float128 constant.
 100     constexpr std::size_t alignment =
 101         std::max(alignof(Type), alignof(std::uint64_t));
 102     alignas(alignment) char data[sizeof(Type)];
 103
 104     // First, verify that our interpretation of __float128 format is correct,
 105     // e.g. by checking at least one known constant.
 106     *reinterpret_cast<Type *>(data) = Type(1.0);
 107     if (*reinterpret_cast<std::uint64_t *>(data) != 0 ||
 108         *(reinterpret_cast<std::uint64_t *>(data) + 1) != 0x3FFF000000000000) {
 109       Terminator terminator{__FILE__, __LINE__};
 110       terminator.Crash("not yet implemented: no full support for __float128");
 111     }
 112
 113     // Recreate FLT128_MIN.
 114     *reinterpret_cast<std::uint64_t *>(data) = 0;
 115     *(reinterpret_cast<std::uint64_t *>(data) + 1) = 0x1000000000000;
 116     return *reinterpret_cast<Type *>(data);
 117   }
 118 };
 119 #endif // HAS_FLOAT128
 120
 121 template <typename T> struct ABSTy {
 122   static constexpr RT_API_ATTRS T compute(T x) { return std::abs(x); }
 123 };
 124
 125 // Suppress the warnings about calling __host__-only
 126 // 'long double' std::frexp, from __device__ code.
 127 RT_DIAG_PUSH
 128 RT_DIAG_DISABLE_CALL_HOST_FROM_DEVICE_WARN
 129
 130 template <typename T> struct FREXPTy {
 131   static constexpr RT_API_ATTRS T compute(T x, int *e) {
 132     return std::frexp(x, e);
 133   }
 134 };
 135
 136 RT_DIAG_POP
 137
 138 template <typename T> struct ILOGBTy {
 139   static constexpr RT_API_ATTRS int compute(T x) { return std::ilogb(x); }
 140 };
 141
 142 template <typename T> struct ISINFTy {
 143   static constexpr RT_API_ATTRS bool compute(T x) { return std::isinf(x); }
 144 };
 145
 146 template <typename T> struct ISNANTy {
 147   static constexpr RT_API_ATTRS bool compute(T x) { return std::isnan(x); }
 148 };
 149
 150 template <typename T> struct LDEXPTy {
 151   template <typename ET> static constexpr RT_API_ATTRS T compute(T x, ET e) {
 152     return std::ldexp(x, e);
 153   }
 154 };
 155
 156 template <typename T> struct MAXTy {
 157   static constexpr RT_API_ATTRS T compute() {
 158     return std::numeric_limits<T>::max();
 159   }
 160 };
 161
 162 #if LDBL_MANT_DIG == 113 || HAS_FLOAT128
 163 template <> struct MAXTy<CppTypeFor<TypeCategory::Real, 16>> {
 164   static CppTypeFor<TypeCategory::Real, 16> compute() {
 165     return MaxOrMinIdentity<TypeCategory::Real, 16, true>::Value();
 166   }
 167 };
 168 #endif
 169
 170 template <typename T> struct MINTy {
 171   static constexpr RT_API_ATTRS T compute() { return MinValue<T>::get(); }
 172 };
 173
 174 template <typename T> struct QNANTy {
 175   static constexpr RT_API_ATTRS T compute() {
 176     return std::numeric_limits<T>::quiet_NaN();
 177   }
 178 };
 179
 180 template <typename T> struct SQRTTy {
 181   static constexpr RT_API_ATTRS T compute(T x) { return std::sqrt(x); }
 182 };
 183
 184 // EXPONENT (16.9.75)
 185 template <typename RESULT, typename ARG>
 186 inline RT_API_ATTRS RESULT Exponent(ARG x) {
 187   if (ISINFTy<ARG>::compute(x) || ISNANTy<ARG>::compute(x)) {
 188     return MAXTy<RESULT>::compute(); // +/-Inf, NaN -> HUGE(0)
 189   } else if (x == 0) {
 190     return 0; // 0 -> 0
 191   } else {
 192     return ILOGBTy<ARG>::compute(x) + 1;
 193   }
 194 }
 195
 196 // FRACTION (16.9.80)
 197 template <typename T> inline RT_API_ATTRS T Fraction(T x) {
 198   if (ISNANTy<T>::compute(x)) {
 199     return x; // NaN -> same NaN
 200   } else if (ISINFTy<T>::compute(x)) {
 201     return QNANTy<T>::compute(); // +/-Inf -> NaN
 202   } else if (x == 0) {
 203     return x; // 0 -> same 0
 204   } else {
 205     int ignoredExp;
 206     return FREXPTy<T>::compute(x, &ignoredExp);
 207   }
 208 }
 209
 210 // SET_EXPONENT (16.9.171)
 211 template <typename T> inline RT_API_ATTRS T SetExponent(T x, std::int64_t p) {
 212   if (ISNANTy<T>::compute(x)) {
 213     return x; // NaN -> same NaN
 214   } else if (ISINFTy<T>::compute(x)) {
 215     return QNANTy<T>::compute(); // +/-Inf -> NaN
 216   } else if (x == 0) {
 217     return x; // return negative zero if x is negative zero
 218   } else {
 219     int expo{ILOGBTy<T>::compute(x) + 1};
 220     auto ip{static_cast<int>(p - expo)};
 221     if (ip != p - expo) {
 222       ip = p < 0 ? std::numeric_limits<int>::min()
 223                  : std::numeric_limits<int>::max();
 224     }
 225     return LDEXPTy<T>::compute(x, ip); // x*2**(p-e)
 226   }
 227 }
 228
 229 // MOD & MODULO (16.9.135, .136)
 230 template <bool IS_MODULO, typename T>
 231 inline RT_API_ATTRS T RealMod(
 232     T a, T p, const char *sourceFile, int sourceLine) {
 233   if (p == 0) {
 234     Terminator{sourceFile, sourceLine}.Crash(
 235         IS_MODULO ? "MODULO with P==0" : "MOD with P==0");
 236   }
 237   if (ISNANTy<T>::compute(a) || ISNANTy<T>::compute(p) ||
 238       ISINFTy<T>::compute(a)) {
 239     return QNANTy<T>::compute();
 240   } else if (IS_MODULO && ISINFTy<T>::compute(p)) {
 241     // Other compilers behave consistently for MOD(x, +/-INF)
 242     // and always return x. This is probably related to
 243     // implementation of std::fmod(). Stick to this behavior
 244     // for MOD, but return NaN for MODULO(x, +/-INF).
 245     return QNANTy<T>::compute();
 246   }
 247   T aAbs{ABSTy<T>::compute(a)};
 248   T pAbs{ABSTy<T>::compute(p)};
 249   if (aAbs <= static_cast<T>(std::numeric_limits<std::int64_t>::max()) &&
 250       pAbs <= static_cast<T>(std::numeric_limits<std::int64_t>::max())) {
 251     if (auto aInt{static_cast<std::int64_t>(a)}; a == aInt) {
 252       if (auto pInt{static_cast<std::int64_t>(p)}; p == pInt) {
 253         // Fast exact case for integer operands
 254         auto mod{aInt - (aInt / pInt) * pInt};
 255         if constexpr (IS_MODULO) {
 256           if (mod == 0) {
 257             // Return properly signed zero.
 258             return pInt > 0 ? T{0} : -T{0};
 259           }
 260           if ((aInt > 0) != (pInt > 0)) {
 261             mod += pInt;
 262           }
 263         } else {
 264           if (mod == 0) {
 265             // Return properly signed zero.
 266             return aInt > 0 ? T{0} : -T{0};
 267           }
 268         }
 269         return static_cast<T>(mod);
 270       }
 271     }
 272   }
 273   if constexpr (std::is_same_v<T, float> || std::is_same_v<T, double> ||
 274       std::is_same_v<T, long double>) {
 275     // std::fmod() semantics on signed operands seems to match
 276     // the requirements of MOD().  MODULO() needs adjustment.
 277     T result{std::fmod(a, p)};
 278     if constexpr (IS_MODULO) {
 279       if ((a < 0) != (p < 0)) {
 280         if (result == 0.) {
 281           result = -result;
 282         } else {
 283           result += p;
 284         }
 285       }
 286     }
 287     return result;
 288   } else {
 289     // The standard defines MOD(a,p)=a-AINT(a/p)*p and
 290     // MODULO(a,p)=a-FLOOR(a/p)*p, but those definitions lose
 291     // precision badly due to cancellation when ABS(a) is
 292     // much larger than ABS(p).
 293     // Insights:
 294     //  - MOD(a,p)=MOD(a-n*p,p) when a>0, p>0, integer n>0, and a>=n*p
 295     //  - when n is a power of two, n*p is exact
 296     //  - as a>=n*p, a-n*p does not round.
 297     // So repeatedly reduce a by all n*p in decreasing order of n;
 298     // what's left is the desired remainder.  This is basically
 299     // the same algorithm as arbitrary precision binary long division,
 300     // discarding the quotient.
 301     T tmp{aAbs};
 302     for (T adj{SetExponent(pAbs, Exponent<int>(aAbs))}; tmp >= pAbs; adj /= 2) {
 303       if (tmp >= adj) {
 304         tmp -= adj;
 305         if (tmp == 0) {
 306           break;
 307         }
 308       }
 309     }
 310     if (a < 0) {
 311       tmp = -tmp;
 312     }
 313     if constexpr (IS_MODULO) {
 314       if ((a < 0) != (p < 0)) {
 315         if (tmp == 0.) {
 316           tmp = -tmp;
 317         } else {
 318           tmp += p;
 319         }
 320       }
 321     }
 322     return tmp;
 323   }
 324 }
 325
 326 // RRSPACING (16.9.164)
 327 template <int PREC, typename T> inline RT_API_ATTRS T RRSpacing(T x) {
 328   if (ISNANTy<T>::compute(x)) {
 329     return x; // NaN -> same NaN
 330   } else if (ISINFTy<T>::compute(x)) {
 331     return QNANTy<T>::compute(); // +/-Inf -> NaN
 332   } else if (x == 0) {
 333     return 0; // 0 -> 0
 334   } else {
 335     return LDEXPTy<T>::compute(
 336         ABSTy<T>::compute(x), PREC - (ILOGBTy<T>::compute(x) + 1));
 337   }
 338 }
 339
 340 // SPACING (16.9.180)
 341 template <int PREC, typename T> inline RT_API_ATTRS T Spacing(T x) {
 342   if (ISNANTy<T>::compute(x)) {
 343     return x; // NaN -> same NaN
 344   } else if (ISINFTy<T>::compute(x)) {
 345     return QNANTy<T>::compute(); // +/-Inf -> NaN
 346   } else if (x == 0) {
 347     // The standard-mandated behavior seems broken, since TINY() can't be
 348     // subnormal.
 349     return MINTy<T>::compute(); // 0 -> TINY(x)
 350   } else {
 351     T result{LDEXPTy<T>::compute(
 352         static_cast<T>(1.0), ILOGBTy<T>::compute(x) + 1 - PREC)}; // 2**(e-p)
 353     return result == 0 ? /*TINY(x)*/ MINTy<T>::compute() : result;
 354   }
 355 }
 356
 357 // ERFC_SCALED (16.9.71)
 358 template <typename T> inline RT_API_ATTRS T ErfcScaled(T arg) {
 359   // Coefficients for approximation to erfc in the first interval.
 360   static const T a[5] = {3.16112374387056560e00, 1.13864154151050156e02,
 361       3.77485237685302021e02, 3.20937758913846947e03, 1.85777706184603153e-1};
 362   static const T b[4] = {2.36012909523441209e01, 2.44024637934444173e02,
 363       1.28261652607737228e03, 2.84423683343917062e03};
 364
 365   // Coefficients for approximation to erfc in the second interval.
 366   static const T c[9] = {5.64188496988670089e-1, 8.88314979438837594e00,
 367       6.61191906371416295e01, 2.98635138197400131e02, 8.81952221241769090e02,
 368       1.71204761263407058e03, 2.05107837782607147e03, 1.23033935479799725e03,
 369       2.15311535474403846e-8};
 370   static const T d[8] = {1.57449261107098347e01, 1.17693950891312499e02,
 371       5.37181101862009858e02, 1.62138957456669019e03, 3.29079923573345963e03,
 372       4.36261909014324716e03, 3.43936767414372164e03, 1.23033935480374942e03};
 373
 374   // Coefficients for approximation to erfc in the third interval.
 375   static const T p[6] = {3.05326634961232344e-1, 3.60344899949804439e-1,
 376       1.25781726111229246e-1, 1.60837851487422766e-2, 6.58749161529837803e-4,
 377       1.63153871373020978e-2};
 378   static const T q[5] = {2.56852019228982242e00, 1.87295284992346047e00,
 379       5.27905102951428412e-1, 6.05183413124413191e-2, 2.33520497626869185e-3};
 380
 381   constexpr T sqrtpi{1.7724538509078120380404576221783883301349L};
 382   constexpr T rsqrtpi{0.5641895835477562869480794515607725858440L};
 383   constexpr T epsilonby2{std::numeric_limits<T>::epsilon() * 0.5};
 384   constexpr T xneg{-26.628e0};
 385   constexpr T xhuge{6.71e7};
 386   constexpr T thresh{0.46875e0};
 387   constexpr T zero{0.0};
 388   constexpr T one{1.0};
 389   constexpr T four{4.0};
 390   constexpr T sixteen{16.0};
 391   constexpr T xmax{1.0 / (sqrtpi * std::numeric_limits<T>::min())};
 392   static_assert(xmax > xhuge, "xmax must be greater than xhuge");
 393
 394   T ysq;
 395   T xnum;
 396   T xden;
 397   T del;
 398   T result;
 399
 400   auto x{arg};
 401   auto y{std::fabs(x)};
 402
 403   if (y <= thresh) {
 404     // evaluate erf for  |x| <= 0.46875
 405     ysq = zero;
 406     if (y > epsilonby2) {
 407       ysq = y * y;
 408     }
 409     xnum = a[4] * ysq;
 410     xden = ysq;
 411     for (int i{0}; i < 3; i++) {
 412       xnum = (xnum + a[i]) * ysq;
 413       xden = (xden + b[i]) * ysq;
 414     }
 415     result = x * (xnum + a[3]) / (xden + b[3]);
 416     result = one - result;
 417     result = std::exp(ysq) * result;
 418     return result;
 419   } else if (y <= four) {
 420     //  evaluate erfc for 0.46875 < |x| <= 4.0
 421     xnum = c[8] * y;
 422     xden = y;
 423     for (int i{0}; i < 7; ++i) {
 424       xnum = (xnum + c[i]) * y;
 425       xden = (xden + d[i]) * y;
 426     }
 427     result = (xnum + c[7]) / (xden + d[7]);
 428   } else {
 429     //  evaluate erfc for |x| > 4.0
 430     result = zero;
 431     if (y >= xhuge) {
 432       if (y < xmax) {
 433         result = rsqrtpi / y;
 434       }
 435     } else {
 436       ysq = one / (y * y);
 437       xnum = p[5] * ysq;
 438       xden = ysq;
 439       for (int i{0}; i < 4; ++i) {
 440         xnum = (xnum + p[i]) * ysq;
 441         xden = (xden + q[i]) * ysq;
 442       }
 443       result = ysq * (xnum + p[4]) / (xden + q[4]);
 444       result = (rsqrtpi - result) / y;
 445     }
 446   }
 447   //  fix up for negative argument, erf, etc.
 448   if (x < zero) {
 449     if (x < xneg) {
 450       result = std::numeric_limits<T>::max();
 451     } else {
 452       ysq = trunc(x * sixteen) / sixteen;
 453       del = (x - ysq) * (x + ysq);
 454       y = std::exp((ysq * ysq)) * std::exp((del));
 455       result = (y + y) - result;
 456     }
 457   }
 458   return result;
 459 }
 460
 461 } // namespace Fortran::runtime
 462
 463 #endif // FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_