From 670555039265d83945b0d3e205aefb020425585b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Patrik=20M=C3=BCller?= Date: Wed, 6 Apr 2022 08:00:09 +0200 Subject: Start definition.tex --- buch/papers/laguerre/scripts/quadrature_gama.py | 178 ++++++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 buch/papers/laguerre/scripts/quadrature_gama.py (limited to 'buch/papers/laguerre/scripts/quadrature_gama.py') diff --git a/buch/papers/laguerre/scripts/quadrature_gama.py b/buch/papers/laguerre/scripts/quadrature_gama.py new file mode 100644 index 0000000..37a9cd8 --- /dev/null +++ b/buch/papers/laguerre/scripts/quadrature_gama.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +"""Use Gauss-Laguerre quadrature to calculate gamma function.""" +# import sympy +from cmath import exp, pi, sin, sqrt + +import matplotlib.pyplot as plt +import numpy as np +import scipy.special as ss + +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) + + +zeros = np.array( + [ + 3.22547689619392312e-1, + 1.74576110115834658e0, + 4.53662029692112798e0, + 9.39507091230113313e0, + ], + np.longdouble, +) +weights = np.array( + [ + 6.03154104341633602e-1, + 3.57418692437799687e-1, + 3.88879085150053843e-2, + 5.39294705561327450e-4, + ], + np.longdouble, +) + +zeros = np.array( + [ + 1.70279632305101000e-1, + 9.03701776799379912e-1, + 2.25108662986613069e0, + 4.26670017028765879e0, + 7.04590540239346570e0, + 1.07585160101809952e1, + 1.57406786412780046e1, + 2.28631317368892641e1, + ], + np.longdouble, +) + +weights = np.array( + [ + 3.69188589341637530e-1, + 4.18786780814342956e-1, + 1.75794986637171806e-1, + 3.33434922612156515e-2, + 2.79453623522567252e-3, + 9.07650877335821310e-5, + 8.48574671627253154e-7, + 1.04800117487151038e-9, + ], + np.longdouble, +) + + +def calc_gamma(z, n, x, w): + res = 0.0 + z = complex(z) + for xi, wi in zip(x, w): + res += xi ** (z + n - 1) * wi + for i in range(int(n)): + res /= z + i + res = drop_imag(res) + return res + +small = 1e-3 +Z = np.linspace(small, 1-small, 101) + +# Z = [-3/2, -1/2, 1/2, 3/2] +# target = +# targets = np.array([gamma(z) for z in Z]) +targets1 = ss.gamma(Z) +targets2 = np.array([gamma(z) for z in Z]) +approxs = np.array([calc_gamma(z, 11, zeros, weights) for z in Z]) +rel_error1 = np.abs(targets1 - approxs) / targets1 +rel_error2 = np.abs(targets2 - approxs) / targets2 + +_, axs = plt.subplots(2, num=1, clear=True, constrained_layout=True) +axs[0].plot(Z, rel_error1) +axs[1].semilogy(Z, rel_error1) +axs[0].plot(Z, rel_error2) +axs[1].semilogy(Z, rel_error2) +axs[1].semilogy(Z, np.abs(targets1-targets2)/targets1) +plt.show() +# values = np.array([calc_gamma]) +# _ = [ +# print( +# n, +# [ +# float( +# f"{np.abs((calc_gamma(z, n, zeros, weights) - gamma(z)) / gamma(z)):.3g}" +# ) +# for z in Z +# ], +# ) +# for n in range(21) +# ] + + +# target = ss.gamma(z) +# target = np.sqrt(np.pi) + +# _, ax = plt.subplots(num=1, clear=True, constrained_layout=True) +# for i, degree in enumerate(degrees): +# samples_points, weights = np.polynomial.laguerre.laggauss(degree) +# values = np.sum( +# samples_points[:, None] ** (z + shifts[None] - 1) * weights[:, None], 0 +# ) / ss.poch(z, shifts) +# # print(np.abs(target - values)) +# print(values) +# ax.plot(shifts, values, label=f"N={degree}") +# ax.legend() +# plt.show() + + +# def count_equal_digits(x, y): +# for i in range(1, 13): +# try: +# np.testing.assert_almost_equal(x, y, i) +# except AssertionError: +# break +# return i + + +# Z = np.linspace(1.0, 11.0, 11) +# # degrees = [2, 4, 8, 16, 32, 64, 100] +# d = 100 +# X = np.zeros(len(Z)) +# for i, z in enumerate(Z): +# samples_points, weights = np.polynomial.laguerre.laggauss(d) +# X[i] = np.sum(samples_points ** (z - 1) * weights) +# # X[i] = np.sum(np.sin(z * samples_points) * weights) +# Y = ss.gamma(Z) +# # Y = Z / (Z ** 2 + 1) +# ed = [count_equal_digits(x, y) for x, y in zip(X, Y)] +# for x,y in zip(X,Y): +# print(x,y) + +# _, ax = plt.subplots(num=1, clear=True, constrained_layout=True) +# ax.plot(Z, ed) +# plt.show() -- cgit v1.2.1