Kokkos Core Kernels Package Version of the Day
Loading...
Searching...
No Matches
Kokkos_MathematicalSpecialFunctions.hpp
1//@HEADER
2// ************************************************************************
3//
4// Kokkos v. 4.0
5// Copyright (2022) National Technology & Engineering
6// Solutions of Sandia, LLC (NTESS).
7//
8// Under the terms of Contract DE-NA0003525 with NTESS,
9// the U.S. Government retains certain rights in this software.
10//
11// Part of Kokkos, under the Apache License v2.0 with LLVM Exceptions.
12// See https://kokkos.org/LICENSE for license information.
13// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
14//
15//@HEADER
16
17#ifndef KOKKOS_MATHEMATICAL_SPECIAL_FUNCTIONS_HPP
18#define KOKKOS_MATHEMATICAL_SPECIAL_FUNCTIONS_HPP
19#ifndef KOKKOS_IMPL_PUBLIC_INCLUDE
20#define KOKKOS_IMPL_PUBLIC_INCLUDE
21#define KOKKOS_IMPL_PUBLIC_INCLUDE_NOTDEFINED_MATHSPECFUNCTIONS
22#endif
23
24#include <Kokkos_Macros.hpp>
25#include <cmath>
26#include <algorithm>
27#include <type_traits>
28#include <Kokkos_MathematicalConstants.hpp>
29#include <Kokkos_MathematicalFunctions.hpp>
30#include <Kokkos_NumericTraits.hpp>
31#include <Kokkos_Complex.hpp>
32
33namespace Kokkos {
34namespace Experimental {
35
37template <class RealType>
38KOKKOS_INLINE_FUNCTION RealType expint1(RealType x) {
39 // This function is a conversion of the corresponding Fortran program in
40 // S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).
41 using Kokkos::exp;
42 using Kokkos::fabs;
43 using Kokkos::log;
44 using Kokkos::pow;
45 using Kokkos::Experimental::epsilon;
46 using Kokkos::Experimental::infinity;
47
48 RealType e1;
49
50 if (x < 0) {
51 e1 = -infinity<RealType>::value;
52 } else if (x == 0.0) {
53 e1 = infinity<RealType>::value;
54 } else if (x <= 1.0) {
55 e1 = 1.0;
56 RealType r = 1.0;
57 for (int k = 1; k <= 25; k++) {
58 RealType k_real = static_cast<RealType>(k);
59 r = -r * k_real * x / pow(k_real + 1.0, 2.0);
60 e1 = e1 + r;
61 if (fabs(r) <= fabs(e1) * epsilon<RealType>::value) break;
62 }
63 e1 = -0.5772156649015328 - log(x) + x * e1;
64 } else {
65 int m = 20 + static_cast<int>(80.0 / x);
66 RealType t0 = 0.0;
67 for (int k = m; k >= 1; k--) {
68 RealType k_real = static_cast<RealType>(k);
69 t0 = k_real / (1.0 + k_real / (x + t0));
70 }
71 e1 = exp(-x) * (1.0 / (x + t0));
72 }
73 return e1;
74}
75
77template <class RealType>
78KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> erf(
80 // This function is a conversion of the corresponding Fortran program written
81 // by D.E. Amos, May,1974. D.E. Amos' revisions of Jan 86 incorporated by
82 // Ken Damrau on 27-Jan-1986 14:37:13
83 //
84 // Reference: NBS HANDBOOK OF MATHEMATICAL FUNCTIONS, AMS 55, By
85 // M. ABRAMOWITZ AND I.A. STEGUN, December,1955.
86 // Summary:
87 // If x < 0, z is replaced by -z and all computation is done in the right
88 // half lane, except for z inside the circle abs(z)<=2, since
89 // erf(-z)=-erf(z). The regions for computation are divided as follows
90 // (1) abs(z)<=2 - Power series, NBS Handbook, p. 298
91 // (2) abs(z)>2 and x>1 - continued fraction, NBS Handbook, p. 298
92 // (3) abs(z)>2 and 0<=x<=1 and abs(y)<6 - series, NBS Handbook, p. 299
93 // (4) abs(z)>2 and 0<=x<=1 and abs(y)>=6 - asymptotic expansion
94 // Error condition: abs(z^2) > 670 is a fatal overflow error
95 using Kokkos::cos;
96 using Kokkos::exp;
97 using Kokkos::fabs;
98 using Kokkos::sin;
99 using Kokkos::Experimental::epsilon_v;
100 using Kokkos::Experimental::infinity_v;
101 using Kokkos::numbers::pi_v;
102
103 using CmplxType = Kokkos::complex<RealType>;
104
105 constexpr auto inf = infinity_v<RealType>;
106 constexpr auto tol = epsilon_v<RealType>;
107
108 const RealType fnorm = 1.12837916709551;
109 const RealType gnorm = 0.564189583547756;
110 const RealType eh = 0.606530659712633;
111 const RealType ef = 0.778800783071405;
112 // const RealType tol = 1.0e-13;
113 constexpr auto pi = pi_v<RealType>;
114
115 CmplxType cans;
116
117 RealType az = Kokkos::abs(z);
118 if (az <= 2.0) { // Series for abs(z)<=2.0
119 CmplxType cz = z * z;
120 CmplxType accum = CmplxType(1.0, 0.0);
121 CmplxType term = accum;
122 RealType ak = 1.5;
123 for (int i = 1; i <= 35; i++) {
124 term = term * cz / ak;
125 accum = accum + term;
126 if (Kokkos::abs(term) <= tol) break;
127 ak = ak + 1.0;
128 }
129 cz = -cz;
130 RealType er = cz.real();
131 RealType ei = cz.imag();
132 accum = accum * z * fnorm;
133 cz = exp(er) * CmplxType(cos(ei), sin(ei));
134 cans = accum * cz;
135 } // end (az <= 2.0)
136 else { //(az > 2.0)
137 CmplxType zp = z;
138 if (z.real() < 0.0) zp = -z;
139 CmplxType cz = zp * zp;
140 RealType xp = zp.real();
141 RealType yp = zp.imag();
142 if (xp > 1.0) {
143 // continued fraction for erfc(z), abs(Z)>2
144 int n = static_cast<int>(100.0 / az + 5.0);
145 int fn = n;
146 CmplxType term = cz;
147 for (int i = 1; i <= n; i++) {
148 RealType fnh = fn - 0.5;
149 term = cz + (fnh * term) / (fn + term);
150 fn = fn - 1;
151 }
152 if (Kokkos::abs(cz) > 670.0) return CmplxType(inf, inf);
153 cz = -cz;
154 RealType er = cz.real();
155 RealType ei = cz.imag();
156 cz = exp(er) * CmplxType(cos(ei), sin(ei));
157 CmplxType accum = zp * gnorm * cz;
158 cans = 1.0 - accum / term;
159 if (z.real() < 0.0) cans = -cans;
160 } // end (xp > 1.0)
161 else { //(xp <= 1.0)
162 if (fabs(yp) <
163 6.0) { // Series (3) for abs(z)>2 and 0<=xp<=1 and abs(yp)<6
164 RealType s1 = 0.0;
165 RealType s2 = 0.0;
166 RealType x2 = xp * xp;
167 RealType fx2 = 4.0 * x2;
168 RealType tx = xp + xp;
169 RealType xy = xp * yp;
170 RealType sxyh = sin(xy);
171 RealType sxy = sin(xy + xy);
172 RealType cxy = cos(xy + xy);
173 RealType fn = 1.0;
174 RealType fnh = 0.5;
175 RealType ey = exp(yp);
176 RealType en = ey;
177 RealType ehn = eh;
178 RealType un = ef;
179 RealType vn = 1.0;
180 for (int i = 1; i <= 50; i++) {
181 RealType ren = 1.0 / en;
182 RealType csh = en + ren;
183 RealType tm = xp * csh;
184 RealType ssh = en - ren;
185 RealType tmp = fnh * ssh;
186 RealType rn = tx - tm * cxy + tmp * sxy;
187 RealType ain = tm * sxy + tmp * cxy;
188 RealType cf = un / (vn + fx2);
189 rn = cf * rn;
190 ain = cf * ain;
191 s1 = s1 + rn;
192 s2 = s2 + ain;
193 if ((fabs(rn) + fabs(ain)) < tol * (fabs(s1) + fabs(s2))) break;
194 un = un * ehn * ef;
195 ehn = ehn * eh;
196 en = en * ey;
197 vn = vn + fn + fn + 1.0;
198 fnh = fnh + 0.5;
199 fn = fn + 1.0;
200 }
201 s1 = s1 + s1;
202 s2 = s2 + s2;
203 if (z.real() == 0.0)
204 s2 = s2 + yp;
205 else {
206 s1 = s1 + sxyh * sxyh / xp;
207 s2 = s2 + sxy / tx;
208 }
209 // Power series for erf(xp), 0<=xp<=1
210 RealType w = 1.0;
211 RealType ak = 1.5;
212 RealType tm = 1.0;
213 for (int i = 1; i <= 17; i++) {
214 tm = tm * x2 / ak;
215 w = w + tm;
216 if (tm <= tol) break;
217 ak = ak + 1.0;
218 }
219 RealType ex = exp(-x2);
220 w = w * xp * fnorm * ex;
221 RealType cf = ex / pi;
222 s1 = cf * s1 + w;
223 s2 = cf * s2;
224 cans = CmplxType(s1, s2);
225 if (z.real() < 0.0) cans = -cans;
226 } // end (abs(yp) < 6.0)
227 else { //(abs(YP)>=6.0)
228 // Asymptotic expansion for 0<=xp<=1 and abs(yp)>=6
229 CmplxType rcz = 0.5 / cz;
230 CmplxType accum = CmplxType(1.0, 0.0);
231 CmplxType term = accum;
232 RealType ak = 1.0;
233 for (int i = 1; i <= 35; i++) {
234 term = -term * ak * rcz;
235 accum = accum + term;
236 if (Kokkos::abs(term) / Kokkos::abs(accum) <= tol) break;
237 ak = ak + 2.0;
238 }
239 accum = accum * gnorm / zp;
240 cz = -cz;
241 RealType er = cz.real();
242 if (fabs(er) > 670.0) return CmplxType(inf, inf);
243 RealType ei = cz.imag();
244 cz = exp(er) * CmplxType(cos(ei), sin(ei));
245 cans = 1.0 - accum * cz;
246 if (z.real() < 0.0) cans = -cans;
247 } // end (abs(YP)>=6.0)
248 } // end (xp <= 1.0)
249 } // end (az > 2.0)
250 return cans;
251}
252
255template <class RealType>
256KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> erfcx(
257 const Kokkos::complex<RealType>& z) {
258 // This function is a conversion of the corresponding Fortran program written
259 // by D.E. Amos, May,1974. D.E. Amos' revisions of Jan 86 incorporated by
260 // Ken Damrau on 27-Jan-1986 14:37:13
261 //
262 // Reference: NBS HANDBOOK OF MATHEMATICAL FUNCTIONS, AMS 55, By
263 // M. ABRAMOWITZ AND I.A. STEGUN, December,1955.
264 // Summary:
265 // If x < 0, z is replaced by -z and all computation is done in the right
266 // half lane, except for z inside the circle abs(z)<=2, since
267 // erfc(-z)=2-erfc(z). The regions for computation are divided as follows
268 // (1) abs(z)<=2 - Power series, NBS Handbook, p. 298
269 // (2) abs(z)>2 and x>1 - continued fraction, NBS Handbook, p. 298
270 // (3) abs(z)>2 and 0<=x<=1 and abs(y)<6 - series, NBS Handbook, p. 299
271 // (4) abs(z)>2 and 0<=x<=1 and abs(y)>=6 - asymptotic expansion
272 // Error condition: abs(z^2) > 670 is a fatal overflow error when x<0
273 using Kokkos::cos;
274 using Kokkos::exp;
275 using Kokkos::fabs;
276 using Kokkos::isinf;
277 using Kokkos::sin;
278 using Kokkos::Experimental::epsilon_v;
279 using Kokkos::Experimental::infinity_v;
280 using Kokkos::numbers::inv_sqrtpi_v;
281 using Kokkos::numbers::pi_v;
282
283 using CmplxType = Kokkos::complex<RealType>;
284
285 constexpr auto inf = infinity_v<RealType>;
286 constexpr auto tol = epsilon_v<RealType>;
287
288 const RealType fnorm = 1.12837916709551;
289 constexpr auto gnorm = inv_sqrtpi_v<RealType>;
290 const RealType eh = 0.606530659712633;
291 const RealType ef = 0.778800783071405;
292 // const RealType tol = 1.0e-13;
293 constexpr auto pi = pi_v<RealType>;
294
295 CmplxType cans;
296
297 if ((isinf(z.real())) && (z.real() > 0)) {
298 cans = CmplxType(0.0, 0.0);
299 return cans;
300 }
301 if ((isinf(z.real())) && (z.real() < 0)) {
302 cans = CmplxType(inf, inf);
303 return cans;
304 }
305
306 RealType az = Kokkos::abs(z);
307 if (az <= 2.0) { // Series for abs(z)<=2.0
308 CmplxType cz = z * z;
309 CmplxType accum = CmplxType(1.0, 0.0);
310 CmplxType term = accum;
311 RealType ak = 1.5;
312 for (int i = 1; i <= 35; i++) {
313 term = term * cz / ak;
314 accum = accum + term;
315 if (Kokkos::abs(term) <= tol) break;
316 ak = ak + 1.0;
317 }
318 cz = -cz;
319 RealType er = cz.real();
320 RealType ei = cz.imag();
321 accum = accum * z * fnorm;
322 cz = exp(er) * CmplxType(cos(ei), sin(ei));
323 cans = 1.0 / cz - accum;
324 } // end (az <= 2.0)
325 else { //(az > 2.0)
326 CmplxType zp = z;
327 if (z.real() < 0.0) zp = -z;
328 CmplxType cz = zp * zp;
329 RealType xp = zp.real();
330 RealType yp = zp.imag();
331 if (xp > 1.0) {
332 // continued fraction for erfc(z), abs(z)>2
333 int n = static_cast<int>(100.0 / az + 5.0);
334 int fn = n;
335 CmplxType term = cz;
336 for (int i = 1; i <= n; i++) {
337 RealType fnh = fn - 0.5;
338 term = cz + (fnh * term) / (fn + term);
339 fn = fn - 1;
340 }
341 cans = zp * gnorm / term;
342 if (z.real() >= 0.0) return cans;
343 if (Kokkos::abs(cz) > 670.0) return CmplxType(inf, inf);
344 ;
345 cz = -cz;
346 RealType er = cz.real();
347 RealType ei = cz.imag();
348 cz = exp(er) * CmplxType(cos(ei), sin(ei));
349 cz = 1.0 / cz;
350 cans = cz + cz - cans;
351 } // end (xp > 1.0)
352 else { //(xp <= 1.0)
353 if (fabs(yp) <
354 6.0) { // Series (3) for abs(z)>2 and 0<=xp<=1 and abs(yp)<6
355 RealType s1 = 0.0;
356 RealType s2 = 0.0;
357 RealType x2 = xp * xp;
358 RealType fx2 = 4.0 * x2;
359 RealType tx = xp + xp;
360 RealType xy = xp * yp;
361 RealType sxyh = sin(xy);
362 RealType sxy = sin(xy + xy);
363 RealType cxy = cos(xy + xy);
364 RealType fn = 1.0;
365 RealType fnh = 0.5;
366 RealType ey = exp(yp);
367 RealType en = ey;
368 RealType ehn = eh;
369 RealType un = ef;
370 RealType vn = 1.0;
371 for (int i = 1; i <= 50; i++) {
372 RealType ren = 1.0 / en;
373 RealType csh = en + ren;
374 RealType tm = xp * csh;
375 RealType ssh = en - ren;
376 RealType tmp = fnh * ssh;
377 RealType rn = tx - tm * cxy + tmp * sxy;
378 RealType ain = tm * sxy + tmp * cxy;
379 RealType cf = un / (vn + fx2);
380 rn = cf * rn;
381 ain = cf * ain;
382 s1 = s1 + rn;
383 s2 = s2 + ain;
384 if ((fabs(rn) + fabs(ain)) < tol * (fabs(s1) + fabs(s2))) break;
385 un = un * ehn * ef;
386 ehn = ehn * eh;
387 en = en * ey;
388 vn = vn + fn + fn + 1.0;
389 fnh = fnh + 0.5;
390 fn = fn + 1.0;
391 }
392 s1 = s1 + s1;
393 s2 = s2 + s2;
394 if (z.real() == 0.0)
395 s2 = s2 + yp;
396 else {
397 s1 = s1 + sxyh * sxyh / xp;
398 s2 = s2 + sxy / tx;
399 }
400 // Power series for erf(xp), 0<=xp<=1
401 RealType w = 1.0;
402 RealType ak = 1.5;
403 RealType tm = 1.0;
404 for (int i = 1; i <= 17; i++) {
405 tm = tm * x2 / ak;
406 w = w + tm;
407 if (tm <= tol) break;
408 ak = ak + 1.0;
409 }
410 RealType ex = exp(-x2);
411 w = w * xp * fnorm * ex;
412 CmplxType rcz = CmplxType(cxy, sxy);
413 RealType y2 = yp * yp;
414 cz = exp(x2 - y2) * rcz;
415 rcz = exp(-y2) * rcz;
416 if (z.real() >= 0.0)
417 cans = cz * (1.0 - w) - rcz * CmplxType(s1, s2) / pi;
418 else
419 cans = cz * (1.0 + w) + rcz * CmplxType(s1, s2) / pi;
420 } // end (abs(yp) < 6.0)
421 else { //(abs(YP)>=6.0)
422 // Asymptotic expansion for 0<=xp<=1 and abs(yp)>=6
423 CmplxType rcz = 0.5 / cz;
424 CmplxType accum = CmplxType(1.0, 0.0);
425 CmplxType term = accum;
426 RealType ak = 1.0;
427 for (int i = 1; i <= 35; i++) {
428 term = -term * ak * rcz;
429 accum = accum + term;
430 if (Kokkos::abs(term) / Kokkos::abs(accum) <= tol) break;
431 ak = ak + 2.0;
432 }
433 accum = accum * gnorm / zp;
434 if (z.real() < 0.0) accum = -accum;
435 cans = accum;
436 } // end (abs(YP)>=6.0)
437 } // end (xp <= 1.0)
438 } // end (az > 2.0)
439 return cans;
440}
441
444template <class RealType>
445KOKKOS_INLINE_FUNCTION RealType erfcx(RealType x) {
446 using CmplxType = Kokkos::complex<RealType>;
447 // Note: using erfcx(complex) for now
448 // TODO: replace with an implementation of erfcx(real)
449 CmplxType zin = CmplxType(x, 0.0);
450 CmplxType zout = erfcx(zin);
451 return zout.real();
452}
453
456template <class CmplxType, class RealType, class IntType>
457KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_j0(const CmplxType& z,
458 const RealType& joint_val = 25,
459 const IntType& bw_start = 70) {
460 // This function is converted and modified from the corresponding Fortran
461 // program CJYNB in S. Zhang & J. Jin "Computation of Special Functions"
462 //(Wiley, 1996).
463 // Input : z --- Complex argument
464 // joint_val --- Joint point of abs(z) separating small and large
465 // argument regions
466 // bw_start --- Starting point for backward recurrence
467 // Output: cbj0 --- J0(z)
468 using Kokkos::fabs;
469 using Kokkos::pow;
470 using Kokkos::numbers::pi_v;
471
472 CmplxType cbj0;
473 constexpr auto pi = pi_v<RealType>;
474 const RealType a[12] = {
475 -0.703125e-01, 0.112152099609375e+00, -0.5725014209747314e+00,
476 0.6074042001273483e+01, -0.1100171402692467e+03, 0.3038090510922384e+04,
477 -0.1188384262567832e+06, 0.6252951493434797e+07, -0.4259392165047669e+09,
478 0.3646840080706556e+11, -0.3833534661393944e+13, 0.4854014686852901e+15};
479 const RealType b[12] = {0.732421875e-01, -0.2271080017089844e+00,
480 0.1727727502584457e+01, -0.2438052969955606e+02,
481 0.5513358961220206e+03, -0.1825775547429318e+05,
482 0.8328593040162893e+06, -0.5006958953198893e+08,
483 0.3836255180230433e+10, -0.3649010818849833e+12,
484 0.4218971570284096e+14, -0.5827244631566907e+16};
485
486 RealType r2p = 2.0 / pi;
487 RealType a0 = Kokkos::abs(z);
488 RealType y0 = fabs(z.imag());
489 CmplxType z1 = z;
490
491 if (a0 < 1e-100) { // Treat z=0 as a special case
492 cbj0 = CmplxType(1.0, 0.0);
493 } else {
494 if (z.real() < 0.0) z1 = -z;
495 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
496 // (default:25)
497 CmplxType cbs = CmplxType(0.0, 0.0);
498 CmplxType csu = CmplxType(0.0, 0.0);
499 CmplxType csv = CmplxType(0.0, 0.0);
500 CmplxType cf2 = CmplxType(0.0, 0.0);
501 CmplxType cf1 = CmplxType(1e-100, 0.0);
502 CmplxType cf, cs0;
503 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
504 // 70)
505 cf = 2.0 * (k + 1.0) / z * cf1 - cf2;
506 RealType tmp_exponent = static_cast<RealType>(k / 2);
507 if (k == 0) cbj0 = cf;
508 if ((k == 2 * (k / 2)) && (k != 0)) {
509 if (y0 <= 1.0)
510 cbs = cbs + 2.0 * cf;
511 else
512 cbs = cbs + pow(-1.0, tmp_exponent) * 2.0 * cf;
513 csu = csu + pow(-1.0, tmp_exponent) * cf / k;
514 } else if (k > 1) {
515 csv = csv + pow(-1.0, tmp_exponent) * k / (k * k - 1.0) * cf;
516 }
517 cf2 = cf1;
518 cf1 = cf;
519 }
520 if (y0 <= 1.0)
521 cs0 = cbs + cf;
522 else
523 cs0 = (cbs + cf) / Kokkos::cos(z);
524 cbj0 = cbj0 / cs0;
525 } else { // Using asymptotic expansion (5.2.5) for |z|>joint_val
526 // (default:25)
527 CmplxType ct1 = z1 - 0.25 * pi;
528 CmplxType cp0 = CmplxType(1.0, 0.0);
529 for (int k = 1; k <= 12; k++) { // Calculate (5.2.9)
530 cp0 = cp0 + a[k - 1] * Kokkos::pow(z1, -2.0 * k);
531 }
532 CmplxType cq0 = -0.125 / z1;
533 for (int k = 1; k <= 12; k++) { // Calculate (5.2.10)
534 cq0 = cq0 + b[k - 1] * Kokkos::pow(z1, -2.0 * k - 1);
535 }
536 CmplxType cu = Kokkos::sqrt(r2p / z1);
537 cbj0 = cu * (cp0 * Kokkos::cos(ct1) - cq0 * Kokkos::sin(ct1));
538 }
539 }
540 return cbj0;
541}
542
545template <class CmplxType, class RealType, class IntType>
546KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_y0(const CmplxType& z,
547 const RealType& joint_val = 25,
548 const IntType& bw_start = 70) {
549 // This function is converted and modified from the corresponding Fortran
550 // program CJYNB in S. Zhang & J. Jin "Computation of Special Functions"
551 //(Wiley, 1996).
552 // Input : z --- Complex argument
553 // joint_val --- Joint point of abs(z) separating small and large
554 // argument regions
555 // bw_start --- Starting point for backward recurrence
556 // Output: cby0 --- Y0(z)
557 using Kokkos::fabs;
558 using Kokkos::pow;
559 using Kokkos::Experimental::infinity_v;
560 using Kokkos::numbers::egamma_v;
561 using Kokkos::numbers::pi_v;
562
563 constexpr auto inf = infinity_v<RealType>;
564
565 CmplxType cby0, cbj0;
566 constexpr auto pi = pi_v<RealType>;
567 constexpr auto el = egamma_v<RealType>;
568 const RealType a[12] = {
569 -0.703125e-01, 0.112152099609375e+00, -0.5725014209747314e+00,
570 0.6074042001273483e+01, -0.1100171402692467e+03, 0.3038090510922384e+04,
571 -0.1188384262567832e+06, 0.6252951493434797e+07, -0.4259392165047669e+09,
572 0.3646840080706556e+11, -0.3833534661393944e+13, 0.4854014686852901e+15};
573 const RealType b[12] = {0.732421875e-01, -0.2271080017089844e+00,
574 0.1727727502584457e+01, -0.2438052969955606e+02,
575 0.5513358961220206e+03, -0.1825775547429318e+05,
576 0.8328593040162893e+06, -0.5006958953198893e+08,
577 0.3836255180230433e+10, -0.3649010818849833e+12,
578 0.4218971570284096e+14, -0.5827244631566907e+16};
579
580 RealType r2p = 2.0 / pi;
581 RealType a0 = Kokkos::abs(z);
582 RealType y0 = fabs(z.imag());
583 CmplxType ci = CmplxType(0.0, 1.0);
584 CmplxType z1 = z;
585
586 if (a0 < 1e-100) { // Treat z=0 as a special case
587 cby0 = -CmplxType(inf, 0.0);
588 } else {
589 if (z.real() < 0.0) z1 = -z;
590 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
591 // (default:25)
592 CmplxType cbs = CmplxType(0.0, 0.0);
593 CmplxType csu = CmplxType(0.0, 0.0);
594 CmplxType csv = CmplxType(0.0, 0.0);
595 CmplxType cf2 = CmplxType(0.0, 0.0);
596 CmplxType cf1 = CmplxType(1e-100, 0.0);
597 CmplxType cf, cs0, ce;
598 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
599 // 70)
600 cf = 2.0 * (k + 1.0) / z * cf1 - cf2;
601 RealType tmp_exponent = static_cast<RealType>(k / 2);
602 if (k == 0) cbj0 = cf;
603 if ((k == 2 * (k / 2)) && (k != 0)) {
604 if (y0 <= 1.0)
605 cbs = cbs + 2.0 * cf;
606 else
607 cbs = cbs + pow(-1.0, tmp_exponent) * 2.0 * cf;
608 csu = csu + pow(-1.0, tmp_exponent) * cf / k;
609 } else if (k > 1) {
610 csv = csv + pow(-1.0, tmp_exponent) * k / (k * k - 1.0) * cf;
611 }
612 cf2 = cf1;
613 cf1 = cf;
614 }
615 if (y0 <= 1.0)
616 cs0 = cbs + cf;
617 else
618 cs0 = (cbs + cf) / Kokkos::cos(z);
619 cbj0 = cbj0 / cs0;
620 ce = Kokkos::log(z / 2.0) + el;
621 cby0 = r2p * (ce * cbj0 - 4.0 * csu / cs0);
622 } else { // Using asymptotic expansion (5.2.6) for |z|>joint_val
623 // (default:25)
624 CmplxType ct1 = z1 - 0.25 * pi;
625 CmplxType cp0 = CmplxType(1.0, 0.0);
626 for (int k = 1; k <= 12; k++) { // Calculate (5.2.9)
627 cp0 = cp0 + a[k - 1] * Kokkos::pow(z1, -2.0 * k);
628 }
629 CmplxType cq0 = -0.125 / z1;
630 for (int k = 1; k <= 12; k++) { // Calculate (5.2.10)
631 cq0 = cq0 + b[k - 1] * Kokkos::pow(z1, -2.0 * k - 1);
632 }
633 CmplxType cu = Kokkos::sqrt(r2p / z1);
634 cbj0 = cu * (cp0 * Kokkos::cos(ct1) - cq0 * Kokkos::sin(ct1));
635 cby0 = cu * (cp0 * Kokkos::sin(ct1) + cq0 * Kokkos::cos(ct1));
636
637 if (z.real() < 0.0) { // Apply (5.4.2)
638 if (z.imag() < 0.0) cby0 = cby0 - 2.0 * ci * cbj0;
639 if (z.imag() >= 0.0) cby0 = cby0 + 2.0 * ci * cbj0;
640 }
641 }
642 }
643 return cby0;
644}
645
648template <class CmplxType, class RealType, class IntType>
649KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_j1(const CmplxType& z,
650 const RealType& joint_val = 25,
651 const IntType& bw_start = 70) {
652 // This function is converted and modified from the corresponding Fortran
653 // program CJYNB in S. Zhang & J. Jin "Computation of Special Functions"
654 //(Wiley, 1996).
655 // Input : z --- Complex argument
656 // joint_val --- Joint point of abs(z) separating small and large
657 // argument regions
658 // bw_start --- Starting point for backward recurrence
659 // Output: cbj1 --- J1(z)
660 using Kokkos::fabs;
661 using Kokkos::pow;
662 using Kokkos::numbers::pi_v;
663
664 CmplxType cbj1;
665 constexpr auto pi = pi_v<RealType>;
666 const RealType a1[12] = {0.1171875e+00, -0.144195556640625e+00,
667 0.6765925884246826e+00, -0.6883914268109947e+01,
668 0.1215978918765359e+03, -0.3302272294480852e+04,
669 0.1276412726461746e+06, -0.6656367718817688e+07,
670 0.4502786003050393e+09, -0.3833857520742790e+11,
671 0.4011838599133198e+13, -0.5060568503314727e+15};
672 const RealType b1[12] = {
673 -0.1025390625e+00, 0.2775764465332031e+00, -0.1993531733751297e+01,
674 0.2724882731126854e+02, -0.6038440767050702e+03, 0.1971837591223663e+05,
675 -0.8902978767070678e+06, 0.5310411010968522e+08, -0.4043620325107754e+10,
676 0.3827011346598605e+12, -0.4406481417852278e+14, 0.6065091351222699e+16};
677
678 RealType r2p = 2.0 / pi;
679 RealType a0 = Kokkos::abs(z);
680 RealType y0 = fabs(z.imag());
681 CmplxType z1 = z;
682
683 if (a0 < 1e-100) { // Treat z=0 as a special case
684 cbj1 = CmplxType(0.0, 0.0);
685 } else {
686 if (z.real() < 0.0) z1 = -z;
687 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
688 // (default:25)
689 CmplxType cbs = CmplxType(0.0, 0.0);
690 CmplxType csu = CmplxType(0.0, 0.0);
691 CmplxType csv = CmplxType(0.0, 0.0);
692 CmplxType cf2 = CmplxType(0.0, 0.0);
693 CmplxType cf1 = CmplxType(1e-100, 0.0);
694 CmplxType cf, cs0;
695 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
696 // 70)
697 cf = 2.0 * (k + 1.0) / z * cf1 - cf2;
698 RealType tmp_exponent = static_cast<RealType>(k / 2);
699 if (k == 1) cbj1 = cf;
700 if ((k == 2 * (k / 2)) && (k != 0)) {
701 if (y0 <= 1.0)
702 cbs = cbs + 2.0 * cf;
703 else
704 cbs = cbs + pow(-1.0, tmp_exponent) * 2.0 * cf;
705 csu = csu + pow(-1.0, tmp_exponent) * cf / k;
706 } else if (k > 1) {
707 csv = csv + pow(-1.0, tmp_exponent) * k / (k * k - 1.0) * cf;
708 }
709 cf2 = cf1;
710 cf1 = cf;
711 }
712 if (y0 <= 1.0)
713 cs0 = cbs + cf;
714 else
715 cs0 = (cbs + cf) / Kokkos::cos(z);
716 cbj1 = cbj1 / cs0;
717 } else { // Using asymptotic expansion (5.2.5) for |z|>joint_val
718 // (default:25)
719 CmplxType ct2 = z1 - 0.75 * pi;
720 CmplxType cp1 = CmplxType(1.0, 0.0);
721 for (int k = 1; k <= 12; k++) { // Calculate (5.2.11)
722 cp1 = cp1 + a1[k - 1] * Kokkos::pow(z1, -2.0 * k);
723 }
724 CmplxType cq1 = 0.375 / z1;
725 for (int k = 1; k <= 12; k++) { // Calculate (5.2.12)
726 cq1 = cq1 + b1[k - 1] * Kokkos::pow(z1, -2.0 * k - 1);
727 }
728 CmplxType cu = Kokkos::sqrt(r2p / z1);
729 cbj1 = cu * (cp1 * Kokkos::cos(ct2) - cq1 * Kokkos::sin(ct2));
730
731 if (real(z) < 0.0) { // Apply (5.4.2)
732 cbj1 = -cbj1;
733 }
734 }
735 }
736 return cbj1;
737}
738
741template <class CmplxType, class RealType, class IntType>
742KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_y1(const CmplxType& z,
743 const RealType& joint_val = 25,
744 const IntType& bw_start = 70) {
745 // This function is converted and modified from the corresponding Fortran
746 // program CJYNB in S. Zhang & J. Jin "Computation of Special Functions"
747 //(Wiley, 1996).
748 // Input : z --- Complex argument
749 // joint_val --- Joint point of abs(z) separating small and large
750 // argument regions
751 // bw_start --- Starting point for backward recurrence
752 // Output: cby1 --- Y1(z)
753 using Kokkos::fabs;
754 using Kokkos::pow;
755 using Kokkos::Experimental::infinity_v;
756 using Kokkos::numbers::egamma_v;
757 using Kokkos::numbers::pi_v;
758
759 constexpr auto inf = infinity_v<RealType>;
760
761 CmplxType cby1, cbj0, cbj1, cby0;
762 constexpr auto pi = pi_v<RealType>;
763 constexpr auto el = egamma_v<RealType>;
764 const RealType a1[12] = {0.1171875e+00, -0.144195556640625e+00,
765 0.6765925884246826e+00, -0.6883914268109947e+01,
766 0.1215978918765359e+03, -0.3302272294480852e+04,
767 0.1276412726461746e+06, -0.6656367718817688e+07,
768 0.4502786003050393e+09, -0.3833857520742790e+11,
769 0.4011838599133198e+13, -0.5060568503314727e+15};
770 const RealType b1[12] = {
771 -0.1025390625e+00, 0.2775764465332031e+00, -0.1993531733751297e+01,
772 0.2724882731126854e+02, -0.6038440767050702e+03, 0.1971837591223663e+05,
773 -0.8902978767070678e+06, 0.5310411010968522e+08, -0.4043620325107754e+10,
774 0.3827011346598605e+12, -0.4406481417852278e+14, 0.6065091351222699e+16};
775
776 RealType r2p = 2.0 / pi;
777 RealType a0 = Kokkos::abs(z);
778 RealType y0 = fabs(z.imag());
779 CmplxType ci = CmplxType(0.0, 1.0);
780 CmplxType z1 = z;
781
782 if (a0 < 1e-100) { // Treat z=0 as a special case
783 cby1 = -CmplxType(inf, 0.0);
784 } else {
785 if (z.real() < 0.0) z1 = -z;
786 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
787 // (default:25)
788 CmplxType cbs = CmplxType(0.0, 0.0);
789 CmplxType csu = CmplxType(0.0, 0.0);
790 CmplxType csv = CmplxType(0.0, 0.0);
791 CmplxType cf2 = CmplxType(0.0, 0.0);
792 CmplxType cf1 = CmplxType(1e-100, 0.0);
793 CmplxType cf, cs0, ce;
794 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
795 // 70)
796 cf = 2.0 * (k + 1.0) / z * cf1 - cf2;
797 RealType tmp_exponent = static_cast<RealType>(k / 2);
798 if (k == 1) cbj1 = cf;
799 if (k == 0) cbj0 = cf;
800 if ((k == 2 * (k / 2)) && (k != 0)) {
801 if (y0 <= 1.0)
802 cbs = cbs + 2.0 * cf;
803 else
804 cbs = cbs + pow(-1.0, tmp_exponent) * 2.0 * cf;
805 csu = csu + pow(-1.0, tmp_exponent) * cf / k;
806 } else if (k > 1) {
807 csv = csv + pow(-1.0, tmp_exponent) * k / (k * k - 1.0) * cf;
808 }
809 cf2 = cf1;
810 cf1 = cf;
811 }
812 if (y0 <= 1.0)
813 cs0 = cbs + cf;
814 else
815 cs0 = (cbs + cf) / Kokkos::cos(z);
816 cbj0 = cbj0 / cs0;
817 ce = Kokkos::log(z / 2.0) + el;
818 cby0 = r2p * (ce * cbj0 - 4.0 * csu / cs0);
819 cbj1 = cbj1 / cs0;
820 cby1 = (cbj1 * cby0 - 2.0 / (pi * z)) / cbj0;
821 } else { // Using asymptotic expansion (5.2.5) for |z|>joint_val
822 // (default:25)
823 CmplxType ct2 = z1 - 0.75 * pi;
824 CmplxType cp1 = CmplxType(1.0, 0.0);
825 for (int k = 1; k <= 12; k++) { // Calculate (5.2.11)
826 cp1 = cp1 + a1[k - 1] * Kokkos::pow(z1, -2.0 * k);
827 }
828 CmplxType cq1 = 0.375 / z1;
829 for (int k = 1; k <= 12; k++) { // Calculate (5.2.12)
830 cq1 = cq1 + b1[k - 1] * Kokkos::pow(z1, -2.0 * k - 1);
831 }
832 CmplxType cu = Kokkos::sqrt(r2p / z1);
833 cbj1 = cu * (cp1 * Kokkos::cos(ct2) - cq1 * Kokkos::sin(ct2));
834 cby1 = cu * (cp1 * Kokkos::sin(ct2) + cq1 * Kokkos::cos(ct2));
835
836 if (z.real() < 0.0) { // Apply (5.4.2)
837 if (z.imag() < 0.0) cby1 = -(cby1 - 2.0 * ci * cbj1);
838 if (z.imag() >= 0.0) cby1 = -(cby1 + 2.0 * ci * cbj1);
839 }
840 }
841 }
842 return cby1;
843}
844
847template <class CmplxType, class RealType, class IntType>
848KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_i0(const CmplxType& z,
849 const RealType& joint_val = 25,
850 const IntType& bw_start = 70) {
851 // This function is converted and modified from the corresponding Fortran
852 // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
853 // Functions" (Wiley, 1996).
854 // Input : z --- Complex argument
855 // joint_val --- Joint point of abs(z) separating small and large
856 // argument regions
857 // bw_start --- Starting point for backward recurrence
858 // Output: cbi0 --- I0(z)
859 using Kokkos::numbers::pi_v;
860
861 CmplxType cbi0;
862 constexpr auto pi = pi_v<RealType>;
863 const RealType a[12] = {0.125,
864 7.03125e-2,
865 7.32421875e-2,
866 1.1215209960938e-1,
867 2.2710800170898e-1,
868 5.7250142097473e-1,
869 1.7277275025845e0,
870 6.0740420012735e0,
871 2.4380529699556e1,
872 1.1001714026925e2,
873 5.5133589612202e2,
874 3.0380905109224e3};
875
876 RealType a0 = Kokkos::abs(z);
877 CmplxType z1 = z;
878
879 if (a0 < 1e-100) { // Treat z=0 as a special case
880 cbi0 = CmplxType(1.0, 0.0);
881 } else {
882 if (z.real() < 0.0) z1 = -z;
883 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
884 // (default:25)
885 CmplxType cbs = CmplxType(0.0, 0.0);
886 // CmplxType csk0 = CmplxType(0.0,0.0);
887 CmplxType cf0 = CmplxType(0.0, 0.0);
888 CmplxType cf1 = CmplxType(1e-100, 0.0);
889 CmplxType cf, cs0;
890 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
891 // 70)
892 cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
893 if (k == 0) cbi0 = cf;
894 // if ((k == 2*(k/2)) && (k != 0)) {
895 // csk0 = csk0+4.0*cf/static_cast<RealType>(k);
896 //}
897 cbs = cbs + 2.0 * cf;
898 cf0 = cf1;
899 cf1 = cf;
900 }
901 cs0 = Kokkos::exp(z1) / (cbs - cf);
902 cbi0 = cbi0 * cs0;
903 } else { // Using asymptotic expansion (6.2.1) for |z|>joint_val
904 // (default:25)
905 CmplxType ca = Kokkos::exp(z1) / Kokkos::sqrt(2.0 * pi * z1);
906 cbi0 = CmplxType(1.0, 0.0);
907 CmplxType zr = 1.0 / z1;
908 for (int k = 1; k <= 12; k++) {
909 cbi0 = cbi0 + a[k - 1] * Kokkos::pow(zr, 1.0 * k);
910 }
911 cbi0 = ca * cbi0;
912 }
913 }
914 return cbi0;
915}
916
919template <class CmplxType, class RealType, class IntType>
920KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_k0(const CmplxType& z,
921 const RealType& joint_val = 9,
922 const IntType& bw_start = 30) {
923 // This function is converted and modified from the corresponding Fortran
924 // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
925 // Functions" (Wiley, 1996).
926 // Purpose: Compute modified Bessel function K0(z) of the second kind of
927 // order zero for a complex argument
928 // Input : z --- Complex argument
929 // joint_val --- Joint point of abs(z) separating small and large
930 // argument regions
931 // bw_start --- Starting point for backward recurrence
932 // Output: cbk0 --- K0(z)
933 using Kokkos::pow;
934 using Kokkos::Experimental::infinity_v;
935 using Kokkos::numbers::egamma_v;
936 using Kokkos::numbers::pi_v;
937
938 constexpr auto inf = infinity_v<RealType>;
939
940 CmplxType cbk0, cbi0;
941 constexpr auto pi = pi_v<RealType>;
942 constexpr auto el = egamma_v<RealType>;
943
944 RealType a0 = Kokkos::abs(z);
945 CmplxType ci = CmplxType(0.0, 1.0);
946 CmplxType z1 = z;
947
948 if (a0 < 1e-100) { // Treat z=0 as a special case
949 cbk0 = CmplxType(inf, 0.0);
950 } else {
951 if (z.real() < 0.0) z1 = -z;
952 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
953 // (default:9)
954 CmplxType cbs = CmplxType(0.0, 0.0);
955 CmplxType csk0 = CmplxType(0.0, 0.0);
956 CmplxType cf0 = CmplxType(0.0, 0.0);
957 CmplxType cf1 = CmplxType(1e-100, 0.0);
958 CmplxType cf, cs0;
959 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
960 // 30)
961 cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
962 if (k == 0) cbi0 = cf;
963 if ((k == 2 * (k / 2)) && (k != 0)) {
964 csk0 = csk0 + 4.0 * cf / static_cast<RealType>(k);
965 }
966 cbs = cbs + 2.0 * cf;
967 cf0 = cf1;
968 cf1 = cf;
969 }
970 cs0 = Kokkos::exp(z1) / (cbs - cf);
971 cbi0 = cbi0 * cs0;
972 cbk0 = -(Kokkos::log(0.5 * z1) + el) * cbi0 + cs0 * csk0;
973 } else { // Using asymptotic expansion (6.2.2) for |z|>joint_val
974 // (default:9)
975 CmplxType ca0 = Kokkos::sqrt(pi / (2.0 * z1)) * Kokkos::exp(-z1);
976 CmplxType cbkl = CmplxType(1.0, 0.0);
977 CmplxType cr = CmplxType(1.0, 0.0);
978 for (int k = 1; k <= 30; k++) {
979 cr = 0.125 * cr * (0.0 - pow(2.0 * k - 1.0, 2.0)) / (k * z1);
980 cbkl = cbkl + cr;
981 }
982 cbk0 = ca0 * cbkl;
983 }
984 if (z.real() < 0.0) { // Apply (6.4.4)
985 if (z.imag() < 0.0)
986 cbk0 = cbk0 + ci * pi * cyl_bessel_i0<CmplxType, RealType, IntType>(z);
987 if (z.imag() >= 0.0)
988 cbk0 = cbk0 - ci * pi * cyl_bessel_i0<CmplxType, RealType, IntType>(z);
989 }
990 }
991 return cbk0;
992}
993
996template <class CmplxType, class RealType, class IntType>
997KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_i1(const CmplxType& z,
998 const RealType& joint_val = 25,
999 const IntType& bw_start = 70) {
1000 // This function is converted and modified from the corresponding Fortran
1001 // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
1002 // Functions" (Wiley, 1996).
1003 // Input : z --- Complex argument
1004 // joint_val --- Joint point of abs(z) separating small and large
1005 // argument regions
1006 // bw_start --- Starting point for backward recurrence
1007 // Output: cbi1 --- I1(z)
1008 using Kokkos::numbers::pi_v;
1009
1010 CmplxType cbi1;
1011 constexpr auto pi = pi_v<RealType>;
1012 const RealType b[12] = {-0.375,
1013 -1.171875e-1,
1014 -1.025390625e-1,
1015 -1.4419555664063e-1,
1016 -2.7757644653320e-1,
1017 -6.7659258842468e-1,
1018 -1.9935317337513,
1019 -6.8839142681099,
1020 -2.7248827311269e1,
1021 -1.2159789187654e2,
1022 -6.0384407670507e2,
1023 -3.3022722944809e3};
1024
1025 RealType a0 = Kokkos::abs(z);
1026 CmplxType z1 = z;
1027
1028 if (a0 < 1e-100) { // Treat z=0 as a special case
1029 cbi1 = CmplxType(0.0, 0.0);
1030 } else {
1031 if (z.real() < 0.0) z1 = -z;
1032 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
1033 // (default:25)
1034 CmplxType cbs = CmplxType(0.0, 0.0);
1035 // CmplxType csk0 = CmplxType(0.0,0.0);
1036 CmplxType cf0 = CmplxType(0.0, 0.0);
1037 CmplxType cf1 = CmplxType(1e-100, 0.0);
1038 CmplxType cf, cs0;
1039 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
1040 // 70)
1041 cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
1042 if (k == 1) cbi1 = cf;
1043 // if ((k == 2*(k/2)) && (k != 0)) {
1044 // csk0 = csk0+4.0*cf/static_cast<RealType>(k);
1045 //}
1046 cbs = cbs + 2.0 * cf;
1047 cf0 = cf1;
1048 cf1 = cf;
1049 }
1050 cs0 = Kokkos::exp(z1) / (cbs - cf);
1051 cbi1 = cbi1 * cs0;
1052 } else { // Using asymptotic expansion (6.2.1) for |z|>joint_val
1053 // (default:25)
1054 CmplxType ca = Kokkos::exp(z1) / Kokkos::sqrt(2.0 * pi * z1);
1055 cbi1 = CmplxType(1.0, 0.0);
1056 CmplxType zr = 1.0 / z1;
1057 for (int k = 1; k <= 12; k++) {
1058 cbi1 = cbi1 + b[k - 1] * Kokkos::pow(zr, 1.0 * k);
1059 }
1060 cbi1 = ca * cbi1;
1061 }
1062 if (z.real() < 0.0) { // Apply (6.4.4)
1063 cbi1 = -cbi1;
1064 }
1065 }
1066 return cbi1;
1067}
1068
1071template <class CmplxType, class RealType, class IntType>
1072KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_k1(const CmplxType& z,
1073 const RealType& joint_val = 9,
1074 const IntType& bw_start = 30) {
1075 // This function is converted and modified from the corresponding Fortran
1076 // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
1077 // Functions" (Wiley, 1996).
1078 // Input : z --- Complex argument
1079 // joint_val --- Joint point of abs(z) separating small and large
1080 // argument regions
1081 // bw_start --- Starting point for backward recurrence
1082 // Output: cbk1 --- K1(z)
1083 using Kokkos::pow;
1084 using Kokkos::Experimental::infinity_v;
1085 using Kokkos::numbers::egamma_v;
1086 using Kokkos::numbers::pi_v;
1087
1088 constexpr auto inf = infinity_v<RealType>;
1089
1090 CmplxType cbk0, cbi0, cbk1, cbi1;
1091 constexpr auto pi = pi_v<RealType>;
1092 constexpr auto el = egamma_v<RealType>;
1093
1094 RealType a0 = Kokkos::abs(z);
1095 CmplxType ci = CmplxType(0.0, 1.0);
1096 CmplxType z1 = z;
1097
1098 if (a0 < 1e-100) { // Treat z=0 as a special case
1099 cbk1 = CmplxType(inf, 0.0);
1100 } else {
1101 if (z.real() < 0.0) z1 = -z;
1102 if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
1103 // (default:9)
1104 CmplxType cbs = CmplxType(0.0, 0.0);
1105 CmplxType csk0 = CmplxType(0.0, 0.0);
1106 CmplxType cf0 = CmplxType(0.0, 0.0);
1107 CmplxType cf1 = CmplxType(1e-100, 0.0);
1108 CmplxType cf, cs0;
1109 for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
1110 // 30)
1111 cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
1112 if (k == 1) cbi1 = cf;
1113 if (k == 0) cbi0 = cf;
1114 if ((k == 2 * (k / 2)) && (k != 0)) {
1115 csk0 = csk0 + 4.0 * cf / static_cast<RealType>(k);
1116 }
1117 cbs = cbs + 2.0 * cf;
1118 cf0 = cf1;
1119 cf1 = cf;
1120 }
1121 cs0 = Kokkos::exp(z1) / (cbs - cf);
1122 cbi0 = cbi0 * cs0;
1123 cbi1 = cbi1 * cs0;
1124 cbk0 = -(Kokkos::log(0.5 * z1) + el) * cbi0 + cs0 * csk0;
1125 cbk1 = (1.0 / z1 - cbi1 * cbk0) / cbi0;
1126 } else { // Using asymptotic expansion (6.2.2) for |z|>joint_val
1127 // (default:9)
1128 CmplxType ca0 = Kokkos::sqrt(pi / (2.0 * z1)) * Kokkos::exp(-z1);
1129 CmplxType cbkl = CmplxType(1.0, 0.0);
1130 CmplxType cr = CmplxType(1.0, 0.0);
1131 for (int k = 1; k <= 30; k++) {
1132 cr = 0.125 * cr * (4.0 - pow(2.0 * k - 1.0, 2.0)) / (k * z1);
1133 cbkl = cbkl + cr;
1134 }
1135 cbk1 = ca0 * cbkl;
1136 }
1137 if (z.real() < 0.0) { // Apply (6.4.4)
1138 if (z.imag() < 0.0)
1139 cbk1 = -cbk1 - ci * pi * cyl_bessel_i1<CmplxType, RealType, IntType>(z);
1140 if (z.imag() >= 0.0)
1141 cbk1 = -cbk1 + ci * pi * cyl_bessel_i1<CmplxType, RealType, IntType>(z);
1142 }
1143 }
1144 return cbk1;
1145}
1146
1149template <class CmplxType>
1150KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h10(const CmplxType& z) {
1151 // This function is converted and modified from the corresponding Fortran
1152 // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1153 //(Wiley, 1996).
1154 using RealType = typename CmplxType::value_type;
1155 using Kokkos::Experimental::infinity_v;
1156 using Kokkos::numbers::pi_v;
1157
1158 constexpr auto inf = infinity_v<RealType>;
1159
1160 CmplxType ch10, cbk0, cbj0, cby0;
1161 constexpr auto pi = pi_v<RealType>;
1162 CmplxType ci = CmplxType(0.0, 1.0);
1163
1164 if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1165 ch10 = CmplxType(1.0, -inf);
1166 } else if (z.imag() <= 0.0) {
1167 cbj0 = cyl_bessel_j0<CmplxType, RealType, int>(z);
1168 cby0 = cyl_bessel_y0<CmplxType, RealType, int>(z);
1169 ch10 = cbj0 + ci * cby0;
1170 } else { //(z.imag() > 0.0)
1171 cbk0 = cyl_bessel_k0<CmplxType, RealType, int>(-ci * z, 18.0, 70);
1172 ch10 = 2.0 / (pi * ci) * cbk0;
1173 }
1174
1175 return ch10;
1176}
1177
1180template <class CmplxType>
1181KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h11(const CmplxType& z) {
1182 // This function is converted and modified from the corresponding Fortran
1183 // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1184 //(Wiley, 1996).
1185 using RealType = typename CmplxType::value_type;
1186 using Kokkos::Experimental::infinity_v;
1187 using Kokkos::numbers::pi_v;
1188
1189 constexpr auto inf = infinity_v<RealType>;
1190
1191 CmplxType ch11, cbk1, cbj1, cby1;
1192 constexpr auto pi = pi_v<RealType>;
1193 CmplxType ci = CmplxType(0.0, 1.0);
1194
1195 if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1196 ch11 = CmplxType(0.0, -inf);
1197 } else if (z.imag() <= 0.0) {
1198 cbj1 = cyl_bessel_j1<CmplxType, RealType, int>(z);
1199 cby1 = cyl_bessel_y1<CmplxType, RealType, int>(z);
1200 ch11 = cbj1 + ci * cby1;
1201 } else { //(z.imag() > 0.0)
1202 cbk1 = cyl_bessel_k1<CmplxType, RealType, int>(-ci * z, 18.0, 70);
1203 ch11 = -2.0 / pi * cbk1;
1204 }
1205
1206 return ch11;
1207}
1208
1211template <class CmplxType>
1212KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h20(const CmplxType& z) {
1213 // This function is converted and modified from the corresponding Fortran
1214 // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1215 //(Wiley, 1996).
1216 using RealType = typename CmplxType::value_type;
1217 using Kokkos::Experimental::infinity_v;
1218 using Kokkos::numbers::pi_v;
1219
1220 constexpr auto inf = infinity_v<RealType>;
1221
1222 CmplxType ch20, cbk0, cbj0, cby0;
1223 constexpr auto pi = pi_v<RealType>;
1224 CmplxType ci = CmplxType(0.0, 1.0);
1225
1226 if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1227 ch20 = CmplxType(1.0, inf);
1228 } else if (z.imag() >= 0.0) {
1229 cbj0 = cyl_bessel_j0<CmplxType, RealType, int>(z);
1230 cby0 = cyl_bessel_y0<CmplxType, RealType, int>(z);
1231 ch20 = cbj0 - ci * cby0;
1232 } else { //(z.imag() < 0.0)
1233 cbk0 = cyl_bessel_k0<CmplxType, RealType, int>(ci * z, 18.0, 70);
1234 ch20 = 2.0 / pi * ci * cbk0;
1235 }
1236
1237 return ch20;
1238}
1239
1242template <class CmplxType>
1243KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h21(const CmplxType& z) {
1244 // This function is converted and modified from the corresponding Fortran
1245 // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1246 //(Wiley, 1996).
1247 using RealType = typename CmplxType::value_type;
1248 using Kokkos::Experimental::infinity_v;
1249 using Kokkos::numbers::pi_v;
1250
1251 constexpr auto inf = infinity_v<RealType>;
1252
1253 CmplxType ch21, cbk1, cbj1, cby1;
1254 constexpr auto pi = pi_v<RealType>;
1255 CmplxType ci = CmplxType(0.0, 1.0);
1256
1257 if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1258 ch21 = CmplxType(0.0, inf);
1259 } else if (z.imag() >= 0.0) {
1260 cbj1 = cyl_bessel_j1<CmplxType, RealType, int>(z);
1261 cby1 = cyl_bessel_y1<CmplxType, RealType, int>(z);
1262 ch21 = cbj1 - ci * cby1;
1263 } else { //(z.imag() < 0.0)
1264 cbk1 = cyl_bessel_k1<CmplxType, RealType, int>(ci * z, 18.0, 70);
1265 ch21 = -2.0 / pi * cbk1;
1266 }
1267
1268 return ch21;
1269}
1270
1271} // namespace Experimental
1272} // namespace Kokkos
1273
1274#ifdef KOKKOS_IMPL_PUBLIC_INCLUDE_NOTDEFINED_MATHSPECFUNCTIONS
1275#undef KOKKOS_IMPL_PUBLIC_INCLUDE
1276#undef KOKKOS_IMPL_PUBLIC_INCLUDE_NOTDEFINED_MATHSPECFUNCTIONS
1277#endif
1278#endif
A thread safe view to a bitset.