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
|
from pathlib import Path
import numpy as np
import scipy.special
EPSILON = 1e-7
root = str(Path(__file__).parent)
img_path = f"{root}/../images"
fontsize = "medium"
cmap = "plasma"
def _prep_zeros_and_weights(x, w, n):
if x is None or w is None:
return np.polynomial.laguerre.laggauss(n)
return x, w
def drop_imag(z):
if abs(z.imag) <= EPSILON:
z = z.real
return z
def pochhammer(z, n):
return np.prod(z + np.arange(n))
def find_optimal_shift(z, n):
mhat = 1.34093 * n + 0.854093
steps = int(np.floor(mhat - np.real(z)))
return steps
def get_shifting_factor(z, steps):
factor = 1.0
if steps > 0:
factor = 1 / pochhammer(z, steps)
elif steps < 0:
factor = pochhammer(z + steps, -steps)
return factor
def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
x, w = _prep_zeros_and_weights(x, w, n)
n = len(x)
z += 0j
steps = int(np.floor(target - np.real(z)))
z_shifted = z + steps
correction_factor = get_shifting_factor(z, steps)
res = np.sum(x ** (z_shifted - 1) * w)
res *= correction_factor
res = drop_imag(res)
return res
def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8):
if z == 0.0:
return np.infty
x, w = _prep_zeros_and_weights(x, w, n)
n = len(x)
z += 0j
opt_shift = find_optimal_shift(z, n)
correction_factor = get_shifting_factor(z, opt_shift)
z_shifted = z + opt_shift
res = np.sum(x ** (z_shifted - 1) * w)
res *= correction_factor
res = drop_imag(res)
return res
def laguerre_gamma_simple(z, x=None, w=None, n=8):
if z == 0.0:
return np.infty
x, w = _prep_zeros_and_weights(x, w, n)
z += 0j
res = np.sum(x ** (z - 1) * w)
res = drop_imag(res)
return res
def laguerre_gamma_mirror(z, x=None, w=None, n=8):
if z == 0.0:
return np.infty
x, w = _prep_zeros_and_weights(x, w, n)
z += 0j
if z.real < 1e-3:
return np.pi / (
np.sin(np.pi * z) * laguerre_gamma_simple(1 - z, x, w)
) # Reflection formula
return laguerre_gamma_simple(z, x, w)
def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
x, w = _prep_zeros_and_weights(x, w, n)
if func == "simple":
f = laguerre_gamma_simple
elif func == "mirror":
f = laguerre_gamma_mirror
elif func == "optimal_shifted":
f = laguerre_gamma_opt_shifted
else:
f = laguerre_gamma_shifted
return np.array([f(zi, x, w, n, **kwargs) for zi in z])
def calc_rel_error(x, y):
mask = np.abs(x) != np.infty
rel_error = np.zeros_like(y)
rel_error[mask] = (y[mask] - x[mask]) / x[mask]
rel_error[~mask] = 0.0
return rel_error
|
