Block-Structured AMR Software Framework
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
AMReX_Math.H
Go to the documentation of this file.
1#ifndef AMREX_MATH_H_
2#define AMREX_MATH_H_
3#include <AMReX_Config.H>
4
6#include <AMReX_Extension.H>
7#include <AMReX_INT.H>
8#include <AMReX_REAL.H>
9#include <cmath>
10#include <cstdlib>
11#include <type_traits>
12#include <utility>
13
14#ifdef AMREX_USE_SYCL
15# include <sycl/sycl.hpp>
16#endif
17
18namespace amrex { inline namespace disabled {
19 // If it is inside namespace amrex, or amrex namespace is imported with using namespace amrex or
20 // amrex::disabled, unqualified abs functions are disabled with a compile time error such as,
21 // call of overload abs(int&) is ambiguous, or a link time error such as, undefined reference to
22 // `amrex::disabled::abs(double)'. To fix it, one can use `std::abs` or `amrex::Math::abs`.
23 // The latter works in both host and device functions, whereas `std::abs` does not currently
24 // work on device with HIP and SYCL.
25 AMREX_GPU_HOST_DEVICE double abs (double);
27 AMREX_GPU_HOST_DEVICE long double abs (long double);
30 AMREX_GPU_HOST_DEVICE long long abs (long long);
31}}
32
33namespace amrex::Math {
34
35// Since Intel's SYCL compiler now supports the following std functions on device,
36// one no longer needs to use amrex::Math::abs, etc. They are kept here for
37// backward compatibility.
38
39using std::abs;
40using std::ceil;
41using std::copysign;
42using std::floor;
43using std::round;
44
45// However, since Intel's SYCL compiler is very aggressive with fast floating
46// point optimisations, the following must be kept, as using the std functions
47// always evaluates to false (even at -O1).
48
49#ifdef AMREX_USE_SYCL
50
51using sycl::isfinite;
52using sycl::isinf;
53
54#else
55
56using std::isfinite;
57using std::isinf;
58
59#endif
60
61template <typename T>
62constexpr std::enable_if_t<std::is_floating_point_v<T>,T> pi ()
63{
64 return T(3.1415926535897932384626433832795029L);
65}
66
69double cospi (double x)
70{
71#if defined(AMREX_USE_SYCL)
72 return sycl::cospi(x);
73#else
74 AMREX_IF_ON_DEVICE(( return ::cospi(x); ))
75 AMREX_IF_ON_HOST(( return std::cos(pi<double>()*x); ))
76#endif
77}
78
81float cospi (float x)
82{
83#if defined(AMREX_USE_SYCL)
84 return sycl::cospi(x);
85#else
86 AMREX_IF_ON_DEVICE(( return ::cospif(x); ))
87 AMREX_IF_ON_HOST(( return std::cos(pi<float>()*x); ))
88#endif
89}
90
93double sinpi (double x)
94{
95#if defined(AMREX_USE_SYCL)
96 return sycl::sinpi(x);
97#else
98 AMREX_IF_ON_DEVICE(( return ::sinpi(x); ))
99 AMREX_IF_ON_HOST(( return std::sin(pi<double>()*x); ))
100#endif
101}
102
105float sinpi (float x)
106{
107#if defined(AMREX_USE_SYCL)
108 return sycl::sinpi(x);
109#else
110 AMREX_IF_ON_DEVICE(( return ::sinpif(x); ))
111 AMREX_IF_ON_HOST(( return std::sin(pi<float>()*x); ))
112#endif
113}
114
115namespace detail {
116 AMREX_FORCE_INLINE void sincos (double x, double* sinx, double* cosx) {
117#if defined(_GNU_SOURCE) && !defined(__APPLE__)
118 ::sincos(x, sinx, cosx);
119#else
120 *sinx = std::sin(x);
121 *cosx = std::cos(x);
122#endif
123 }
124
125 AMREX_FORCE_INLINE void sincosf (float x, float* sinx, float* cosx) {
126#if defined(_GNU_SOURCE) && !defined(__APPLE__)
127 ::sincosf(x, sinx, cosx);
128#else
129 *sinx = std::sin(x);
130 *cosx = std::cos(x);
131#endif
132 }
133}
134
137std::pair<double,double> sincos (double x)
138{
139 std::pair<double,double> r;
140#if defined(AMREX_USE_SYCL)
141 r.first = sycl::sincos(x, sycl::private_ptr<double>(&r.second));
142#else
143 AMREX_IF_ON_DEVICE(( ::sincos(x, &r.first, &r.second); ))
144 AMREX_IF_ON_HOST(( detail::sincos(x, &r.first, &r.second); ))
145#endif
146 return r;
147}
148
151std::pair<float,float> sincos (float x)
152{
153 std::pair<float,float> r;
154#if defined(AMREX_USE_SYCL)
155 r.first = sycl::sincos(x, sycl::private_ptr<float>(&r.second));
156#else
157 AMREX_IF_ON_DEVICE(( ::sincosf(x, &r.first, &r.second); ))
158 AMREX_IF_ON_HOST(( detail::sincosf(x, &r.first, &r.second); ))
159#endif
160 return r;
161}
162
165std::pair<double,double> sincospi (double x)
166{
167 std::pair<double,double> r;
168#if defined(AMREX_USE_SYCL)
169 r = sincos(pi<double>()*x);
170#else
171 AMREX_IF_ON_DEVICE(( ::sincospi(x, &r.first, &r.second); ))
172 AMREX_IF_ON_HOST(( r = sincos(pi<double>()*x); ))
173#endif
174 return r;
175}
176
179std::pair<float,float> sincospi (float x)
180{
181 std::pair<float,float> r;
182#if defined(AMREX_USE_SYCL)
183 r = sincos(pi<float>()*x);
184#else
185 AMREX_IF_ON_DEVICE(( ::sincospif(x, &r.first, &r.second); ))
186 AMREX_IF_ON_HOST(( r = sincos(pi<float>()*x); ))
187#endif
188 return r;
189}
190
192template <int Power, typename T,
193 typename = std::enable_if_t<!std::is_integral<T>() || Power>=0>>
195constexpr T powi (T x) noexcept
196{
197 if constexpr (Power < 0) {
198 return T(1)/powi<-Power>(x);
199 } else if constexpr (Power == 0) {
200 //note: 0^0 is implementation-defined, but most compilers return 1
201 return T(1);
202 } else if constexpr (Power == 1) {
203 return x;
204 } else if constexpr (Power == 2) {
205 return x*x;
206 } else if constexpr (Power%2 == 0) {
207 return powi<2>(powi<Power/2>(x));
208 } else {
209 return x*powi<Power-1>(x);
210 }
211}
212
213#if defined(AMREX_INT128_SUPPORTED)
215std::uint64_t umulhi (std::uint64_t a, std::uint64_t b)
216{
217#if defined(AMREX_USE_SYCL)
218 return sycl::mul_hi(a,b);
219#else
220 AMREX_IF_ON_DEVICE(( return __umul64hi(a, b); ))
222 auto tmp = amrex::UInt128_t(a) * amrex::UInt128_t(b);
223 return std::uint64_t(tmp >> 64);
224 ))
225#endif
226}
227#endif
228
229template <typename T>
232{
233 // Computing K based on DLMF
234 // https://dlmf.nist.gov/19.8
235 T tol = std::numeric_limits<T>::epsilon();
236
237 T a0 = 1.0;
238 T g0 = std::sqrt(1.0 - k*k);
239 T a = a0;
240 T g = g0;
241
242 // Find Arithmetic Geometric mean
243 while(std::abs(a0 - g0) > tol) {
244 a = 0.5*(a0 + g0);
245 g = std::sqrt(a0 * g0);
246
247 a0 = a;
248 g0 = g;
249 }
250
251 return 0.5*pi<T>()/a;
252}
253
254template <typename T>
257{
258 // Computing E based on DLMF
259 // https://dlmf.nist.gov/19.8
260 T Kcomp = amrex::Math::comp_ellint_1<T>(k);
261 T tol = std::numeric_limits<T>::epsilon();
262
263 // Step Zero
264 T a0 = 1.0;
265 T g0 = std::sqrt(1.0 - k*k);
266 T cn = std::sqrt(a0*a0 - g0*g0);
267
268 // Step 1
269 int n = 1;
270 T a = 0.5 * (a0 + g0);
271 T g = std::sqrt(a0*g0);
272 cn = 0.25*cn*cn/a;
273
274 T sum_val = a*a;
275 a0 = a;
276 g0 = g;
277
278 while(std::abs(cn*cn) > tol) {
279 // Compute coefficients for this iteration
280 a = 0.5 * (a0 + g0);
281 g = std::sqrt(a0*g0);
282 cn = 0.25*cn*cn/a;
283
284 n++;
285 sum_val -= std::pow(2,n-1)*cn*cn;
286
287 // Save a and g for next iteration
288 a0 = a;
289 g0 = g;
290 }
291
292 return Kcomp*sum_val;
293}
294
295/***************************************************************************************************
296 * Copyright (c) 2017 - 2024 NVIDIA CORPORATION & AFFILIATES. All rights reserved.
297 * SPDX-License-Identifier: BSD-3-Clause
298 *
299 * Redistribution and use in source and binary forms, with or without
300 * modification, are permitted provided that the following conditions are met:
301 *
302 * 1. Redistributions of source code must retain the above copyright notice, this
303 * list of conditions and the following disclaimer.
304 *
305 * 2. Redistributions in binary form must reproduce the above copyright notice,
306 * this list of conditions and the following disclaimer in the documentation
307 * and/or other materials provided with the distribution.
308 *
309 * 3. Neither the name of the copyright holder nor the names of its
310 * contributors may be used to endorse or promote products derived from
311 * this software without specific prior written permission.
312 *
313 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
314 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
315 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
316 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
317 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
318 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
319 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
320 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
321 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
322 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
323 *
324 **************************************************************************************************/
325
327
343{
344 std::uint64_t divisor = 0;
345
346#ifdef AMREX_INT128_SUPPORTED
347 std::uint64_t multiplier = 1U;
348 unsigned int shift_right = 0;
349 unsigned int round_up = 0;
350
351 //
352 // Static methods
353 //
354
356 static std::uint32_t integer_log2 (std::uint64_t x)
357 {
358 std::uint32_t n = 0;
359 while (x >>= 1) {
360 ++n;
361 }
362 return n;
363 }
364
368 FastDivmodU64 (std::uint64_t divisor_)
369 : divisor(divisor_)
370 {
371 if (divisor) {
372 shift_right = integer_log2(divisor);
373
374 if ((divisor & (divisor - 1)) == 0) {
375 multiplier = 0;
376 }
377 else {
378 std::uint64_t power_of_two = (std::uint64_t(1) << shift_right);
379 auto n = amrex::UInt128_t(power_of_two) << 64;
380 std::uint64_t multiplier_lo = n / divisor;
381 n += power_of_two;
382 multiplier = n / divisor;
383 round_up = (multiplier_lo == multiplier ? 1 : 0);
384 }
385 }
386 }
387
388#else
389
390 FastDivmodU64 (std::uint64_t divisor_) : divisor(divisor_) {}
391
392#endif
393
395 FastDivmodU64 () = default;
396
398 [[nodiscard]] AMREX_GPU_HOST_DEVICE
399 std::uint64_t divide (std::uint64_t dividend) const
400 {
401#if defined(AMREX_INT128_SUPPORTED)
402 auto x = dividend;
403 if (multiplier) {
404 x = amrex::Math::umulhi(dividend + round_up, multiplier);
405 }
406 return (x >> shift_right);
407#else
408 return dividend / divisor;
409#endif
410 }
411
413 [[nodiscard]] AMREX_GPU_HOST_DEVICE
414 std::uint64_t modulus (std::uint64_t quotient, std::uint64_t dividend) const
415 {
416 return dividend - quotient * divisor;
417 }
418
420 [[nodiscard]] AMREX_GPU_HOST_DEVICE
421 std::uint64_t divmod (std::uint64_t &remainder, std::uint64_t dividend) const
422 {
423 auto quotient = divide(dividend);
424 remainder = modulus(quotient, dividend);
425 return quotient;
426 }
427
431 void operator() (std::uint64_t &quotient, std::uint64_t &remainder, std::uint64_t dividend) const
432 {
433 quotient = divmod(remainder, dividend);
434 }
435};
436
437}
438
439#endif
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_IF_ON_DEVICE(CODE)
Definition AMReX_GpuQualifiers.H:56
#define AMREX_IF_ON_HOST(CODE)
Definition AMReX_GpuQualifiers.H:58
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
AMREX_FORCE_INLINE void sincos(double x, double *sinx, double *cosx)
Definition AMReX_Math.H:116
AMREX_FORCE_INLINE void sincosf(float x, float *sinx, float *cosx)
Definition AMReX_Math.H:125
Definition AMReX_Math.H:33
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T comp_ellint_1(T k)
Definition AMReX_Math.H:231
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE double sinpi(double x)
Return sin(x*pi) given x.
Definition AMReX_Math.H:93
constexpr std::enable_if_t< std::is_floating_point_v< T >, T > pi()
Definition AMReX_Math.H:62
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T comp_ellint_2(T k)
Definition AMReX_Math.H:256
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::pair< double, double > sincospi(double x)
Return sin(pi*x) and cos(pi*x) given x.
Definition AMReX_Math.H:165
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE double cospi(double x)
Return cos(x*pi) given x.
Definition AMReX_Math.H:69
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::pair< double, double > sincos(double x)
Return sine and cosine of given number.
Definition AMReX_Math.H:137
AMREX_FORCE_INLINE constexpr T powi(T x) noexcept
Return pow(x, Power), where Power is an integer known at compile time.
Definition AMReX_Math.H:195
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition AMReX_GpuComplex.H:356
Definition AMReX_FabArrayCommI.H:896
Definition AMReX_Math.H:343
AMREX_GPU_HOST_DEVICE std::uint64_t divmod(std::uint64_t &remainder, std::uint64_t dividend) const
Returns the quotient of floor(dividend / divisor) and computes the remainder.
Definition AMReX_Math.H:421
AMREX_GPU_HOST_DEVICE std::uint64_t modulus(std::uint64_t quotient, std::uint64_t dividend) const
Computes the remainder given a computed quotient and dividend.
Definition AMReX_Math.H:414
std::uint64_t divisor
Definition AMReX_Math.H:344
AMREX_GPU_HOST_DEVICE void operator()(std::uint64_t &quotient, std::uint64_t &remainder, std::uint64_t dividend) const
Definition AMReX_Math.H:431
AMREX_GPU_HOST_DEVICE std::uint64_t divide(std::uint64_t dividend) const
Returns the quotient of floor(dividend / divisor)
Definition AMReX_Math.H:399
FastDivmodU64(std::uint64_t divisor_)
Definition AMReX_Math.H:390
FastDivmodU64()=default
Default construct an invalid FastDivmodU64.