#ifndef BESSEL_H
#define BESSEL_H

namespace Distribution {

	/** helper class */
	template <typename T> class Bessel {

	public:

		/**
		 * calculation for I_v(z)
		 * https://en.wikipedia.org/wiki/Bessel_function
		 * http://www.boost.org/doc/libs/1_35_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/bessel/mbessel.html
		 */
		T getModified(const int v, const T z) const {

			// running factorials (updated within for-loops)
			uint64_t runFac1 = 1;
			uint64_t runFac2 = factorial(v);	// Delta(k+v+1) = (k+v+1-1)! = (k+v)! [k starts at 0] -> (v)!

			// running potential
			T runPot = 1;
			const T pot = T(0.25) * z * z;

			// start for-loop at k=1 instead of k=0. allows for running factorial without using if (k==0)
			T sum = 0;
			for (int k = 0; k < 16; ) {
				sum += runPot / runFac1 / runFac2;
				++k;
				runFac1 *= k;
				runFac2 *= k;
				runPot *= pot;
			}

			// done
			return std::pow( (T(0.5) * z), v) * sum;

		}

	private:





		static uint64_t factorial(const int n) {
			uint64_t res = 1;
			for (int i = 2; i <= n; ++i) {res *= i;}
			return res;
		}

	};

}

#endif // BESSEL_H