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
28namespace Fortran::runtime {
29
30// MAX/MIN/LOWEST values for different data types.
31
32// MaxOrMinIdentity returns MAX or LOWEST value of the given type.
33template <TypeCategory CAT, int KIND, bool IS_MAXVAL, typename Enable = void>
34struct 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
43template <bool IS_MAXVAL>
44struct 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.
60template <bool IS_MAXVAL>
61struct 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.
91template <typename T> struct MinValue {
92 static RT_API_ATTRS T get() { return std::numeric_limits<T>::min(); }
93};
94
95#if HAS_FLOAT128
96template <> 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
121template <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.
127RT_DIAG_PUSH
128RT_DIAG_DISABLE_CALL_HOST_FROM_DEVICE_WARN
129
130template <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
136RT_DIAG_POP
137
138template <typename T> struct ILOGBTy {
139 static constexpr RT_API_ATTRS int compute(T x) { return std::ilogb(x); }
140};
141
142template <typename T> struct ISINFTy {
143 static constexpr RT_API_ATTRS bool compute(T x) { return std::isinf(x); }
144};
145
146template <typename T> struct ISNANTy {
147 static constexpr RT_API_ATTRS bool compute(T x) { return std::isnan(x); }
148};
149
150template <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
156template <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
163template <> 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
170template <typename T> struct MINTy {
171 static constexpr RT_API_ATTRS T compute() { return MinValue<T>::get(); }
172};
173
174template <typename T> struct QNANTy {
175 static constexpr RT_API_ATTRS T compute() {
176 return std::numeric_limits<T>::quiet_NaN();
177 }
178};
179
180template <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)
185template <typename RESULT, typename ARG>
186inline 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)
197template <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)
211template <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)
230template <bool IS_MODULO, typename T>
231inline 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)
327template <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)
341template <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} // namespace Fortran::runtime
358
359#endif // FORTRAN_RUNTIME_NUMERIC_TEMPLATES_H_
360

source code of flang/runtime/numeric-templates.h