Bessel_K: Adapted from Boost by David Ossim

CSharp · Console

/Program.cs

using System;
					
public class Program
{
	public static void Main()
	{
		Console.WriteLine(bessel_k(1,2));
	}
	
	public static double bessel_k(double v, double x)
        {
            if (x <= 0)
            {
                throw new ArgumentOutOfRangeException("x", "Value must be greater than 0");
            }
            return bessel_ik(v, x);
        }

        private static double bessel_ik(double v, double x)
        {
            double u, Kv, Kv1, Ku, Ku1, W, current, prev, next;
            bool reflect = false;
            uint n, k;
            if (v < 0)
            {
                reflect = true;
                v = -v;
            }
            n = (uint)Math.Round(v);
            u = v - n;
            W = 1 / x;
            CF2_ik(u, x, out Ku, out Ku1);
            prev = Ku;
            current = Ku1;
            double scale = 1;
            for (k = 1; k <= n; k++)
            {
                double fact = 2 * (u + k) / x;
                if ((Double.MaxValue - Math.Abs(prev)) / fact < Math.Abs(current))
                {
                    prev /= current;
                    scale /= current;
                    current = 1;
                }
                next = fact * current + prev;
                prev = current;
                current = next;
            }
            Kv = prev;
            Kv1 = current;
            if (reflect)
            {
                double z = (u + n % 2);
                double fact = (2 / Math.PI) * (sin_pi(z) * Kv);
                if (Double.MaxValue * scale < fact)
                {
                    throw new ArgumentOutOfRangeException("fact", "Error. Calculation exceeds Double.MaxValue");
                }
            }
            if (Double.MaxValue * scale < Kv)
            {
                throw new ArgumentOutOfRangeException("fact", "Error. Calculation exceeds Double.MaxValue");
            }
            return Kv / scale;
        }

        static double sin_pi(double x)
        {
            if (x < 0)
            {
                return -sin_pi(-x);
            }
            bool invert;
            if (x < 0.5)
            {
                return Math.Sin(Math.PI * x);
            }
            if (x < 1)
            {
                invert = true;
                x = -x;
            }
            else
            {
                invert = false;
            }
            double rem = Math.Floor(x);
            if (itrunc(rem) == 1)
            {
                invert = !invert;
            }
            rem = x - rem;
            if (rem > 0.5f)
            {
                rem = 1 - rem;
            }
            if (rem == 0.5f)
            {
                return (double)(invert ? -1 : 1);
            }
            rem = Math.Sin(Math.PI * rem);
            return invert ? -rem : rem;
        }

        static int itrunc(double v)
        {
            double r = Math.Truncate(v);
            if ((Math.Abs(r)) > (Int64.MaxValue))
            {
                throw new ArgumentOutOfRangeException("r", "Error. Calculation exceeds Int64.MaxValue");
            }
            return (int)(r);
        }

        static int CF2_ik(double v, double x, out double Kv, out  double Kv1)
        {
            double S, C, Q, D, f, a, b, q, delta, tolerance, current, prev;
            ulong k;
            tolerance = Double.Epsilon;
            a = v * v - 0.25f;
            b = 2 * (x + 1);
            D = 1 / b;
            f = delta = D;
            prev = 0;
            current = 1;
            Q = C = -a;
            S = 1 + Q * delta;
            for (k = 2; k < 50; k++)
            {
                a -= 2 * (k - 1);
                b += 2;
                D = 1 / (b + a * D);
                delta *= b * D - 1;
                f += delta;
                q = (prev - (b - 2) * current) / a;
                prev = current;
                current = q;
                C *= -a / k;
                Q += C * q;
                S += Q * delta;
                if (q < Double.Epsilon)
                {
                    C *= q;
                    prev /= q;
                    current /= q;
                    q = 1;
                }
                if (Math.Abs(Q * delta) < Math.Abs(S) * tolerance)
                {
                    break;
                }
            }
            if (x >= Math.Log(Double.MaxValue))
            {
                Kv = Math.Exp(0.5f * Math.Log(Math.PI / (2 * x)) - x - Math.Log(S));
            }
            else
            {
                Kv = Math.Sqrt(Math.PI / (2 * x)) * Math.Exp(-x) / S;
            }
            Kv1 = Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
            return 0;
        }
}