Kokkos Core Kernels Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 
33 namespace Kokkos {
34 namespace Experimental {
35 
37 template <class RealType>
38 KOKKOS_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 
77 template <class RealType>
78 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> erf(
79  const Kokkos::complex<RealType>& z) {
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 
255 template <class RealType>
256 KOKKOS_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 
444 template <class RealType>
445 KOKKOS_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 
456 template <class CmplxType, class RealType, class IntType>
457 KOKKOS_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 
545 template <class CmplxType, class RealType, class IntType>
546 KOKKOS_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 
648 template <class CmplxType, class RealType, class IntType>
649 KOKKOS_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 
741 template <class CmplxType, class RealType, class IntType>
742 KOKKOS_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 
847 template <class CmplxType, class RealType, class IntType>
848 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_i0(const CmplxType& z,
849  const RealType& joint_val = 18,
850  const IntType& n_terms = 50) {
851  // This function is converted and modified from the corresponding Fortran
852  // programs 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  // n_terms --- Numbers of terms used in the power series
858  // Output: cbi0 --- I0(z)
859 
860  CmplxType cbi0(1.0, 0.0);
861  RealType a0 = Kokkos::abs(z);
862  CmplxType z1 = z;
863 
864  if (a0 > 1e-100) {
865  if (z.real() < 0.0) z1 = -z;
866  if (a0 <= joint_val) {
867  // Using power series definition for |z|<=joint_val (default:18)
868  CmplxType cr = CmplxType(1.0e+00, 0.0e+00);
869  CmplxType z2 = z * z;
870  for (int k = 1; k < n_terms; ++k) {
871  cr = RealType(.25) * cr * z2 / CmplxType(k * k);
872  cbi0 += cr;
873  if (Kokkos::abs(cr / cbi0) < RealType(1.e-15)) continue;
874  }
875  } else {
876  // Using asymptotic expansion (6.2.1) for |z|>joint_val (default:18)
877  const RealType a[12] = {0.125,
878  7.03125e-2,
879  7.32421875e-2,
880  1.1215209960938e-1,
881  2.2710800170898e-1,
882  5.7250142097473e-1,
883  1.7277275025845e0,
884  6.0740420012735e0,
885  2.4380529699556e1,
886  1.1001714026925e2,
887  5.5133589612202e2,
888  3.0380905109224e3};
889 
890  for (int k = 1; k <= 12; k++) {
891  cbi0 += a[k - 1] * Kokkos::pow(z1, -k);
892  }
893  cbi0 *= Kokkos::exp(z1) /
894  Kokkos::sqrt(2.0 * Kokkos::numbers::pi_v<RealType> * z1);
895  }
896  }
897  return cbi0;
898 }
899 
902 template <class CmplxType, class RealType, class IntType>
903 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_k0(const CmplxType& z,
904  const RealType& joint_val = 9,
905  const IntType& bw_start = 30) {
906  // This function is converted and modified from the corresponding Fortran
907  // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
908  // Functions" (Wiley, 1996).
909  // Purpose: Compute modified Bessel function K0(z) of the second kind of
910  // order zero for a complex argument
911  // Input : z --- Complex argument
912  // joint_val --- Joint point of abs(z) separating small and large
913  // argument regions
914  // bw_start --- Starting point for backward recurrence
915  // Output: cbk0 --- K0(z)
916  using Kokkos::pow;
917  using Kokkos::Experimental::infinity_v;
918  using Kokkos::numbers::egamma_v;
919  using Kokkos::numbers::pi_v;
920 
921  constexpr auto inf = infinity_v<RealType>;
922 
923  CmplxType cbk0, cbi0;
924  constexpr auto pi = pi_v<RealType>;
925  constexpr auto el = egamma_v<RealType>;
926 
927  RealType a0 = Kokkos::abs(z);
928  CmplxType ci = CmplxType(0.0, 1.0);
929  CmplxType z1 = z;
930 
931  if (a0 < 1e-100) { // Treat z=0 as a special case
932  cbk0 = CmplxType(inf, 0.0);
933  } else {
934  if (z.real() < 0.0) z1 = -z;
935  if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
936  // (default:9)
937  CmplxType cbs = CmplxType(0.0, 0.0);
938  CmplxType csk0 = CmplxType(0.0, 0.0);
939  CmplxType cf0 = CmplxType(0.0, 0.0);
940  CmplxType cf1 = CmplxType(1e-100, 0.0);
941  CmplxType cf, cs0;
942  for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
943  // 30)
944  cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
945  if (k == 0) cbi0 = cf;
946  if ((k == 2 * (k / 2)) && (k != 0)) {
947  csk0 = csk0 + 4.0 * cf / static_cast<RealType>(k);
948  }
949  cbs = cbs + 2.0 * cf;
950  cf0 = cf1;
951  cf1 = cf;
952  }
953  cs0 = Kokkos::exp(z1) / (cbs - cf);
954  cbi0 = cbi0 * cs0;
955  cbk0 = -(Kokkos::log(0.5 * z1) + el) * cbi0 + cs0 * csk0;
956  } else { // Using asymptotic expansion (6.2.2) for |z|>joint_val
957  // (default:9)
958  CmplxType ca0 = Kokkos::sqrt(pi / (2.0 * z1)) * Kokkos::exp(-z1);
959  CmplxType cbkl = CmplxType(1.0, 0.0);
960  CmplxType cr = CmplxType(1.0, 0.0);
961  for (int k = 1; k <= 30; k++) {
962  cr = 0.125 * cr * (0.0 - pow(2.0 * k - 1.0, 2.0)) / (k * z1);
963  cbkl = cbkl + cr;
964  }
965  cbk0 = ca0 * cbkl;
966  }
967  if (z.real() < 0.0) { // Apply (6.4.4)
968  if (z.imag() < 0.0)
969  cbk0 = cbk0 + ci * pi * cyl_bessel_i0<CmplxType, RealType, IntType>(z);
970  if (z.imag() >= 0.0)
971  cbk0 = cbk0 - ci * pi * cyl_bessel_i0<CmplxType, RealType, IntType>(z);
972  }
973  }
974  return cbk0;
975 }
976 
979 template <class CmplxType, class RealType, class IntType>
980 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_i1(const CmplxType& z,
981  const RealType& joint_val = 25,
982  const IntType& bw_start = 70) {
983  // This function is converted and modified from the corresponding Fortran
984  // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
985  // Functions" (Wiley, 1996).
986  // Input : z --- Complex argument
987  // joint_val --- Joint point of abs(z) separating small and large
988  // argument regions
989  // bw_start --- Starting point for backward recurrence
990  // Output: cbi1 --- I1(z)
991  using Kokkos::numbers::pi_v;
992 
993  CmplxType cbi1;
994  constexpr auto pi = pi_v<RealType>;
995  const RealType b[12] = {-0.375,
996  -1.171875e-1,
997  -1.025390625e-1,
998  -1.4419555664063e-1,
999  -2.7757644653320e-1,
1000  -6.7659258842468e-1,
1001  -1.9935317337513,
1002  -6.8839142681099,
1003  -2.7248827311269e1,
1004  -1.2159789187654e2,
1005  -6.0384407670507e2,
1006  -3.3022722944809e3};
1007 
1008  RealType a0 = Kokkos::abs(z);
1009  CmplxType z1 = z;
1010 
1011  if (a0 < 1e-100) { // Treat z=0 as a special case
1012  cbi1 = CmplxType(0.0, 0.0);
1013  } else {
1014  if (z.real() < 0.0) z1 = -z;
1015  if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
1016  // (default:25)
1017  CmplxType cbs = CmplxType(0.0, 0.0);
1018  // CmplxType csk0 = CmplxType(0.0,0.0);
1019  CmplxType cf0 = CmplxType(0.0, 0.0);
1020  CmplxType cf1 = CmplxType(1e-100, 0.0);
1021  CmplxType cf, cs0;
1022  for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
1023  // 70)
1024  cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
1025  if (k == 1) cbi1 = cf;
1026  // if ((k == 2*(k/2)) && (k != 0)) {
1027  // csk0 = csk0+4.0*cf/static_cast<RealType>(k);
1028  //}
1029  cbs = cbs + 2.0 * cf;
1030  cf0 = cf1;
1031  cf1 = cf;
1032  }
1033  cs0 = Kokkos::exp(z1) / (cbs - cf);
1034  cbi1 = cbi1 * cs0;
1035  } else { // Using asymptotic expansion (6.2.1) for |z|>joint_val
1036  // (default:25)
1037  CmplxType ca = Kokkos::exp(z1) / Kokkos::sqrt(2.0 * pi * z1);
1038  cbi1 = CmplxType(1.0, 0.0);
1039  CmplxType zr = 1.0 / z1;
1040  for (int k = 1; k <= 12; k++) {
1041  cbi1 = cbi1 + b[k - 1] * Kokkos::pow(zr, 1.0 * k);
1042  }
1043  cbi1 = ca * cbi1;
1044  }
1045  if (z.real() < 0.0) { // Apply (6.4.4)
1046  cbi1 = -cbi1;
1047  }
1048  }
1049  return cbi1;
1050 }
1051 
1054 template <class CmplxType, class RealType, class IntType>
1055 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_k1(const CmplxType& z,
1056  const RealType& joint_val = 9,
1057  const IntType& bw_start = 30) {
1058  // This function is converted and modified from the corresponding Fortran
1059  // programs CIKNB and CIK01 in S. Zhang & J. Jin "Computation of Special
1060  // Functions" (Wiley, 1996).
1061  // Input : z --- Complex argument
1062  // joint_val --- Joint point of abs(z) separating small and large
1063  // argument regions
1064  // bw_start --- Starting point for backward recurrence
1065  // Output: cbk1 --- K1(z)
1066  using Kokkos::pow;
1067  using Kokkos::Experimental::infinity_v;
1068  using Kokkos::numbers::egamma_v;
1069  using Kokkos::numbers::pi_v;
1070 
1071  constexpr auto inf = infinity_v<RealType>;
1072 
1073  CmplxType cbk0, cbi0, cbk1, cbi1;
1074  constexpr auto pi = pi_v<RealType>;
1075  constexpr auto el = egamma_v<RealType>;
1076 
1077  RealType a0 = Kokkos::abs(z);
1078  CmplxType ci = CmplxType(0.0, 1.0);
1079  CmplxType z1 = z;
1080 
1081  if (a0 < 1e-100) { // Treat z=0 as a special case
1082  cbk1 = CmplxType(inf, 0.0);
1083  } else {
1084  if (z.real() < 0.0) z1 = -z;
1085  if (a0 <= joint_val) { // Using backward recurrence for |z|<=joint_val
1086  // (default:9)
1087  CmplxType cbs = CmplxType(0.0, 0.0);
1088  CmplxType csk0 = CmplxType(0.0, 0.0);
1089  CmplxType cf0 = CmplxType(0.0, 0.0);
1090  CmplxType cf1 = CmplxType(1e-100, 0.0);
1091  CmplxType cf, cs0;
1092  for (int k = bw_start; k >= 0; k--) { // Backward recurrence (default:
1093  // 30)
1094  cf = 2.0 * (k + 1.0) * cf1 / z1 + cf0;
1095  if (k == 1) cbi1 = cf;
1096  if (k == 0) cbi0 = cf;
1097  if ((k == 2 * (k / 2)) && (k != 0)) {
1098  csk0 = csk0 + 4.0 * cf / static_cast<RealType>(k);
1099  }
1100  cbs = cbs + 2.0 * cf;
1101  cf0 = cf1;
1102  cf1 = cf;
1103  }
1104  cs0 = Kokkos::exp(z1) / (cbs - cf);
1105  cbi0 = cbi0 * cs0;
1106  cbi1 = cbi1 * cs0;
1107  cbk0 = -(Kokkos::log(0.5 * z1) + el) * cbi0 + cs0 * csk0;
1108  cbk1 = (1.0 / z1 - cbi1 * cbk0) / cbi0;
1109  } else { // Using asymptotic expansion (6.2.2) for |z|>joint_val
1110  // (default:9)
1111  CmplxType ca0 = Kokkos::sqrt(pi / (2.0 * z1)) * Kokkos::exp(-z1);
1112  CmplxType cbkl = CmplxType(1.0, 0.0);
1113  CmplxType cr = CmplxType(1.0, 0.0);
1114  for (int k = 1; k <= 30; k++) {
1115  cr = 0.125 * cr * (4.0 - pow(2.0 * k - 1.0, 2.0)) / (k * z1);
1116  cbkl = cbkl + cr;
1117  }
1118  cbk1 = ca0 * cbkl;
1119  }
1120  if (z.real() < 0.0) { // Apply (6.4.4)
1121  if (z.imag() < 0.0)
1122  cbk1 = -cbk1 - ci * pi * cyl_bessel_i1<CmplxType, RealType, IntType>(z);
1123  if (z.imag() >= 0.0)
1124  cbk1 = -cbk1 + ci * pi * cyl_bessel_i1<CmplxType, RealType, IntType>(z);
1125  }
1126  }
1127  return cbk1;
1128 }
1129 
1132 template <class CmplxType>
1133 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h10(const CmplxType& z) {
1134  // This function is converted and modified from the corresponding Fortran
1135  // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1136  //(Wiley, 1996).
1137  using RealType = typename CmplxType::value_type;
1138  using Kokkos::Experimental::infinity_v;
1139  using Kokkos::numbers::pi_v;
1140 
1141  constexpr auto inf = infinity_v<RealType>;
1142 
1143  CmplxType ch10, cbk0, cbj0, cby0;
1144  constexpr auto pi = pi_v<RealType>;
1145  CmplxType ci = CmplxType(0.0, 1.0);
1146 
1147  if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1148  ch10 = CmplxType(1.0, -inf);
1149  } else if (z.imag() <= 0.0) {
1150  cbj0 = cyl_bessel_j0<CmplxType, RealType, int>(z);
1151  cby0 = cyl_bessel_y0<CmplxType, RealType, int>(z);
1152  ch10 = cbj0 + ci * cby0;
1153  } else { //(z.imag() > 0.0)
1154  cbk0 = cyl_bessel_k0<CmplxType, RealType, int>(-ci * z, 18.0, 70);
1155  ch10 = 2.0 / (pi * ci) * cbk0;
1156  }
1157 
1158  return ch10;
1159 }
1160 
1163 template <class CmplxType>
1164 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h11(const CmplxType& z) {
1165  // This function is converted and modified from the corresponding Fortran
1166  // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1167  //(Wiley, 1996).
1168  using RealType = typename CmplxType::value_type;
1169  using Kokkos::Experimental::infinity_v;
1170  using Kokkos::numbers::pi_v;
1171 
1172  constexpr auto inf = infinity_v<RealType>;
1173 
1174  CmplxType ch11, cbk1, cbj1, cby1;
1175  constexpr auto pi = pi_v<RealType>;
1176  CmplxType ci = CmplxType(0.0, 1.0);
1177 
1178  if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1179  ch11 = CmplxType(0.0, -inf);
1180  } else if (z.imag() <= 0.0) {
1181  cbj1 = cyl_bessel_j1<CmplxType, RealType, int>(z);
1182  cby1 = cyl_bessel_y1<CmplxType, RealType, int>(z);
1183  ch11 = cbj1 + ci * cby1;
1184  } else { //(z.imag() > 0.0)
1185  cbk1 = cyl_bessel_k1<CmplxType, RealType, int>(-ci * z, 18.0, 70);
1186  ch11 = -2.0 / pi * cbk1;
1187  }
1188 
1189  return ch11;
1190 }
1191 
1194 template <class CmplxType>
1195 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h20(const CmplxType& z) {
1196  // This function is converted and modified from the corresponding Fortran
1197  // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1198  //(Wiley, 1996).
1199  using RealType = typename CmplxType::value_type;
1200  using Kokkos::Experimental::infinity_v;
1201  using Kokkos::numbers::pi_v;
1202 
1203  constexpr auto inf = infinity_v<RealType>;
1204 
1205  CmplxType ch20, cbk0, cbj0, cby0;
1206  constexpr auto pi = pi_v<RealType>;
1207  CmplxType ci = CmplxType(0.0, 1.0);
1208 
1209  if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1210  ch20 = CmplxType(1.0, inf);
1211  } else if (z.imag() >= 0.0) {
1212  cbj0 = cyl_bessel_j0<CmplxType, RealType, int>(z);
1213  cby0 = cyl_bessel_y0<CmplxType, RealType, int>(z);
1214  ch20 = cbj0 - ci * cby0;
1215  } else { //(z.imag() < 0.0)
1216  cbk0 = cyl_bessel_k0<CmplxType, RealType, int>(ci * z, 18.0, 70);
1217  ch20 = 2.0 / pi * ci * cbk0;
1218  }
1219 
1220  return ch20;
1221 }
1222 
1225 template <class CmplxType>
1226 KOKKOS_INLINE_FUNCTION CmplxType cyl_bessel_h21(const CmplxType& z) {
1227  // This function is converted and modified from the corresponding Fortran
1228  // programs CH12N in S. Zhang & J. Jin "Computation of Special Functions"
1229  //(Wiley, 1996).
1230  using RealType = typename CmplxType::value_type;
1231  using Kokkos::Experimental::infinity_v;
1232  using Kokkos::numbers::pi_v;
1233 
1234  constexpr auto inf = infinity_v<RealType>;
1235 
1236  CmplxType ch21, cbk1, cbj1, cby1;
1237  constexpr auto pi = pi_v<RealType>;
1238  CmplxType ci = CmplxType(0.0, 1.0);
1239 
1240  if ((z.real() == 0.0) && (z.imag() == 0.0)) {
1241  ch21 = CmplxType(0.0, inf);
1242  } else if (z.imag() >= 0.0) {
1243  cbj1 = cyl_bessel_j1<CmplxType, RealType, int>(z);
1244  cby1 = cyl_bessel_y1<CmplxType, RealType, int>(z);
1245  ch21 = cbj1 - ci * cby1;
1246  } else { //(z.imag() < 0.0)
1247  cbk1 = cyl_bessel_k1<CmplxType, RealType, int>(ci * z, 18.0, 70);
1248  ch21 = -2.0 / pi * cbk1;
1249  }
1250 
1251  return ch21;
1252 }
1253 
1254 } // namespace Experimental
1255 } // namespace Kokkos
1256 
1257 #ifdef KOKKOS_IMPL_PUBLIC_INCLUDE_NOTDEFINED_MATHSPECFUNCTIONS
1258 #undef KOKKOS_IMPL_PUBLIC_INCLUDE
1259 #undef KOKKOS_IMPL_PUBLIC_INCLUDE_NOTDEFINED_MATHSPECFUNCTIONS
1260 #endif
1261 #endif
Kokkos::complex
Partial reimplementation of std::complex that works as the result of a Kokkos::parallel_reduce.
Definition: Kokkos_Complex.hpp:42
Kokkos::complex::real
KOKKOS_INLINE_FUNCTION constexpr RealType & real() noexcept
The real part of this complex number.
Definition: Kokkos_Complex.hpp:139
Generated on Fri Apr 26 2024 09:20:01 for Kokkos Core Kernels Package by   doxygen 1.8.5