From 2d4e1886c35e003cb7dc4e4f447ef8b31c2a53cb Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Sat, 3 Mar 2018 19:49:56 +0100
Subject: [PATCH] simplified i0() and added some LaTeX code

---
 window.hh | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/window.hh b/window.hh
index 396bf46..bff9d2b 100644
--- a/window.hh
+++ b/window.hh
@@ -89,6 +89,14 @@ template <int TAPS, typename TYPE>
 class Kaiser
 {
 	TYPE w[TAPS];
+	/*
+	i0() implements the zero-th order modified Bessel function of the first kind:
+	https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1
+	$I_\alpha(x) = i^{-\alpha} J_\alpha(ix) = \sum_{m=0}^\infty \frac{1}{m!\, \Gamma(m+\alpha+1)}\left(\frac{x}{2}\right)^{2m+\alpha}$
+	$I_0(x) = J_0(ix) = \sum_{m=0}^\infty \frac{1}{m!\, \Gamma(m+1)}\left(\frac{x}{2}\right)^{2m} = \sum_{m=0}^\infty \left(\frac{x^m}{2^m\,m!}\right)^{2}$
+	We obviously can't use the factorial here, so let's get rid of it:
+	$= 1 + \left(\frac{x}{2 \cdot 1}\right)^2 + \left(\frac{x}{2 \cdot 1}\cdot \frac{x}{2 \cdot 2}\right)^2 + \left(\frac{x}{2 \cdot 1}\cdot \frac{x}{2 \cdot 2}\cdot \frac{x}{2 \cdot 3}\right)^2 + .. = 1 + \sum_{m=1}^\infty \left(\prod_{n=1}^m \frac{x}{2n}\right)^2$
+	*/
 	TYPE i0(TYPE x)
 	{
 		Kahan<TYPE> sum(1.0);
@@ -97,8 +105,8 @@ class Kaiser
 		// float: 25 iterations
 		// double: 35 iterations
 		for (int n = 1; n < 35; ++n) {
-			TYPE tmp = x / TYPE(2 * n);
-			if (sum.same(val *= tmp * tmp))
+			val *= x / TYPE(2 * n);
+			if (sum.same(val * val))
 				return sum();
 		}
 		return sum();