diff options
author | Patrik Müller <patrik.mueller@ost.ch> | 2022-04-06 08:00:09 +0200 |
---|---|---|
committer | Patrik Müller <patrik.mueller@ost.ch> | 2022-04-06 08:00:09 +0200 |
commit | 670555039265d83945b0d3e205aefb020425585b (patch) | |
tree | 21038cb909c97aff834de131cf0a8c5e276cb756 /buch/papers/laguerre/scripts/lanczos_approximation.py | |
parent | add notes for MSE session 4 (diff) | |
download | SeminarSpezielleFunktionen-670555039265d83945b0d3e205aefb020425585b.tar.gz SeminarSpezielleFunktionen-670555039265d83945b0d3e205aefb020425585b.zip |
Start definition.tex
Diffstat (limited to 'buch/papers/laguerre/scripts/lanczos_approximation.py')
-rw-r--r-- | buch/papers/laguerre/scripts/lanczos_approximation.py | 47 |
1 files changed, 47 insertions, 0 deletions
diff --git a/buch/papers/laguerre/scripts/lanczos_approximation.py b/buch/papers/laguerre/scripts/lanczos_approximation.py new file mode 100644 index 0000000..3c48266 --- /dev/null +++ b/buch/papers/laguerre/scripts/lanczos_approximation.py @@ -0,0 +1,47 @@ +from cmath import exp, pi, sin, sqrt + +p = [ + 676.5203681218851, + -1259.1392167224028, + 771.32342877765313, + -176.61502916214059, + 12.507343278686905, + -0.13857109526572012, + 9.9843695780195716e-6, + 1.5056327351493116e-7, +] + +EPSILON = 1e-07 + + +def drop_imag(z): + if abs(z.imag) <= EPSILON: + z = z.real + return z + + +def gamma(z): + z = complex(z) + if z.real < 0.5: + y = pi / (sin(pi * z) * gamma(1 - z)) # Reflection formula + else: + z -= 1 + x = 0.99999999999980993 + for (i, pval) in enumerate(p): + x += pval / (z + i + 1) + t = z + len(p) - 0.5 + y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x + return drop_imag(y) + + +""" +The above use of the reflection (thus the if-else structure) is necessary, even though +it may look strange, as it allows to extend the approximation to values of z where +Re(z) < 0.5, where the Lanczos method is not valid. +""" + +print(gamma(1)) +print(gamma(5)) +print(gamma(0.5)) +print(gamma(0.5* (1 + 1j))) +print(gamma(-0.5)) |