1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
from libc.math cimport pow, sqrt, floor, log, log1p, exp, M_PI, fabs
from numpy.math cimport NAN, isinf
cimport numpy as np
from _xlogy cimport xlogy
from _complexstuff cimport (
zsqrt, zpow, zabs, npy_cdouble_from_double_complex,
double_complex_from_npy_cdouble)
cdef extern from "float.h":
double DBL_MAX, DBL_MIN
cdef extern from "cephes.h":
double iv(double v, double x) nogil
double jv(double n, double x) nogil
double Gamma(double x) nogil
double lgam(double x) nogil
cdef extern from "c_misc/misc.h":
double gammasgn(double x) nogil
cdef extern from "amos_wrappers.h":
np.npy_cdouble cbesi_wrap(double v, np.npy_cdouble z) nogil
np.npy_cdouble cbesj_wrap(double v, np.npy_cdouble z) nogil
double sin_pi(double x) nogil
#
# Real-valued kernel
#
cdef inline double _hyp0f1_real(double v, double z) nogil:
cdef double arg, v1, arg_exp, bess_val
# handle poles, zeros
if v <= 0.0 and v == floor(v):
return NAN
if z == 0.0 and v != 0.0:
return 1.0
# both v and z small: truncate the Taylor series at O(z**2)
if fabs(z) < 1e-6*(1.0 + fabs(v)):
return 1.0 + z/v + z*z/(2.0*v*(v+1.0))
if z > 0:
arg = sqrt(z)
arg_exp = xlogy(1.0-v, arg) + lgam(v)
bess_val = iv(v-1, 2.0*arg)
if (arg_exp > log(DBL_MAX) or bess_val == 0 or # overflow
arg_exp < log(DBL_MIN) or isinf(bess_val)): # underflow
return _hyp0f1_asy(v, z)
else:
return exp(arg_exp) * gammasgn(v) * bess_val
else:
arg = sqrt(-z)
return pow(arg, 1.0 - v) * Gamma(v) * jv(v - 1, 2*arg)
cdef inline double _hyp0f1_asy(double v, double z) nogil:
r"""Asymptotic expansion for I_{v-1}(2*sqrt(z)) * Gamma(v)
for real $z > 0$ and $v\to +\infty$.
Based off DLMF 10.41
"""
cdef:
double arg = sqrt(z)
double v1 = fabs(v - 1)
double x = 2.0 * arg / v1
double p1 = sqrt(1.0 + x*x)
double eta = p1 + log(x) - log1p(p1)
double arg_exp_i, arg_exp_k
double pp, p2, p4, p6, u1, u2, u3, u_corr_i, u_corr_k
double result, gs
arg_exp_i = -0.5*log(p1)
arg_exp_i -= 0.5*log(2.0*M_PI*v1)
arg_exp_i += lgam(v)
gs = gammasgn(v)
arg_exp_k = arg_exp_i
arg_exp_i += v1 * eta
arg_exp_k -= v1 * eta
# large-v asymptotic correction, DLMF 10.41.10
pp = 1.0/p1
p2 = pp*pp
p4 = p2*p2
p6 = p4*p2
u1 = (3.0 - 5.0*p2) * pp / 24.0
u2 = (81.0 - 462.0*p2 + 385.0*p4) * p2 / 1152.0
u3 = (30375.0 - 369603.0*p2 + 765765.0*p4 - 425425.0*p6) * pp * p2 / 414720.0
u_corr_i = 1.0 + u1/v1 + u2/(v1*v1) + u3/(v1*v1*v1)
result = exp(arg_exp_i - xlogy(v1, arg)) * gs * u_corr_i
if v - 1 < 0:
# DLMF 10.27.2: I_{-v} = I_{v} + (2/pi) sin(pi*v) K_v
u_corr_k = 1.0 - u1/v1 + u2/(v1*v1) - u3/(v1*v1*v1)
result += exp(arg_exp_k + xlogy(v1, arg)) * gs * 2.0 * sin_pi(v1) * u_corr_k
return result
#
# Complex valued kernel
#
cdef inline double complex _hyp0f1_cmplx(double v, double complex z) nogil:
cdef:
np.npy_cdouble zz = npy_cdouble_from_double_complex(z)
np.npy_cdouble r
double complex arg, s
double complex t1, t2
# handle poles, zeros
if v <= 0.0 and v == floor(v):
return NAN
if z.real == 0.0 and z.imag == 0.0 and v != 0.0:
return 1.0
# both v and z small: truncate the Taylor series at O(z**2)
if zabs(z) < 1e-6*(1.0 + zabs(v)):
t1 = 1.0 + z/v
t2 = z*z / (2.0*v*(v+1.0))
return t1 + t2
if zz.real > 0:
arg = zsqrt(z)
s = 2.0 * arg
r = cbesi_wrap(v-1.0, npy_cdouble_from_double_complex(s))
else:
arg = zsqrt(-z)
s = 2.0 * arg
r = cbesj_wrap(v-1.0, npy_cdouble_from_double_complex(s))
return double_complex_from_npy_cdouble(r) * Gamma(v) * zpow(arg, 1.0 - v)
|
