Skip to content

Commit

Permalink
adding equation numbering / derivations of gammanu
Browse files Browse the repository at this point in the history
  • Loading branch information
@hatemhelal
hatemhelal committed Oct 9, 2023
1 parent c571c6d commit 9e53f74
Show file tree
Hide file tree
Showing 2 changed files with 82 additions and 15 deletions.
95 changes: 81 additions & 14 deletions notebooks/gammanu.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,10 @@
"The evaluation of integrals involving a $\\frac{1}{r}$ term require computing the following:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{1}\n",
"F(\\nu, t) = \\int_0^1 u^{2 \\nu} e^{-t u^2} du\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Taketa et al. call this the auxiliary function $F$ and we can show its relationship to the $\\gamma$ function with SymPy.\n",
Expand All @@ -25,7 +28,7 @@
},
{
"cell_type": "code",
"execution_count": 1,
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -40,7 +43,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 3,
"metadata": {},
"outputs": [
{
Expand All @@ -52,7 +55,7 @@
"Integral(u**(2*nu)*exp(-t*u**2), (u, 0, 1))"
]
},
"execution_count": 2,
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
Expand Down Expand Up @@ -381,11 +384,11 @@
"From this, we propose the following recurrence:\n",
"\n",
"$$\n",
"F(0, t) = \\frac{1}{2} \\text{erf}(\\sqrt{t}) \\sqrt{\\frac{\\pi}{t}}\n",
"$$\n",
"\n",
"$$\n",
"F(\\nu, t) = \\frac{1}{2 t} \\left( (2\\nu - 1) F(\\nu - 1, t) - e^{-t} \\right)\n",
"\\begin{aligned}\n",
"\\tag{2}\n",
"F(0, t) &= \\frac{1}{2} \\text{erf}(\\sqrt{t}) \\sqrt{\\frac{\\pi}{t}} \\\\\n",
"F(\\nu, t) &= \\frac{1}{2 t} \\left( (2\\nu - 1) F(\\nu - 1, t) - e^{-t} \\right)\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"We test this recurrence out by comparing to the earlier lambda implementation derived by SymPy"
Expand Down Expand Up @@ -493,36 +496,57 @@
"\n",
"1. The lower incomplete gamma function of $a$ and $z$\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{3}\n",
"\\gamma(a, z) = \\int_0^z t^{a-1} e^{-t} dt\n",
"\\end{equation}\n",
"$$ \n",
"\n",
"2. normalised lower incomplete gamma function of $a$ and $z$. This definition is consistent with the [scipy.special.gammainc](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gammainc.html) implementation\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{4}\n",
"P(a, z) = \\frac{\\gamma(a, z)}{\\Gamma(a)} = \\frac{1}{\\Gamma(a)} \\int_0^z t^{a-1} e^{-t} dt\n",
"\\end{equation}\n",
"$$\n",
"\n",
"3. upper incomplete gamma function of $a$ and $z$.\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{5}\n",
"\\Gamma(a, z) = \\int_z^\\infty t^{a-1} e^{-t} dt\n",
"\\end{equation}\n",
"$$\n",
"\n",
"4. normalised upper incomplete gamma function of $a$ and $z$. This definition is consisent with the [scipy.special.gammaincc](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gammaincc.html) implementation\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{6}\n",
"Q(a, z) = \\frac{\\Gamma(a, z)}{\\Gamma(a)} = \\frac{1}{\\Gamma(a)} \\int_z^\\infty t^{a-1} e^{-t} dt\n",
"\\end{equation}\n",
"$$\n",
"\n",
"An additional identiy that follows from the integral definition:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{7}\n",
"\\gamma(a, z) + \\Gamma(a, z) = \\Gamma(a)\n",
"\\end{equation}\n",
"$$\n",
"and for the normalised functions we have:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{8}\n",
"P(a, z) + Q(a, z) = 1\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Back to the boost libary, see [#10](https://beta.boost.org/doc/libs/1_82_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html#math_toolkit.sf_gamma.igamma.implementation) of their implementation notes which states for half-integers $a \\in [\\frac{1}{2}, 30]$ that the following series expansion should be used to evaluate the upper incomplete regularised gamma function:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{9}\n",
"Q(a, x) = \\text{erfc}(\\sqrt{x}) + \\frac{e^{-x}}{\\sqrt{\\pi x}} \\sum_{n=1}^{a - \\frac{1}{2}} \\frac{x^n}{(1 - \\frac{1}{2}) \\cdots (n - \\frac{1}{2})}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"\n"
Expand Down Expand Up @@ -609,19 +633,31 @@
"\n",
"The special case above appears to have excessive errors as `a` increases. We can instead try the standard series expansion method:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{10}\n",
"\\gamma(a, z) = z^a e^{-z} \\sum_{k=0}^{\\infty} \\frac{\\Gamma(a)}{\\Gamma(a + k + 1)} z^k\n",
"\\end{equation}\n",
"$$\n",
"See [#4](https://beta.boost.org/doc/libs/1_82_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html#math_toolkit.sf_gamma.igamma.implementation) on the boost documentation which also shows how this series can be simplfied to:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{11}\n",
"\\gamma(a, z) = z^a e^{-z} \\sum_{k=0}^{\\infty} \\frac{z^k}{a^{\\overline{k + 1}}}\n",
"\\end{equation}\n",
"$$\n",
"by using the recursion relation:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{12}\n",
"\\Gamma(z+1) = z \\Gamma(z)\n",
"\\end{equation}\n",
"$$\n",
"where we use $a^{\\overline{k + 1}}$ as the [Pochhammer symbol](https://en.wikipedia.org/wiki/Falling_and_rising_factorials) evaluated as:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{13}\n",
"a^{\\overline{k + 1}} = a(a + 1)(a+2) \\cdots (a + k)\n",
"\\end{equation}\n",
"$$\n",
"\n",
"This can be used to build a log-space implementation of the lower incomplete gamma function."
Expand Down Expand Up @@ -909,34 +945,65 @@
"\n",
"Remember from earlier that goal is to evaluate:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{14}\n",
"F(\\nu, t) = \\int_0^1 u^{2 \\nu} e^{-t u^2} du\n",
"\\end{equation}\n",
"$$\n",
"which we used SymPy to show that\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{15}\n",
"F(\\nu, t) = \\frac{1}{2} t^{-\\nu - \\frac{1}{2}} \\gamma(\\nu + \\frac{1}{2}, t)\n",
"\\end{equation}\n",
"$$\n",
"now combining this result with the series representation for the lower incomplete gamma function:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{16}\n",
"\\gamma(\\nu + \\frac{1}{2}, t) = t^{\\nu + \\frac{1}{2}} e^{-t} \\sum_{k=0}^{\\infty} \\frac{t^k}{(\\nu + \\frac{1}{2})^{\\overline{k+1}}}\n",
"\\end{equation}\n",
"$$\n",
"we arrive at:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{17}\n",
"F(\\nu, t) = \\frac{e^{-t}}{2} \\sum_{k=0}^{\\infty} \\frac{t^k}{(\\nu + \\frac{1}{2})^{\\overline{k+1}}}\n",
"$$"
"\\end{equation}\n",
"$$\n",
"recall that $a^{\\overline{k + 1}}$ is the [Pochhammer symbol](https://en.wikipedia.org/wiki/Falling_and_rising_factorials) evaluated as:\n",
"$$\n",
"\\begin{equation}\n",
"\\tag{18}\n",
"a^{\\overline{k + 1}} = a(a + 1)(a+2) \\cdots (a + k)\n",
"\\end{equation}\n",
"$$\n",
"which leads to a simple recurrence that can be used to evaluate the series:\n",
"$$\n",
"\\begin{aligned}\n",
"\\tag{19}\n",
"a &\\equiv \\nu + \\frac{1}{2} \\\\ \n",
"f_0 &= \\frac{1}{a}\\\\\n",
"f_k &= \\frac{t}{a + 1} f_{k-1} \\\\\n",
"F(\\nu, t) &= \\frac{e^{-t}}{2} \\sum_{k=0}^{\\infty} f_k \n",
"\\end{aligned}\n",
"$$\n",
"This series will converge as $k \\rightarrow \\infty$ since the denominator terms will grow faster then the numerator. \n",
"In practice the number of terms needed to converge should be expected to be similar to the study above for the convergence of the lower incomplete gamma function in log space.\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0267741724495078e-07"
"3.331779296901205e-13"
]
},
"execution_count": 17,
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -958,8 +1025,8 @@
" return np.exp(-t) / 2 * total\n",
"\n",
"\n",
"nu = 20.0\n",
"t = 20.0\n",
"nu = 0.0\n",
"t = 0.0\n",
"\n",
"abs((gammanu(nu, t) - Fgamma(nu, t)) / Fgamma(nu, t))"
]
Expand Down
2 changes: 1 addition & 1 deletion pyscf_ipu/experimental/special.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ def binom_lookup(x: IntN, y: IntN, nmax: int = LMAX) -> IntN:

def gammanu(nu: IntN, t: FloatN, num_terms: int = 128) -> FloatN:
"""
eq 2.11 from THO but simplified as derived in gammanu.ipynb
eq 2.11 from THO but simplified as derived in equation 19 of gammanu.ipynb
"""
an = nu + 0.5
tn = 1 / an
Expand Down

0 comments on commit 9e53f74

Please sign in to comment.