```1 /**
2 Authors: [Ilya Yaroshenko](http://9il.github.io)
3
5
7 */
8 module atmosphere.math;
9
10 //import core.stdc.tgmath;
11 import std.range;
12 import std.traits;
13 import std.typecons;
14 import std.math;
15
16
17 version (LDC)
18 {
20 		T llvm_fmuladd(T)(T vala, T valb, T valc) @safe pure nothrow @nogc;
21 }
22
23 /**
24 Computes accurate sum of binary logarithms of input range \$(D r).
25  */
26 ElementType!Range sumOfLog2s(Range)(Range r)
27 	if (isInputRange!Range && isFloatingPoint!(ElementType!Range))
28 {
29 	import std.math: frexp;
30 	long exp = 0;
31 	Unqual!(typeof(return)) x = 1;
32 	foreach (e; r)
33 	{
34 		if (e < 0)
35 			return typeof(return).nan;
36 		int lexp = void;
37 		x *= frexp(e, lexp);
38 		exp += lexp;
39 		if (x < 0.5)
40 		{
41 			x *= 2;
42 			exp--;
43 		}
44 	}
45 	return exp + log2(x);
46 }
47
48 //Issue 14231
49 package T findRoot(T, DF, DT)(scope DF f, in T a, in T b,
50 	scope DT tolerance) //= (T a, T b) => false)
51 	if(
52 		isFloatingPoint!T &&
53 		is(typeof(tolerance(T.init, T.init)) : bool) &&
54 		is(typeof(f(T.init)) == R, R) && isFloatingPoint!R
55 	)
56 {
57 	import std.numeric: findRoot;
58 	import std.functional;
59 	immutable fa = f(a);
60 	if(fa == 0)
61 		return a;
62 	immutable fb = f(b);
63 	if(fb == 0)
64 		return b;
65 	immutable r = findRoot!(Unqual!T, Unqual!T)(f.toDelegate, a, b, fa, fb, tolerance.toDelegate);
66 	// Return the first value if it is smaller or NaN
67 	return !(fabs(r[2]) > fabs(r[3])) ? r[0] : r[1];
68 }
69
70 //Issue 14231
71 package T findRoot(T, DF)(scope DF f, in T a, in T b)
72 {
73 	return findRoot(f, a, b, (T a, T b) => false);
74 }
75
76 /** Inverse of the Log Minus Digamma function
77  *
78  *   Returns x such \$(D log(x) - digamma(x) == y).
79  *
80  * References:
81  *   1. Abramowitz, M., and Stegun, I. A. (1970).
82  *	 Handbook of mathematical functions. Dover, New York,
83  *	 pages 258-259, equation 6.3.18.
84  *
85  * Authors: Ilya Yaroshenko
86  */
87 T logmdigammaInverse(T)(const T y)
88 {
89 	//Issue 14231
90 	//import std.numeric: findRoot;
91 	enum maxY = 1 / T.min_normal;
92 	static assert(maxY > 0 && maxY <= T.max);
93 	if (y >= maxY)
94 		return 1 / y; //lim x->0 (log(x)-digamma(x))*x == 1
95 	if (y < 0)
96 		return T.nan;
97 	if (y < 2*T.min_normal)
98 		return 0.5 / y; //6.3.18
99 	if (y > 0)
100 		// x/2 <= logmdigamma(1 / x) <= x, x > 0
101 		// calls logmdigamma ~6 times
102 		return 1 / findRoot((T x) => logmdigamma(1 / x) - y, y,  2*y); // ~6 times
103 	return y; //NaN
104 }
105
106
107 unittest {
108 	alias T = double;
109 	import std.range : iota;
110 	import std.math : approxEqual, nextDown;
111 	//import std.mathspecial: logmdigamma;
112 	foreach (x; iota(1.3L, 2L, 0.1L))
113 	{
114 		assert(logmdigammaInverse(logmdigamma(x)).approxEqual(x));
115 		assert(logmdigammaInverse(logmdigamma(1/x)).approxEqual(1/x));
116 	}
117 	//WolframAlpha, 22.02.2015
118 	immutable T[2][5] testData = [
119 		[1.0L, 0.615556766479594378978099158335549201923L],
120 		[1.0L/8, 4.15937801516894947161054974029150730555L],
121 		[1.0L/1024, 512.166612384991507850643277924243523243L],
122 		[0.000500083333325000003968249801594877323784632117L, 1000.0L],
123 		[1017.644138623741168814449776695062817947092468536L, 1.0L/1024],
124 	];
125 	foreach(test; testData)
126 	{
127 		assert(approxEqual(logmdigammaInverse(test[0]), test[1], 1e-10, 0));
128 	}
129 	assert(approxEqual(logmdigamma(logmdigammaInverse(1.0)), 1, 1e-15, 0));
130 	assert(approxEqual(logmdigamma(logmdigammaInverse(T.min_normal)), T.min_normal, 1e-15, 0));
131 	assert(approxEqual(logmdigamma(logmdigammaInverse(T.max/2)), T.max/2, 1e-15, 0));
132 	assert(approxEqual(logmdigammaInverse(logmdigamma(1.0)), 1, 1e-10, 0));
133 	assert(approxEqual(logmdigammaInverse(logmdigamma(T.min_normal)), T.min_normal, 1e-15, 0));
134 	assert(approxEqual(logmdigammaInverse(logmdigamma(T.max/2)), T.max/2, 1e-15, 0));
135 }
136
137 T logmdigamma(T)(const T x)
138 {
139 	static immutable T [7] Bn_n  = [
140 	1.0L/(6*2), -1.0L/(30*4), 1.0L/(42*6), -1.0L/(30*8),
141 	5.0L/(66*10), -691.0L/(2730*12), 7.0L/(6*14) ];
142
143 	if (x <= 0.0)
144 	{
145 		if (x == 0.0)
146 		{
147 		   return T.infinity;
148 	}
149 		return T.nan;
150 	}
151
152 	T s = x;
153 	T w = 0f;
154 	while ( s < 10f ) {
155 		w += 1f/s;
156 		s += 1f;
157 	}
158 	T y;
159 	if ( s < 1.0e17 ) {
160 		immutable T z = 1.0/(s * s);
161 		version(LDC)
162 		{
163 			y = z *
170 				Bn_n[6] ,
171 				Bn_n[5]),
172 				Bn_n[4]),
173 				Bn_n[3]),
174 				Bn_n[2]),
175 				Bn_n[1]),
176 				Bn_n[0]);
177 		}
178 		else
179 		{
180 			y =   z * (Bn_n[0]
181 				+ z * (Bn_n[1]
182 				+ z * (Bn_n[2]
183 				+ z * (Bn_n[3]
184 				+ z * (Bn_n[4]
185 				+ z * (Bn_n[5]
186 				+ z *  Bn_n[6]))))));
187 		}
188 	} else
189 		y = 0f;
190
191 	return x == s ? y + 0.5f/s : (log(x/s) + 0.5f/s + y + w);
192 }
193
194 unittest {
195 	import std.typetuple;
196 	import std.math: isNaN, isIdentical, NaN, approxEqual;
197 	import std.mathspecial: digamma;
198 	foreach(T; TypeTuple!(real, double, float))
199 	{
200 		alias lmd = logmdigamma!T;
201 		assert(lmd(-5.0).isNaN());
202 		assert(isIdentical(lmd(NaN(0xABC)), NaN(0xABC)));
203 		assert(lmd(0.0) == T.infinity);
204 		for(auto x = 0.01; x < 1.0; x += 0.1)
205 			assert(approxEqual(digamma(x), log(x) - lmd(x)));
206 		for(auto x = 1.0; x < 15.0; x += 1.0)
207 			assert(approxEqual(digamma(x), log(x) - lmd(x)));
208 	}
209 }
210
211
212 /++
213 Find a real minimum of a real function \$(D f(x)) via bracketing.
214 Given a function \$(D f) and a range \$(D (ax..bx)),
215 returns the value of \$(D x) in the range which is closest to a minimum of \$(D f(x)).
216 \$(D f) is never evaluted at the endpoints of \$(D ax) and \$(D bx).
217 If \$(D f(x)) has more than one minimum in the range, one will be chosen arbitrarily.
218 If \$(D f(x)) returns NaN or -Infinity, \$(D (x, f(x), NaN)) will be returned;
219 otherwise, this algorithm is guaranteed to succeed.
220
221 Params:
222 	f = Function to be analyzed
223 	ax = Left bound of initial range of f known to contain the minimum.
224 	bx = Right bound of initial range of f known to contain the minimum.
225 	relTolerance = Relative tolerance.
226 	absTolerance = Absolute tolerance.
227
228 Preconditions:
229 	\$(D ax) and \$(D bx) shall be finite reals. \$(BR)
230 	\$(D relTolerance) shall be normal positive real. \$(BR)
231 	\$(D absTolerance) shall be normal positive real no less then \$(D T.epsilon*2).
232
233 Returns:
234 	A tuple consisting of \$(D x), \$(D y = f(x)) and \$(D error = 3 * (absTolerance * fabs(x) + relTolerance)).
235
236 	The method used is a combination of golden section search and
237 successive parabolic interpolation. Convergence is never much slower
238 than that for a Fibonacci search.
239
240 References:
241 	"Algorithms for Minimization without Derivatives", Richard Brent, Prentice-Hall, Inc. (1973)
242
243 See_Also: \$(LREF findRoot), \$(XREF math, isNormal)
244 +/
245
246 Tuple!(T, "x", Unqual!(ReturnType!DF), "y", T, "error")
247 findLocalMin(T, DF)(
248 		scope DF f,
249 		in T ax,
250 		in T bx,
251 		in T relTolerance = T.epsilon^^0.25,
252 		in T absTolerance = sqrt(T.epsilon),
253 		)
254 	if(isFloatingPoint!T)
255 in
256 {
257 	assert(isFinite(ax) && isFinite(bx), "ax and bx shall be finite reals");
258 	assert(isNormal(absTolerance) && absTolerance >= T.epsilon*2,
259 		"absTolerance shall be normal positive real no less then \$(D T.epsilon*2)");
260 	assert(isNormal(relTolerance) && relTolerance > 0,
261 		"relTolerance shall be normal positive real");
262 }
263 out(result)
264 {
265 	assert(isFinite(result.x));
266 }
267 body
268 {
269 	alias R = Unqual!(CommonType!(ReturnType!DF, T));
270
271 	// c is the squared inverse of the golden ratio
272 	// (3 - sqrt(5))/2
273 	// Value obtained from Wolfram Alpha.
274 	immutable T c = 0x0.61c8864680b583ea0c633f9fa31237p+0L;
275 	immutable T cm1 = 0x0.9e3779b97f4a7c15f39cc0605cedc8p+0L;
276
277 	R tolerance;
278
279 	T a = ax > bx ? bx : ax;
280 	T b = ax > bx ? ax : bx;
281
282 	//sequence of declarations suitable for SIMD instructions
283
284 	T  v = a * cm1 + b * c;
285 	assert(isFinite(v));
286
287 	R fv = f(v);
288
289 	//R fv = f(v);
290 	if(isNaN(fv) || fv == -T.infinity)
291 		return typeof(return)(v, fv, T.init);
292
293 	T  w = v;
294 	R fw = fv;
295
296 	T  x = v;
297 	R fx = fv;
298
299 	for(R d = 0, e = 0;;)
300 	{
301 		T m = (a + b) / 2;
302
303 		// This fix is not part of original algorithm
304 		if(!isFinite(m)) // fix infinity loop. Issue can be reproduced in R.
305 		{
306 			m = a / 2 + b / 2;
307 			if(!isFinite(m)) // fast-math compiler switch is enabled
308 			{
309 				//SIMD instructions can be used by compiler, do not reduce declarations
310 				int a_exp = void;
311 				int b_exp = void;
312 				immutable an = frexp(a, a_exp);
313 				immutable bn = frexp(b, b_exp);
314 				immutable am = ldexp(an, a_exp-1);
315 				immutable bm = ldexp(bn, b_exp-1);
316 				m = am + bm;
317 				if(!isFinite(m)) // wrong input: constraints is disabled in release mode
318 				{
319 					return typeof(return).init;
320 				}
321 			}
322 		}
323 		assert(isFinite(m));
324 		tolerance = absTolerance * fabs(x) + relTolerance;
325 		immutable t2 = tolerance * 2;
326 		//assert(x + tolerance < b);
327 		//assert(x - tolerance > a);
328
329
330 		// check stopping criterion
331 		if (!(fabs(x - m) > t2 - (b - a) / 2))
332 			break;
333
334 		R p = 0;
335 		R q = 0;
336 		R r = 0;
337
338 		// fit parabola
339 		if (fabs(e) > tolerance)
340 		{
341 			immutable  xw =  x -  w;
342 			immutable fxw = fx - fw;
343 			immutable  xv =  x -  v;
344 			immutable fxv = fx - fv;
345 			immutable xwfxv = xw * fxv;
346 			immutable xvfxw = xv * fxw;
347 			p = xv * xvfxw - xw * xwfxv;
348 			q = (xvfxw - xwfxv) * 2;
349 			if (q > 0)
350 				p = -p;
351 			else
352 				q = -q;
353 			r = e;
354 			e = d;
355 		}
356
357 		T u;
358 		// a parabolic-interpolation step
359 		if (fabs(p) < fabs(q * r / 2) && p > q * (a - x) && p < q * (b - x))
360 		{
361 			d = p / q;
362 			u = x + d;
363 			// f must not be evaluated too close to a or b
364 			if (u - a < t2 || b - u < t2)
365 				d = x < m ? tolerance : -tolerance;
366 		}
367 		// a golden-section step
368 		else
369 		{
370 			e = (x < m ? b : a) - x;
371 			d = c * e;
372 		}
373
374 		// f must not be evaluated too close to x
375 		u = x + (fabs(d) >= tolerance ? d : d > 0 ? tolerance : -tolerance);
376 		immutable fu = f(u);
377 		if(isNaN(fu) || fu == -T.infinity)
378 			return typeof(return)(u, fu, T.init);
379
380 		//  update  a, b, v, w, and x
381 		if (fu <= fx)
382 		{
383 			u < x ? b : a = x;
384 			v = w; fv = fw;
385 			w = x; fw = fx;
386 			x = u; fx = fu;
387 		}
388 		else
389 		{
390 			u < x ? a : b = u;
391 			if (fu <= fw || w == x)
392 			{
393 				v = w; fv = fw;
394 				w = u; fw = fu;
395 			}
396 			else if (fu <= fv || v == x || v == w)
397 			{ // do not remove this braces
398 				v = u; fv = fu;
399 			}
400 		}
401 	}
402 	return typeof(return)(x, fx, tolerance * 3);
403 }
404
405 ///
406 unittest
407 {
408 	auto ret = findLocalMin((double x) => (x-4)^^2, -1e9, 1e9);
409 	assert(ret.x.approxEqual(4.0));
410 	assert(ret.y.approxEqual(0.0));
411 }
412
413 unittest
414 {
415 	int i;
416 	double f(double x)
417 	{
418 		i++;
419 		return fabs(x-1);
420 	}
421
422 	//slow
423 	auto minD = findLocalMin(&f, -double.max/2, double.max/2, double.min_normal, 2*double.epsilon);
424 	with(minD)
425 	{
426 		assert(approxEqual(x, 1));
427 		assert(approxEqual(y, 0));
428 		assert(error <= 10 * double.epsilon);
429 		//assert(minD.error == 0);
430 	}
431 }
432
433 unittest
434 {
435 	auto ret = findLocalMin((double x) => double.init, 0.0, 1.0, double.min_normal, 2*double.epsilon);
436 	assert(!ret.x.isNaN);
437 	assert(ret.y.isNaN);
438 	assert(ret.error.isNaN);
439 }
440
441 unittest
442 {
443 	auto ret = findLocalMin((double x) => log(x), 0.0, 1.0, double.min_normal, 2*double.epsilon);
444 	assert(ret.error < 3.00001 * ((2*double.epsilon)*fabs(ret.x)+ double.min_normal));
445 	assert(ret.x >= 0 && ret.x <= ret.error);
446 }
447
448
449 unittest
450 {
451 	size_t i;
452 	auto ret = findLocalMin((double x) {i++; return log(x);}, 0.0, double.max, double.min_normal, 2*double.epsilon);
453 	assert(ret.y < -18);
454 	assert(ret.error < 5e-08);
455 	assert(ret.x >= 0 && ret.x <= ret.error);
456 }
457
458
459 unittest
460 {
461 	size_t i;
462 	auto ret = findLocalMin((double x) {i++; return -fabs(x);}, -1.0, 1.0, double.min_normal, 2*double.epsilon);
463 	assert(ret.x.fabs.approxEqual(1));
464 	assert(ret.y.fabs.approxEqual(1));
465 	assert(ret.error.approxEqual(0));
466 }
467
468 unittest
469 {
470 	size_t i;
471 	auto ret = findLocalMin((double x) {i++; return -fabs(x);}, -1.0, 1.0, double.min_normal, 2*double.epsilon);
472 	assert(ret.x.fabs.approxEqual(1));
473 	assert(ret.y.fabs.approxEqual(1));
474 	assert(ret.error.approxEqual(0));
475 }
476
477
478 version(none)
479 // check speed regressions and infinity loops
480 unittest
481 {
482 	import std.typetuple;
483 	import std.math;
484 	size_t i = 0;
485 	R f00(T, R)(T x) { i++; return 0; }
486 	R f01(T, R)(T x) { i++; return 1; }
487 	R f02(T, R)(T x) { i++; return R.min_normal/2; } //subnormal
488 	R f03(T, R)(T x) { i++; return R(PI); }
489 	R f04(T, R)(T x) { i++; return log2(cast(R)x); }
490 	R f05(T, R)(T x) { i++; return log(cast(R)x); }
491 	R f06(T, R)(T x) { i++; return exp2(cast(R)x); }
492 	R f07(T, R)(T x) { i++; return exp(cast(R)x); }
493 	R f08(T, R)(T x) { i++; return sqrt(cast(R)x); }
494 	R f09(T, R)(T x) { i++; return cbrt(cast(R)x); }
495 	R f10(T, R)(T x) { i++; return x; }
496 	R f11(T, R)(T x) { i++; return cast(R)(x)^^2; }
497 	R f12(T, R)(T x) { i++; return cast(R)(x)^^PI; }
498 	R f13(T, R)(T x) { i++; return cast(R)(x)^^(1/PI); }
499 	R f14(T, R)(T x) { i++; return sin(cast(R)x); }
500 	R f15(T, R)(T x) { i++; return cos(cast(R)x); }
501 	R f16(T, R)(T x) { i++; return sqrt(abs(cast(R)(x)^^2 - 1)); }
502 	R f17(T, R)(T x) { i++; return sqrt(abs(cast(R)(x)^^2 - 1)); }
503 	R f18(T, R)(T x) { i++; return floor(exp(cast(R)x)); } //multiminimum
504 	R f19(T, R)(T x) { i++; return floor(log(cast(R)x)); } //multiminimum
505 	R f20(T, R)(T x) { i++; return floor(cast(R)x); }      //multiminimum
506 	// vars for global checks
507 	int s1, s2;
508 	int c1, c2, c;
509 	foreach(T; TypeTuple!(
510 		float,
511 		double,
512 		real,
513 		)) // 1
514 	foreach(R; TypeTuple!(
515 		float,
516 		double,
517 		real,
518 		)) // 2
519 	{
520 		immutable ar1 = [
521 			&f00!(T, R),
522 			&f01!(T, R),
523 			&f02!(T, R),
524 			&f03!(T, R),
525 			&f04!(T, R),
526 			&f05!(T, R),
527 			&f06!(T, R),
528 			&f07!(T, R),
529 			&f08!(T, R),
530 			&f09!(T, R),
531 			&f10!(T, R),
532 			&f11!(T, R),
533 			&f12!(T, R),
534 			&f13!(T, R),
535 			&f14!(T, R),
536 			&f15!(T, R),
537 			&f16!(T, R),
538 			&f17!(T, R),
539 			&f18!(T, R),
540 			&f19!(T, R),
541 			&f20!(T, R),
542 		];
543
544 		foreach(relTol; [
545 			T.min_normal,
546 			T.epsilon,
547 			sqrt(T.epsilon),
548 			T.epsilon^^0.25,
549 			]) // 3
550 		foreach(absTol; [
551 			2*T.epsilon,
552 			sqrt(T.epsilon),
553 			]) // 4
554 		{
555 			foreach(sign; [-1, 1]) // 5
556 			foreach(rBound; [1, T.max]) // 6
557 			foreach(shift; [0, -1, 1]) // 7
558 			foreach(j, f; ar1) // 8
559 			{
560 				auto m2 = findLocalMin!T((T x) => sign * f(x-shift), shift, rBound+shift, relTol, absTol);
561 			}
562 		}
563 	}
564 }
565
566
567 Unqual!(CommonType!(T1, T2))
568 chebevReversed(T1, T2)(in T1[] c, in T2 x)
569 {
570 	typeof(return) d = 0;
571 	typeof(return) dd = 0;
572 	foreach(e; c[0..\$-1])
573 	{
574 		immutable sv = d;
575 		d = 2 * x * d - dd + e;
576 		dd = sv;
577 	}
578     return x * d - dd + c[\$-1] / 2;
579 }
580
581
582
583 T modifiedBesselCF2(T)(in T nu, in T x) @nogc @safe pure nothrow
584 	if(isFloatingPoint!T)
585 in {
586 	assert(fabs(nu) <= 0.5f);
587 	assert(x >= 2);
588 	assert(isFinite(x));
589 }
590 body {
591 	T bc = 2 * (1 + x);
592 	if(isInfinity(bc))
593 		return 1;
594 	T d = 1 / bc;
595 	T h = d;
596 	T delh = d;
597 	T ac = nu^^2 - 0.25f;
598 	int i = 1;
599 	do
600 	{
601 		immutable ip1 = i + 1;
602 		immutable a = ac - ip1 * i;
603 		immutable b = bc + i * 2;
604 		i = ip1;
605 		d = 1 / (b + a * d);
606 		delh = (b * d - 1) * delh;
607 		h += delh;
608 		debug assert(!isNaN(delh / h));
609 	}
610 	while(fabs(delh / h) > T.epsilon);
611 	return (nu + 0.5f + x + ac * h) / x;
612 }
613
614
615
616 T[2] modifiedBesselCF2Full(T)(in T nu, in T x) @nogc @safe pure nothrow
617 	if(isFloatingPoint!T)
618 in {
619 	assert(fabs(nu) <= 0.5f);
620 	assert(x >= 2);
621 	assert(isFinite(x));
622 }
623 body {
624 	T bc = 2 * (1 + x);
625 	if(isInfinity(bc))
626 		return [nu < 0.5f ? 0 : 1.2533141373155003, 1];
627 	T d = 1 / bc;
628 	T h = d;
629 	T delh = d;
630 	T ac = nu^^2 - 0.25f;
631 	T q1 = 0;
632 	T q2 = 1;
633 	T q = -ac;
634 	T c = -ac;
635 	T s = 1 + q*delh;
636 	T dels = void;
637 	uint i = 1;
638 	T b = bc;
639 	do
640 	{
641 		immutable ip1 = i + 1;
642 		immutable a = ac - ip1 * i;
643 		c = -a * c / ip1;
644 		immutable qnew = (q1 - b * q2) / a;
645 		q1 = q2;
646 		q2 = qnew;
647 		q += c * qnew;
648 		b = bc + i * 2;
649 		d = 1 / (b + a * d);
650 		delh = (b * d - 1) * delh;
651 		dels = q * delh;
652 		h += delh;
653 		s += dels;
654 		i = ip1;
655 	}
656 	//Need only test convergence of sum since CF2 itself converges more quickly.
657 	while(fabs(dels / s) > T.epsilon);
658
659 	enum T SQRTPI_2 = sqrt(PI / 2);
660 	return [SQRTPI_2 / s, (nu + 0.5 + x + ac * h) / x];
661 }
662
663
664 unittest
665 {
666 	assert(approxEqual(modifiedBesselCF2(0.4, 2.00001), 1.44212, 0.0, 1e-4));
667 	assert(approxEqual(modifiedBesselCF2(0.5, 2.00001), 1.5, 0.0, 1e-4));
668 	assert(approxEqual(modifiedBesselCF2(0.0, 2.00001), 1.22804, 0.0, 1e-4));
669 	assert(approxEqual(modifiedBesselCF2(0.4, 1000.0), 1.00100, 0.0, 1e-3));
670 	assert(approxEqual(modifiedBesselCF2(0.499, 1000.0), 1.00090, 0.0, 1e-4));
671 	assert(approxEqual(modifiedBesselCF2(0.0, 1000), 1.000499875124805092705355767307116656100795175108335495396004, 0.0, 1e-14));
672 }
673
674 /++
675 Returns:
676 	[K(x, nu), 2 / x * K(x, nu+1)]
677 +/
678 T[2] modifiedBesselTemmeSeriesImpl(T)(in T nu, in T x) @nogc @safe pure nothrow
679 in {
680 	assert(fabs(nu) <= 0.5f);
681 	assert(x <= 2);
682 }
683 body {
684 	//TODO: recalculate c1 and c2 for different FP formats
685 	static if(is(T == float))
686 	{
687 		static immutable float[4] c1 = [
688 			-3.4706269649e-6,
689 			3.087090173086e-4,
690 			6.5165112670737e-3,
691 			-1.142022680371168e0,
692 			];
693 		static immutable float[8] c2 = [
694 			-4.9717367042e-6,
695 			1.2719271366546e-3,
696 			-7.68528408447867e-2,
697 			1.843740587300905e0,
698 			];
699 	}
700 	else
701 	{
702 		static immutable double[7] c1 = [
703 			-1.356e-13,
704 			3.67795e-11,
705 			6.9437664e-9,
706 			-3.4706269649e-6,
707 			3.087090173086e-4,
708 			6.5165112670737e-3,
709 			-1.142022680371168e0,
710 			];
711 		static immutable double[8] c2 = [
712 			-1.49e-15,
713 			-1.702e-13,
714 			2.423096e-10,
715 			-3.31261198e-8,
716 			-4.9717367042e-6,
717 			1.2719271366546e-3,
718 			-7.68528408447867e-2,
719 			1.843740587300905e0,
720 			];
721 	}
722
723 	immutable T x2 = x / 2;
724 	immutable T pimu = cast(T) PI * nu;
725 	immutable T nu2 = nu^^2;
726 	immutable T fact = abs(pimu) < T.epsilon ? 1 : pimu / sin(pimu);
727 	T d = -log(x2);
728 	T e = nu * d;
729 	immutable T fact2 = abs(e) < T.epsilon ? 1 : sinh(e) / e;
730 	immutable T xx = 8 * nu^^2 - 1;
731 	immutable T gam1 = chebevReversed(c1, xx);
732 	immutable T gam2 = chebevReversed(c2, xx);
733 	immutable T gampl = gam2 - nu * gam1;
734 	immutable T gammi = gam2 + nu * gam1;
735 	T ff = fact * (gam1 * cosh(e) + gam2 * fact2 * d);
736 	T sum = ff;
737 	e = exp(e);
738 	T p = 0.5f * e / gampl;
739 	T q = 0.5f / (e * gammi);
740 	T c = 1;
741 	d = x2^^2;
742 	T sum1 = p;
743 	int i;
744 	T del = void;
745 	do {
746 		i++;
747 		ff = (i * ff + p + q) / (i * i - nu2);
748 		c *= d / i;
749 		p /= i - nu;
750 		q /= i + nu;
751 		del = c * ff;
752 		sum += del;
753 		immutable del1 = c * (p - i * ff);
754 		sum1 += del1;
755 	}
756 	while(abs(del) >= abs(sum) * T.epsilon);
757 	return [sum, sum1];
758 }
759
760 /++
761 Returns:
762 	K(x, nu+1) / K(x, nu)
763 +/
764 T modifiedBesselTemmeSeries(T)(in T nu, in T x)
765 in {
766 	assert(fabs(nu) <= 0.5f);
767 	assert(x <= 2);
768 }
769 body {
770 	immutable sums = modifiedBesselTemmeSeriesImpl(nu, x);
771 	return sums[1] / sums[0] / x * 2;
772 }
773
774 unittest
775 {
776 	assert(approxEqual(modifiedBesselTemmeSeries(0.4, 1.0), 1.87491, 0.0, 1e-4));
777 	assert(approxEqual(modifiedBesselTemmeSeries(0.499, 1.0), 1.9987, 0.0, 1e-4));
778 	assert(approxEqual(modifiedBesselTemmeSeries(0.0, 1.0), 1.429625398260401758028108023450365630080017018192662983459034, 0.0, 1e-14));
779 }
780
781
782 T besselKD(T)(in T nu, in T x)
783 	if(isFloatingPoint!T)
784 out(r) {
785 	import std.conv;
786 	assert(r.isNaN || r >= 1, text("x = ", x, " nu = ", nu, " r = ", r));
787 }
788 in {
789 	assert(x >= 0);
790 	assert(x < T.infinity);
791 }
792 body {
793 	version(LDC)
794 	{
795 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
796 	}
797 	immutable anu = fabs(nu);
798 	int nl = cast(int)floor(anu+0.5);
799 	T mu=anu-nl;
800 	//import std.stdio;
801 	////writeln(nu);
802 	assert(fabs(mu) <= 0.5f);
803 	T ret;
804 	if(x >= 2)
805 	{
806 		T r = modifiedBesselCF2(mu, x);
807 		mu *= 2;
808 		if(nl)
809 		{
810 			T d;
811 			do
812 			{
813 				d = r;
814 				mu += 2;
815 				r = mu / x + 1 / r;
816 			}
817 			while(--nl);
818 			ret = r / d;
819 		}
820 		else
821 		{
822 			ret = r * (r - mu / x);
823 		}
824 	}
825 	else
826 	{
827 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
828 		T r = sums[1] / sums[0];
829 		if(nl == 0)
830 		{
831 			if(mu < 0)
832 			{
833 				if(r >= mu)
834 					r = mu.nextDown;
835 			}
836 			else if(mu > 0)
837 			{
838 				if(r <= mu)
839 					r = mu.nextUp;
840 			}
841 			else
842 			{
843 				if(r == 0)
844 					return x ? 1 / (x * log(x))^^2 : T.infinity;
845 			}
846 		}
847 		//import std.stdio;
848 		//writeln("r = ", r);
849 		//writeln("sums[1] = ", sums[1]);
850 		//writeln("sums[0] = ", sums[0]);
851 		if(nl)
852 		{
853 			immutable x2_4 = x * x / 4;
854 			T d;
855 			do
856 			{
857 				d = r;
858 				mu += 1;
859 				r = mu + x2_4 / r;
860 			}
861 			while(--nl);
862 			ret = r / d;
863 		}
864 		else
865 		{
866 			//
867 			//writeln("r - mu = ", r - mu);
868 			ret = r / x * (r - mu) / x * 4;
869 		}
870 	}
871 	//import std.stdio;
872 	//writeln("ret = ", ret);
873 	//writeln("x = ", x);
874 	return fmax(ret, cast(T)1);
875 }
876
877 unittest {
878 	assert(approxEqual(besselKD(0.46, 1.0),  2.006492004938732508413453526510952827555303123690945926246942, 0.0, 1e-14));
879 	assert(approxEqual(besselKD(0.46, 1e10),  1.000000000100000000000000000000019199999996160000000869529599, 0.0, 1e-14));
880 	assert(approxEqual(besselKD(0.46, 1e5),  1.000010000000000019199616008695062658208936100866055206637995, 0.0, 1e-14));
881 	//import std.stdio;
882 	//writeln(besselKD(0.23, 5.20964e-116 ));
883 	assert(approxEqual(besselKD(0.46, 1e-10), 5.2419685330917511351588283276691731294601340736654433005e+10, 1, 1e-14));
884 	assert(approxEqual(besselKD(0.0, 1.0),  2.043828779351212337989476332008573915738881135574432645409779, 0.0, 1e-14));
885 	assert(approxEqual(besselKD(0.0, 2.0), 1.508074700999049512999886217349628893603326842760291910336296, 0.0, 1e-14));
886 	assert(approxEqual(besselKD(0.0, 0.5), 3.210807086475177172355017867637452321623564307416114645664507, 0.0, 1e-14));
887 	assert(approxEqual(besselKD(1.0, 0.5), 2.543749846614157209231745521424755087195088446583780182048739, 0.0, 1e-14));
888 	assert(approxEqual(besselKD(1.0, 1.0), 1.888245647341297344338366806469818100073176922119421727713106, 0.0, 1e-14));
889 	assert(approxEqual(besselKD(1.0, 2.0), 1.477404884746695468820498957141015915645967793158828983150281, 0.0, 1e-14));
890 	assert(approxEqual(besselKD(0.4, 2.0), 1.502873827552600466534737445302356963995715395265386806359790, 0.0, 1e-14));
891 	assert(approxEqual(besselKD(0.4, 1.0), 2.015349430323277001160289680464116945013807398639704169522157, 0.0, 1e-14));
892 	assert(approxEqual(besselKD(0.4, 0.5), 3.071599175777718550825656123800662591905836277372864541044359, 0.0, 1e-13));
893 	assert(approxEqual(besselKD(11.1, 1.11e+10), 1.000000000090090090090090090045136443975807419926316735329998, 0.0, 1e-14));
894 	assert(approxEqual(besselKD(11.1, 1.11e-10), 1.099009900990099009900983462621096186432918126639008518272442, 0.0, 1e-14));
895 	assert(approxEqual(besselKD(2.0, 3.0), 1.296650971433360224728795947829976056477719397636727742896, 0.0, 1e-14));
896 	assert(approxEqual(besselKD(2.2, 0.8), 1.624595391806645609233071414803368115353081130470326304255669, 0.0, 1e-14));
897 	assert(approxEqual(besselKD(0.4, 3.0), 1.334349384832101337822320699631685431560403480505548447388101, 0.0, 1e-14));
898 	assert(approxEqual(besselKD(0.4, 0.8), 2.275575715445329347201695021084846610619371526255383391028726, 0.0, 1e-14));
899 	assert(approxEqual(besselKD(1.0, 1.0), 1.888245647341297344338366806469818100073176922119421727713106268163595832389436074069065151380128808, 0.0, 1e-14));
900 	assert(approxEqual(besselKD(1.0, 10.0), 1.099687904173359294358694759872440909233436542663077290833713959244496776538436835904856529604568407, 0.0, 1e-14));
901 	assert(approxEqual(besselKD(10.0, 1.0), 1.110348556190446111957216138204031855721452559868300859853180425692157659153049491662221677392988991, 0.0, 1e-14));
902 	assert(approxEqual(besselKD(10.0, 1.49011e-08), 1.11111, 0.0, 1e-4));
903 	assert(approxEqual(besselKD(1.0, 1.49011e-08), 36.2755, 0.0, 1e-4));
904 	assert(approxEqual(besselKD(1.0, 2.22044e-16), 72.3192, 0.0, 1e-4));
905 	assert(approxEqual(besselKD(1.0, 1.49166e-154), 708.628, 0.0, 1e-3));
906 	assert(approxEqual(besselKD(10.0, 1.49166e-154), 1.111111111111111, 0.0, 1e-4));
907 	assert(approxEqual(besselKD(1.1, double.min_normal), 11.0, 0.0, 1e-4));
908 	assert(isFinite(besselKD(0.1, double.epsilon)));
909 	assert(isFinite(besselKD(0, double.epsilon)));
910 }
911
912
913 T besselKRM(T)(in T nu, in T x)
914 	if(isFloatingPoint!T)
915 out(r){
916 	//import std.stdio;
917 	//writefln("KRM(%s, %s) = %s", nu, x, r);
918 }
919 in {
920
921 }
922 body {
923 	version(LDC)
924 	{
925 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
926 	}
927 	if(nu < 0)
928 		return 1 / besselKRM(-nu, x);
929 	immutable anu = fabs(nu);
930 	int nl = cast(int)floor(anu+0.5f);
931 	T mu=anu-nl;
932 	assert(fabs(mu) <= 0.5f);
933
934 	if(x >= 2)
935 	{
936 		T r = modifiedBesselCF2(mu, x);
937 		mu *= 2;
938 		if(nl)
939 		{
940 			T d;
941 			do
942 			{
943 				d = r;
944 				mu += 2;
945 				r = mu / x + 1 / r;
946 			}
947 			while(--nl);
948 			return sqrt(r) * sqrt(d);
949 		}
950 		else
951 		{
952 			return sqrt(r / (r - mu / x));
953 		}
954 	}
955 	else
956 	{
957 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
958 		T r = sums[1] / sums[0];
959 		//import std.stdio;
960 		//writeln(sums, " ", r);
961 		if(nl)
962 		{
963 			immutable x2_4 = x * x / 4;
964 			T d;
965 			do
966 			{
967 				d = r;
968 				mu += 1;
969 				r = mu + x2_4 / r;
970 			}
971 			while(--nl);
972 			//import std.stdio;
973 			//writeln("sqrt(r) = ", sqrt(r), " sqrt(d) = ", sqrt(d));
974 			return sqrt(r) * sqrt(d) * 2 / x;
975 		}
976 		else
977 		{
978 			//writeln("r = ", r, " mu = ", mu);
979 			T ret = sums[1] / (sums[1] - mu * sums[0]);
980 			if(ret < 0)
981 			{
982 				return T.infinity;
983 			}
984 			//import std.stdio;
985 			//writeln("sums[1] - mu * sums[0] = %s, ", sums[1], sums[1] - mu * sums[0]);
986 			if(ret == T.infinity)
987 			{
988 				//writeln(ret);
989 				ret = T.max;
990 				//return T.max;
991 			}
992 			return sqrt(ret);
993 		}
994 	}
995 }
996
997 unittest {
998 	assert(approxEqual(besselKRM(0.0, 1.0),  1, 0.0, 1e-14));
999 	assert(approxEqual(besselKRM(0.0, 2.0), 1, 0.0, 1e-14));
1000 	assert(approxEqual(besselKRM(0.0, 0.5), 1, 0.0, 1e-14));
1001 	assert(approxEqual(besselKRM(1.0, 0.5), 2.857882088842869132027932428611809106473984497664375438969650, 0.0, 1e-14));
1002 	assert(approxEqual(besselKRM(1.0, 1.0), 1.964497593920848638120029967522901757996791347548236796086390, 0.0, 1e-14));
1003 	assert(approxEqual(besselKRM(1.0, 2.0), 1.492661023078886446121280337304042151124577136818374559084502, 0.0, 1e-14));
1004 	assert(approxEqual(besselKRM(0.4, 2.0), 1.176363556186439775048077295541564293759652878052851708370510, 0.0, 1e-14));
1005 	assert(approxEqual(besselKRM(0.4, 1.0), 1.320700844681798447461890099116285783758321865242462619205560, 0.0, 1e-14));
1006 	assert(approxEqual(besselKRM(0.4, 0.5), 1.555719773577105697258998484956087203635813355802430411799783, 0.0, 1e-14));
1007 	assert(approxEqual(besselKRM(11.1, 1.11e+10), 1.000000001000000000454954954912953493935881553614516153308770, 0.0, 1e-14));
1008 	assert(approxEqual(besselKRM(11.1, 1.11e-10), 1.9077839604210010339594792869174419072661582323995475082e+11, 0.0, 1e-14));
1009 }
1010
1011
1012 T logBesselK(T)(in T nu, in T x)
1013 	if(isFloatingPoint!T)
1014 in {
1015 	assert(x.isFinite);
1016 	assert(x >= T.min_normal);
1017 }
1018 body {
1019 	version(LDC)
1020 	{
1021 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
1022 	}
1023 	immutable anu = fabs(nu);
1024 	int nl = cast(int)floor(anu+0.5f);
1025 	T mu = anu - nl;
1026 	assert(fabs(mu) <= 0.5f);
1027 	if(x >= 2)
1028 	{
1029 		immutable r2 = modifiedBesselCF2Full(mu, x);
1030 		T r = r2[1];
1031 		T ret = log(r2[0]) - log(x) / 2 - x;
1032 		if(nl)
1033 		{
1034 			mu *= 2;
1035 			T l = 1;
1036 			long exs;
1037 			do
1038 			{
1039 				int ex = void;
1040 				l *= frexp(r, ex);
1041 				exs += ex;
1042 				mu += 2;
1043 				r = mu / x + 1 / r;
1044 			}
1045 			while(--nl);
1046 			ret += cast(T)LN2 * (exs + log2(l));
1047 		}
1048 		return ret;
1049 	}
1050 	else
1051 	{
1052 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
1053 		T r = sums[1] / sums[0];
1054 		T ret = log(sums[0]);
1055 		if(nl)
1056 		{
1057 			ret += nl * (cast(T)LN2 - log(x));
1058 			immutable x2_4 = x * x / 4;
1059 			T l = 1;
1060 			long exs;
1061 			do
1062 			{
1063 				int ex = void;
1064 				l *= frexp(r, ex);
1065 				exs += ex;
1066 				mu += 1;
1067 				r = mu + x2_4 / r;
1068 			}
1069 			while(--nl);
1070 			ret += cast(T)LN2 * (exs + log2(l));
1071 		}
1072 		return ret;
1073 	}
1074 }
1075
1076 unittest {
1077 	assert(approxEqual(logBesselK(0.0, 1.0),  -0.86506439890678809679875790803368568022489315035161867839839, 0.0, 1e-14));
1078 	assert(approxEqual(logBesselK(0.0, 2.0), -2.17248820497570993473841333643717923143973636448830725037508, 0.0, 1e-14));
1079 	assert(approxEqual(logBesselK(0.0, 0.5), -0.07858976986908141689523697453802973224104044591707028228710, 0.0, 1e-14));
1080 	assert(approxEqual(logBesselK(1.0, 0.5), 0.504671397304651177308416839874408827783443947120148152746119, 0.0, 1e-14));
1081 	assert(approxEqual(logBesselK(1.0, 1.0), -0.50765194821075233094791485120634298590875979858568077818450, 0.0, 1e-14));
1082 	assert(approxEqual(logBesselK(1.0, 2.0), -1.96707130256051389147686464265533209478058062247688850653003, 0.0, 1e-14));
1083 	assert(approxEqual(logBesselK(0.4, 2.0), -2.13936877477972262972321591624176676086808535490162854194273, 0.0, 1e-14));
1084 	assert(approxEqual(logBesselK(0.4, 1.0), -0.80679541168661951923202239125477303379665660825595791745329, 0.0, 1e-14));
1085 	assert(approxEqual(logBesselK(0.4, 0.5), 0.018456437585950072585426399714417863872587115447191398524780, 0.0, 1e-14));
1086 	assert(approxEqual(logBesselK(11.1, 1.11e+10), -1.110000001133931411444888363310129265544563412394720435e10, 0.0, 1e-14));
1087 	assert(approxEqual(logBesselK(11.1, 1.11e-10), 276.7693978383758755547170249282452519578463074349871817713218, 0.0, 1e-11));
1088 }
1089
1090
1091 T besselK(Flag!"ExponentiallyScaled" expFlag = Flag!"ExponentiallyScaled".no, T)(in T nu, in T x)
1092 	if(isFloatingPoint!T)
1093 in {
1094 	assert(x.isFinite);
1095 	assert(x >= T.min_normal);
1096 }
1097 body {
1098 	version(LDC)
1099 	{
1100 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
1101 	}
1102 	immutable anu = fabs(nu);
1103 	int nl = cast(int)floor(anu+0.5f);
1104 	T mu = anu - nl;
1105 	assert(fabs(mu) <= 0.5f);
1106 	T r = void, ret = void;
1107 	if(x >= 2)
1108 	{
1109 		immutable r2 = modifiedBesselCF2Full(mu, x);
1110 		r = r2[1];
1111 		static if(expFlag)
1112 			ret = r2[0] / sqrt(x);
1113 		else
1114 			ret = r2[0] / sqrt(x) * exp(-x);
1115 	}
1116 	else
1117 	{
1118 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
1119 		r = sums[1] / sums[0] / x * 2;
1120 		static if(expFlag)
1121 			ret = sums[0] * exp(x);
1122 		else
1123 			ret = sums[0];
1124 	}
1125 	if(nl)
1126 	{
1127 		mu *= 2;
1128 		do
1129 		{
1130 			ret *= r;
1131 			mu += 2;
1132 			r = mu / x + 1 / r;
1133 		}
1134 		while(--nl);
1135 	}
1136 	return ret;
1137
1138 }
1139
1140 unittest {
1141 	assert(approxEqual(besselK(0.0, 1.0),  0.421024438240708333335627379212609036136219748226660472298969, 0.0, 1e-14));
1142 	assert(approxEqual(besselK(0.0, 2.0), 0.113893872749533435652719574932481832998326624388808882892529, 0.0, 1e-14));
1143 	assert(approxEqual(besselK(0.0, 0.5), 0.924419071227665861781924167530216989538768311953529684815019, 0.0, 1e-14));
1144 	assert(approxEqual(besselK(1.0, 0.5), 1.656441120003300893696445403174091511534100759464077446055427, 0.0, 1e-14));
1145 	assert(approxEqual(besselK(1.0, 1.0), 0.601907230197234574737540001535617339261586889968106456017767, 0.0, 1e-14));
1146 	assert(approxEqual(besselK(1.0, 2.0), 0.139865881816522427284598807035411023887234584841515530384442, 0.0, 1e-14));
1147 	assert(approxEqual(besselK(0.4, 2.0), 0.117729133170423325690699234439848483526261145336868148793547, 0.0, 1e-14));
1148 	assert(approxEqual(besselK(0.4, 1.0), 0.446285939834668179310023409179745077844437039396299727961180, 0.0, 1e-14));
1149 	assert(approxEqual(besselK(0.4, 0.5), 1.018627810316608462220290727117637366662549451305627759999544, 0.0, 1e-14));
1150 	assert(approxEqual(besselK(11.1, 100000.0), 0.0, 0.0, 1e-14));
1151 	assert(approxEqual(besselK(11.1, 1/100000.0), 1.5940696949048155233993419471665743544335615887236248959e+65, 1e51, 1e-14));
1152
1153 	assert(approxEqual(besselK!(Flag!"ExponentiallyScaled".yes)(0.4, 2.0), 0.869907169474735192472996963053995791191543506954884340734219, 0.0, 1e-14));
1154 	assert(approxEqual(besselK!(Flag!"ExponentiallyScaled".yes)(0.4, 1.0), 1.213130960549345270517674559709580335203213288805317162444177, 0.0, 1e-14));
1155 	assert(approxEqual(besselK!(Flag!"ExponentiallyScaled".yes)(0.4, 0.5), 1.679433337795687807090050796759423855788542163267174953125487, 0.0, 1e-14));
1156 	assert(approxEqual(besselK!(Flag!"ExponentiallyScaled".yes)(11.1, 100000.0), 0.003965764688216290045701642597190265433950245246337093558273, 0.0, 1e-14));
1157
1158 	assert(approxEqual(besselK(1.0, 2.0), 0.139865881816522427284, 0.0, 1e-14));
1159 	assert(approxEqual(besselK(1.0, 20.0), 5.88305796955703817765028e-10, 0.0, 1e-14));
1160 	assert(approxEqual(besselK(1.2, 1.0), 0.701080));
1161 	assert(approxEqual(besselK(10.0, 400.0), 1.359308280937876101049800424530298700140591995078599604e-175, 0.0, 1e-14));
1162 }
1163
1164
1165 T besselKRS(T)(in T nu, in T x)
1166 	if(isFloatingPoint!T)
1167 in {
1168 	assert(x.isFinite);
1169 	assert(x >= T.min_normal);
1170 }
1171 body {
1172 	version(LDC)
1173 	{
1174 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
1175 	}
1176 	immutable anu = fabs(nu);
1177 	int nl = cast(int)floor(anu+0.5f);
1178 	T mu=anu-nl;
1179 	assert(fabs(mu) <= 0.5f);
1180 	if(x >= 2)
1181 	{
1182 		T r = modifiedBesselCF2(mu, x);
1183 		mu *= 2;
1184 		if(nl)
1185 		{
1186 			T d = void;
1187 			do
1188 			{
1189 				d = r;
1190 				mu += 2;
1191 				r = mu / x + 1 / r;
1192 			}
1193 			while(--nl);
1194 			return r + 1 / d;
1195 		}
1196 		else
1197 		{
1198 			return r + (r - mu / x);
1199 		}
1200 	}
1201 	else
1202 	{
1203 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
1204 		T r = sums[1] / sums[0];
1205 		immutable x_2 = x / 2;
1206 		if(nl)
1207 		{
1208 			immutable x2_4 = x_2^^2;
1209 			T d = void;
1210 			do
1211 			{
1212 				d = r;
1213 				mu += 1;
1214 				r = mu + x2_4 / r;
1215 			}
1216 			while(--nl);
1217 			return r / x_2 + x_2 / d;
1218 		}
1219 		else
1220 		{
1221 			if(fabs(r) < fabs(mu))
1222 				r = mu;
1223 			return r / x_2 + (r - mu) / x_2;
1224 		}
1225 	}
1226 }
1227
1228 unittest {
1229 	assert(approxEqual(besselKRS(0.0, 1.0), 2.859250796520803516056216046900731260160034036385325966918068, 0.0, 1e-14));
1230 	assert(approxEqual(besselKRS(0.0, 2.0), 2.456073859637815951485344904163437808478733905321719388934076, 0.0, 1e-14));
1231 	assert(approxEqual(besselKRS(0.0, 0.5), 3.583745016864440467305949390266788044526362750672196202869812, 0.0, 1e-14));
1232 	assert(approxEqual(besselKRS(1.0, 0.5), 5.116150836953170679741316342708831083392598534819729276449587, 0.0, 1e-14));
1233 	assert(approxEqual(besselKRS(1.0, 1.0), 3.398967871187544687785354799542729509747041034806318935543498, 0.0, 1e-14));
1234 	assert(approxEqual(besselKRS(1.0, 2.0), 2.628615517527578979897586948721927750995112315346619061302174, 0.0, 1e-14));
1235 	assert(approxEqual(besselKRS(0.4, 2.0), 2.484249446052147531028619260526756105681146308342088119522102, 0.0, 1e-14));
1236 	assert(approxEqual(besselKRS(0.4, 1.0), 2.949813167184170666766713209833464015196615131065026837642303, 0.0, 1e-14));
1237 	assert(approxEqual(besselKRS(0.4, 0.5), 3.853102218097889263846807474862866521371874509444807669438352, 0.0, 1e-14));
1238 	assert(approxEqual(besselKRS(11.1, 1.11e+10), 2.000000000090090091088061033917072660444297152286161364522101, 0.0, 1e-14));
1239 	assert(approxEqual(besselKRS(11.1, 1.11e-10), 2.0000000000000000000001099009900990099009900953267052034e+11, 0.0, 1e-14));
1240 }
1241
1242
1243 T besselKRX(T)(in T nu, in T x)
1244 	if(isFloatingPoint!T)
1245 in {
1246 	assert(nu >= 0);
1247 	assert(x.isFinite);
1248 	assert(x >= T.min_normal);
1249 }
1250 body {
1251 	version(LDC)
1252 	{
1253 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
1254 	}
1255 	immutable anu = fabs(nu);
1256 	int nl = cast(int)floor(anu+0.5f);
1257 	T mu = anu - nl;
1258 	assert(fabs(mu) <= 0.5f);
1259 	if(x >= 2)
1260 	{
1261 		T r = modifiedBesselCF2(mu, x);
1262 		if(nl)
1263 		{
1264 			mu *= 2;
1265 			do
1266 			{
1267 				mu += 2;
1268 				r = mu / x + 1 / r;
1269 			}
1270 			while(--nl);
1271 		}
1272 		return r * x;
1273 	}
1274 	else
1275 	{
1276 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
1277 		T r = sums[1] / sums[0];
1278 		if(nl)
1279 		{
1280 			immutable x2_4 = x * x / 4;
1281 			do
1282 			{
1283 				mu += 1;
1284 				r = mu + x2_4 / r;
1285 			}
1286 			while(--nl);
1287 		}
1288 		else
1289 		{
1290 			if(fabs(r) < fabs(mu))
1291 				r = mu;
1292 		}
1293 		return r * 2;
1294 	}
1295 }
1296
1297 unittest {
1298 	assert(approxEqual(besselKRX(10.0, double.min_normal), 20.0, 0.0, 1e-14));
1299 	assert(approxEqual(besselKRX(0.0, 1.0), 1.429625398260401758028108023450365630080017018192662983459034, 0.0, 1e-14));
1300 	assert(approxEqual(besselKRX(0.0, 2.0), 2.456073859637815951485344904163437808478733905321719388934076, 0.0, 1e-14));
1301 	assert(approxEqual(besselKRX(0.0, 0.5), 0.895936254216110116826487347566697011131590687668049050717453, 0.0, 1e-14));
1302 	assert(approxEqual(besselKRX(1.0, 0.5), 2.279037709238292669935329085677207770848149633704932319112396, 0.0, 1e-14));
1303 	assert(approxEqual(besselKRX(1.0, 1.0), 2.699483935593772343892677399771364754873520517403159467771749, 0.0, 1e-14));
1304 	assert(approxEqual(besselKRX(1.0, 2.0), 3.628615517527578979897586948721927750995112315346619061302174, 0.0, 1e-14));
1305 	assert(approxEqual(besselKRX(0.4, 2.0), 2.884249446052147531028619260526756105681146308342088119522102, 0.0, 1e-14));
1306 	assert(approxEqual(besselKRX(0.4, 1.0), 1.874906583592085333383356604916732007598307565532513418821151, 0.0, 1e-14));
1307 	assert(approxEqual(besselKRX(0.4, 0.5), 1.363275554524472315961701868715716630342968627361201917359588, 0.0, 1e-14));
1308 	assert(approxEqual(besselKRX(11.1, 11.1), 27.05342065662651512004991741552637606461741839973809374665713, 0.0, 1e-14));
1309 	assert(approxEqual(besselKRX(11, 2.0^^30), 1.000000010710209660444967811966826255158285934025032904639974*2.0^^30, 0.0, 1e-14));
1310 	assert(approxEqual(besselKRX(11, 2.0^^30), 1.000000010710209660444967811966826255158285934025032904639974*2.0^^30, 0.0, 1e-14));
1311 	assert(approxEqual(besselKRX(11, 2.0^^40), 1.000000000010459189070438615765505792843518644341506580233109*2.0^^40, 0.0, 1e-14));
1312 	assert(approxEqual(besselKRX(11, 1e+10), 1.0000000011500000006037499999396249998267992188210915625e+10, 1, 1e-14));
1313 	assert(approxEqual(besselKRX(11.1, 1.11e+10), 1.1100000011600000005538738738239753265465849195188195573e+10, 1, 1e-14));
1314 	assert(approxEqual(besselKRX(11.1, 1.11e-10), 22.20000000000000000000060995049504950495049502906321387905301, 0.0, 1e-14));
1315 }
1316
1317
1318 T besselKR(T)(in T nu, in T x)
1319 	if(isFloatingPoint!T)
1320 in {
1321 }
1322 body {
1323 	version(LDC)
1324 	{
1325 		import ldc.intrinsics: log = llvm_log, log2 = llvm_log2, fabs = llvm_fabs;
1326 	}
1327 	if(nu < 0)
1328 		return 1 / besselKR(-nu-1, x);
1329 	if(x.isNaN() || nu.isNaN() || x < 0)
1330 		return T.init;
1331 	if(x < T.min_normal)
1332 		return T.infinity;
1333 	immutable anu = fabs(nu);
1334 	int nl = cast(int)floor(anu+0.5f);
1335 	T mu = anu - nl;
1336 	assert(fabs(mu) <= 0.5f);
1337 	if(x >= 2)
1338 	{
1339 		T r = modifiedBesselCF2(mu, x);
1340 		if(nl)
1341 		{
1342 			mu *= 2;
1343 			do
1344 			{
1345 				mu += 2;
1346 				r = mu / x + 1 / r;
1347 			}
1348 			while(--nl);
1349 		}
1350 		return r;
1351 	}
1352 	else
1353 	{
1354 		immutable sums = modifiedBesselTemmeSeriesImpl(mu, x);
1355 		T r = sums[1] / sums[0];
1356 		if(nl)
1357 		{
1358 			immutable x2_4 = x * x / 4;
1359 			do
1360 			{
1361 				mu += 1;
1362 				r = mu + x2_4 / r;
1363 			}
1364 			while(--nl);
1365 		}
1366 		return r * 2 / x;
1367 	}
1368 }
1369
1370 unittest {
1371 	assert(approxEqual(besselKR(0.0, 1.0)    , 1.429625398260401758028108023450365630080017018192662983459034, 0.0, 1e-14));
1372 	assert(approxEqual(besselKR(0.0, 2.0) * 2, 2.456073859637815951485344904163437808478733905321719388934076, 0.0, 1e-14));
1373 	assert(approxEqual(besselKR(0.0, 0.5) / 2, 0.895936254216110116826487347566697011131590687668049050717453, 0.0, 1e-14));
1374 	assert(approxEqual(besselKR(1.0, 0.5) / 2, 2.279037709238292669935329085677207770848149633704932319112396, 0.0, 1e-14));
1375 	assert(approxEqual(besselKR(1.0, 1.0)    , 2.699483935593772343892677399771364754873520517403159467771749, 0.0, 1e-14));
1376 	assert(approxEqual(besselKR(1.0, 2.0) * 2, 3.628615517527578979897586948721927750995112315346619061302174, 0.0, 1e-14));
1377 	assert(approxEqual(besselKR(0.4, 2.0) * 2, 2.884249446052147531028619260526756105681146308342088119522102, 0.0, 1e-14));
1378 	assert(approxEqual(besselKR(0.4, 1.0)    , 1.874906583592085333383356604916732007598307565532513418821151, 0.0, 1e-14));
1379 	assert(approxEqual(besselKR(0.4, 0.5) / 2, 1.363275554524472315961701868715716630342968627361201917359588, 0.0, 1e-14));
1380 	assert(approxEqual(besselKR(11.1, 1.11e+10), 1.1100000011600000005538738738239753265465849195188195573e+10 / 1.11e+10, 1, 1e-14));
1381 	assert(approxEqual(besselKR(11.1, 1.11e-10) * 1.11e-10, 22.20000000000000000000060995049504950495049502906321387905301, 0.0, 1e-14));
1382
1383 	assert(approxEqual(besselKR(-1.0, 1.0), 1 / 1.429625398260401758028108023450365630080017018192662983459034, 0.0, 1e-14));
1384 }
1385
1386
1387 T besselKR(T)(in T nu, in T x, uint k)
1388 {
1389 	T l = besselKR(nu, x);
1390 	T r = l;
1391 	T mu = nu * 2;
1392 	while(k--)
1393 	{
1394 		mu += 2;
1395 		r = mu / x + 1 / r;
1396 		l *= r;
1397 	}
1398 	return l;
1399 }
1400
1401 unittest {
1402 	assert(approxEqual(besselKR(0.0, 1.0, 0)    , 1.429625398260401758028108023450365630080017018192662983459034, 0.0, 1e-14));
1403 	assert(approxEqual(besselKR(0.0, 2.0, 0) * 2, 2.456073859637815951485344904163437808478733905321719388934076, 0.0, 1e-14));
1404 	assert(approxEqual(besselKR(0.0, 1.0, 3), 105.0590223025824984495740493132204752844809530187891270737059, 1e-10, 1e-14));
1405 	assert(approxEqual(besselKR(0.0, 2.0, 3), 19.28036929818907975742672452081718904239366952660859694467038, 1e-10, 1e-14));
1406 	assert(approxEqual(besselKR(0.0, 0.5, 3), 813.7490033728880934611898780533576089052725501344392405739625, 1e-10, 1e-14));
1407 	assert(approxEqual(besselKR(1.0, 0.5, 3), 7303.597652823473130198226747312888245046226857159388835630678, 1e-10, 1e-14));
1408 	assert(approxEqual(besselKR(1.0, 1.0, 3), 599.6947228611295581541061895533584099941981855502445314254368, 1e-10, 1e-14));
1409 	assert(approxEqual(besselKR(1.0, 2.0, 3), 67.42923276291368469846380423082891626492668473019928591953262, 1e-10, 1e-14));
1410 	assert(approxEqual(besselKR(0.4, 2.0, 3), 32.55703150637502077170415944139304552928545885327101882556329, 1e-10, 1e-14));
1411 	assert(approxEqual(besselKR(0.4, 1.0, 3), 222.9905656901318819890519502437505989113682776582595951935857, 1e-10, 1e-14));
1412 	assert(approxEqual(besselKR(0.4, 0.5, 3), 2177.389452959348919338879066729351907090043415959389603727847, 1e-10, 1e-14));
1413 }```