1 /**
2  *  Mathlib : A C Library of Special Functions
3  *  Copyright (C) 1998-2014 Ross Ihaka and the R Core team.
4  *  Copyright (C) 2002-3	The R Foundation
5  *  Copyright (C) 2014-3	Ilya Yaroshenko (Ported to D)
6  *
7  *  This program is free software; you can redistribute it and/or modify
8  *  it under the terms of the GNU General Public License as published b
9  *  the Free Software Foundation; either version 2 of the License, or
10  *  (at your option) any later version.
11  *
12  *  This program is distributed in the hope that it will be useful,
13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *  GNU General Public License for more details.
16  *
17  *  You should have received a copy of the GNU General Public License
18  *  along with this program; if not, a copy is available at
19  *  http://www.r-project.org/Licenses/
20  */
21 module bessel;
22 
23 static import std.math;
24 import core.stdc.stdlib;
25 import core.stdc.tgmath;
26 import std.algorithm;
27 import std.exception;
28 import std..string;
29 import std.traits;
30 import std.typecons;
31 
32 public import std.typecons : Flag;
33 
34 /* Documentation that apply to several of the bessel[IJKY] functions */
35 
36 /* *******************************************************************
37 
38  Explanation of machine-dependent constants
39 
40 	beta	  = Radix for the floating-point system
41 	minexp = Smallest representable power of beta
42 	maxexp = Smallest power of beta that overflows
43 	it = p = Number of bits (base-beta digits) in the mantissa
44 		(significand) of a working precision (floating-point) variable
45 	NSIG	  = Decimal significance desired.  Should be set to
46 		INT(LOG10(2)*it+1).	 Setting NSIG lower will result
47 		in decreased accuracy while setting NSIG higher will
48 		increase CPU time without increasing accuracy.  The
49 		truncation error is limited to a relative error of
50 		T=.5*10^(-NSIG).
51 	ENTEN  = 10 ^ K, where K is the largest int such that
52 		ENTEN is machine-representable in working precision
53 	ENSIG  = 10 ^ NSIG
54 	RTNSIG = 10 ^ (-K) for the smallest int K such that K >= NSIG/4
55 	ENMTEN = Smallest ABS(X) such that X/4 does not underflow
56 	XINF	  = Largest positive machine number; approximately beta ^ maxexp
57 		== double.max (defined in  #include <float.h>)
58 	SQXMIN = Square root of beta ^ minexp = sqrt(double.min_normal)
59 
60 	EPS	  = The smallest positive floating-point number such that 1.0+EPS > 1.0
61 	  = beta ^ (-p)	 == double.epsilon
62 
63 
64   For I :
65 
66 	EXPARG = Largest working precision argument that the liaary
67 		EXP routine can handle and upper limit on the
68 		magnitude of X when IZE=1; approximately LOG(beta ^ maxexp)
69 
70   For I and J :
71 
72 	xlrg_IJ = xlrg_BESS_IJ (was = XLARGE). Upper limit on the magnitude of X
73 		(when IZE=2 for I()).  Bear in mind that if floor(abs(x)) =: N, then
74 		at least N iterations of the backward recursion will be executed.
75 		The value of 10 ^ 4 was used till Feb.2009, when it was increased
76 		to 10 ^ 5 (= 1e5).
77 
78   For j :
79 	XMIN_J  = Smallest acceptable argument for RBESY; approximately
80 		max(2*beta ^ minexp, 2/XINF), rounded up
81 
82   For Y :
83 
84 	xlrg_Y =  (was = XLARGE). Upper bound on X;
85 		approximately 1/DEL, because the sine and cosine functions
86 		have lost about half of their precision at that point.
87 
88 	EPS_SINC = Machine number below which sin(x)/x = 1; approximately SQRT(EPS).
89 	THRESH = Lower bound for use of the asymptotic form;
90 		approximately AINT(-LOG10(EPS/2.0))+1.0
91 
92 
93   For K :
94 
95 	xmax_k =  (was = XMAX). Upper limit on the magnitude of X when ize = 1;
96 		i.e. maximal x for UNscaled answer.
97 
98 		Solution to equation:
99 			W(X) * (1 -1/8 X + 9/128 X^2) = beta ^ minexp
100 		where  W(X) = EXP(-X)*SQRT(PI/2X)
101 
102  --------------------------------------------------------------------
103 
104 	 Approximate values for some important machines are:
105 
106 		  beta minexp maxexp it NSIG ENTEN ENSIG RTNSIG ENMTEN	 EXPARG
107  IEEE (IBM/XT,
108 	SUN, etc.) (S.P.)  2	  -126	128  24	  8  1e38	1e8	  1e-2	4.70e-38	 88
109  IEEE	(...) (D.P.)  2	 -1022 1024  53	 16  1e308  1e16  1e-4	8.90e-308	709
110  CRAY-1		  (S.P.)  2	 -8193 8191  48	 15  1e2465 1e15  1e-4	1.84e-2466 5677
111  Cyber 180/855
112 	under NOS  (S.P.)  2	  -975 1070  48	 15  1e322  1e15  1e-4	1.25e-293	741
113  IBM 3033	 (D.P.) 16		-65	 63  14	  5  1e75	1e5	  1e-2	2.16e-78	174
114  VAX		  (S.P.)  2	  -128	127  24	  8  1e38	1e8	  1e-2	1.17e-38	 88
115  VAX D-Format (D.P.)  2	  -128	127  56	 17  1e38	1e17  1e-5	1.17e-38	 88
116  VAX G-Format (D.P.)  2	 -1024 1023  53	 16  1e307  1e16  1e-4	2.22e-308	709
117 
118 
119 And routine specific :
120 
121 			xlrg_IJ xlrg_Y xmax_k EPS_SINC XMIN_J	XINF	THRESH
122  IEEE (IBM/XT,
123 	SUN, etc.) (S.P.)	1e4  1e4	85.337  1e-4	 2.36e-38	3.40e38	8.
124  IEEE	(...) (D.P.)	1e4  1e8  705.342  1e-8	 4.46e-308  1.79e308	16.
125  CRAY-1		  (S.P.)	1e4  2e7 5674.858  5e-8	 3.67e-2466 5.45e2465  15.
126  Cyber 180/855
127 	under NOS  (S.P.)	1e4  2e7  672.788  5e-8	 6.28e-294  1.26e322	15.
128  IBM 3033	 (D.P.)	1e4  1e8  177.852  1e-8	 2.77e-76	7.23e75	17.
129  VAX		  (S.P.)	1e4  1e4	86.715  1e-4	 1.18e-38	1.70e38	8.
130  VAX e-Format (D.P.)	1e4  1e9	86.715  1e-9	 1.18e-38	1.70e38	17.
131  VAX G-Format (D.P.)	1e4  1e8  706.728  1e-8	 2.23e-308  8.98e307	16.
132 
133 */
134 private enum double nsig_BESS = 16;
135 private enum double ensig_BESS = 1e16;
136 private enum double rtnsig_BESS = 1e-4;
137 private enum double enmten_BESS = 8.9e-308;
138 private enum double enten_BESS = 1e308;
139 private enum double exparg_BESS = 709.;
140 private enum double xlrg_BESS_IJ = 1e5;
141 private enum double xlrg_BESS_Y = 1e8;
142 private enum double thresh_BESS_Y = 16.;
143 private enum double xmax_BESS_K = 705.342 /* maximal x for UNscaled answer */ ;
144 /* sqrt(double.min_normal) =	1.491668e-154 */
145 private enum double sqxmin_BESS_K = 1.49e-154;
146 /* x < eps_sinc	 <==>  sin(x)/x == 1 (particularly "==>");
147   Linux (around 2001-02) gives 2.14946906753213e-08
148   Solaris 2.5.1		 gives 2.14911933289084e-08
149 */
150 private enum double M_eps_sinc = 2.149e-8;
151 private enum double M_SQRT_2dPI = 0.797884560802865355879892119869;
152 ///
153 double besselK(double x, double alpha,
154 	Flag!"ExponentiallyScaled" expFlag = Flag!"ExponentiallyScaled".no)
155 {
156 	ptrdiff_t nb;
157 	ptrdiff_t ncalc;
158 	double* b;
159 	/* NaNs propagated correctly */
160 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
161 		return x + alpha;
162 	if (x < 0)
163 		return double.nan;
164 	if (alpha < 0)
165 		alpha = -alpha;
166 	nb = 1 + cast(int) floor(alpha); /* b.length-1 <= |alpha| < b.length */
167 	alpha -= nb - 1;
168 	b = cast(double* ) calloc(nb, double.sizeof).enforce("besselK allocation error");
169 	scope(exit) b.free;
170 	besselK(x, alpha, expFlag, b[0 .. nb], ncalc);
171 	if (ncalc != nb)
172 	{
173 		 /* error input */
174 		if (ncalc < 0)
175 			throw new Exception(format("besselK(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
176 			x, ncalc, nb, alpha));
177 		else throw new Exception(format(
178 			"besselK(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
179 	}
180 	x = b[nb - 1];
181 	return x;
182 }
183 
184 ///
185 unittest
186 {
187 	assert(std.math.approxEqual(besselK(2, 1), 0.139865881816522427284));
188 	assert(std.math.approxEqual(besselK(20, 1), 5.88305796955703817765028e-10));
189 	assert(std.math.approxEqual(besselK(1, 1.2), 0.701080));
190 }
191 
192 /* modified version of besselK that accepts a work array instead of
193 	allocating one. */
194 double besselK(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b)
195 {
196 	ptrdiff_t nb;
197 	ptrdiff_t ncalc;
198 	/* NaNs propagated correctly */
199 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
200 		return x + alpha;
201 	if (x < 0)
202 		return double.nan;
203 	if (alpha < 0)
204 		alpha = -alpha;
205 	nb = 1 + cast(int) floor(alpha); /* b.length-1 <= |alpha| < b.length */
206 	alpha -= nb - 1;
207 	besselK(x, alpha, expFlag, b[0 .. nb], ncalc);
208 	if (ncalc != nb)
209 	{
210 		 /* error input */
211 		if (ncalc < 0)
212 			throw new Exception(format("besselK(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
213 			x, ncalc, nb, alpha));
214 		else throw new Exception(format(
215 			"besselK(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
216 	}
217 	x = b[nb - 1];
218 	return x;
219 }
220 
221 void besselK(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b,
222 	ref ptrdiff_t ncalc)
223 {
224 	/*-------------------------------------------------------------------
225 
226   This routine calculates modified Bessel functions
227   of the third kind, K_(N+ALPHA) (X), for non-negative
228   argument X, and non-negative order N+ALPHA, with or without
229   exponential scaling.
230 
231   Explanation of variables in the calling sequence
232 
233  X	 - Non-negative argument for which
234 	 K's or exponentially scaled K's (K*EXP(X))
235 	 are to be calculated.	If K's are to be calculated,
236 	 X must not be greater than XMAX_BESS_K.
237  ALPHA - Fractional part of order for which
238 	 K's or exponentially scaled K's (K*EXP(X)) are
239 	 to be calculated.  0 <= ALPHA < 1.0.
240  NB	- Number of functions to be calculated, NB > 0.
241 	 The first function calculated is of order ALPHA, and the
242 	 last is of order (NB - 1 + ALPHA).
243  expFlag	- Type.	expFlag = No if unscaled K's are to be calculated,
244 			= Yes if exponentially scaled K's are to be calculated.
245  BK	- Output vector of length NB.	If the
246 	 routine terminates normally (NCALC=NB), the vector BK
247 	 contains the functions K(ALPHA,X), ... , K(NB-1+ALPHA,X),
248 	 or the corresponding exponentially scaled functions.
249 	 If (0 < NCALC < NB), BK(I) contains correct function
250 	 values for I <= NCALC, and contains the ratios
251 	 K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
252  NCALC - Output variable indicating possible errors.
253 	 Before using the vector BK, the user should check that
254 	 NCALC=NB, i.e., all orders have been calculated to
255 	 the desired accuracy.	See error returns below.
256 
257 
258  *******************************************************************
259 
260  Error returns
261 
262   In case of an error, NCALC != NB, and not all K's are
263   calculated to the desired accuracy.
264 
265   NCALC < -1:  An argument is out of range. For example,
266 	NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= XMAX_BESS_K.
267 	In this case, the B-vector is not calculated,
268 	and NCALC is set to MIN0(NB,0)-2	 so that NCALC != NB.
269   NCALC = -1:  Either  K(ALPHA,X) >= XINF  or
270 	K(ALPHA+NB-1,X)/K(ALPHA+NB-2,X) >= XINF.	 In this case,
271 	the B-vector is not calculated.	Note that again
272 	NCALC != NB.
273 
274   0 < NCALC < NB: Not all requested function values could
275 	be calculated accurately.  BK(I) contains correct function
276 	values for I <= NCALC, and contains the ratios
277 	K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
278 
279 
280  Intrinsic functions required are:
281 
282 	 ABS, AINT, EXP, INT, LOG, MAX, MIN, SINH, SQRT
283 
284 
285  Acknowledgement
286 
287 	This program is based on a program written b J. B. Campbell
288 	(2) that computes values of the Bessel functions K of float
289 	argument and float order.  Modifications include the addition
290 	of non-scaled functions, parameterization of machine
291 	dependencies, and the use of more accurate approximations
292 	for SINH and SIN.
293 
294  References: "On Temme's Algorithm for the Modified Bessel
295 		  Functions of the Third Kind," Campbell, J. B.,
296 		  TOMS 6(4), Dec. 1980, pp. 581-586.
297 
298 		 "A FORTRAN IV Suaoutine for the Modified Bessel
299 		  Functions of the Third Kind of Real Order and Real
300 		  Argument," Campbell, J. B., Report NRC/ERB-925,
301 		  National Research Council, Canada.
302 
303   Latest modification: May 30, 1989
304 
305   Modified b: W. J. Cody and L. Stoltz
306 			Applied Mathematics Division
307 			Argonne National Laboratory
308 			Argonne, IL  60439
309 
310  -------------------------------------------------------------------
311 */
312 	
313 	/*---------------------------------------------------------------------
314 	 * Mathematical constants
315 	 *	A = LOG(2) - Euler's constant
316 	 *	D = SQRT(2/PI)
317 	 ---------------------------------------------------------------------*/
318 	static immutable double a = .11593151565841244881;
319 	/*---------------------------------------------------------------------
320 	  P, Q - Approximation for LOG(GAMMA(1+ALPHA))/ALPHA + Euler's constant
321 	  Coefficients converted from hex to decimal and modified
322 	  b W. J. Cody, 2/26/82 */
323 	static immutable double[8] p = [.805629875690432845, 20.4045500205365151,
324 		157.705605106676174, 536.671116469207504, 900.382759291288778,
325 		730.923886650660393, 229.299301509425145, .822467033424113231];
326 	static immutable double[7] q = [29.4601986247850434, 277.577868510221208,
327 		1206.70325591027438, 2762.91444159791519, 3443.74050506564618,
328 		2210.63190113378647, 572.267338359892221];
329 	/* R, S - Approximation for (1-ALPHA*PI/SIN(ALPHA*PI))/(2.D0*ALPHA) */
330 	static immutable double[5] r = [-.48672575865218401848,
331 		13.079485869097804016,
332 		-101.96490580880537526, 347.65409106507813131, 3.495898124521934782e-4];
333 	static immutable double[4] s = [-25.579105509976461286,
334 		212.57260432226544008, -610.69018684944109624, 422.69668805777760407];
335 	/* T	- Approximation for SINH(Y)/Y */
336 	static immutable double[6] t = [1.6125990452916363814e-10,
337 		2.5051878502858255354e-8, 2.7557319615147964774e-6,
338 		1.9841269840928373686e-4, .0083333333333334751799,
339 		.16666666666666666446];
340 	/*---------------------------------------------------------------------*/
341 	static immutable double[6] estm = [52.0583, 5.7607, 2.7782, 14.4303, 185.3004,
342 		9.3715];
343 	static immutable double[7] estf = [41.8341, 7.1075, 6.4306, 42.511, 1.35633,
344 		84.5096, 20.];
345 	/* Local variables */
346 	ptrdiff_t iend, i, j, k, m, ii, mplus1;
347 	double x2b4, twox, c, blpha, ratio, wminf;
348 	double d1, d2, d3, f0, f1, f2, p0, q0, t1, t2, twonu;
349 	double dm, ex, b1, b2, nu;
350 	ii = 0; /* -Wall */
351 	ex = x;
352 	nu = alpha;
353 	ncalc = -2;
354 	if (b.length > 0 && (0. <= nu && nu < 1.))
355 	{
356 		if (ex <= 0 || (expFlag == false  && ex > xmax_BESS_K))
357 		{
358 			if (ex <= 0)
359 			{
360 				enforce(ex == 0, "besselK Range Error");
361 				for (i = 0; i < b.length; i++)
362 					b[i] = double.infinity;
363 			}
364 			else  /* would only have underflow */
365 			for (i = 0; i < b.length; i++)
366 				b[i] = 0;
367 			ncalc = b.length;
368 			return;
369 		}
370 		k = 0;
371 		if (nu < sqxmin_BESS_K)
372 		{
373 			nu = 0;
374 		}
375 		else if (nu > .5)
376 		{
377 			k = 1;
378 			nu -= 1.;
379 		}
380 		twonu = nu + nu;
381 		iend = b.length + k - 1;
382 		c = nu * nu;
383 		d3 = -c;
384 		if (ex <= 1.)
385 		{
386 			/* ------------------------------------------------------------
387 			Calculation of P0 = GAMMA(1+ALPHA) * (2/X)**ALPHA
388 				  Q0 = GAMMA(1-ALPHA) * (X/2)**ALPHA
389 			------------------------------------------------------------ */
390 			d1 = 0;
391 			d2 = p[0];
392 			t1 = 1.;
393 			t2 = q[0];
394 			for (i = 2; i <= 7; i += 2)
395 			{
396 				d1 = c * d1 + p[i - 1];
397 				d2 = c * d2 + p[i];
398 				t1 = c * t1 + q[i - 1];
399 				t2 = c * t2 + q[i];
400 			}
401 			d1 = nu * d1;
402 			t1 = nu * t1;
403 			f1 = log(ex);
404 			f0 = a + nu * (p[7] - nu * (d1 + d2) / (t1 + t2)) - f1;
405 			q0 = exp(-nu * (a - nu * (p[7] + nu * (d1 - d2) / (t1 - t2)) - f1));
406 			f1 = nu * f0;
407 			p0 = exp(f1);
408 			/* -----------------------------------------------------------
409 			Calculation of F0 =
410 			----------------------------------------------------------- */
411 			d1 = r[4];
412 			t1 = 1.;
413 			for (i = 0; i < 4; ++i)
414 			{
415 				d1 = c * d1 + r[i];
416 				t1 = c * t1 + s[i];
417 			}
418 			/* d2 := sinh(f1)/ nu = sinh(f1)/(f1/f0)
419 		 *		= f0 * sinh(f1)/f1 */
420 			if (fabs(f1) <= .5)
421 			{
422 				f1 *= f1;
423 				d2 = 0;
424 				for (i = 0; i < 6; ++i)
425 				{
426 					d2 = f1 * d2 + t[i];
427 				}
428 				d2 = f0 + f0 * f1 * d2;
429 			}
430 			else
431 			{
432 				d2 = sinh(f1) / nu;
433 			}
434 			f0 = d2 - nu * d1 / (t1 * p0);
435 			if (ex <= 1e-10)
436 			{
437 				/* ---------------------------------------------------------
438 			X <= 1.0E-10
439 			Calculation of K(ALPHA,X) and X*K(ALPHA+1,X)/K(ALPHA,X)
440 			--------------------------------------------------------- */
441 				b[0] = f0 + ex * f0;
442 				if (expFlag == false)
443 				{
444 					b[0] -= ex * b[0];
445 				}
446 				ratio = p0 / f0;
447 				c = ex * double.max;
448 				if (k != 0)
449 				{
450 					/* ---------------------------------------------------
451 				Calculation of K(ALPHA,X)
452 				and  X*K(ALPHA+1,X)/K(ALPHA,X),	ALPHA >= 1/2
453 				--------------------------------------------------- */
454 					ncalc = -1;
455 					if (b[0] >= c / ratio)
456 					{
457 						return;
458 					}
459 					b[0] = ratio * b[0] / ex;
460 					twonu += 2.;
461 					ratio = twonu;
462 				}
463 				ncalc = 1;
464 				if (b.length == 1)
465 					return;
466 				/* -----------------------------------------------------
467 			Calculate  K(ALPHA+L,X)/K(ALPHA+L-1,X),
468 			L = 1, 2, ... , NB-1
469 			----------------------------------------------------- */
470 				ncalc = -1;
471 				for (i = 1; i < b.length; ++i)
472 				{
473 					if (ratio >= c)
474 						return;
475 					b[i] = ratio / ex;
476 					twonu += 2.;
477 					ratio = twonu;
478 				}
479 				ncalc = 1;
480 				goto L420;
481 			}
482 			else
483 			{
484 				/* ------------------------------------------------------
485 			10^-10 < X <= 1.0
486 			------------------------------------------------------ */
487 				c = 1.;
488 				x2b4 = ex * ex / 4.;
489 				p0 = .5 * p0;
490 				q0 = .5 * q0;
491 				d1 = -1.;
492 				d2 = 0;
493 				b1 = 0;
494 				b2 = 0;
495 				f1 = f0;
496 				f2 = p0;
497 				do
498 				{
499 					d1 += 2.;
500 					d2 += 1.;
501 					d3 = d1 + d3;
502 					c = x2b4 * c / d2;
503 					f0 = (d2 * f0 + p0 + q0) / d3;
504 					p0 /= d2 - nu;
505 					q0 /= d2 + nu;
506 					t1 = c * f0;
507 					t2 = c * (p0 - d2 * f0);
508 					b1 += t1;
509 					b2 += t2;
510 				}
511 				while (fabs(t1 / (f1 + b1)) > double.epsilon || fabs(t2 / (
512 					f2 + b2)) > double.epsilon);
513 				b1 = f1 + b1;
514 				b2 = 2. * (f2 + b2) / ex;
515 				if (expFlag)
516 				{
517 					d1 = exp(ex);
518 					b1 *= d1;
519 					b2 *= d1;
520 				}
521 				wminf = estf[0] * ex + estf[1];
522 			}
523 		}
524 		else if (double.epsilon * ex > 1.)
525 		{
526 			/* -------------------------------------------------
527 			X > 1./EPS
528 			------------------------------------------------- */
529 			ncalc = b.length;
530 			b1 = 1. / (M_SQRT_2dPI * sqrt(ex));
531 			for (i = 0; i < b.length; ++i)
532 				b[i] = b1;
533 			return;
534 		}
535 		else
536 		{
537 			/* -------------------------------------------------------
538 			X > 1.0
539 			------------------------------------------------------- */
540 			twox = ex + ex;
541 			blpha = 0;
542 			ratio = 0;
543 			if (ex <= 4.)
544 			{
545 			/* ----------------------------------------------------------
546 			Calculation of K(ALPHA+1,X)/K(ALPHA,X),  1.0 <= X <= 4.0
547 			----------------------------------------------------------*/
548 				d2 = trunc(estm[0] / ex + estm[1]);
549 				m = cast(int) d2;
550 				d1 = d2 + d2;
551 				d2 -= .5;
552 				d2 *= d2;
553 				for (i = 2; i <= m; ++i)
554 				{
555 					d1 -= 2.;
556 					d2 -= d1;
557 					ratio = (d3 + d2) / (twox + d1 - ratio);
558 				}
559 				/* -----------------------------------------------------------
560 			Calculation of I(|ALPHA|,X) and I(|ALPHA|+1,X) by backward
561 			recurrence and K(ALPHA,X) from the wronskian
562 			-----------------------------------------------------------*/
563 				d2 = trunc(estm[2] * ex + estm[3]);
564 				m = cast(int) d2;
565 				c = fabs(nu);
566 				d3 = c + c;
567 				d1 = d3 - 1.;
568 				f1 = double.min_normal;
569 				f0 = (2. * (c + d2) / ex + .5 * ex / (c + d2 + 1.)) * double.min_normal;
570 				for (i = 3; i <= m; ++i)
571 				{
572 					d2 -= 1.;
573 					f2 = (d3 + d2 + d2) * f0;
574 					blpha = (1. + d1 / d2) * (f2 + blpha);
575 					f2 = f2 / ex + f1;
576 					f1 = f0;
577 					f0 = f2;
578 				}
579 				f1 = (d3 + 2.) * f0 / ex + f1;
580 				d1 = 0;
581 				t1 = 1.;
582 				for (i = 1; i <= 7; ++i)
583 				{
584 					d1 = c * d1 + p[i - 1];
585 					t1 = c * t1 + q[i - 1];
586 				}
587 				p0 = exp(c * (a + c * (p[7] - c * d1 / t1) - log(ex))) / ex;
588 				f2 = (c + .5 - ratio) * f1 / ex;
589 				b1 = p0 + (d3 * f0 - f2 + f0 + blpha) / (f2 + f1 + f0) * p0;
590 				if (expFlag == false)
591 				{
592 					b1 *= exp(-ex);
593 				}
594 				wminf = estf[2] * ex + estf[3];
595 			}
596 			else
597 			{
598 				/* ---------------------------------------------------------
599 			Calculation of K(ALPHA,X) and K(ALPHA+1,X)/K(ALPHA,X), by
600 			backward recurrence, for  X > 4.0
601 			----------------------------------------------------------*/
602 				dm = trunc(estm[4] / ex + estm[5]);
603 				m = cast(int) dm;
604 				d2 = dm - .5;
605 				d2 *= d2;
606 				d1 = dm + dm;
607 				for (i = 2; i <= m; ++i)
608 				{
609 					dm -= 1.;
610 					d1 -= 2.;
611 					d2 -= d1;
612 					ratio = (d3 + d2) / (twox + d1 - ratio);
613 					blpha = (ratio + ratio * blpha) / dm;
614 				}
615 				b1 = 1. / ((M_SQRT_2dPI + M_SQRT_2dPI * blpha) * sqrt(ex));
616 				if (expFlag == false)
617 					b1 *= exp(-ex);
618 				wminf = estf[4] * (ex - fabs(ex - estf[6])) + estf[5];
619 			}
620 			/* ---------------------------------------------------------
621 			Calculation of K(ALPHA+1,X)
622 			from K(ALPHA,X) and  K(ALPHA+1,X)/K(ALPHA,X)
623 			--------------------------------------------------------- */
624 			b2 = b1 + b1 * (nu + .5 - ratio) / ex;
625 		}
626 		/*--------------------------------------------------------------------
627 	  Calculation of 'NCALC', K(ALPHA+I,X),	I  =  0, 1, ... , NCALC-1,
628 	  &	  K(ALPHA+I,X)/K(ALPHA+I-1,X),	I = NCALC, NCALC+1, ... , NB-1
629 	  -------------------------------------------------------------------*/
630 		ncalc = b.length;
631 		b[0] = b1;
632 		if (iend == 0)
633 			return;
634 		j = 1 - k;
635 		if (j >= 0)
636 			b[j] = b2;
637 		if (iend == 1)
638 			return;
639 		m = min(cast(int)(wminf - nu), iend);
640 		for (i = 2; i <= m; ++i)
641 		{
642 			t1 = b1;
643 			b1 = b2;
644 			twonu += 2.;
645 			if (ex < 1.)
646 			{
647 				if (b1 >= double.max / twonu * ex)
648 					break;
649 			}
650 			else
651 			{
652 				if (b1 / ex >= double.max / twonu)
653 					break;
654 			}
655 			b2 = twonu / ex * b1 + t1;
656 			ii = i;
657 			++j;
658 			if (j >= 0)
659 			{
660 				b[j] = b2;
661 			}
662 		}
663 		m = ii;
664 		if (m == iend)
665 		{
666 			return;
667 		}
668 		ratio = b2 / b1;
669 		mplus1 = m + 1;
670 		ncalc = -1;
671 		for (i = mplus1; i <= iend; ++i)
672 		{
673 			twonu += 2.;
674 			ratio = twonu / ex + 1. / ratio;
675 			++j;
676 			if (j >= 1)
677 			{
678 				b[j] = ratio;
679 			}
680 			else
681 			{
682 				if (b2 >= double.max / ratio)
683 					return;
684 				b2 *= ratio;
685 			}
686 		}
687 		ncalc = max(1, mplus1 - k);
688 		if (ncalc == 1)
689 			b[0] = b2;
690 		if (b.length == 1)
691 			return;
692 		L420: for (i = ncalc; i < b.length; ++i)
693 		{
694 			 /* i == ncalc */
695 			if (b[i - 1] >= double.max / b[i])
696 				return;
697 			b[i] *= b[i - 1];
698 			ncalc++;
699 		}
700 	}
701 }
702 
703 ///
704 double besselI(double x, double alpha,
705 	Flag!"ExponentiallyScaled" expFlag = Flag!"ExponentiallyScaled".no)
706 {
707 	ptrdiff_t nb;
708 	ptrdiff_t ncalc;
709 	double na;
710 	double* b;
711 	/* NaNs propagated correctly */
712 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
713 		return x + alpha;
714 	if (x < 0)
715 		return double.nan;
716 	na = floor(alpha);
717 	if (alpha < 0)
718 	{
719 		/* Using Aaamowitz & Stegun  9.6.2 & 9.6.6
720 	 * this may not be quite optimal (CPU and accuracy wise) */
721 		return (besselI(x, -alpha, expFlag) + ((alpha == na) ?  /* sin(pi * alpha) = 0 */ 0 : besselK(x, -alpha,
722 			expFlag) * (expFlag == true  ? 2. : 2. * exp(-2.0 * x)) / cast(double) std.math.PI * sinPi(-alpha)));
723 	}
724 	nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */
725 	alpha -= nb - 1;
726 	b = cast(double* ) calloc(nb, double.sizeof).enforce("besselI allocation error");
727 	scope(exit) b.free;
728 	besselI(x, alpha, expFlag, b[0 .. nb], ncalc);
729 	if (ncalc != nb)
730 	{
731 		 /* error input */
732 		if (ncalc < 0)
733 			throw new Exception(format("besselI(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
734 			x, ncalc, nb, alpha));
735 		else throw new Exception(format(
736 			"besselI(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
737 	}
738 	x = b[nb - 1];
739 	return x;
740 }
741 
742 ///
743 unittest
744 {
745 	assert(std.math.approxEqual(besselI(2, 1),
746 		1.590636854637329063382254424999666));
747 	assert(std.math.approxEqual(besselI(20, 1),
748 		4.245497338512777018140990665e+7));
749 	assert(std.math.approxEqual(besselI(1, 1.2), 0.441739));
750 }
751 
752 /* modified version of besselI that accepts a work array instead of
753 	allocating one. */
754 double besselI(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b)
755 {
756 	ptrdiff_t nb;
757 	ptrdiff_t ncalc;
758 	double na;
759 	/* NaNs propagated correctly */
760 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
761 		return x + alpha;
762 	if (x < 0)
763 		return double.nan;
764 	na = floor(alpha);
765 	if (alpha < 0)
766 	{
767 		/* Using Aaamowitz & Stegun  9.6.2 & 9.6.6
768 	 * this may not be quite optimal (CPU and accuracy wise) */
769 		return (besselI(x, -alpha, expFlag, b) + ((alpha == na) ? 0 : besselK(x,
770 			-alpha, expFlag, b) * (expFlag == true  ? 2. : 2. * exp(-2.0 * x)) / cast(
771 			double) std.math.PI * sinPi(-alpha)));
772 	}
773 	nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */
774 	alpha -= nb - 1;
775 	besselI(x, alpha, expFlag, b[0 .. nb], ncalc);
776 	if (ncalc != nb)
777 	{
778 		 /* error input */
779 		if (ncalc < 0)
780 			throw new Exception(format("besselI(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
781 			x, ncalc, nb, alpha));
782 		else throw new Exception(format(
783 			"besselI(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
784 	}
785 	x = b[nb - 1];
786 	return x;
787 }
788 
789 void besselI(double x, double alpha, Flag!"ExponentiallyScaled" expFlag, double[] b,
790 	ref ptrdiff_t ncalc)
791 {
792 	/* -------------------------------------------------------------------
793 
794  This routine calculates Bessel functions I_(N+ALPHA) (X)
795  for non-negative argument X, and non-negative order N+ALPHA,
796  with or without exponential scaling.
797 
798 
799  Explanation of variables in the calling sequence
800 
801  X	 - Non-negative argument for which
802 	 I's or exponentially scaled I's (I*EXP(-X))
803 	 are to be calculated.	If I's are to be calculated,
804 	 X must be less than exparg_BESS (IZE=1) or xlrg_BESS_IJ (IZE=2),
805 	 (see bessel.h).
806  ALPHA - Fractional part of order for which
807 	 I's or exponentially scaled I's (I*EXP(-X)) are
808 	 to be calculated.  0 <= ALPHA < 1.0.
809  NB	- Number of functions to be calculated, NB > 0.
810 	 The first function calculated is of order ALPHA, and the
811 	 last is of order (NB - 1 + ALPHA).
812  IZE	- Type.	IZE = 1 if unscaled I's are to be calculated,
813 			= 2 if exponentially scaled I's are to be calculated.
814  BI	- Output vector of length NB.	If the routine
815 	 terminates normally (NCALC=NB), the vector BI contains the
816 	 functions I(ALPHA,X) through I(NB-1+ALPHA,X), or the
817 	 corresponding exponentially scaled functions.
818  NCALC - Output variable indicating possible errors.
819 	 Before using the vector BI, the user should check that
820 	 NCALC=NB, i.e., all orders have been calculated to
821 	 the desired accuracy.	See error returns below.
822 
823 
824  *******************************************************************
825  *******************************************************************
826 
827  Error returns
828 
829   In case of an error,	NCALC != NB, and not all I's are
830   calculated to the desired accuracy.
831 
832   NCALC < 0:  An argument is out of range. For example,
833 	 NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >= EXPARG_BESS.
834 	 In this case, the BI-vector is not calculated, and NCALC is
835 	 set to MIN0(NB,0)-1 so that NCALC != NB.
836 
837   NB > NCALC > 0: Not all requested function values could
838 	 be calculated accurately.	This usually occurs because NB is
839 	 much larger than ABS(X).  In this case, BI[N] is calculated
840 	 to the desired accuracy for N <= NCALC, but precision
841 	 is lost for NCALC < N <= NB.  If BI[N] does not vanish
842 	 for N > NCALC (because it is too small to be represented),
843 	 and BI[N]/BI[NCALC] = 10**(-K), then only the first NSIG-K
844 	 significant figures of BI[N] can be trusted.
845 
846 
847  Intrinsic functions required are:
848 
849 	 DBLE, EXP, gamma_cody, INT, MAX, MIN, REAL, SQRT
850 
851 
852  Acknowledgement
853 
854   This program is based on a program written b David J.
855   Sookne (2) that computes values of the Bessel functions J or
856   I of float argument and long order.  Modifications include
857   the restriction of the computation to the I Bessel function
858   of non-negative float argument, the extension of the computation
859   to arbtrary positive order, the inclusion of optional
860   exponential scaling, and the elimination of most underflow.
861   An earlier version was published in (3).
862 
863  References: "A Note on Backward Recurrence Algorithms," Olver,
864 		  F. W. J., and Sookne, D. J., Math. Comp. 26, 1972,
865 		  pp 941-947.
866 
867 		 "Bessel Functions of Real Argument and Integer Order,"
868 		  Sookne, D. J., NBS Jour. of Res. B. 77B, 1973, pp
869 		  125-132.
870 
871 		 "ALGORITHM 597, Sequence of Modified Bessel Functions
872 		  of the First Kind," Cody, W. J., Trans. Math. Soft.,
873 		  1983, pp. 242-245.
874 
875   Latest modification: May 30, 1989
876 
877   Modified b: W. J. Cody and L. Stoltz
878 			Applied Mathematics Division
879 			Argonne National Laboratory
880 			Argonne, IL  60439
881 */
882 	
883 	/*-------------------------------------------------------------------
884 	  Mathematical constants
885 	  -------------------------------------------------------------------*/
886 	enum double const__ = 1.585;
887 	/* Local variables */
888 	ptrdiff_t nend, intx, nbmx, k, l, n, nstart;
889 	double pold, test, p, em, en, empal, emp2al, halfx, aa, bb, cc, psave, plast, tover,
890 		psavel, sum, nu, twonu;
891 	/*Parameter adjustments */
892 	//	--b;
893 	nu = alpha;
894 	twonu = nu + nu;
895 	/*-------------------------------------------------------------------
896 	  Check for X, NB, OR IZE out of range.
897 	  ------------------------------------------------------------------- */
898 	if (b.length > 0 && x >= 0. && (0. <= nu && nu < 1.))
899 	{
900 		ncalc = b.length;
901 		if (expFlag == false  && x > exparg_BESS)
902 		{
903 			b[] = double.infinity; /* the limit *is* = Inf */
904 			return;
905 		}
906 		if (expFlag && x > xlrg_BESS_IJ)
907 		{
908 			b[] = 0; /* The limit exp(-x) * I_nu(x) --> 0 : */
909 			return;
910 		}
911 		intx = cast(int)(x); /* fine, since x <= xlrg_BESS_IJ <<< LONG_MAX */
912 		if (x >= rtnsig_BESS)
913 		{
914 			 /* "non-small" x ( >= 1e-4 ) */
915 			/* -------------------------------------------------------------------
916 	Initialize the forward sweep, the P-sequence of Olver
917 	------------------------------------------------------------------- */
918 			nbmx = b.length - intx;
919 			n = intx + 1;
920 			en = cast(double)(n + n) + twonu;
921 			plast = 1.;
922 			p = en / x;
923 			/* ------------------------------------------------
924 			Calculate general significance test
925 			------------------------------------------------ */
926 			test = ensig_BESS + ensig_BESS;
927 			if (intx << 1 > nsig_BESS * 5)
928 			{
929 				test = sqrt(test * p);
930 			}
931 			else
932 			{
933 				test /= const__ ^^ intx;
934 			}
935 			if (nbmx >= 3)
936 			{
937 				/* --------------------------------------------------
938 			Calculate P-sequence until N = NB-1
939 			Check for possible overflow.
940 			------------------------------------------------ */
941 				tover = enten_BESS / ensig_BESS;
942 				nstart = intx + 2;
943 				nend = b.length - 1;
944 				for (k = nstart; k <= nend; ++k)
945 				{
946 					n = k;
947 					en += 2.;
948 					pold = plast;
949 					plast = p;
950 					p = en * plast / x + pold;
951 					if (p > tover)
952 					{
953 						/* ------------------------------------------------
954 				To avoid overflow, divide P-sequence b TOVER.
955 				Calculate P-sequence until ABS(P) > 1.
956 				---------------------------------------------- */
957 						tover = enten_BESS;
958 						p /= tover;
959 						plast /= tover;
960 						psave = p;
961 						psavel = plast;
962 						nstart = n + 1;
963 						do
964 						{
965 							++n;
966 							en += 2.;
967 							pold = plast;
968 							plast = p;
969 							p = en * plast / x + pold;
970 						}
971 						while (p <= 1.);
972 						bb = en / x;
973 						/* ------------------------------------------------
974 				Calculate backward test, and find NCALC,
975 				the highest N such that the test is passed.
976 				------------------------------------------------ */
977 						test = pold * plast / ensig_BESS;
978 						test *= .5 - .5 / (bb * bb);
979 						p = plast * tover;
980 						--n;
981 						en -= 2.;
982 						nend = min(b.length, n);
983 						for (l = nstart; l <= nend; ++l)
984 						{
985 							ncalc = l;
986 							pold = psavel;
987 							psavel = psave;
988 							psave = en * psavel / x + pold;
989 							if (psave * psavel > test)
990 							{
991 								goto L90;
992 							}
993 						}
994 						ncalc = nend + 1;
995 						L90: ncalc--;
996 						goto L120;
997 					}
998 				}
999 				n = nend;
1000 				en = cast(double)(n + n) + twonu;
1001 				/*---------------------------------------------------
1002 		  Calculate special significance test for NBMX > 2.
1003 		  --------------------------------------------------- */
1004 				test = max(test, sqrt(plast * ensig_BESS) * sqrt(p + p));
1005 			}
1006 			/* --------------------------------------------------------
1007 			Calculate P-sequence until significance test passed.
1008 			-------------------------------------------------------- */
1009 			do
1010 			{
1011 				++n;
1012 				en += 2.;
1013 				pold = plast;
1014 				plast = p;
1015 				p = en * plast / x + pold;
1016 			}
1017 			while (p < test);
1018 			L120: /* -------------------------------------------------------------------
1019  Initialize the backward recursion and the normalization sum.
1020  ------------------------------------------------------------------- */
1021 			++n;
1022 			en += 2.;
1023 			bb = 0;
1024 			aa = 1. / p;
1025 			em = cast(double) n - 1.;
1026 			empal = em + nu;
1027 			emp2al = em - 1. + twonu;
1028 			sum = aa * empal * emp2al / em;
1029 			nend = n - b.length;
1030 			if (nend < 0)
1031 			{
1032 				/* -----------------------------------------------------
1033 			N < NB, so store BI[N] and set higher orders to 0..
1034 			----------------------------------------------------- */
1035 				b[n - 1] = aa;
1036 				nend = -nend;
1037 				b[n .. n + nend] = 0;
1038 			}
1039 			else
1040 			{
1041 				if (nend > 0)
1042 				{
1043 					/* -----------------------------------------------------
1044 				Recur backward via difference equation,
1045 				calculating (but not storing) BI[N], until N = NB.
1046 				--------------------------------------------------- */
1047 					for (l = 1; l <= nend; ++l)
1048 					{
1049 						--n;
1050 						en -= 2.;
1051 						cc = bb;
1052 						bb = aa;
1053 						/* for x ~= 1500,  sum would overflow to 'inf' here,
1054 			 * and the final b[] /= sum would give 0 wrongly;
1055 			 * RE-normalize (aa, sum) here -- no need to undo */
1056 						if (nend > 100 && aa > 1e200)
1057 						{
1058 							/* multiply b  2^-900 = 1.18e-271 */
1059 							cc = ldexp(cc, -900);
1060 							bb = ldexp(bb, -900);
1061 							sum = ldexp(sum, -900);
1062 						}
1063 						aa = en * bb / x + cc;
1064 						em -= 1.;
1065 						emp2al -= 1.;
1066 						if (n == 1)
1067 						{
1068 							break;
1069 						}
1070 						if (n == 2)
1071 						{
1072 							emp2al = 1.;
1073 						}
1074 						empal -= 1.;
1075 						sum = (sum + aa * empal) * emp2al / em;
1076 					}
1077 				}
1078 				/* ---------------------------------------------------
1079 			Store BI[NB]
1080 			--------------------------------------------------- */
1081 				b[n - 1] = aa;
1082 				if (b.length <= 1)
1083 				{
1084 					sum = sum + sum + aa;
1085 					goto L230;
1086 				}
1087 				/* -------------------------------------------------
1088 			Calculate and Store BI[NB-1]
1089 			------------------------------------------------- */
1090 				--n;
1091 				en -= 2.;
1092 				b[n - 1] = en * aa / x + bb;
1093 				if (n == 1)
1094 				{
1095 					goto L220;
1096 				}
1097 				em -= 1.;
1098 				if (n == 2)
1099 					emp2al = 1.;
1100 				else emp2al -= 1.;
1101 				empal -= 1.;
1102 				sum = (sum + b[n - 1] * empal) * emp2al / em;
1103 			}
1104 			nend = n - 2;
1105 			if (nend > 0)
1106 			{
1107 				/* --------------------------------------------
1108 				Calculate via difference equation
1109 				and store BI[N], until N = 2.
1110 				------------------------------------------ */
1111 				for (l = 1; l <= nend; ++l)
1112 				{
1113 					--n;
1114 					en -= 2.;
1115 					b[n - 1] = en * b[n + 0] / x + b[n + 1];
1116 					em -= 1.;
1117 					if (n == 2)
1118 						emp2al = 1.;
1119 					else emp2al -= 1.;
1120 					empal -= 1.;
1121 					sum = (sum + b[n - 1] * empal) * emp2al / em;
1122 				}
1123 			}
1124 			/* ----------------------------------------------
1125 			Calculate BI[1]
1126 			-------------------------------------------- */
1127 			b[0] = 2. * empal * b[1] / x + b[2];
1128 			L220: sum = sum + sum + b[0];
1129 			L230: /* ---------------------------------------------------------
1130 			Normalize.  Divide all BI[N] b sum.
1131 			--------------------------------------------------------- */
1132 			if (nu != 0.)
1133 				sum *= (gamma_cody(1. + nu) * pow(x * .5, -nu));
1134 			if (expFlag == false)
1135 				sum *= exp(-(x));
1136 			aa = enmten_BESS;
1137 			if (sum > 1.)
1138 				aa *= sum;
1139 			for (n = 1; n <= b.length; ++n)
1140 			{
1141 				if (b[n - 1] < aa)
1142 					b[n - 1] = 0;
1143 				else b[n - 1] /= sum;
1144 			}
1145 			return;
1146 		}
1147 		else
1148 		{
1149 			 /* small x  < 1e-4 */
1150 			/* -----------------------------------------------------------
1151 			Two-term ascending series for small X.
1152 			-----------------------------------------------------------*/
1153 			aa = 1.;
1154 			empal = 1. + nu;
1155 			/* No need to check for underflow */
1156 			halfx = .5 * x;
1157 			if (nu != 0.)
1158 				aa = pow(halfx, nu) / gamma_cody(empal);
1159 			if (expFlag)
1160 				aa *= exp(-(x));
1161 			bb = halfx * halfx;
1162 			b[0] = aa + aa * bb / empal;
1163 			if (x != 0. && b[0] == 0.)
1164 				ncalc = 0;
1165 			if (b.length > 1)
1166 			{
1167 				if (x == 0.)
1168 				{
1169 					b[1 .. $] = 0;
1170 				}
1171 				else
1172 				{
1173 					/* -------------------------------------------------
1174 				Calculate higher-order functions.
1175 				------------------------------------------------- */
1176 					cc = halfx;
1177 					tover = (enmten_BESS + enmten_BESS) / x;
1178 					if (bb != 0.)
1179 						tover = enmten_BESS / bb;
1180 					for (n = 2; n <= b.length; ++n)
1181 					{
1182 						aa /= empal;
1183 						empal += 1.;
1184 						aa *= cc;
1185 						if (aa <= tover * empal)
1186 							b[n - 1] = aa = 0;
1187 						else b[n - 1] = aa + aa * bb / empal;
1188 						if (b[n - 1] == 0. && ncalc > n)
1189 							ncalc = n - 1;
1190 					}
1191 				}
1192 			}
1193 		}
1194 	}
1195 	else
1196 	{
1197 		 /* argument out of range */
1198 		ncalc = min(b.length, 0) - 1;
1199 	}
1200 }
1201 
1202 ///
1203 double besselJ(double x, double alpha)
1204 {
1205 	ptrdiff_t nb, ncalc;
1206 	double na;
1207 	double* b;
1208 	/* NaNs propagated correctly */
1209 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
1210 		return x + alpha;
1211 	if (x < 0)
1212 		return double.nan;
1213 	na = floor(alpha);
1214 	if (alpha < 0)
1215 	{
1216 		/* Using Aaamowitz & Stegun  9.1.2
1217 	 * this may not be quite optimal (CPU and accuracy wise) */
1218 		return (besselJ(x, -alpha) * cosPi(alpha) + ((alpha == na) ? 0 : besselY(x,
1219 			-alpha) * sinPi(alpha)));
1220 	}
1221 	nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */
1222 	alpha -= nb - 1;
1223 	b = cast(double* ) calloc(nb, double.sizeof).enforce("besselJ allocation error");
1224 	scope(exit) b.free;
1225 	besselJ(x, alpha, b[0 .. nb], ncalc);
1226 	if (ncalc != nb)
1227 	{
1228 		 /* error input */
1229 		if (ncalc < 0)
1230 			throw new Exception(format("besselJ(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
1231 			x, ncalc, nb, alpha));
1232 		else throw new Exception(format(
1233 			"besselJ(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
1234 	}
1235 	x = b[nb - 1];
1236 	return x;
1237 }
1238 
1239 ///
1240 unittest
1241 {
1242 	assert(std.math.approxEqual(besselJ(2, 1),
1243 		0.5767248077568733872024482422691));
1244 	assert(std.math.approxEqual(besselJ(20, 1),
1245 		0.066833124175850045578992974193));
1246 	assert(std.math.approxEqual(besselJ(1, 1.2), 0.351884));
1247 }
1248 
1249 /* modified version of besselJ that accepts a work array instead of
1250 	allocating one. */
1251 double besselJ(double x, double alpha, double[] b)
1252 {
1253 	ptrdiff_t nb, ncalc;
1254 	double na;
1255 	/* NaNs propagated correctly */
1256 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
1257 		return x + alpha;
1258 	if (x < 0)
1259 		return double.nan;
1260 	na = floor(alpha);
1261 	if (alpha < 0)
1262 	{
1263 		/* Using Aaamowitz & Stegun  9.1.2
1264 	 * this may not be quite optimal (CPU and accuracy wise) */
1265 		return (besselJ(x, -alpha, b) * cosPi(alpha) + ((alpha == na) ? 0 : besselY(x,
1266 			-alpha, b) * sinPi(alpha)));
1267 	}
1268 	nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */
1269 	alpha -= nb - 1;
1270 	besselJ(x, alpha, b[0 .. nb], ncalc);
1271 	if (ncalc != nb)
1272 	{
1273 		 /* error input */
1274 		if (ncalc < 0)
1275 			throw new Exception(format("besselJ(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
1276 			x, ncalc, nb, alpha));
1277 		else throw new Exception(format(
1278 			"besselJ(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
1279 	}
1280 	x = b[b.length - 1];
1281 	return x;
1282 }
1283 
1284 void besselJ(double x, double alpha, double[] b, ref ptrdiff_t ncalc)
1285 {
1286 	/*
1287  Calculates Bessel functions J_{n+alpha} (x)
1288  for non-negative argument x, and non-negative order n+alpha, n = 0,1,..,b.length-1.
1289 
1290   Explanation of variables in the calling sequence.
1291 
1292  X	 - Non-negative argument for which J's are to be calculated.
1293  ALPHA - Fractional part of order for which
1294 	 J's are to be calculated.  0 <= ALPHA < 1.
1295  NB	- Number of functions to be calculated, NB >= 1.
1296 	 The first function calculated is of order ALPHA, and the
1297 	 last is of order (NB - 1 + ALPHA).
1298  B	 - Output vector of length NB.  If RJBESL
1299 	 terminates normally (NCALC=NB), the vector B contains the
1300 	 functions J/ALPHA/(X) through J/NB-1+ALPHA/(X).
1301  NCALC - Output variable indicating possible errors.
1302 	 Before using the vector B, the user should check that
1303 	 NCALC=NB, i.e., all orders have been calculated to
1304 	 the desired accuracy.	See the following
1305 
1306 	 ****************************************************************
1307 
1308  Error return codes
1309 
1310 	In case of an error,  NCALC != NB, and not all J's are
1311 	calculated to the desired accuracy.
1312 
1313 	NCALC < 0:	An argument is out of range. For example,
1314 		NBES <= 0, ALPHA < 0 or > 1, or X is too large.
1315 		In this case, b[1] is set to zero, the remainder of the
1316 		B-vector is not calculated, and NCALC is set to
1317 		MIN(NB,0)-1 so that NCALC != NB.
1318 
1319 	NB > NCALC > 0: Not all requested function values could
1320 		be calculated accurately.  This usually occurs because NB is
1321 		much larger than ABS(X).	 In this case, b[N] is calculated
1322 		to the desired accuracy for N <= NCALC, but precision
1323 		is lost for NCALC < N <= NB.  If b[N] does not vanish
1324 		for N > NCALC (because it is too small to be represented),
1325 		and b[N]/b[NCALC] = 10^(-K), then only the first NSIG - K
1326 		significant figures of b[N] can be trusted.
1327 
1328 
1329   Acknowledgement
1330 
1331 	This program is based on a program written b David J. Sookne
1332 	(2) that computes values of the Bessel functions J or I of float
1333 	argument and long order.  Modifications include the restriction
1334 	of the computation to the J Bessel function of non-negative float
1335 	argument, the extension of the computation to arbtrary positive
1336 	order, and the elimination of most underflow.
1337 
1338   References:
1339 
1340 	Olver, F.W.J., and Sookne, D.J. (1972)
1341 	"A Note on Backward Recurrence Algorithms";
1342 	Math. Comp. 26, 941-947.
1343 
1344 	Sookne, D.J. (1973)
1345 	"Bessel Functions of Real Argument and Integer Order";
1346 	NBS Jour. of Res. B. 77B, 125-132.
1347 
1348   Latest modification: March 19, 1990
1349 
1350   Author: W. J. Cody
1351 	  Applied Mathematics Division
1352 	  Argonne National Laboratory
1353 	  Argonne, IL  60439
1354  *******************************************************************
1355  */
1356 	
1357 	/* ---------------------------------------------------------------------
1358   Mathematical constants
1359 
1360 	PI2	  = 2 / PI
1361 	TWOPI1 = first few significant digits of 2 * PI
1362 	TWOPI2 = (2*PI - TWOPI1) to working precision, i.e.,
1363 		TWOPI1 + TWOPI2 = 2 * PI to extra precision.
1364  --------------------------------------------------------------------- */
1365 	static immutable double pi2 = .636619772367581343075535;
1366 	static immutable double twopi1 = 6.28125;
1367 	static immutable double twopi2 = .001935307179586476925286767;
1368 	/*---------------------------------------------------------------------
1369  *  Factorial(N)
1370  *--------------------------------------------------------------------- */
1371 	static immutable double[25] fact = [1., 1., 2., 6., 24., 120., 720., 5040., 40320.,
1372 		362880., 3628800., 39916800., 479001600., 6227020800., 87178291200.,
1373 		1.307674368e12, 2.0922789888e13, 3.55687428096e14, 6.402373705728e15,
1374 		1.21645100408832e17, 2.43290200817664e18, 5.109094217170944e19,
1375 		1.12400072777760768e21,
1376 		2.585201673888497664e22, 6.2044840173323943936e23];
1377 	/* Local variables */
1378 	ptrdiff_t nend, intx, nbmx, i, j, k, l, m, n, nstart;
1379 	double nu, twonu, capp, capq, pold, vcos, test, vsin;
1380 	double p, s, t, z, alpem, halfx, aa, bb, cc, psave, plast;
1381 	double tover, t1, alp2em, em, en, xc, xk, xm, psavel, gnu, xin, sum;
1382 	/* Parameter adjustment */
1383 	//	--b;
1384 	nu = alpha;
1385 	twonu = nu + nu;
1386 	/*-------------------------------------------------------------------
1387 	  Check for out of range arguments.
1388 	  -------------------------------------------------------------------*/
1389 	if (b.length > 0 && x >= 0. && 0. <= nu && nu < 1.)
1390 	{
1391 		ncalc = b.length;
1392 		if (x > xlrg_BESS_IJ)
1393 		{
1394 			/* indeed, the limit is 0,
1395 		 * but the cutoff happens too early */
1396 			b[] = double.infinity;
1397 			return;
1398 		}
1399 		intx = cast(int)(x);
1400 		/* Initialize result array to zero. */
1401 		b[] = 0;
1402 		/*===================================================================
1403 	  Branch into  3 cases :
1404 	  1) use 2-term ascending series for small X
1405 	  2) use asymptotic form for large X when NB is not too large
1406 	  3) use recursion otherwise
1407 	  ===================================================================*/
1408 		if (x < rtnsig_BESS)
1409 		{
1410 			/* ---------------------------------------------------------------
1411 		 Two-term ascending series for small X.
1412 		 --------------------------------------------------------------- */
1413 			alpem = 1. + nu;
1414 			halfx = (x > enmten_BESS) ? .5 * x : 0;
1415 			aa = (nu != 0.) ? pow(halfx, nu) / (nu * gamma_cody(nu)) : 1.;
1416 			bb = (x + 1. > 1.) ? -halfx * halfx : 0;
1417 			b[0] = aa + aa * bb / alpem;
1418 			if (x != 0. && b[0] == 0.)
1419 				ncalc = 0;
1420 			if (b.length != 1)
1421 			{
1422 				if (x <= 0.)
1423 				{
1424 					for (n = 2; n <= b.length; ++n)
1425 						b[n - 1] = 0;
1426 				}
1427 				else
1428 				{
1429 					/* ----------------------------------------------
1430 				Calculate higher order functions.
1431 				---------------------------------------------- */
1432 					if (bb == 0.)
1433 						tover = (enmten_BESS + enmten_BESS) / x;
1434 					else tover = enmten_BESS / bb;
1435 					cc = halfx;
1436 					for (n = 2; n <= b.length; ++n)
1437 					{
1438 						aa /= alpem;
1439 						alpem += 1.;
1440 						aa *= cc;
1441 						if (aa <= tover * alpem)
1442 							aa = 0;
1443 						b[n - 1] = aa + aa * bb / alpem;
1444 						if (b[n - 1] == 0. && ncalc > n)
1445 							ncalc = n - 1;
1446 					}
1447 				}
1448 			}
1449 		}
1450 		else if (x > 25. && b.length <= intx + 1)
1451 		{
1452 			/* ------------------------------------------------------------
1453 			Asymptotic series for X > 25 (and not too large b.length)
1454 			------------------------------------------------------------ */
1455 			xc = sqrt(pi2 / x);
1456 			xin = 1 / (64 * x * x);
1457 			if (x >= 130.)
1458 				m = 4;
1459 			else if (x >= 35.)
1460 				m = 8;
1461 			else m = 11;
1462 			xm = 4. * cast(double) m;
1463 			/* ------------------------------------------------
1464 			Argument reduction for SIN and COS routines.
1465 			------------------------------------------------ */
1466 			t = trunc(x / (twopi1 + twopi2) + .5);
1467 			z = (x - t * twopi1) - t * twopi2 - (nu + .5) / pi2;
1468 			vsin = sin(z);
1469 			vcos = cos(z);
1470 			gnu = twonu;
1471 			for (i = 1; i <= 2; ++i)
1472 			{
1473 				s = (xm - 1. - gnu) * (xm - 1. + gnu) * xin * .5;
1474 				t = (gnu - (xm - 3.)) * (gnu + (xm - 3.));
1475 				t1 = (gnu - (xm + 1.)) * (gnu + (xm + 1.));
1476 				k = m + m;
1477 				capp = s * t / fact[k];
1478 				capq = s * t1 / fact[k + 1];
1479 				xk = xm;
1480 				for (; k >= 4; k -= 2)
1481 				{
1482 					 /* k + 2(j-2) == 2m */
1483 					xk -= 4.;
1484 					s = (xk - 1. - gnu) * (xk - 1. + gnu);
1485 					t1 = t;
1486 					t = (gnu - (xk - 3.)) * (gnu + (xk - 3.));
1487 					capp = (capp + 1. / fact[k - 2]) * s * t * xin;
1488 					capq = (capq + 1. / fact[k - 1]) * s * t1 * xin;
1489 				}
1490 				capp += 1.;
1491 				capq = (capq + 1.) * (gnu * gnu - 1.) * (.125 / x);
1492 				b[i - 1] = xc * (capp * vcos - capq * vsin);
1493 				if (b.length == 1)
1494 					return;
1495 				/* vsin <--> vcos */ t = vsin;
1496 				vsin = -vcos;
1497 				vcos = t;
1498 				gnu += 2.;
1499 			}
1500 			/* -----------------------------------------------
1501 			If  NB > 2, compute J(X,ORDER+I)	for I = 2, NB-1
1502 			----------------------------------------------- */
1503 			if (b.length > 2)
1504 				for (gnu = twonu + 2., j = 2; j <= b.length; j++, gnu += 2.)
1505 				b[j] = gnu * b[j - 1] / x - b[j - 2];
1506 		}
1507 		else
1508 		{
1509 			/* rtnsig_BESS <= x && ( x <= 25 || intx+1 < b.length ) :
1510 			--------------------------------------------------------
1511 			Use recurrence to generate results.
1512 			First initialize the calculation of P*S.
1513 			-------------------------------------------------------- */
1514 			nbmx = b.length - intx;
1515 			n = intx + 1;
1516 			en = cast(double)(n + n) + twonu;
1517 			plast = 1.;
1518 			p = en / x;
1519 			/* ---------------------------------------------------
1520 			Calculate general significance test.
1521 			--------------------------------------------------- */
1522 			test = ensig_BESS + ensig_BESS;
1523 			if (nbmx >= 3)
1524 			{
1525 				/* ------------------------------------------------------------
1526 			Calculate P*S until N = NB-1.  Check for possible overflow.
1527 			---------------------------------------------------------- */
1528 				tover = enten_BESS / ensig_BESS;
1529 				nstart = intx + 2;
1530 				nend = b.length - 1;
1531 				en = cast(double)(nstart + nstart) - 2. + twonu;
1532 				for (k = nstart; k <= nend; ++k)
1533 				{
1534 					n = k;
1535 					en += 2.;
1536 					pold = plast;
1537 					plast = p;
1538 					p = en * plast / x - pold;
1539 					if (p > tover)
1540 					{
1541 						/* -------------------------------------------
1542 				To avoid overflow, divide P*S b TOVER.
1543 				Calculate P*S until ABS(P) > 1.
1544 				-------------------------------------------*/
1545 						tover = enten_BESS;
1546 						p /= tover;
1547 						plast /= tover;
1548 						psave = p;
1549 						psavel = plast;
1550 						nstart = n + 1;
1551 						do
1552 						{
1553 							++n;
1554 							en += 2.;
1555 							pold = plast;
1556 							plast = p;
1557 							p = en * plast / x - pold;
1558 						}
1559 						while (p <= 1.);
1560 						bb = en / x;
1561 						/* -----------------------------------------------
1562 				Calculate backward test and find NCALC,
1563 				the highest N such that the test is passed.
1564 				----------------------------------------------- */
1565 						test = pold * plast * (.5 - .5 / (bb * bb));
1566 						test /= ensig_BESS;
1567 						p = plast * tover;
1568 						--n;
1569 						en -= 2.;
1570 						nend = min(b.length, n);
1571 						for (l = nstart; l <= nend; ++l)
1572 						{
1573 							pold = psavel;
1574 							psavel = psave;
1575 							psave = en * psavel / x - pold;
1576 							if (psave * psavel > test)
1577 							{
1578 								ncalc = l - 1;
1579 								goto L190;
1580 							}
1581 						}
1582 						ncalc = nend;
1583 						goto L190;
1584 					}
1585 				}
1586 				n = nend;
1587 				en = cast(double)(n + n) + twonu;
1588 				/* -----------------------------------------------------
1589 			Calculate special significance test for NBMX > 2.
1590 			-----------------------------------------------------*/
1591 				test = max(test, sqrt(plast * ensig_BESS) * sqrt(p + p));
1592 			}
1593 			/* ------------------------------------------------
1594 			Calculate P*S until significance test passes. */
1595 			do
1596 			{
1597 				++n;
1598 				en += 2.;
1599 				pold = plast;
1600 				plast = p;
1601 				p = en * plast / x - pold;
1602 			}
1603 			while (p < test);
1604 			L190: /*---------------------------------------------------------------
1605 		  Initialize the backward recursion and the normalization sum.
1606 		  --------------------------------------------------------------- */
1607 			++n;
1608 			en += 2.;
1609 			bb = 0;
1610 			aa = 1. / p;
1611 			m = n / 2;
1612 			em = cast(double) m;
1613 			m = (n << 1) - (m << 2); /* = 2 n - 4 (n/2)
1614 						= 0 for even, 2 for odd n */
1615 			if (m == 0)
1616 				sum = 0;
1617 			else
1618 			{
1619 				alpem = em - 1. + nu;
1620 				alp2em = em + em + nu;
1621 				sum = aa * alpem * alp2em / em;
1622 			}
1623 			nend = n - b.length;
1624 			/* if (nend > 0) */
1625 			/* --------------------------------------------------------
1626 			Recur backward via difference equation, calculating
1627 			(but not storing) b[N], until N = NB.
1628 			-------------------------------------------------------- */
1629 			for (l = 1; l <= nend; ++l)
1630 			{
1631 				--n;
1632 				en -= 2.;
1633 				cc = bb;
1634 				bb = aa;
1635 				aa = en * bb / x - cc;
1636 				m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
1637 				if (m != 0)
1638 				{
1639 					em -= 1.;
1640 					alp2em = em + em + nu;
1641 					if (n == 1)
1642 						break;
1643 					alpem = em - 1. + nu;
1644 					if (alpem == 0.)
1645 						alpem = 1.;
1646 					sum = (sum + aa * alp2em) * alpem / em;
1647 				}
1648 			}
1649 			/*--------------------------------------------------
1650 		  Store b[NB].
1651 		  --------------------------------------------------*/
1652 			b[n - 1] = aa;
1653 			if (nend >= 0)
1654 			{
1655 				if (b.length <= 1)
1656 				{
1657 					if (nu + 1. == 1.)
1658 						alp2em = 1.;
1659 					else alp2em = nu;
1660 					sum += b[0] * alp2em;
1661 					goto L250;
1662 				}
1663 				else
1664 				{
1665 					 /*-- b.length >= 2 : ---------------------------
1666 			Calculate and store b[NB-1].
1667 			----------------------------------------*/
1668 					--n;
1669 					en -= 2.;
1670 					b[n - 1] = en * aa / x - bb;
1671 					if (n == 1)
1672 						goto L240;
1673 					m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
1674 					if (m != 0)
1675 					{
1676 						em -= 1.;
1677 						alp2em = em + em + nu;
1678 						alpem = em - 1. + nu;
1679 						if (alpem == 0.)
1680 							alpem = 1.;
1681 						sum = (sum + b[n - 1] * alp2em) * alpem / em;
1682 					}
1683 				}
1684 			}
1685 			/* if (n - 2 != 0) */
1686 			/* --------------------------------------------------------
1687 			Calculate via difference equation and store b[N],
1688 			until N = 2.
1689 			-------------------------------------------------------- */
1690 			for (n = n - 1; n >= 2; n--)
1691 			{
1692 				en -= 2.;
1693 				b[n - 1] = en * b[n + 0] / x - b[n + 1];
1694 				m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
1695 				if (m != 0)
1696 				{
1697 					em -= 1.;
1698 					alp2em = em + em + nu;
1699 					alpem = em - 1. + nu;
1700 					if (alpem == 0.)
1701 						alpem = 1.;
1702 					sum = (sum + b[n - 1] * alp2em) * alpem / em;
1703 				}
1704 			}
1705 			/* ---------------------------------------
1706 			Calculate b[1].
1707 			-----------------------------------------*/
1708 			b[0] = 2. * (nu + 1.) * b[1] / x - b[2];
1709 			L240: em -= 1.;
1710 			alp2em = em + em + nu;
1711 			if (alp2em == 0.)
1712 				alp2em = 1.;
1713 			sum += b[0] * alp2em;
1714 			L250: /* ---------------------------------------------------
1715 			Normalize.  Divide all b[N] b sum.
1716 			---------------------------------------------------*/
1717 			
1718 			/*		if (nu + 1. != 1.) poor test */
1719 			if (fabs(nu) > 1e-15)
1720 				sum *= (gamma_cody(nu) * pow(.5 * x, -nu));
1721 			aa = enmten_BESS;
1722 			if (sum > 1.)
1723 				aa *= sum;
1724 			for (n = 1; n <= b.length; ++n)
1725 			{
1726 				if (fabs(b[n - 1]) < aa)
1727 					b[n - 1] = 0;
1728 				else b[n - 1] /= sum;
1729 			}
1730 		}
1731 	}
1732 	else
1733 	{
1734 		/* Error return -- X, NB, or ALPHA is out of range : */
1735 		b[0] = 0;
1736 		ncalc = min(b.length, 0) - 1;
1737 	}
1738 }
1739 
1740 ///
1741 double besselY(double x, double alpha)
1742 {
1743 	ptrdiff_t nb, ncalc;
1744 	double na;
1745 	double* b;
1746 	/* NaNs propagated correctly */
1747 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
1748 		return x + alpha;
1749 	if (x < 0)
1750 		return double.nan;
1751 	na = floor(alpha);
1752 	if (alpha < 0)
1753 	{
1754 		/* Using Aaamowitz & Stegun  9.1.2
1755 	 * this may not be quite optimal (CPU and accuracy wise) */
1756 		return (besselY(x, -alpha) * cosPi(alpha) - ((alpha == na) ? 0 : besselJ(x,
1757 			-alpha) * sinPi(alpha)));
1758 	}
1759 	nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */
1760 	alpha -= nb - 1;
1761 	b = cast(double* ) calloc(nb, double.sizeof).enforce("besselY allocation error");
1762 	scope(exit) b.free;
1763 	besselY(x, alpha, b[0 .. nb], ncalc);
1764 	if (ncalc != nb)
1765 	{
1766 		 /* error input */
1767 		if (ncalc == -1)
1768 			return double.infinity;
1769 		else if (ncalc < -1)
1770 			throw new Exception(format("besselY(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
1771 			x, ncalc, nb, alpha));
1772 		else  /* ncalc >= 0 */
1773 		throw new Exception(format(
1774 			"besselY(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
1775 	}
1776 	x = b[nb - 1];
1777 	return x;
1778 }
1779 
1780 ///
1781 unittest
1782 {
1783 	assert(std.math.approxEqual(besselY(2, 1),
1784 		-0.1070324315409375468883707722774));
1785 	assert(std.math.approxEqual(besselY(20, 1),
1786 		-0.165511614362521295863976023243));
1787 	assert(std.math.approxEqual(besselY(1, 1.2), -0.901215));
1788 }
1789 
1790 /* modified version of besselY that accepts a work array instead of
1791 	allocating one. */
1792 double besselY(double x, double alpha, double[] b)
1793 {
1794 	ptrdiff_t nb, ncalc;
1795 	double na;
1796 	/* NaNs propagated correctly */
1797 	if (std.math.isNaN(x) || std.math.isNaN(alpha))
1798 		return x + alpha;
1799 	if (x < 0)
1800 		return double.nan;
1801 	na = floor(alpha);
1802 	if (alpha < 0)
1803 	{
1804 		/* Using Aaamowitz & Stegun  9.1.2
1805 	 * this may not be quite optimal (CPU and accuracy wise) */
1806 		return (besselY(x, -alpha, b) * cosPi(alpha) - ((alpha == na) ? 0 : besselJ(x,
1807 			-alpha, b) * sinPi(alpha)));
1808 	}
1809 	nb = 1 + cast(int) na; /* b.length-1 <= alpha < b.length */
1810 	alpha -= nb - 1;
1811 	besselY(x, alpha, b[0 .. nb], ncalc);
1812 	if (ncalc != nb)
1813 	{
1814 		 /* error input */
1815 		if (ncalc == -1)
1816 			return double.infinity;
1817 		else if (ncalc < -1)
1818 			throw new Exception(format("besselY(%g): ncalc (=%d) != b.length (=%d); alpha=%g. Arg. out of range?\n",
1819 			x, ncalc, nb, alpha));
1820 		else  /* ncalc >= 0 */
1821 		throw new Exception(format(
1822 			"besselY(%g,nu=%g): precision lost in result\n", x, alpha + (nb - 1)));
1823 	}
1824 	x = b[nb - 1];
1825 	return x;
1826 }
1827 
1828 void besselY(double x, double alpha, double[] b, ref ptrdiff_t ncalc)
1829 {
1830 	/* ----------------------------------------------------------------------
1831 
1832  This routine calculates Bessel functions Y_(N+ALPHA) (X)
1833 v for non-negative argument X, and non-negative order N+ALPHA.
1834 
1835 
1836  Explanation of variables in the calling sequence
1837 
1838  X	 - Non-negative argument for which
1839 	 Y's are to be calculated.
1840  ALPHA - Fractional part of order for which
1841 	 Y's are to be calculated.  0 <= ALPHA < 1.0.
1842  NB	- Number of functions to be calculated, NB > 0.
1843 	 The first function calculated is of order ALPHA, and the
1844 	 last is of order (NB - 1 + ALPHA).
1845  BY	- Output vector of length NB.	If the
1846 	 routine terminates normally (NCALC=NB), the vector BY
1847 	 contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X),
1848 	 If (0 < NCALC < NB), BY(I) contains correct function
1849 	 values for I <= NCALC, and contains the ratios
1850 	 Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array.
1851  NCALC - Output variable indicating possible errors.
1852 	 Before using the vector BY, the user should check that
1853 	 NCALC=NB, i.e., all orders have been calculated to
1854 	 the desired accuracy.	See error returns below.
1855 
1856 
1857  *******************************************************************
1858 
1859  Error returns
1860 
1861   In case of an error, NCALC != NB, and not all Y's are
1862   calculated to the desired accuracy.
1863 
1864   NCALC < -1:  An argument is out of range. For example,
1865 	NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >=
1866 	XMAX.  In this case, BY[0] = 0.0, the remainder of the
1867 	BY-vector is not calculated, and NCALC is set to
1868 	MIN0(NB,0)-2  so that NCALC != NB.
1869   NCALC = -1:  Y(ALPHA,X) >= XINF.  The requested function
1870 	values are set to 0.0.
1871   1 < NCALC < NB: Not all requested function values could
1872 	be calculated accurately.  BY(I) contains correct function
1873 	values for I <= NCALC, and and the remaining NB-NCALC
1874 	array elements contain 0.0.
1875 
1876 
1877  Intrinsic functions required are:
1878 
1879 	 DBLE, EXP, INT, MAX, MIN, REAL, SQRT
1880 
1881 
1882  Acknowledgement
1883 
1884 	This program draws heavily on Temme's Algol program for Y(a,x)
1885 	and Y(a+1,x) and on Campbell's programs for Y_nu(x).	Temme's
1886 	scheme is used for  x < THRESH, and Campbell's scheme is used
1887 	in the asymptotic region.  Segments of code from both sources
1888 	have been translated into Fortran 77, merged, and heavily modified.
1889 	Modifications include parameterization of machine dependencies,
1890 	use of a new approximation for ln(gamma(x)), and built-in
1891 	protection against over/underflow.
1892 
1893  References: "Bessel functions J_nu(x) and Y_nu(x) of float
1894 		  order and float argument," Campbell, J. B.,
1895 		  Comp. Phy. Comm. 18, 1979, pp. 133-142.
1896 
1897 		 "On the numerical evaluation of the ordinary
1898 		  Bessel function of the second kind," Temme,
1899 		  N. M., J. Comput. Phys. 21, 1976, pp. 343-350.
1900 
1901   Latest modification: March 19, 1990
1902 
1903   Modified b: W. J. Cody
1904 			Applied Mathematics Division
1905 			Argonne National Laboratory
1906 			Argonne, IL  60439
1907  ----------------------------------------------------------------------*/
1908 	
1909 	/* ----------------------------------------------------------------------
1910   Mathematical constants
1911 	FIVPI = 5*PI
1912 	PIM5 = 5*PI - 15
1913  ----------------------------------------------------------------------*/
1914 	static immutable double fivpi = 15.707963267948966192;
1915 	static immutable double pim5 = .70796326794896619231;
1916 	/*----------------------------------------------------------------------
1917 	  Coefficients for Chebshev polynomial expansion of
1918 	  1/gamma(1-x), abs(x) <= .5
1919 	  ----------------------------------------------------------------------*/
1920 	static immutable double[21] ch = [-6.7735241822398840964e-24,
1921 		-6.1455180116049879894e-23, 2.9017595056104745456e-21,
1922 		1.3639417919073099464e-19, 2.3826220476859635824e-18,
1923 		-9.0642907957550702534e-18,
1924 		-1.4943667065169001769e-15,
1925 		-3.3919078305362211264e-14,
1926 		-1.7023776642512729175e-13, 9.1609750938768647911e-12,
1927 		2.4230957900482704055e-10, 1.7451364971382984243e-9,
1928 		-3.3126119768180852711e-8, -8.6592079961391259661e-7,
1929 		-4.9717367041957398581e-6,
1930 		7.6309597585908126618e-5,
1931 		.0012719271366545622927,
1932 		.0017063050710955562222,
1933 		-.07685284084478667369,
1934 		-.28387654227602353814, .92187029365045265648];
1935 	/* Local variables */
1936 	ptrdiff_t i, k, na;
1937 	double alfa, div, ddiv, even, gamma, term, cosmu, sinmu, a, c, d, e, f, g, h, p, q, r,
1938 		s, d1, d2, q0, pa, pa1, qa, qa1, en, en1, nu, ex, ya, ya1, twobx, den, odd, aye,
1939 		dmu, x2, xna;
1940 	en1 = ya = ya1 = 0; /* -Wall */
1941 	ex = x;
1942 	nu = alpha;
1943 	if (b.length > 0 && 0. <= nu && nu < 1.)
1944 	{
1945 		if (ex < double.min_normal || ex > xlrg_BESS_Y)
1946 		{
1947 			/* Warning is not really appropriate, give
1948 		 * proper limit:
1949 		 * ML_ERROR(ME_RANGE, "besselY"); */
1950 			ncalc = b.length;
1951 			if (ex > xlrg_BESS_Y)
1952 				b[0] = 0; /*was double.infinity */
1953 			else if (ex < double.min_normal)
1954 				b[0] = -double.infinity;
1955 			for (i = 0; i < b.length; i++)
1956 				b[i] = b[0];
1957 			return;
1958 		}
1959 		xna = trunc(nu + .5);
1960 		na = cast(int) xna;
1961 		if (na == 1)
1962 		{
1963 			 /* <==>  .5 <= alpha < 1	 <==>  -5. <= nu < 0 */
1964 			nu -= xna;
1965 		}
1966 		if (nu == -.5)
1967 		{
1968 			p = M_SQRT_2dPI / sqrt(ex);
1969 			ya = p * sin(ex);
1970 			ya1 = -p * cos(ex);
1971 		}
1972 		else if (ex < 3.)
1973 		{
1974 			/* -------------------------------------------------------------
1975 			Use Temme's scheme for small X
1976 			------------------------------------------------------------- */
1977 			a = ex * .5;
1978 			d = -log(a);
1979 			f = nu * d;
1980 			e = pow(a, -nu);
1981 			if (fabs(nu) < M_eps_sinc)
1982 				c = std.math.M_1_PI;
1983 			else c = nu / sinPi(nu);
1984 			/* ------------------------------------------------------------
1985 			Computation of sinh(f)/f
1986 			------------------------------------------------------------ */
1987 			if (fabs(f) < 1.)
1988 			{
1989 				x2 = f * f;
1990 				en = 19.;
1991 				s = 1.;
1992 				for (i = 1; i <= 9; ++i)
1993 				{
1994 					s = s * x2 / en / (en - 1.) + 1.;
1995 					en -= 2.;
1996 				}
1997 			}
1998 			else
1999 			{
2000 				s = (e - 1. / e) * .5 / f;
2001 			}
2002 			/* --------------------------------------------------------
2003 			Computation of 1/gamma(1-a) using Chebshev polynomials */
2004 			x2 = nu * nu * 8.;
2005 			aye = ch[0];
2006 			even = 0;
2007 			alfa = ch[1];
2008 			odd = 0;
2009 			for (i = 3; i <= 19; i += 2)
2010 			{
2011 				even = -(aye + aye + even);
2012 				aye = -even * x2 - aye + ch[i - 1];
2013 				odd = -(alfa + alfa + odd);
2014 				alfa = -odd * x2 - alfa + ch[i];
2015 			}
2016 			even = (even * .5 + aye) * x2 - aye + ch[20];
2017 			odd = (odd + alfa) * 2.;
2018 			gamma = odd * nu + even;
2019 			/* End of computation of 1/gamma(1-a)
2020 			----------------------------------------------------------- */
2021 			g = e * gamma;
2022 			e = (e + 1. / e) * .5;
2023 			f = 2. * c * (odd * e + even * s * d);
2024 			e = nu * nu;
2025 			p = g * c;
2026 			q = std.math.M_1_PI / g;
2027 			c = nu * std.math.PI_2;
2028 			if (fabs(c) < M_eps_sinc)
2029 				r = 1.;
2030 			else r = sinPi(nu / 2) / c;
2031 			r = cast(double) std.math.PI * c * r * r;
2032 			c = 1.;
2033 			d = -(a * a);
2034 			h = 0;
2035 			ya = f + r * q;
2036 			ya1 = p;
2037 			en = 1.;
2038 			while (fabs(g / (1. + fabs(ya))) + fabs(h / (1. + fabs(ya1))) > double.epsilon)
2039 			{
2040 				f = (f * en + p + q) / (en * en - e);
2041 				c *= (d / en);
2042 				p /= en - nu;
2043 				q /= en + nu;
2044 				g = c * (f + r * q);
2045 				h = c * p - en * g;
2046 				ya += g;
2047 				ya1 += h;
2048 				en += 1.;
2049 			}
2050 			ya = -ya;
2051 			ya1 = -ya1 / a;
2052 		}
2053 		else if (ex < thresh_BESS_Y)
2054 		{
2055 			/* --------------------------------------------------------------
2056 			Use Temme's scheme for moderate X :  3 <= x < 16
2057 			-------------------------------------------------------------- */
2058 			c = (.5 - nu) * (.5 + nu);
2059 			a = ex + ex;
2060 			e = ex * std.math.M_1_PI * cosPi(nu) / double.epsilon;
2061 			e *= e;
2062 			p = 1.;
2063 			q = -ex;
2064 			r = 1. + ex * ex;
2065 			s = r;
2066 			en = 2.;
2067 			while (r * en * en < e)
2068 			{
2069 				en1 = en + 1.;
2070 				d = (en - 1. + c / en) / s;
2071 				p = (en + en - p * d) / en1;
2072 				q = (-a + q * d) / en1;
2073 				s = p * p + q * q;
2074 				r *= s;
2075 				en = en1;
2076 			}
2077 			f = p / s;
2078 			p = f;
2079 			g = -q / s;
2080 			q = g;
2081 			L220: en -= 1.;
2082 			if (en > 0.)
2083 			{
2084 				r = en1 * (2. - p) - 2.;
2085 				s = a + en1 * q;
2086 				d = (en - 1. + c / en) / (r * r + s * s);
2087 				p = d * r;
2088 				q = d * s;
2089 				e = f + 1.;
2090 				f = p * e - g * q;
2091 				g = q * e + p * g;
2092 				en1 = en;
2093 				goto L220;
2094 			}
2095 			f = 1. + f;
2096 			d = f * f + g * g;
2097 			pa = f / d;
2098 			qa = -g / d;
2099 			d = nu + .5 - p;
2100 			q += ex;
2101 			pa1 = (pa * q - qa * d) / ex;
2102 			qa1 = (qa * q + pa * d) / ex;
2103 			a = ex - std.math.PI_2 * (nu + .5);
2104 			c = cos(a);
2105 			s = sin(a);
2106 			d = M_SQRT_2dPI / sqrt(ex);
2107 			ya = d * (pa * s + qa * c);
2108 			ya1 = d * (qa1 * s - pa1 * c);
2109 		}
2110 		else
2111 		{
2112 			 /* x > thresh_BESS_Y */
2113 			/* ----------------------------------------------------------
2114 			Use Campbell's asymptotic scheme.
2115 			---------------------------------------------------------- */
2116 			na = 0;
2117 			d1 = trunc(ex / fivpi);
2118 			i = cast(int) d1;
2119 			dmu = ex - 15. * d1 - d1 * pim5 - (alpha + .5) * std.math.PI_2;
2120 			if (i - (i / 2 << 1) == 0)
2121 			{
2122 				cosmu = cos(dmu);
2123 				sinmu = sin(dmu);
2124 			}
2125 			else
2126 			{
2127 				cosmu = -cos(dmu);
2128 				sinmu = -sin(dmu);
2129 			}
2130 			ddiv = 8. * ex;
2131 			dmu = alpha;
2132 			den = sqrt(ex);
2133 			for (k = 1; k <= 2; ++k)
2134 			{
2135 				p = cosmu;
2136 				cosmu = sinmu;
2137 				sinmu = -p;
2138 				d1 = (2. * dmu - 1.) * (2. * dmu + 1.);
2139 				d2 = 0;
2140 				div = ddiv;
2141 				p = 0;
2142 				q = 0;
2143 				q0 = d1 / div;
2144 				term = q0;
2145 				for (i = 2; i <= 20; ++i)
2146 				{
2147 					d2 += 8.;
2148 					d1 -= d2;
2149 					div += ddiv;
2150 					term = -term * d1 / div;
2151 					p += term;
2152 					d2 += 8.;
2153 					d1 -= d2;
2154 					div += ddiv;
2155 					term *= (d1 / div);
2156 					q += term;
2157 					if (fabs(term) <= double.epsilon)
2158 					{
2159 						break;
2160 					}
2161 				}
2162 				p += 1.;
2163 				q += q0;
2164 				if (k == 1)
2165 					ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
2166 				else ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
2167 				dmu += 1.;
2168 			}
2169 		}
2170 		if (na == 1)
2171 		{
2172 			h = 2. * (nu + 1.) / ex;
2173 			if (h > 1.)
2174 			{
2175 				if (fabs(ya1) > double.max / h)
2176 				{
2177 					h = 0;
2178 					ya = 0;
2179 				}
2180 			}
2181 			h = h * ya1 - ya;
2182 			ya = ya1;
2183 			ya1 = h;
2184 		}
2185 		/* ---------------------------------------------------------------
2186 		Now have first one or two Y's
2187 		--------------------------------------------------------------- */
2188 		b[0] = ya;
2189 		ncalc = 1;
2190 		if (b.length > 1)
2191 		{
2192 			b[1] = ya1;
2193 			if (ya1 != 0.)
2194 			{
2195 				aye = 1. + alpha;
2196 				twobx = 2. / ex;
2197 				ncalc = 2;
2198 				for (i = 2; i < b.length; ++i)
2199 				{
2200 					if (twobx < 1.)
2201 					{
2202 						if (fabs(b[i - 1]) * twobx >= double.max / aye)
2203 							goto L450;
2204 					}
2205 					else
2206 					{
2207 						if (fabs(b[i - 1]) >= double.max / aye / twobx)
2208 							goto L450;
2209 					}
2210 					b[i] = twobx * aye * b[i - 1] - b[i - 2];
2211 					aye += 1.;
2212 					ncalc++;
2213 				}
2214 			}
2215 		}
2216 		L450: for (i = ncalc; i < b.length; ++i)
2217 			b[i] = -double.infinity; /* was 0 */
2218 	}
2219 	else
2220 	{
2221 		b[0] = 0;
2222 		ncalc = min(b.length, 0) - 1;
2223 	}
2224 }
2225 
2226 private:
2227 double gamma_cody(double x)
2228 {
2229 	/* ----------------------------------------------------------------------
2230 
2231 	This routine calculates the GAMMA function for a float argument X.
2232 	Computation is based on an algorithm outlined in reference [1].
2233 	The program uses rational functions that approximate the GAMMA
2234 	function to at least 20 significant decimal digits.	Coefficients
2235 	for the approximation over the interval (1,2) are unpublished.
2236 	Those for the approximation for X >= 12 are from reference [2].
2237 	The accuracy achieved depends on the arithmetic system, the
2238 	compiler, the intrinsic functions, and proper selection of the
2239 	machine-dependent constants.
2240 
2241 	*******************************************************************
2242 
2243 	Error returns
2244 
2245 	The program returns the value XINF for singularities or
2246 	when overflow would occur.	 The computation is believed
2247 	to be free of underflow and overflow.
2248 
2249 	Intrinsic functions required are:
2250 
2251 	INT, DBLE, EXP, LOG, REAL, SIN
2252 
2253 
2254 	References:
2255 	[1]  "An Overview of Software Development for Special Functions",
2256 	W. J. Cody, Lecture Notes in Mathematics, 506,
2257 	Numerical Analysis Dundee, 1975, G. A. Watson (ed.),
2258 	Springer Verlag, Berlin, 1976.
2259 
2260 	[2]  Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968.
2261 
2262 	Latest modification: October 12, 1989
2263 
2264 	Authors: W. J. Cody and L. Stoltz
2265 	Applied Mathematics Division
2266 	Argonne National Laboratory
2267 	Argonne, IL 60439
2268 	----------------------------------------------------------------------*/
2269 	
2270 	/* ----------------------------------------------------------------------
2271 	Mathematical constants
2272 	----------------------------------------------------------------------*/
2273 	enum double sqrtpi = .9189385332046727417803297; /* == ??? */
2274 	
2275 	/* *******************************************************************
2276 
2277 	Explanation of machine-dependent constants
2278 
2279 	beta	- radix for the floating-point representation
2280 	maxexp - the smallest positive power of beta that overflows
2281 	XBIG	- the largest argument for which GAMMA(X) is representable
2282 	in the machine, i.e., the solution to the equation
2283 	GAMMA(XBIG) = beta**maxexp
2284 	XINF	- the largest machine representable floating-point number;
2285 	approximately beta**maxexp
2286 	EPS	- the smallest positive floating-point number such that  1.0+EPS > 1.0
2287 	XMININ - the smallest positive floating-point number such that
2288 	1/XMININ is machine representable
2289 
2290 	Approximate values for some important machines are:
2291 
2292 	beta		  maxexp		 XBIG
2293 
2294 	CRAY-1		(S.P.)		  2		8191		966.961
2295 	Cyber 180/855
2296 	under NOS	(S.P.)		  2		1070		177.803
2297 	IEEE (IBM/XT,
2298 	SUN, etc.)	(S.P.)		  2		 128		35.040
2299 	IEEE (IBM/XT,
2300 	SUN, etc.)	(D.P.)		  2		1024		171.624
2301 	IBM 3033	(D.P.)		 16		  63		57.574
2302 	VAX D-Format	(D.P.)		  2		 127		34.844
2303 	VAX G-Format	(D.P.)		  2		1023		171.489
2304 
2305 	XINF	 EPS		XMININ
2306 
2307 	CRAY-1		(S.P.)	 5.45E+2465	7.11E-15	  1.84E-2466
2308 	Cyber 180/855
2309 	under NOS	(S.P.)	 1.26E+322	3.55E-15	  3.14E-294
2310 	IEEE (IBM/XT,
2311 	SUN, etc.)	(S.P.)	 3.40E+38	 1.19E-7	  1.18E-38
2312 	IEEE (IBM/XT,
2313 	SUN, etc.)	(D.P.)	 1.79D+308	2.22D-16	  2.23D-308
2314 	IBM 3033	(D.P.)	 7.23D+75	 2.22D-16	  1.39D-76
2315 	VAX D-Format	(D.P.)	 1.70D+38	 1.39D-17	  5.88D-39
2316 	VAX G-Format	(D.P.)	 8.98D+307	1.11D-16	  1.12D-308
2317 
2318 	*******************************************************************
2319 
2320 	----------------------------------------------------------------------
2321 	Machine dependent parameters
2322 	----------------------------------------------------------------------
2323 	*/
2324 	enum double xbg = 171.624;
2325 	/* ML_POSINF ==	const double xinf = 1.79e308;*/
2326 	/* DBL_EPSILON = const double eps = 2.22e-16;*/
2327 	/* DBL_MIN ==	const double xminin = 2.23e-308;*/
2328 	
2329 	/*----------------------------------------------------------------------
2330 	  Numerator and denominator coefficients for rational minimax
2331 	  approximation over (1,2).
2332 	  ----------------------------------------------------------------------*/
2333 	static immutable double[8] p = [-1.71618513886549492533811,
2334 		24.7656508055759199108314, -379.804256470945635097577,
2335 		629.331155312818442661052, 866.966202790413211295064,
2336 		-31451.2729688483675254357,
2337 		-36144.4134186911729807069,
2338 		66456.1438202405440627855];
2339 	static immutable double[8] q = [-30.8402300119738975254353,
2340 		315.350626979604161529144, -1015.15636749021914166146,
2341 		-3107.77167157231109440444, 22538.1184209801510330112,
2342 		4755.84627752788110767815,
2343 		-134659.959864969306392456,
2344 		-115132.259675553483497211];
2345 	/*----------------------------------------------------------------------
2346 	  Coefficients for minimax approximation over (12, INF).
2347 	  ----------------------------------------------------------------------*/
2348 	static immutable double[7] c = [-.001910444077728, 8.4171387781295e-4,
2349 		-5.952379913043012e-4, 7.93650793500350248e-4, -.002777777777777681622553,
2350 		.08333333333333333331554247, .0057083835261];
2351 	/* Local variables */
2352 	int i, n;
2353 	int parity; /*logical*/
2354 	double fact, xden, xnum, y, z, yi, res, sum, ysq;
2355 	parity = (0);
2356 	fact = 1.;
2357 	n = 0;
2358 	y = x;
2359 	if (y <= 0.)
2360 	{
2361 		/* -------------------------------------------------------------
2362 		Argument is negative
2363 		------------------------------------------------------------- */
2364 		y = -x;
2365 		yi = trunc(y);
2366 		res = y - yi;
2367 		if (res != 0.)
2368 		{
2369 			if (yi != trunc(yi * .5) * 2.)
2370 				parity = (1);
2371 			fact = -std.math.PI / sinPi(res);
2372 			y += 1.;
2373 		}
2374 		else
2375 		{
2376 			return (double.infinity);
2377 		}
2378 	}
2379 	/* -----------------------------------------------------------------
2380 		Argument is positive
2381 		-----------------------------------------------------------------*/
2382 	if (y < double.epsilon)
2383 	{
2384 		/* --------------------------------------------------------------
2385 		Argument < EPS
2386 		-------------------------------------------------------------- */
2387 		if (y >= double.min_normal)
2388 		{
2389 			res = 1. / y;
2390 		}
2391 		else
2392 		{
2393 			return (double.infinity);
2394 		}
2395 	}
2396 	else if (y < 12.)
2397 	{
2398 		yi = y;
2399 		if (y < 1.)
2400 		{
2401 			/* ---------------------------------------------------------
2402 			EPS < argument < 1
2403 			--------------------------------------------------------- */
2404 			z = y;
2405 			y += 1.;
2406 		}
2407 		else
2408 		{
2409 			/* -----------------------------------------------------------
2410 			1 <= argument < 12, reduce argument if necessary
2411 			----------------------------------------------------------- */
2412 			n = cast(int) y - 1;
2413 			y -= cast(double) n;
2414 			z = y - 1.;
2415 		}
2416 		/* ---------------------------------------------------------
2417 		Evaluate approximation for 1. < argument < 2.
2418 		---------------------------------------------------------*/
2419 		xnum = 0;
2420 		xden = 1.;
2421 		for (i = 0; i < 8; ++i)
2422 		{
2423 			xnum = (xnum + p[i]) * z;
2424 			xden = xden * z + q[i];
2425 		}
2426 		res = xnum / xden + 1.;
2427 		if (yi < y)
2428 		{
2429 			/* --------------------------------------------------------
2430 			Adjust result for case  0. < argument < 1.
2431 			-------------------------------------------------------- */
2432 			res /= yi;
2433 		}
2434 		else if (yi > y)
2435 		{
2436 			/* ----------------------------------------------------------
2437 			Adjust result for case  2. < argument < 12.
2438 			---------------------------------------------------------- */
2439 			for (i = 0; i < n; ++i)
2440 			{
2441 				res *= y;
2442 				y += 1.;
2443 			}
2444 		}
2445 	}
2446 	else
2447 	{
2448 		/* -------------------------------------------------------------
2449 		Evaluate for argument >= 12.,
2450 		------------------------------------------------------------- */
2451 		if (y <= xbg)
2452 		{
2453 			ysq = y * y;
2454 			sum = c[6];
2455 			for (i = 0; i < 6; ++i)
2456 			{
2457 				sum = sum / ysq + c[i];
2458 			}
2459 			sum = sum / y - y + sqrtpi;
2460 			sum += (y - .5) * log(y);
2461 			res = exp(sum);
2462 		}
2463 		else
2464 		{
2465 			return double.infinity;
2466 		}
2467 	}
2468 	/* ----------------------------------------------------------------------
2469 		Final adjustments and return
2470 		----------------------------------------------------------------------*/
2471 	if (parity)
2472 		res = -res;
2473 	if (fact != 1.)
2474 		res = fact / res;
2475 	return res;
2476 }
2477 
2478 /// sin of pix:
2479 T sinPi(T)(T x) if (isFloatingPoint!T)
2480 {
2481 	if (x < 0)
2482 		return -sinPi(-x);
2483 	// sin of pix:
2484 	bool invert;
2485 	if (x < cast(T)0.5)
2486 		return sin(T(std.math.PI) * x);
2487 	if (x < 1)
2488 	{
2489 		invert = true;
2490 		x = -x;
2491 	}
2492 	else invert = false;
2493 	T rem = floor(x);
2494 	if (cast(int)rem & 1)
2495 		invert = !invert;
2496 	rem = x - rem;
2497 	if (rem > cast(T)0.5)
2498 		rem = 1 - rem;
2499 	if (rem == cast(T)0.5)
2500 		return invert ? -1 : 1;
2501 	rem = sin(T(std.math.PI) * rem);
2502 	return invert ? T(-rem) : rem;
2503 }
2504 
2505 unittest
2506 {
2507 	import std.math : approxEqual;
2508 
2509 	assert(sinPi(0.0) == 0);
2510 	assert(sinPi(1.0) == 0);
2511 	assert(sinPi(+1.0 / 6).approxEqual(+0.5));
2512 	assert(sinPi(-1.0 / 6).approxEqual(-0.5));
2513 	assert(sinPi(+5.0 / 6).approxEqual(+0.5));
2514 	assert(sinPi(-5.0 / 6).approxEqual(-0.5));
2515 	assert(sinPi(+2.0) == 0);
2516 	assert(sinPi(-2.0) == 0);
2517 	assert(sinPi(+0.5) == +1);
2518 	assert(sinPi(-0.5) == -1);
2519 }
2520 
2521 /// cos of pix:
2522 T cosPi(T)(T x) if (isFloatingPoint!T)
2523 {
2524 	// cos of pix:
2525 	bool invert = false;
2526 	if (fabs(x) < 0.5)
2527 		return cos(T(std.math.PI) * x);
2528 	if (x < cast(T)1)
2529 	{
2530 		x = -x;
2531 	}
2532 	T rem = floor(x);
2533 	if (cast(int)rem & 1)
2534 		invert = !invert;
2535 	rem = x - rem;
2536 	if (rem > cast(T)0.5)
2537 	{
2538 		rem = 1 - rem;
2539 		invert = !invert;
2540 	}
2541 	if (rem == cast(T)0.5)
2542 		return 0;
2543 	rem = cos(T(std.math.PI) * rem);
2544 	return invert ? T(-rem) : rem;
2545 }
2546 
2547 unittest
2548 {
2549 	import std.math : approxEqual;
2550 
2551 	assert(cosPi(0.0) == 1);
2552 	assert(cosPi(1.0) == -1);
2553 	assert(cosPi(+1.0 / 3).approxEqual(+0.5));
2554 	assert(cosPi(-1.0 / 3).approxEqual(+0.5));
2555 	assert(cosPi(+2.0 / 3).approxEqual(-0.5));
2556 	assert(cosPi(-2.0 / 3).approxEqual(-0.5));
2557 	assert(cosPi(2.0) == 1);
2558 	assert(cosPi(0.5) == 0);
2559 }