Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.1.0/ALGLIB-3.1.0/integration.cs @ 6345

Visit:

Last change on this file since 6345 was 4977, checked in by gkronber, 14 years ago
Added new alglib classes (version 3.1.0) #875
File size: 131.3 KB

Rev	Line
[4977]	1	/*************************************************************************
	2	Copyright (c) Sergey Bochkanov (ALGLIB project).
	3
	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
	9
	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
	14
	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
	17	>>> END OF LICENSE >>>
	18	*************************************************************************/
	19	#pragma warning disable 162
	20	#pragma warning disable 219
	21	using System;
	22
	23	public partial class alglib
	24	{
	25
	26
	27	/*************************************************************************
	28	Integration report:
	29	* TerminationType = completetion code:
	30	* -5 non-convergence of Gauss-Kronrod nodes
	31	calculation subroutine.
	32	* -1 incorrect parameters were specified
	33	* 1 OK
	34	* Rep.NFEV countains number of function calculations
	35	* Rep.NIntervals contains number of intervals [a,b]
	36	was partitioned into.
	37	*************************************************************************/
	38	public class autogkreport
	39	{
	40	//
	41	// Public declarations
	42	//
	43	public int terminationtype { get { return _innerobj.terminationtype; } set { _innerobj.terminationtype = value; } }
	44	public int nfev { get { return _innerobj.nfev; } set { _innerobj.nfev = value; } }
	45	public int nintervals { get { return _innerobj.nintervals; } set { _innerobj.nintervals = value; } }
	46
	47	public autogkreport()
	48	{
	49	_innerobj = new autogk.autogkreport();
	50	}
	51
	52	//
	53	// Although some of declarations below are public, you should not use them
	54	// They are intended for internal use only
	55	//
	56	private autogk.autogkreport _innerobj;
	57	public autogk.autogkreport innerobj { get { return _innerobj; } }
	58	public autogkreport(autogk.autogkreport obj)
	59	{
	60	_innerobj = obj;
	61	}
	62	}
	63
	64
	65	/*************************************************************************
	66	This structure stores state of the integration algorithm.
	67
	68	Although this class has public fields, they are not intended for external
	69	use. You should use ALGLIB functions to work with this class:
	70	* autogksmooth()/AutoGKSmoothW()/... to create objects
	71	* autogkintegrate() to begin integration
	72	* autogkresults() to get results
	73	*************************************************************************/
	74	public class autogkstate
	75	{
	76	//
	77	// Public declarations
	78	//
	79	public bool needf { get { return _innerobj.needf; } set { _innerobj.needf = value; } }
	80	public double x { get { return _innerobj.x; } set { _innerobj.x = value; } }
	81	public double xminusa { get { return _innerobj.xminusa; } set { _innerobj.xminusa = value; } }
	82	public double bminusx { get { return _innerobj.bminusx; } set { _innerobj.bminusx = value; } }
	83	public double f { get { return _innerobj.f; } set { _innerobj.f = value; } }
	84
	85	public autogkstate()
	86	{
	87	_innerobj = new autogk.autogkstate();
	88	}
	89
	90	//
	91	// Although some of declarations below are public, you should not use them
	92	// They are intended for internal use only
	93	//
	94	private autogk.autogkstate _innerobj;
	95	public autogk.autogkstate innerobj { get { return _innerobj; } }
	96	public autogkstate(autogk.autogkstate obj)
	97	{
	98	_innerobj = obj;
	99	}
	100	}
	101
	102	/*************************************************************************
	103	Integration of a smooth function F(x) on a finite interval [a,b].
	104
	105	Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
	106	is calculated with accuracy close to the machine precision.
	107
	108	Algorithm works well only with smooth integrands. It may be used with
	109	continuous non-smooth integrands, but with less performance.
	110
	111	It should never be used with integrands which have integrable singularities
	112	at lower or upper limits - algorithm may crash. Use AutoGKSingular in such
	113	cases.
	114
	115	INPUT PARAMETERS:
	116	A, B - interval boundaries (A<B, A=B or A>B)
	117
	118	OUTPUT PARAMETERS
	119	State - structure which stores algorithm state
	120
	121	SEE ALSO
	122	AutoGKSmoothW, AutoGKSingular, AutoGKResults.
	123
	124
	125	-- ALGLIB --
	126	Copyright 06.05.2009 by Bochkanov Sergey
	127	*************************************************************************/
	128	public static void autogksmooth(double a, double b, out autogkstate state)
	129	{
	130	state = new autogkstate();
	131	autogk.autogksmooth(a, b, state.innerobj);
	132	return;
	133	}
	134
	135	/*************************************************************************
	136	Integration of a smooth function F(x) on a finite interval [a,b].
	137
	138	This subroutine is same as AutoGKSmooth(), but it guarantees that interval
	139	[a,b] is partitioned into subintervals which have width at most XWidth.
	140
	141	Subroutine can be used when integrating nearly-constant function with
	142	narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth
	143	subroutine can overlook them.
	144
	145	INPUT PARAMETERS:
	146	A, B - interval boundaries (A<B, A=B or A>B)
	147
	148	OUTPUT PARAMETERS
	149	State - structure which stores algorithm state
	150
	151	SEE ALSO
	152	AutoGKSmooth, AutoGKSingular, AutoGKResults.
	153
	154
	155	-- ALGLIB --
	156	Copyright 06.05.2009 by Bochkanov Sergey
	157	*************************************************************************/
	158	public static void autogksmoothw(double a, double b, double xwidth, out autogkstate state)
	159	{
	160	state = new autogkstate();
	161	autogk.autogksmoothw(a, b, xwidth, state.innerobj);
	162	return;
	163	}
	164
	165	/*************************************************************************
	166	Integration on a finite interval [A,B].
	167	Integrand have integrable singularities at A/B.
	168
	169	F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known
	170	alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates
	171	from below can be used (but these estimates should be greater than -1 too).
	172
	173	One of alpha/beta variables (or even both alpha/beta) may be equal to 0,
	174	which means than function F(x) is non-singular at A/B. Anyway (singular at
	175	bounds or not), function F(x) is supposed to be continuous on (A,B).
	176
	177	Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
	178	is calculated with accuracy close to the machine precision.
	179
	180	INPUT PARAMETERS:
	181	A, B - interval boundaries (A<B, A=B or A>B)
	182	Alpha - power-law coefficient of the F(x) at A,
	183	Alpha>-1
	184	Beta - power-law coefficient of the F(x) at B,
	185	Beta>-1
	186
	187	OUTPUT PARAMETERS
	188	State - structure which stores algorithm state
	189
	190	SEE ALSO
	191	AutoGKSmooth, AutoGKSmoothW, AutoGKResults.
	192
	193
	194	-- ALGLIB --
	195	Copyright 06.05.2009 by Bochkanov Sergey
	196	*************************************************************************/
	197	public static void autogksingular(double a, double b, double alpha, double beta, out autogkstate state)
	198	{
	199	state = new autogkstate();
	200	autogk.autogksingular(a, b, alpha, beta, state.innerobj);
	201	return;
	202	}
	203
	204	/*************************************************************************
	205	This function provides reverse communication interface
	206	Reverse communication interface is not documented or recommended to use.
	207	See below for functions which provide better documented API
	208	*************************************************************************/
	209	public static bool autogkiteration(autogkstate state)
	210	{
	211
	212	bool result = autogk.autogkiteration(state.innerobj);
	213	return result;
	214	}
	215
	216
	217	/*************************************************************************
	218	This function is used to launcn iterations of ODE solver
	219
	220	It accepts following parameters:
	221	diff - callback which calculates dy/dx for given y and x
	222	obj - optional object which is passed to diff; can be NULL
	223
	224
	225	-- ALGLIB --
	226	Copyright 07.05.2009 by Bochkanov Sergey
	227
	228	*************************************************************************/
	229	public static void autogkintegrate(autogkstate state, integrator1_func func, object obj)
	230	{
	231	if( func==null )
	232	throw new alglibexception("ALGLIB: error in 'autogkintegrate()' (func is null)");
	233	while( alglib.autogkiteration(state) )
	234	{
	235	if( state.needf )
	236	{
	237	func(state.innerobj.x, state.innerobj.xminusa, state.innerobj.bminusx, ref state.innerobj.f, obj);
	238	continue;
	239	}
	240	throw new alglibexception("ALGLIB: unexpected error in 'autogksolve'");
	241	}
	242	}
	243
	244	/*************************************************************************
	245	Adaptive integration results
	246
	247	Called after AutoGKIteration returned False.
	248
	249	Input parameters:
	250	State - algorithm state (used by AutoGKIteration).
	251
	252	Output parameters:
	253	V - integral(f(x)dx,a,b)
	254	Rep - optimization report (see AutoGKReport description)
	255
	256	-- ALGLIB --
	257	Copyright 14.11.2007 by Bochkanov Sergey
	258	*************************************************************************/
	259	public static void autogkresults(autogkstate state, out double v, out autogkreport rep)
	260	{
	261	v = 0;
	262	rep = new autogkreport();
	263	autogk.autogkresults(state.innerobj, ref v, rep.innerobj);
	264	return;
	265	}
	266
	267	}
	268	public partial class alglib
	269	{
	270
	271
	272	/*************************************************************************
	273	Computation of nodes and weights for a Gauss quadrature formula
	274
	275	The algorithm generates the N-point Gauss quadrature formula with weight
	276	function given by coefficients alpha and beta of a recurrence relation
	277	which generates a system of orthogonal polynomials:
	278
	279	P-1(x) = 0
	280	P0(x) = 1
	281	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	282
	283	and zeroth moment Mu0
	284
	285	Mu0 = integral(W(x)dx,a,b)
	286
	287	INPUT PARAMETERS:
	288	Alpha array[0..N-1], alpha coefficients
	289	Beta array[0..N-1], beta coefficients
	290	Zero-indexed element is not used and may be arbitrary.
	291	Beta[I]>0.
	292	Mu0 zeroth moment of the weight function.
	293	N number of nodes of the quadrature formula, N>=1
	294
	295	OUTPUT PARAMETERS:
	296	Info - error code:
	297	* -3 internal eigenproblem solver hasn't converged
	298	* -2 Beta[i]<=0
	299	* -1 incorrect N was passed
	300	* 1 OK
	301	X - array[0..N-1] - array of quadrature nodes,
	302	in ascending order.
	303	W - array[0..N-1] - array of quadrature weights.
	304
	305	-- ALGLIB --
	306	Copyright 2005-2009 by Bochkanov Sergey
	307	*************************************************************************/
	308	public static void gqgeneraterec(double[] alpha, double[] beta, double mu0, int n, out int info, out double[] x, out double[] w)
	309	{
	310	info = 0;
	311	x = new double[0];
	312	w = new double[0];
	313	gq.gqgeneraterec(alpha, beta, mu0, n, ref info, ref x, ref w);
	314	return;
	315	}
	316
	317	/*************************************************************************
	318	Computation of nodes and weights for a Gauss-Lobatto quadrature formula
	319
	320	The algorithm generates the N-point Gauss-Lobatto quadrature formula with
	321	weight function given by coefficients alpha and beta of a recurrence which
	322	generates a system of orthogonal polynomials.
	323
	324	P-1(x) = 0
	325	P0(x) = 1
	326	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	327
	328	and zeroth moment Mu0
	329
	330	Mu0 = integral(W(x)dx,a,b)
	331
	332	INPUT PARAMETERS:
	333	Alpha array[0..N-2], alpha coefficients
	334	Beta array[0..N-2], beta coefficients.
	335	Zero-indexed element is not used, may be arbitrary.
	336	Beta[I]>0
	337	Mu0 zeroth moment of the weighting function.
	338	A left boundary of the integration interval.
	339	B right boundary of the integration interval.
	340	N number of nodes of the quadrature formula, N>=3
	341	(including the left and right boundary nodes).
	342
	343	OUTPUT PARAMETERS:
	344	Info - error code:
	345	* -3 internal eigenproblem solver hasn't converged
	346	* -2 Beta[i]<=0
	347	* -1 incorrect N was passed
	348	* 1 OK
	349	X - array[0..N-1] - array of quadrature nodes,
	350	in ascending order.
	351	W - array[0..N-1] - array of quadrature weights.
	352
	353	-- ALGLIB --
	354	Copyright 2005-2009 by Bochkanov Sergey
	355	*************************************************************************/
	356	public static void gqgenerategausslobattorec(double[] alpha, double[] beta, double mu0, double a, double b, int n, out int info, out double[] x, out double[] w)
	357	{
	358	info = 0;
	359	x = new double[0];
	360	w = new double[0];
	361	gq.gqgenerategausslobattorec(alpha, beta, mu0, a, b, n, ref info, ref x, ref w);
	362	return;
	363	}
	364
	365	/*************************************************************************
	366	Computation of nodes and weights for a Gauss-Radau quadrature formula
	367
	368	The algorithm generates the N-point Gauss-Radau quadrature formula with
	369	weight function given by the coefficients alpha and beta of a recurrence
	370	which generates a system of orthogonal polynomials.
	371
	372	P-1(x) = 0
	373	P0(x) = 1
	374	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	375
	376	and zeroth moment Mu0
	377
	378	Mu0 = integral(W(x)dx,a,b)
	379
	380	INPUT PARAMETERS:
	381	Alpha array[0..N-2], alpha coefficients.
	382	Beta array[0..N-1], beta coefficients
	383	Zero-indexed element is not used.
	384	Beta[I]>0
	385	Mu0 zeroth moment of the weighting function.
	386	A left boundary of the integration interval.
	387	N number of nodes of the quadrature formula, N>=2
	388	(including the left boundary node).
	389
	390	OUTPUT PARAMETERS:
	391	Info - error code:
	392	* -3 internal eigenproblem solver hasn't converged
	393	* -2 Beta[i]<=0
	394	* -1 incorrect N was passed
	395	* 1 OK
	396	X - array[0..N-1] - array of quadrature nodes,
	397	in ascending order.
	398	W - array[0..N-1] - array of quadrature weights.
	399
	400
	401	-- ALGLIB --
	402	Copyright 2005-2009 by Bochkanov Sergey
	403	*************************************************************************/
	404	public static void gqgenerategaussradaurec(double[] alpha, double[] beta, double mu0, double a, int n, out int info, out double[] x, out double[] w)
	405	{
	406	info = 0;
	407	x = new double[0];
	408	w = new double[0];
	409	gq.gqgenerategaussradaurec(alpha, beta, mu0, a, n, ref info, ref x, ref w);
	410	return;
	411	}
	412
	413	/*************************************************************************
	414	Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
	415	nodes.
	416
	417	INPUT PARAMETERS:
	418	N - number of nodes, >=1
	419
	420	OUTPUT PARAMETERS:
	421	Info - error code:
	422	* -4 an error was detected when calculating
	423	weights/nodes. N is too large to obtain
	424	weights/nodes with high enough accuracy.
	425	Try to use multiple precision version.
	426	* -3 internal eigenproblem solver hasn't converged
	427	* -1 incorrect N was passed
	428	* +1 OK
	429	X - array[0..N-1] - array of quadrature nodes,
	430	in ascending order.
	431	W - array[0..N-1] - array of quadrature weights.
	432
	433
	434	-- ALGLIB --
	435	Copyright 12.05.2009 by Bochkanov Sergey
	436	*************************************************************************/
	437	public static void gqgenerategausslegendre(int n, out int info, out double[] x, out double[] w)
	438	{
	439	info = 0;
	440	x = new double[0];
	441	w = new double[0];
	442	gq.gqgenerategausslegendre(n, ref info, ref x, ref w);
	443	return;
	444	}
	445
	446	/*************************************************************************
	447	Returns nodes/weights for Gauss-Jacobi quadrature on [-1,1] with weight
	448	function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
	449
	450	INPUT PARAMETERS:
	451	N - number of nodes, >=1
	452	Alpha - power-law coefficient, Alpha>-1
	453	Beta - power-law coefficient, Beta>-1
	454
	455	OUTPUT PARAMETERS:
	456	Info - error code:
	457	* -4 an error was detected when calculating
	458	weights/nodes. Alpha or Beta are too close
	459	to -1 to obtain weights/nodes with high enough
	460	accuracy, or, may be, N is too large. Try to
	461	use multiple precision version.
	462	* -3 internal eigenproblem solver hasn't converged
	463	* -1 incorrect N/Alpha/Beta was passed
	464	* +1 OK
	465	X - array[0..N-1] - array of quadrature nodes,
	466	in ascending order.
	467	W - array[0..N-1] - array of quadrature weights.
	468
	469
	470	-- ALGLIB --
	471	Copyright 12.05.2009 by Bochkanov Sergey
	472	*************************************************************************/
	473	public static void gqgenerategaussjacobi(int n, double alpha, double beta, out int info, out double[] x, out double[] w)
	474	{
	475	info = 0;
	476	x = new double[0];
	477	w = new double[0];
	478	gq.gqgenerategaussjacobi(n, alpha, beta, ref info, ref x, ref w);
	479	return;
	480	}
	481
	482	/*************************************************************************
	483	Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf) with
	484	weight function W(x)=Power(x,Alpha)*Exp(-x)
	485
	486	INPUT PARAMETERS:
	487	N - number of nodes, >=1
	488	Alpha - power-law coefficient, Alpha>-1
	489
	490	OUTPUT PARAMETERS:
	491	Info - error code:
	492	* -4 an error was detected when calculating
	493	weights/nodes. Alpha is too close to -1 to
	494	obtain weights/nodes with high enough accuracy
	495	or, may be, N is too large. Try to use
	496	multiple precision version.
	497	* -3 internal eigenproblem solver hasn't converged
	498	* -1 incorrect N/Alpha was passed
	499	* +1 OK
	500	X - array[0..N-1] - array of quadrature nodes,
	501	in ascending order.
	502	W - array[0..N-1] - array of quadrature weights.
	503
	504
	505	-- ALGLIB --
	506	Copyright 12.05.2009 by Bochkanov Sergey
	507	*************************************************************************/
	508	public static void gqgenerategausslaguerre(int n, double alpha, out int info, out double[] x, out double[] w)
	509	{
	510	info = 0;
	511	x = new double[0];
	512	w = new double[0];
	513	gq.gqgenerategausslaguerre(n, alpha, ref info, ref x, ref w);
	514	return;
	515	}
	516
	517	/*************************************************************************
	518	Returns nodes/weights for Gauss-Hermite quadrature on (-inf,+inf) with
	519	weight function W(x)=Exp(-x*x)
	520
	521	INPUT PARAMETERS:
	522	N - number of nodes, >=1
	523
	524	OUTPUT PARAMETERS:
	525	Info - error code:
	526	* -4 an error was detected when calculating
	527	weights/nodes. May be, N is too large. Try to
	528	use multiple precision version.
	529	* -3 internal eigenproblem solver hasn't converged
	530	* -1 incorrect N/Alpha was passed
	531	* +1 OK
	532	X - array[0..N-1] - array of quadrature nodes,
	533	in ascending order.
	534	W - array[0..N-1] - array of quadrature weights.
	535
	536
	537	-- ALGLIB --
	538	Copyright 12.05.2009 by Bochkanov Sergey
	539	*************************************************************************/
	540	public static void gqgenerategausshermite(int n, out int info, out double[] x, out double[] w)
	541	{
	542	info = 0;
	543	x = new double[0];
	544	w = new double[0];
	545	gq.gqgenerategausshermite(n, ref info, ref x, ref w);
	546	return;
	547	}
	548
	549	}
	550	public partial class alglib
	551	{
	552
	553
	554	/*************************************************************************
	555	Computation of nodes and weights of a Gauss-Kronrod quadrature formula
	556
	557	The algorithm generates the N-point Gauss-Kronrod quadrature formula with
	558	weight function given by coefficients alpha and beta of a recurrence
	559	relation which generates a system of orthogonal polynomials:
	560
	561	P-1(x) = 0
	562	P0(x) = 1
	563	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	564
	565	and zero moment Mu0
	566
	567	Mu0 = integral(W(x)dx,a,b)
	568
	569
	570	INPUT PARAMETERS:
	571	Alpha alpha coefficients, array[0..floor(3*K/2)].
	572	Beta beta coefficients, array[0..ceil(3*K/2)].
	573	Beta[0] is not used and may be arbitrary.
	574	Beta[I]>0.
	575	Mu0 zeroth moment of the weight function.
	576	N number of nodes of the Gauss-Kronrod quadrature formula,
	577	N >= 3,
	578	N = 2*K+1.
	579
	580	OUTPUT PARAMETERS:
	581	Info - error code:
	582	* -5 no real and positive Gauss-Kronrod formula can
	583	be created for such a weight function with a
	584	given number of nodes.
	585	* -4 N is too large, task may be ill conditioned -
	586	x[i]=x[i+1] found.
	587	* -3 internal eigenproblem solver hasn't converged
	588	* -2 Beta[i]<=0
	589	* -1 incorrect N was passed
	590	* +1 OK
	591	X - array[0..N-1] - array of quadrature nodes,
	592	in ascending order.
	593	WKronrod - array[0..N-1] - Kronrod weights
	594	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	595	corresponding to extended Kronrod nodes).
	596
	597	-- ALGLIB --
	598	Copyright 08.05.2009 by Bochkanov Sergey
	599	*************************************************************************/
	600	public static void gkqgeneraterec(double[] alpha, double[] beta, double mu0, int n, out int info, out double[] x, out double[] wkronrod, out double[] wgauss)
	601	{
	602	info = 0;
	603	x = new double[0];
	604	wkronrod = new double[0];
	605	wgauss = new double[0];
	606	gkq.gkqgeneraterec(alpha, beta, mu0, n, ref info, ref x, ref wkronrod, ref wgauss);
	607	return;
	608	}
	609
	610	/*************************************************************************
	611	Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre
	612	quadrature with N points.
	613
	614	GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed table) is
	615	used depending on machine precision and number of nodes.
	616
	617	INPUT PARAMETERS:
	618	N - number of Kronrod nodes, must be odd number, >=3.
	619
	620	OUTPUT PARAMETERS:
	621	Info - error code:
	622	* -4 an error was detected when calculating
	623	weights/nodes. N is too large to obtain
	624	weights/nodes with high enough accuracy.
	625	Try to use multiple precision version.
	626	* -3 internal eigenproblem solver hasn't converged
	627	* -1 incorrect N was passed
	628	* +1 OK
	629	X - array[0..N-1] - array of quadrature nodes, ordered in
	630	ascending order.
	631	WKronrod - array[0..N-1] - Kronrod weights
	632	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	633	corresponding to extended Kronrod nodes).
	634
	635
	636	-- ALGLIB --
	637	Copyright 12.05.2009 by Bochkanov Sergey
	638	*************************************************************************/
	639	public static void gkqgenerategausslegendre(int n, out int info, out double[] x, out double[] wkronrod, out double[] wgauss)
	640	{
	641	info = 0;
	642	x = new double[0];
	643	wkronrod = new double[0];
	644	wgauss = new double[0];
	645	gkq.gkqgenerategausslegendre(n, ref info, ref x, ref wkronrod, ref wgauss);
	646	return;
	647	}
	648
	649	/*************************************************************************
	650	Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi
	651	quadrature on [-1,1] with weight function
	652
	653	W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
	654
	655	INPUT PARAMETERS:
	656	N - number of Kronrod nodes, must be odd number, >=3.
	657	Alpha - power-law coefficient, Alpha>-1
	658	Beta - power-law coefficient, Beta>-1
	659
	660	OUTPUT PARAMETERS:
	661	Info - error code:
	662	* -5 no real and positive Gauss-Kronrod formula can
	663	be created for such a weight function with a
	664	given number of nodes.
	665	* -4 an error was detected when calculating
	666	weights/nodes. Alpha or Beta are too close
	667	to -1 to obtain weights/nodes with high enough
	668	accuracy, or, may be, N is too large. Try to
	669	use multiple precision version.
	670	* -3 internal eigenproblem solver hasn't converged
	671	* -1 incorrect N was passed
	672	* +1 OK
	673	* +2 OK, but quadrature rule have exterior nodes,
	674	x[0]<-1 or x[n-1]>+1
	675	X - array[0..N-1] - array of quadrature nodes, ordered in
	676	ascending order.
	677	WKronrod - array[0..N-1] - Kronrod weights
	678	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	679	corresponding to extended Kronrod nodes).
	680
	681
	682	-- ALGLIB --
	683	Copyright 12.05.2009 by Bochkanov Sergey
	684	*************************************************************************/
	685	public static void gkqgenerategaussjacobi(int n, double alpha, double beta, out int info, out double[] x, out double[] wkronrod, out double[] wgauss)
	686	{
	687	info = 0;
	688	x = new double[0];
	689	wkronrod = new double[0];
	690	wgauss = new double[0];
	691	gkq.gkqgenerategaussjacobi(n, alpha, beta, ref info, ref x, ref wkronrod, ref wgauss);
	692	return;
	693	}
	694
	695	/*************************************************************************
	696	Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
	697
	698	Reduction to tridiagonal eigenproblem is used.
	699
	700	INPUT PARAMETERS:
	701	N - number of Kronrod nodes, must be odd number, >=3.
	702
	703	OUTPUT PARAMETERS:
	704	Info - error code:
	705	* -4 an error was detected when calculating
	706	weights/nodes. N is too large to obtain
	707	weights/nodes with high enough accuracy.
	708	Try to use multiple precision version.
	709	* -3 internal eigenproblem solver hasn't converged
	710	* -1 incorrect N was passed
	711	* +1 OK
	712	X - array[0..N-1] - array of quadrature nodes, ordered in
	713	ascending order.
	714	WKronrod - array[0..N-1] - Kronrod weights
	715	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	716	corresponding to extended Kronrod nodes).
	717
	718	-- ALGLIB --
	719	Copyright 12.05.2009 by Bochkanov Sergey
	720	*************************************************************************/
	721	public static void gkqlegendrecalc(int n, out int info, out double[] x, out double[] wkronrod, out double[] wgauss)
	722	{
	723	info = 0;
	724	x = new double[0];
	725	wkronrod = new double[0];
	726	wgauss = new double[0];
	727	gkq.gkqlegendrecalc(n, ref info, ref x, ref wkronrod, ref wgauss);
	728	return;
	729	}
	730
	731	/*************************************************************************
	732	Returns Gauss and Gauss-Kronrod nodes for quadrature with N points using
	733	pre-calculated table. Nodes/weights were computed with accuracy up to
	734	1.0E-32 (if MPFR version of ALGLIB is used). In standard double precision
	735	accuracy reduces to something about 2.0E-16 (depending on your compiler's
	736	handling of long floating point constants).
	737
	738	INPUT PARAMETERS:
	739	N - number of Kronrod nodes.
	740	N can be 15, 21, 31, 41, 51, 61.
	741
	742	OUTPUT PARAMETERS:
	743	X - array[0..N-1] - array of quadrature nodes, ordered in
	744	ascending order.
	745	WKronrod - array[0..N-1] - Kronrod weights
	746	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	747	corresponding to extended Kronrod nodes).
	748
	749
	750	-- ALGLIB --
	751	Copyright 12.05.2009 by Bochkanov Sergey
	752	*************************************************************************/
	753	public static void gkqlegendretbl(int n, out double[] x, out double[] wkronrod, out double[] wgauss, out double eps)
	754	{
	755	x = new double[0];
	756	wkronrod = new double[0];
	757	wgauss = new double[0];
	758	eps = 0;
	759	gkq.gkqlegendretbl(n, ref x, ref wkronrod, ref wgauss, ref eps);
	760	return;
	761	}
	762
	763	}
	764	public partial class alglib
	765	{
	766	public class autogk
	767	{
	768	/*************************************************************************
	769	Integration report:
	770	* TerminationType = completetion code:
	771	* -5 non-convergence of Gauss-Kronrod nodes
	772	calculation subroutine.
	773	* -1 incorrect parameters were specified
	774	* 1 OK
	775	* Rep.NFEV countains number of function calculations
	776	* Rep.NIntervals contains number of intervals [a,b]
	777	was partitioned into.
	778	*************************************************************************/
	779	public class autogkreport
	780	{
	781	public int terminationtype;
	782	public int nfev;
	783	public int nintervals;
	784	};
	785
	786
	787	public class autogkinternalstate
	788	{
	789	public double a;
	790	public double b;
	791	public double eps;
	792	public double xwidth;
	793	public double x;
	794	public double f;
	795	public int info;
	796	public double r;
	797	public double[,] heap;
	798	public int heapsize;
	799	public int heapwidth;
	800	public int heapused;
	801	public double sumerr;
	802	public double sumabs;
	803	public double[] qn;
	804	public double[] wg;
	805	public double[] wk;
	806	public double[] wr;
	807	public int n;
	808	public rcommstate rstate;
	809	public autogkinternalstate()
	810	{
	811	heap = new double[0,0];
	812	qn = new double[0];
	813	wg = new double[0];
	814	wk = new double[0];
	815	wr = new double[0];
	816	rstate = new rcommstate();
	817	}
	818	};
	819
	820
	821	/*************************************************************************
	822	This structure stores state of the integration algorithm.
	823
	824	Although this class has public fields, they are not intended for external
	825	use. You should use ALGLIB functions to work with this class:
	826	* autogksmooth()/AutoGKSmoothW()/... to create objects
	827	* autogkintegrate() to begin integration
	828	* autogkresults() to get results
	829	*************************************************************************/
	830	public class autogkstate
	831	{
	832	public double a;
	833	public double b;
	834	public double alpha;
	835	public double beta;
	836	public double xwidth;
	837	public double x;
	838	public double xminusa;
	839	public double bminusx;
	840	public bool needf;
	841	public double f;
	842	public int wrappermode;
	843	public autogkinternalstate internalstate;
	844	public rcommstate rstate;
	845	public double v;
	846	public int terminationtype;
	847	public int nfev;
	848	public int nintervals;
	849	public autogkstate()
	850	{
	851	internalstate = new autogkinternalstate();
	852	rstate = new rcommstate();
	853	}
	854	};
	855
	856
	857
	858
	859	/*************************************************************************
	860	Integration of a smooth function F(x) on a finite interval [a,b].
	861
	862	Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
	863	is calculated with accuracy close to the machine precision.
	864
	865	Algorithm works well only with smooth integrands. It may be used with
	866	continuous non-smooth integrands, but with less performance.
	867
	868	It should never be used with integrands which have integrable singularities
	869	at lower or upper limits - algorithm may crash. Use AutoGKSingular in such
	870	cases.
	871
	872	INPUT PARAMETERS:
	873	A, B - interval boundaries (A<B, A=B or A>B)
	874
	875	OUTPUT PARAMETERS
	876	State - structure which stores algorithm state
	877
	878	SEE ALSO
	879	AutoGKSmoothW, AutoGKSingular, AutoGKResults.
	880
	881
	882	-- ALGLIB --
	883	Copyright 06.05.2009 by Bochkanov Sergey
	884	*************************************************************************/
	885	public static void autogksmooth(double a,
	886	double b,
	887	autogkstate state)
	888	{
	889	ap.assert(math.isfinite(a), "AutoGKSmooth: A is not finite!");
	890	ap.assert(math.isfinite(b), "AutoGKSmooth: B is not finite!");
	891	autogksmoothw(a, b, 0.0, state);
	892	}
	893
	894
	895	/*************************************************************************
	896	Integration of a smooth function F(x) on a finite interval [a,b].
	897
	898	This subroutine is same as AutoGKSmooth(), but it guarantees that interval
	899	[a,b] is partitioned into subintervals which have width at most XWidth.
	900
	901	Subroutine can be used when integrating nearly-constant function with
	902	narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth
	903	subroutine can overlook them.
	904
	905	INPUT PARAMETERS:
	906	A, B - interval boundaries (A<B, A=B or A>B)
	907
	908	OUTPUT PARAMETERS
	909	State - structure which stores algorithm state
	910
	911	SEE ALSO
	912	AutoGKSmooth, AutoGKSingular, AutoGKResults.
	913
	914
	915	-- ALGLIB --
	916	Copyright 06.05.2009 by Bochkanov Sergey
	917	*************************************************************************/
	918	public static void autogksmoothw(double a,
	919	double b,
	920	double xwidth,
	921	autogkstate state)
	922	{
	923	ap.assert(math.isfinite(a), "AutoGKSmoothW: A is not finite!");
	924	ap.assert(math.isfinite(b), "AutoGKSmoothW: B is not finite!");
	925	ap.assert(math.isfinite(xwidth), "AutoGKSmoothW: XWidth is not finite!");
	926	state.wrappermode = 0;
	927	state.a = a;
	928	state.b = b;
	929	state.xwidth = xwidth;
	930	state.needf = false;
	931	state.rstate.ra = new double[10+1];
	932	state.rstate.stage = -1;
	933	}
	934
	935
	936	/*************************************************************************
	937	Integration on a finite interval [A,B].
	938	Integrand have integrable singularities at A/B.
	939
	940	F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known
	941	alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates
	942	from below can be used (but these estimates should be greater than -1 too).
	943
	944	One of alpha/beta variables (or even both alpha/beta) may be equal to 0,
	945	which means than function F(x) is non-singular at A/B. Anyway (singular at
	946	bounds or not), function F(x) is supposed to be continuous on (A,B).
	947
	948	Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
	949	is calculated with accuracy close to the machine precision.
	950
	951	INPUT PARAMETERS:
	952	A, B - interval boundaries (A<B, A=B or A>B)
	953	Alpha - power-law coefficient of the F(x) at A,
	954	Alpha>-1
	955	Beta - power-law coefficient of the F(x) at B,
	956	Beta>-1
	957
	958	OUTPUT PARAMETERS
	959	State - structure which stores algorithm state
	960
	961	SEE ALSO
	962	AutoGKSmooth, AutoGKSmoothW, AutoGKResults.
	963
	964
	965	-- ALGLIB --
	966	Copyright 06.05.2009 by Bochkanov Sergey
	967	*************************************************************************/
	968	public static void autogksingular(double a,
	969	double b,
	970	double alpha,
	971	double beta,
	972	autogkstate state)
	973	{
	974	ap.assert(math.isfinite(a), "AutoGKSingular: A is not finite!");
	975	ap.assert(math.isfinite(b), "AutoGKSingular: B is not finite!");
	976	ap.assert(math.isfinite(alpha), "AutoGKSingular: Alpha is not finite!");
	977	ap.assert(math.isfinite(beta), "AutoGKSingular: Beta is not finite!");
	978	state.wrappermode = 1;
	979	state.a = a;
	980	state.b = b;
	981	state.alpha = alpha;
	982	state.beta = beta;
	983	state.xwidth = 0.0;
	984	state.needf = false;
	985	state.rstate.ra = new double[10+1];
	986	state.rstate.stage = -1;
	987	}
	988
	989
	990	/*************************************************************************
	991
	992	-- ALGLIB --
	993	Copyright 07.05.2009 by Bochkanov Sergey
	994	*************************************************************************/
	995	public static bool autogkiteration(autogkstate state)
	996	{
	997	bool result = new bool();
	998	double s = 0;
	999	double tmp = 0;
	1000	double eps = 0;
	1001	double a = 0;
	1002	double b = 0;
	1003	double x = 0;
	1004	double t = 0;
	1005	double alpha = 0;
	1006	double beta = 0;
	1007	double v1 = 0;
	1008	double v2 = 0;
	1009
	1010
	1011	//
	1012	// Reverse communication preparations
	1013	// I know it looks ugly, but it works the same way
	1014	// anywhere from C++ to Python.
	1015	//
	1016	// This code initializes locals by:
	1017	// * random values determined during code
	1018	// generation - on first subroutine call
	1019	// * values from previous call - on subsequent calls
	1020	//
	1021	if( state.rstate.stage>=0 )
	1022	{
	1023	s = state.rstate.ra[0];
	1024	tmp = state.rstate.ra[1];
	1025	eps = state.rstate.ra[2];
	1026	a = state.rstate.ra[3];
	1027	b = state.rstate.ra[4];
	1028	x = state.rstate.ra[5];
	1029	t = state.rstate.ra[6];
	1030	alpha = state.rstate.ra[7];
	1031	beta = state.rstate.ra[8];
	1032	v1 = state.rstate.ra[9];
	1033	v2 = state.rstate.ra[10];
	1034	}
	1035	else
	1036	{
	1037	s = -983;
	1038	tmp = -989;
	1039	eps = -834;
	1040	a = 900;
	1041	b = -287;
	1042	x = 364;
	1043	t = 214;
	1044	alpha = -338;
	1045	beta = -686;
	1046	v1 = 912;
	1047	v2 = 585;
	1048	}
	1049	if( state.rstate.stage==0 )
	1050	{
	1051	goto lbl_0;
	1052	}
	1053	if( state.rstate.stage==1 )
	1054	{
	1055	goto lbl_1;
	1056	}
	1057	if( state.rstate.stage==2 )
	1058	{
	1059	goto lbl_2;
	1060	}
	1061
	1062	//
	1063	// Routine body
	1064	//
	1065	eps = 0;
	1066	a = state.a;
	1067	b = state.b;
	1068	alpha = state.alpha;
	1069	beta = state.beta;
	1070	state.terminationtype = -1;
	1071	state.nfev = 0;
	1072	state.nintervals = 0;
	1073
	1074	//
	1075	// smooth function at a finite interval
	1076	//
	1077	if( state.wrappermode!=0 )
	1078	{
	1079	goto lbl_3;
	1080	}
	1081
	1082	//
	1083	// special case
	1084	//
	1085	if( (double)(a)==(double)(b) )
	1086	{
	1087	state.terminationtype = 1;
	1088	state.v = 0;
	1089	result = false;
	1090	return result;
	1091	}
	1092
	1093	//
	1094	// general case
	1095	//
	1096	autogkinternalprepare(a, b, eps, state.xwidth, state.internalstate);
	1097	lbl_5:
	1098	if( !autogkinternaliteration(state.internalstate) )
	1099	{
	1100	goto lbl_6;
	1101	}
	1102	x = state.internalstate.x;
	1103	state.x = x;
	1104	state.xminusa = x-a;
	1105	state.bminusx = b-x;
	1106	state.needf = true;
	1107	state.rstate.stage = 0;
	1108	goto lbl_rcomm;
	1109	lbl_0:
	1110	state.needf = false;
	1111	state.nfev = state.nfev+1;
	1112	state.internalstate.f = state.f;
	1113	goto lbl_5;
	1114	lbl_6:
	1115	state.v = state.internalstate.r;
	1116	state.terminationtype = state.internalstate.info;
	1117	state.nintervals = state.internalstate.heapused;
	1118	result = false;
	1119	return result;
	1120	lbl_3:
	1121
	1122	//
	1123	// function with power-law singularities at the ends of a finite interval
	1124	//
	1125	if( state.wrappermode!=1 )
	1126	{
	1127	goto lbl_7;
	1128	}
	1129
	1130	//
	1131	// test coefficients
	1132	//
	1133	if( (double)(alpha)<=(double)(-1) \| (double)(beta)<=(double)(-1) )
	1134	{
	1135	state.terminationtype = -1;
	1136	state.v = 0;
	1137	result = false;
	1138	return result;
	1139	}
	1140
	1141	//
	1142	// special cases
	1143	//
	1144	if( (double)(a)==(double)(b) )
	1145	{
	1146	state.terminationtype = 1;
	1147	state.v = 0;
	1148	result = false;
	1149	return result;
	1150	}
	1151
	1152	//
	1153	// reduction to general form
	1154	//
	1155	if( (double)(a)<(double)(b) )
	1156	{
	1157	s = 1;
	1158	}
	1159	else
	1160	{
	1161	s = -1;
	1162	tmp = a;
	1163	a = b;
	1164	b = tmp;
	1165	tmp = alpha;
	1166	alpha = beta;
	1167	beta = tmp;
	1168	}
	1169	alpha = Math.Min(alpha, 0);
	1170	beta = Math.Min(beta, 0);
	1171
	1172	//
	1173	// first, integrate left half of [a,b]:
	1174	// integral(f(x)dx, a, (b+a)/2) =
	1175	// = 1/(1+alpha) * integral(t^(-alpha/(1+alpha))f(a+t^(1/(1+alpha)))dt, 0, (0.5(b-a))^(1+alpha))
	1176	//
	1177	autogkinternalprepare(0, Math.Pow(0.5*(b-a), 1+alpha), eps, state.xwidth, state.internalstate);
	1178	lbl_9:
	1179	if( !autogkinternaliteration(state.internalstate) )
	1180	{
	1181	goto lbl_10;
	1182	}
	1183
	1184	//
	1185	// Fill State.X, State.XMinusA, State.BMinusX.
	1186	// Latter two are filled correctly even if B<A.
	1187	//
	1188	x = state.internalstate.x;
	1189	t = Math.Pow(x, 1/(1+alpha));
	1190	state.x = a+t;
	1191	if( (double)(s)>(double)(0) )
	1192	{
	1193	state.xminusa = t;
	1194	state.bminusx = b-(a+t);
	1195	}
	1196	else
	1197	{
	1198	state.xminusa = a+t-b;
	1199	state.bminusx = -t;
	1200	}
	1201	state.needf = true;
	1202	state.rstate.stage = 1;
	1203	goto lbl_rcomm;
	1204	lbl_1:
	1205	state.needf = false;
	1206	if( (double)(alpha)!=(double)(0) )
	1207	{
	1208	state.internalstate.f = state.f*Math.Pow(x, -(alpha/(1+alpha)))/(1+alpha);
	1209	}
	1210	else
	1211	{
	1212	state.internalstate.f = state.f;
	1213	}
	1214	state.nfev = state.nfev+1;
	1215	goto lbl_9;
	1216	lbl_10:
	1217	v1 = state.internalstate.r;
	1218	state.nintervals = state.nintervals+state.internalstate.heapused;
	1219
	1220	//
	1221	// then, integrate right half of [a,b]:
	1222	// integral(f(x)dx, (b+a)/2, b) =
	1223	// = 1/(1+beta) * integral(t^(-beta/(1+beta))f(b-t^(1/(1+beta)))dt, 0, (0.5(b-a))^(1+beta))
	1224	//
	1225	autogkinternalprepare(0, Math.Pow(0.5*(b-a), 1+beta), eps, state.xwidth, state.internalstate);
	1226	lbl_11:
	1227	if( !autogkinternaliteration(state.internalstate) )
	1228	{
	1229	goto lbl_12;
	1230	}
	1231
	1232	//
	1233	// Fill State.X, State.XMinusA, State.BMinusX.
	1234	// Latter two are filled correctly (X-A, B-X) even if B<A.
	1235	//
	1236	x = state.internalstate.x;
	1237	t = Math.Pow(x, 1/(1+beta));
	1238	state.x = b-t;
	1239	if( (double)(s)>(double)(0) )
	1240	{
	1241	state.xminusa = b-t-a;
	1242	state.bminusx = t;
	1243	}
	1244	else
	1245	{
	1246	state.xminusa = -t;
	1247	state.bminusx = a-(b-t);
	1248	}
	1249	state.needf = true;
	1250	state.rstate.stage = 2;
	1251	goto lbl_rcomm;
	1252	lbl_2:
	1253	state.needf = false;
	1254	if( (double)(beta)!=(double)(0) )
	1255	{
	1256	state.internalstate.f = state.f*Math.Pow(x, -(beta/(1+beta)))/(1+beta);
	1257	}
	1258	else
	1259	{
	1260	state.internalstate.f = state.f;
	1261	}
	1262	state.nfev = state.nfev+1;
	1263	goto lbl_11;
	1264	lbl_12:
	1265	v2 = state.internalstate.r;
	1266	state.nintervals = state.nintervals+state.internalstate.heapused;
	1267
	1268	//
	1269	// final result
	1270	//
	1271	state.v = s*(v1+v2);
	1272	state.terminationtype = 1;
	1273	result = false;
	1274	return result;
	1275	lbl_7:
	1276	result = false;
	1277	return result;
	1278
	1279	//
	1280	// Saving state
	1281	//
	1282	lbl_rcomm:
	1283	result = true;
	1284	state.rstate.ra[0] = s;
	1285	state.rstate.ra[1] = tmp;
	1286	state.rstate.ra[2] = eps;
	1287	state.rstate.ra[3] = a;
	1288	state.rstate.ra[4] = b;
	1289	state.rstate.ra[5] = x;
	1290	state.rstate.ra[6] = t;
	1291	state.rstate.ra[7] = alpha;
	1292	state.rstate.ra[8] = beta;
	1293	state.rstate.ra[9] = v1;
	1294	state.rstate.ra[10] = v2;
	1295	return result;
	1296	}
	1297
	1298
	1299	/*************************************************************************
	1300	Adaptive integration results
	1301
	1302	Called after AutoGKIteration returned False.
	1303
	1304	Input parameters:
	1305	State - algorithm state (used by AutoGKIteration).
	1306
	1307	Output parameters:
	1308	V - integral(f(x)dx,a,b)
	1309	Rep - optimization report (see AutoGKReport description)
	1310
	1311	-- ALGLIB --
	1312	Copyright 14.11.2007 by Bochkanov Sergey
	1313	*************************************************************************/
	1314	public static void autogkresults(autogkstate state,
	1315	ref double v,
	1316	autogkreport rep)
	1317	{
	1318	v = 0;
	1319
	1320	v = state.v;
	1321	rep.terminationtype = state.terminationtype;
	1322	rep.nfev = state.nfev;
	1323	rep.nintervals = state.nintervals;
	1324	}
	1325
	1326
	1327	/*************************************************************************
	1328	Internal AutoGK subroutine
	1329	eps<0 - error
	1330	eps=0 - automatic eps selection
	1331
	1332	width<0 - error
	1333	width=0 - no width requirements
	1334	*************************************************************************/
	1335	private static void autogkinternalprepare(double a,
	1336	double b,
	1337	double eps,
	1338	double xwidth,
	1339	autogkinternalstate state)
	1340	{
	1341
	1342	//
	1343	// Save settings
	1344	//
	1345	state.a = a;
	1346	state.b = b;
	1347	state.eps = eps;
	1348	state.xwidth = xwidth;
	1349
	1350	//
	1351	// Prepare RComm structure
	1352	//
	1353	state.rstate.ia = new int[3+1];
	1354	state.rstate.ra = new double[8+1];
	1355	state.rstate.stage = -1;
	1356	}
	1357
	1358
	1359	/*************************************************************************
	1360	Internal AutoGK subroutine
	1361	*************************************************************************/
	1362	private static bool autogkinternaliteration(autogkinternalstate state)
	1363	{
	1364	bool result = new bool();
	1365	double c1 = 0;
	1366	double c2 = 0;
	1367	int i = 0;
	1368	int j = 0;
	1369	double intg = 0;
	1370	double intk = 0;
	1371	double inta = 0;
	1372	double v = 0;
	1373	double ta = 0;
	1374	double tb = 0;
	1375	int ns = 0;
	1376	double qeps = 0;
	1377	int info = 0;
	1378
	1379
	1380	//
	1381	// Reverse communication preparations
	1382	// I know it looks ugly, but it works the same way
	1383	// anywhere from C++ to Python.
	1384	//
	1385	// This code initializes locals by:
	1386	// * random values determined during code
	1387	// generation - on first subroutine call
	1388	// * values from previous call - on subsequent calls
	1389	//
	1390	if( state.rstate.stage>=0 )
	1391	{
	1392	i = state.rstate.ia[0];
	1393	j = state.rstate.ia[1];
	1394	ns = state.rstate.ia[2];
	1395	info = state.rstate.ia[3];
	1396	c1 = state.rstate.ra[0];
	1397	c2 = state.rstate.ra[1];
	1398	intg = state.rstate.ra[2];
	1399	intk = state.rstate.ra[3];
	1400	inta = state.rstate.ra[4];
	1401	v = state.rstate.ra[5];
	1402	ta = state.rstate.ra[6];
	1403	tb = state.rstate.ra[7];
	1404	qeps = state.rstate.ra[8];
	1405	}
	1406	else
	1407	{
	1408	i = 497;
	1409	j = -271;
	1410	ns = -581;
	1411	info = 745;
	1412	c1 = -533;
	1413	c2 = -77;
	1414	intg = 678;
	1415	intk = -293;
	1416	inta = 316;
	1417	v = 647;
	1418	ta = -756;
	1419	tb = 830;
	1420	qeps = -871;
	1421	}
	1422	if( state.rstate.stage==0 )
	1423	{
	1424	goto lbl_0;
	1425	}
	1426	if( state.rstate.stage==1 )
	1427	{
	1428	goto lbl_1;
	1429	}
	1430	if( state.rstate.stage==2 )
	1431	{
	1432	goto lbl_2;
	1433	}
	1434
	1435	//
	1436	// Routine body
	1437	//
	1438
	1439	//
	1440	// initialize quadratures.
	1441	// use 15-point Gauss-Kronrod formula.
	1442	//
	1443	state.n = 15;
	1444	gkq.gkqgenerategausslegendre(state.n, ref info, ref state.qn, ref state.wk, ref state.wg);
	1445	if( info<0 )
	1446	{
	1447	state.info = -5;
	1448	state.r = 0;
	1449	result = false;
	1450	return result;
	1451	}
	1452	state.wr = new double[state.n];
	1453	for(i=0; i<=state.n-1; i++)
	1454	{
	1455	if( i==0 )
	1456	{
	1457	state.wr[i] = 0.5*Math.Abs(state.qn[1]-state.qn[0]);
	1458	continue;
	1459	}
	1460	if( i==state.n-1 )
	1461	{
	1462	state.wr[state.n-1] = 0.5*Math.Abs(state.qn[state.n-1]-state.qn[state.n-2]);
	1463	continue;
	1464	}
	1465	state.wr[i] = 0.5*Math.Abs(state.qn[i-1]-state.qn[i+1]);
	1466	}
	1467
	1468	//
	1469	// special case
	1470	//
	1471	if( (double)(state.a)==(double)(state.b) )
	1472	{
	1473	state.info = 1;
	1474	state.r = 0;
	1475	result = false;
	1476	return result;
	1477	}
	1478
	1479	//
	1480	// test parameters
	1481	//
	1482	if( (double)(state.eps)<(double)(0) \| (double)(state.xwidth)<(double)(0) )
	1483	{
	1484	state.info = -1;
	1485	state.r = 0;
	1486	result = false;
	1487	return result;
	1488	}
	1489	state.info = 1;
	1490	if( (double)(state.eps)==(double)(0) )
	1491	{
	1492	state.eps = 1000*math.machineepsilon;
	1493	}
	1494
	1495	//
	1496	// First, prepare heap
	1497	// * column 0 - absolute error
	1498	// * column 1 - integral of a F(x) (calculated using Kronrod extension nodes)
	1499	// * column 2 - integral of a \|F(x)\| (calculated using modified rect. method)
	1500	// * column 3 - left boundary of a subinterval
	1501	// * column 4 - right boundary of a subinterval
	1502	//
	1503	if( (double)(state.xwidth)!=(double)(0) )
	1504	{
	1505	goto lbl_3;
	1506	}
	1507
	1508	//
	1509	// no maximum width requirements
	1510	// start from one big subinterval
	1511	//
	1512	state.heapwidth = 5;
	1513	state.heapsize = 1;
	1514	state.heapused = 1;
	1515	state.heap = new double[state.heapsize, state.heapwidth];
	1516	c1 = 0.5*(state.b-state.a);
	1517	c2 = 0.5*(state.b+state.a);
	1518	intg = 0;
	1519	intk = 0;
	1520	inta = 0;
	1521	i = 0;
	1522	lbl_5:
	1523	if( i>state.n-1 )
	1524	{
	1525	goto lbl_7;
	1526	}
	1527
	1528	//
	1529	// obtain F
	1530	//
	1531	state.x = c1*state.qn[i]+c2;
	1532	state.rstate.stage = 0;
	1533	goto lbl_rcomm;
	1534	lbl_0:
	1535	v = state.f;
	1536
	1537	//
	1538	// Gauss-Kronrod formula
	1539	//
	1540	intk = intk+v*state.wk[i];
	1541	if( i%2==1 )
	1542	{
	1543	intg = intg+v*state.wg[i];
	1544	}
	1545
	1546	//
	1547	// Integral \|F(x)\|
	1548	// Use rectangles method
	1549	//
	1550	inta = inta+Math.Abs(v)*state.wr[i];
	1551	i = i+1;
	1552	goto lbl_5;
	1553	lbl_7:
	1554	intk = intk(state.b-state.a)0.5;
	1555	intg = intg(state.b-state.a)0.5;
	1556	inta = inta(state.b-state.a)0.5;
	1557	state.heap[0,0] = Math.Abs(intg-intk);
	1558	state.heap[0,1] = intk;
	1559	state.heap[0,2] = inta;
	1560	state.heap[0,3] = state.a;
	1561	state.heap[0,4] = state.b;
	1562	state.sumerr = state.heap[0,0];
	1563	state.sumabs = Math.Abs(inta);
	1564	goto lbl_4;
	1565	lbl_3:
	1566
	1567	//
	1568	// maximum subinterval should be no more than XWidth.
	1569	// so we create Ceil((B-A)/XWidth)+1 small subintervals
	1570	//
	1571	ns = (int)Math.Ceiling(Math.Abs(state.b-state.a)/state.xwidth)+1;
	1572	state.heapsize = ns;
	1573	state.heapused = ns;
	1574	state.heapwidth = 5;
	1575	state.heap = new double[state.heapsize, state.heapwidth];
	1576	state.sumerr = 0;
	1577	state.sumabs = 0;
	1578	j = 0;
	1579	lbl_8:
	1580	if( j>ns-1 )
	1581	{
	1582	goto lbl_10;
	1583	}
	1584	ta = state.a+j*(state.b-state.a)/ns;
	1585	tb = state.a+(j+1)*(state.b-state.a)/ns;
	1586	c1 = 0.5*(tb-ta);
	1587	c2 = 0.5*(tb+ta);
	1588	intg = 0;
	1589	intk = 0;
	1590	inta = 0;
	1591	i = 0;
	1592	lbl_11:
	1593	if( i>state.n-1 )
	1594	{
	1595	goto lbl_13;
	1596	}
	1597
	1598	//
	1599	// obtain F
	1600	//
	1601	state.x = c1*state.qn[i]+c2;
	1602	state.rstate.stage = 1;
	1603	goto lbl_rcomm;
	1604	lbl_1:
	1605	v = state.f;
	1606
	1607	//
	1608	// Gauss-Kronrod formula
	1609	//
	1610	intk = intk+v*state.wk[i];
	1611	if( i%2==1 )
	1612	{
	1613	intg = intg+v*state.wg[i];
	1614	}
	1615
	1616	//
	1617	// Integral \|F(x)\|
	1618	// Use rectangles method
	1619	//
	1620	inta = inta+Math.Abs(v)*state.wr[i];
	1621	i = i+1;
	1622	goto lbl_11;
	1623	lbl_13:
	1624	intk = intk(tb-ta)0.5;
	1625	intg = intg(tb-ta)0.5;
	1626	inta = inta(tb-ta)0.5;
	1627	state.heap[j,0] = Math.Abs(intg-intk);
	1628	state.heap[j,1] = intk;
	1629	state.heap[j,2] = inta;
	1630	state.heap[j,3] = ta;
	1631	state.heap[j,4] = tb;
	1632	state.sumerr = state.sumerr+state.heap[j,0];
	1633	state.sumabs = state.sumabs+Math.Abs(inta);
	1634	j = j+1;
	1635	goto lbl_8;
	1636	lbl_10:
	1637	lbl_4:
	1638
	1639	//
	1640	// method iterations
	1641	//
	1642	lbl_14:
	1643	if( false )
	1644	{
	1645	goto lbl_15;
	1646	}
	1647
	1648	//
	1649	// additional memory if needed
	1650	//
	1651	if( state.heapused==state.heapsize )
	1652	{
	1653	mheapresize(ref state.heap, ref state.heapsize, 4*state.heapsize, state.heapwidth);
	1654	}
	1655
	1656	//
	1657	// TODO: every 20 iterations recalculate errors/sums
	1658	// TODO: one more criterion to prevent infinite loops with too strict Eps
	1659	//
	1660	if( (double)(state.sumerr)<=(double)(state.eps*state.sumabs) )
	1661	{
	1662	state.r = 0;
	1663	for(j=0; j<=state.heapused-1; j++)
	1664	{
	1665	state.r = state.r+state.heap[j,1];
	1666	}
	1667	result = false;
	1668	return result;
	1669	}
	1670
	1671	//
	1672	// Exclude interval with maximum absolute error
	1673	//
	1674	mheappop(ref state.heap, state.heapused, state.heapwidth);
	1675	state.sumerr = state.sumerr-state.heap[state.heapused-1,0];
	1676	state.sumabs = state.sumabs-state.heap[state.heapused-1,2];
	1677
	1678	//
	1679	// Divide interval, create subintervals
	1680	//
	1681	ta = state.heap[state.heapused-1,3];
	1682	tb = state.heap[state.heapused-1,4];
	1683	state.heap[state.heapused-1,3] = ta;
	1684	state.heap[state.heapused-1,4] = 0.5*(ta+tb);
	1685	state.heap[state.heapused,3] = 0.5*(ta+tb);
	1686	state.heap[state.heapused,4] = tb;
	1687	j = state.heapused-1;
	1688	lbl_16:
	1689	if( j>state.heapused )
	1690	{
	1691	goto lbl_18;
	1692	}
	1693	c1 = 0.5*(state.heap[j,4]-state.heap[j,3]);
	1694	c2 = 0.5*(state.heap[j,4]+state.heap[j,3]);
	1695	intg = 0;
	1696	intk = 0;
	1697	inta = 0;
	1698	i = 0;
	1699	lbl_19:
	1700	if( i>state.n-1 )
	1701	{
	1702	goto lbl_21;
	1703	}
	1704
	1705	//
	1706	// F(x)
	1707	//
	1708	state.x = c1*state.qn[i]+c2;
	1709	state.rstate.stage = 2;
	1710	goto lbl_rcomm;
	1711	lbl_2:
	1712	v = state.f;
	1713
	1714	//
	1715	// Gauss-Kronrod formula
	1716	//
	1717	intk = intk+v*state.wk[i];
	1718	if( i%2==1 )
	1719	{
	1720	intg = intg+v*state.wg[i];
	1721	}
	1722
	1723	//
	1724	// Integral \|F(x)\|
	1725	// Use rectangles method
	1726	//
	1727	inta = inta+Math.Abs(v)*state.wr[i];
	1728	i = i+1;
	1729	goto lbl_19;
	1730	lbl_21:
	1731	intk = intk(state.heap[j,4]-state.heap[j,3])0.5;
	1732	intg = intg(state.heap[j,4]-state.heap[j,3])0.5;
	1733	inta = inta(state.heap[j,4]-state.heap[j,3])0.5;
	1734	state.heap[j,0] = Math.Abs(intg-intk);
	1735	state.heap[j,1] = intk;
	1736	state.heap[j,2] = inta;
	1737	state.sumerr = state.sumerr+state.heap[j,0];
	1738	state.sumabs = state.sumabs+state.heap[j,2];
	1739	j = j+1;
	1740	goto lbl_16;
	1741	lbl_18:
	1742	mheappush(ref state.heap, state.heapused-1, state.heapwidth);
	1743	mheappush(ref state.heap, state.heapused, state.heapwidth);
	1744	state.heapused = state.heapused+1;
	1745	goto lbl_14;
	1746	lbl_15:
	1747	result = false;
	1748	return result;
	1749
	1750	//
	1751	// Saving state
	1752	//
	1753	lbl_rcomm:
	1754	result = true;
	1755	state.rstate.ia[0] = i;
	1756	state.rstate.ia[1] = j;
	1757	state.rstate.ia[2] = ns;
	1758	state.rstate.ia[3] = info;
	1759	state.rstate.ra[0] = c1;
	1760	state.rstate.ra[1] = c2;
	1761	state.rstate.ra[2] = intg;
	1762	state.rstate.ra[3] = intk;
	1763	state.rstate.ra[4] = inta;
	1764	state.rstate.ra[5] = v;
	1765	state.rstate.ra[6] = ta;
	1766	state.rstate.ra[7] = tb;
	1767	state.rstate.ra[8] = qeps;
	1768	return result;
	1769	}
	1770
	1771
	1772	private static void mheappop(ref double[,] heap,
	1773	int heapsize,
	1774	int heapwidth)
	1775	{
	1776	int i = 0;
	1777	int p = 0;
	1778	double t = 0;
	1779	int maxcp = 0;
	1780
	1781	if( heapsize==1 )
	1782	{
	1783	return;
	1784	}
	1785	for(i=0; i<=heapwidth-1; i++)
	1786	{
	1787	t = heap[heapsize-1,i];
	1788	heap[heapsize-1,i] = heap[0,i];
	1789	heap[0,i] = t;
	1790	}
	1791	p = 0;
	1792	while( 2*p+1<heapsize-1 )
	1793	{
	1794	maxcp = 2*p+1;
	1795	if( 2*p+2<heapsize-1 )
	1796	{
	1797	if( (double)(heap[2p+2,0])>(double)(heap[2p+1,0]) )
	1798	{
	1799	maxcp = 2*p+2;
	1800	}
	1801	}
	1802	if( (double)(heap[p,0])<(double)(heap[maxcp,0]) )
	1803	{
	1804	for(i=0; i<=heapwidth-1; i++)
	1805	{
	1806	t = heap[p,i];
	1807	heap[p,i] = heap[maxcp,i];
	1808	heap[maxcp,i] = t;
	1809	}
	1810	p = maxcp;
	1811	}
	1812	else
	1813	{
	1814	break;
	1815	}
	1816	}
	1817	}
	1818
	1819
	1820	private static void mheappush(ref double[,] heap,
	1821	int heapsize,
	1822	int heapwidth)
	1823	{
	1824	int i = 0;
	1825	int p = 0;
	1826	double t = 0;
	1827	int parent = 0;
	1828
	1829	if( heapsize==0 )
	1830	{
	1831	return;
	1832	}
	1833	p = heapsize;
	1834	while( p!=0 )
	1835	{
	1836	parent = (p-1)/2;
	1837	if( (double)(heap[p,0])>(double)(heap[parent,0]) )
	1838	{
	1839	for(i=0; i<=heapwidth-1; i++)
	1840	{
	1841	t = heap[p,i];
	1842	heap[p,i] = heap[parent,i];
	1843	heap[parent,i] = t;
	1844	}
	1845	p = parent;
	1846	}
	1847	else
	1848	{
	1849	break;
	1850	}
	1851	}
	1852	}
	1853
	1854
	1855	private static void mheapresize(ref double[,] heap,
	1856	ref int heapsize,
	1857	int newheapsize,
	1858	int heapwidth)
	1859	{
	1860	double[,] tmp = new double[0,0];
	1861	int i = 0;
	1862	int i_ = 0;
	1863
	1864	tmp = new double[heapsize, heapwidth];
	1865	for(i=0; i<=heapsize-1; i++)
	1866	{
	1867	for(i_=0; i_<=heapwidth-1;i_++)
	1868	{
	1869	tmp[i,i_] = heap[i,i_];
	1870	}
	1871	}
	1872	heap = new double[newheapsize, heapwidth];
	1873	for(i=0; i<=heapsize-1; i++)
	1874	{
	1875	for(i_=0; i_<=heapwidth-1;i_++)
	1876	{
	1877	heap[i,i_] = tmp[i,i_];
	1878	}
	1879	}
	1880	heapsize = newheapsize;
	1881	}
	1882
	1883
	1884	}
	1885	public class gq
	1886	{
	1887	/*************************************************************************
	1888	Computation of nodes and weights for a Gauss quadrature formula
	1889
	1890	The algorithm generates the N-point Gauss quadrature formula with weight
	1891	function given by coefficients alpha and beta of a recurrence relation
	1892	which generates a system of orthogonal polynomials:
	1893
	1894	P-1(x) = 0
	1895	P0(x) = 1
	1896	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	1897
	1898	and zeroth moment Mu0
	1899
	1900	Mu0 = integral(W(x)dx,a,b)
	1901
	1902	INPUT PARAMETERS:
	1903	Alpha array[0..N-1], alpha coefficients
	1904	Beta array[0..N-1], beta coefficients
	1905	Zero-indexed element is not used and may be arbitrary.
	1906	Beta[I]>0.
	1907	Mu0 zeroth moment of the weight function.
	1908	N number of nodes of the quadrature formula, N>=1
	1909
	1910	OUTPUT PARAMETERS:
	1911	Info - error code:
	1912	* -3 internal eigenproblem solver hasn't converged
	1913	* -2 Beta[i]<=0
	1914	* -1 incorrect N was passed
	1915	* 1 OK
	1916	X - array[0..N-1] - array of quadrature nodes,
	1917	in ascending order.
	1918	W - array[0..N-1] - array of quadrature weights.
	1919
	1920	-- ALGLIB --
	1921	Copyright 2005-2009 by Bochkanov Sergey
	1922	*************************************************************************/
	1923	public static void gqgeneraterec(double[] alpha,
	1924	double[] beta,
	1925	double mu0,
	1926	int n,
	1927	ref int info,
	1928	ref double[] x,
	1929	ref double[] w)
	1930	{
	1931	int i = 0;
	1932	double[] d = new double[0];
	1933	double[] e = new double[0];
	1934	double[,] z = new double[0,0];
	1935
	1936	info = 0;
	1937	x = new double[0];
	1938	w = new double[0];
	1939
	1940	if( n<1 )
	1941	{
	1942	info = -1;
	1943	return;
	1944	}
	1945	info = 1;
	1946
	1947	//
	1948	// Initialize
	1949	//
	1950	d = new double[n];
	1951	e = new double[n];
	1952	for(i=1; i<=n-1; i++)
	1953	{
	1954	d[i-1] = alpha[i-1];
	1955	if( (double)(beta[i])<=(double)(0) )
	1956	{
	1957	info = -2;
	1958	return;
	1959	}
	1960	e[i-1] = Math.Sqrt(beta[i]);
	1961	}
	1962	d[n-1] = alpha[n-1];
	1963
	1964	//
	1965	// EVD
	1966	//
	1967	if( !evd.smatrixtdevd(ref d, e, n, 3, ref z) )
	1968	{
	1969	info = -3;
	1970	return;
	1971	}
	1972
	1973	//
	1974	// Generate
	1975	//
	1976	x = new double[n];
	1977	w = new double[n];
	1978	for(i=1; i<=n; i++)
	1979	{
	1980	x[i-1] = d[i-1];
	1981	w[i-1] = mu0*math.sqr(z[0,i-1]);
	1982	}
	1983	}
	1984
	1985
	1986	/*************************************************************************
	1987	Computation of nodes and weights for a Gauss-Lobatto quadrature formula
	1988
	1989	The algorithm generates the N-point Gauss-Lobatto quadrature formula with
	1990	weight function given by coefficients alpha and beta of a recurrence which
	1991	generates a system of orthogonal polynomials.
	1992
	1993	P-1(x) = 0
	1994	P0(x) = 1
	1995	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	1996
	1997	and zeroth moment Mu0
	1998
	1999	Mu0 = integral(W(x)dx,a,b)
	2000
	2001	INPUT PARAMETERS:
	2002	Alpha array[0..N-2], alpha coefficients
	2003	Beta array[0..N-2], beta coefficients.
	2004	Zero-indexed element is not used, may be arbitrary.
	2005	Beta[I]>0
	2006	Mu0 zeroth moment of the weighting function.
	2007	A left boundary of the integration interval.
	2008	B right boundary of the integration interval.
	2009	N number of nodes of the quadrature formula, N>=3
	2010	(including the left and right boundary nodes).
	2011
	2012	OUTPUT PARAMETERS:
	2013	Info - error code:
	2014	* -3 internal eigenproblem solver hasn't converged
	2015	* -2 Beta[i]<=0
	2016	* -1 incorrect N was passed
	2017	* 1 OK
	2018	X - array[0..N-1] - array of quadrature nodes,
	2019	in ascending order.
	2020	W - array[0..N-1] - array of quadrature weights.
	2021
	2022	-- ALGLIB --
	2023	Copyright 2005-2009 by Bochkanov Sergey
	2024	*************************************************************************/
	2025	public static void gqgenerategausslobattorec(double[] alpha,
	2026	double[] beta,
	2027	double mu0,
	2028	double a,
	2029	double b,
	2030	int n,
	2031	ref int info,
	2032	ref double[] x,
	2033	ref double[] w)
	2034	{
	2035	int i = 0;
	2036	double[] d = new double[0];
	2037	double[] e = new double[0];
	2038	double[,] z = new double[0,0];
	2039	double pim1a = 0;
	2040	double pia = 0;
	2041	double pim1b = 0;
	2042	double pib = 0;
	2043	double t = 0;
	2044	double a11 = 0;
	2045	double a12 = 0;
	2046	double a21 = 0;
	2047	double a22 = 0;
	2048	double b1 = 0;
	2049	double b2 = 0;
	2050	double alph = 0;
	2051	double bet = 0;
	2052
	2053	alpha = (double[])alpha.Clone();
	2054	beta = (double[])beta.Clone();
	2055	info = 0;
	2056	x = new double[0];
	2057	w = new double[0];
	2058
	2059	if( n<=2 )
	2060	{
	2061	info = -1;
	2062	return;
	2063	}
	2064	info = 1;
	2065
	2066	//
	2067	// Initialize, D[1:N+1], E[1:N]
	2068	//
	2069	n = n-2;
	2070	d = new double[n+2];
	2071	e = new double[n+1];
	2072	for(i=1; i<=n+1; i++)
	2073	{
	2074	d[i-1] = alpha[i-1];
	2075	}
	2076	for(i=1; i<=n; i++)
	2077	{
	2078	if( (double)(beta[i])<=(double)(0) )
	2079	{
	2080	info = -2;
	2081	return;
	2082	}
	2083	e[i-1] = Math.Sqrt(beta[i]);
	2084	}
	2085
	2086	//
	2087	// Caclulate Pn(a), Pn+1(a), Pn(b), Pn+1(b)
	2088	//
	2089	beta[0] = 0;
	2090	pim1a = 0;
	2091	pia = 1;
	2092	pim1b = 0;
	2093	pib = 1;
	2094	for(i=1; i<=n+1; i++)
	2095	{
	2096
	2097	//
	2098	// Pi(a)
	2099	//
	2100	t = (a-alpha[i-1])pia-beta[i-1]pim1a;
	2101	pim1a = pia;
	2102	pia = t;
	2103
	2104	//
	2105	// Pi(b)
	2106	//
	2107	t = (b-alpha[i-1])pib-beta[i-1]pim1b;
	2108	pim1b = pib;
	2109	pib = t;
	2110	}
	2111
	2112	//
	2113	// Calculate alpha'(n+1), beta'(n+1)
	2114	//
	2115	a11 = pia;
	2116	a12 = pim1a;
	2117	a21 = pib;
	2118	a22 = pim1b;
	2119	b1 = a*pia;
	2120	b2 = b*pib;
	2121	if( (double)(Math.Abs(a11))>(double)(Math.Abs(a21)) )
	2122	{
	2123	a22 = a22-a12*a21/a11;
	2124	b2 = b2-b1*a21/a11;
	2125	bet = b2/a22;
	2126	alph = (b1-bet*a12)/a11;
	2127	}
	2128	else
	2129	{
	2130	a12 = a12-a22*a11/a21;
	2131	b1 = b1-b2*a11/a21;
	2132	bet = b1/a12;
	2133	alph = (b2-bet*a22)/a21;
	2134	}
	2135	if( (double)(bet)<(double)(0) )
	2136	{
	2137	info = -3;
	2138	return;
	2139	}
	2140	d[n+1] = alph;
	2141	e[n] = Math.Sqrt(bet);
	2142
	2143	//
	2144	// EVD
	2145	//
	2146	if( !evd.smatrixtdevd(ref d, e, n+2, 3, ref z) )
	2147	{
	2148	info = -3;
	2149	return;
	2150	}
	2151
	2152	//
	2153	// Generate
	2154	//
	2155	x = new double[n+2];
	2156	w = new double[n+2];
	2157	for(i=1; i<=n+2; i++)
	2158	{
	2159	x[i-1] = d[i-1];
	2160	w[i-1] = mu0*math.sqr(z[0,i-1]);
	2161	}
	2162	}
	2163
	2164
	2165	/*************************************************************************
	2166	Computation of nodes and weights for a Gauss-Radau quadrature formula
	2167
	2168	The algorithm generates the N-point Gauss-Radau quadrature formula with
	2169	weight function given by the coefficients alpha and beta of a recurrence
	2170	which generates a system of orthogonal polynomials.
	2171
	2172	P-1(x) = 0
	2173	P0(x) = 1
	2174	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	2175
	2176	and zeroth moment Mu0
	2177
	2178	Mu0 = integral(W(x)dx,a,b)
	2179
	2180	INPUT PARAMETERS:
	2181	Alpha array[0..N-2], alpha coefficients.
	2182	Beta array[0..N-1], beta coefficients
	2183	Zero-indexed element is not used.
	2184	Beta[I]>0
	2185	Mu0 zeroth moment of the weighting function.
	2186	A left boundary of the integration interval.
	2187	N number of nodes of the quadrature formula, N>=2
	2188	(including the left boundary node).
	2189
	2190	OUTPUT PARAMETERS:
	2191	Info - error code:
	2192	* -3 internal eigenproblem solver hasn't converged
	2193	* -2 Beta[i]<=0
	2194	* -1 incorrect N was passed
	2195	* 1 OK
	2196	X - array[0..N-1] - array of quadrature nodes,
	2197	in ascending order.
	2198	W - array[0..N-1] - array of quadrature weights.
	2199
	2200
	2201	-- ALGLIB --
	2202	Copyright 2005-2009 by Bochkanov Sergey
	2203	*************************************************************************/
	2204	public static void gqgenerategaussradaurec(double[] alpha,
	2205	double[] beta,
	2206	double mu0,
	2207	double a,
	2208	int n,
	2209	ref int info,
	2210	ref double[] x,
	2211	ref double[] w)
	2212	{
	2213	int i = 0;
	2214	double[] d = new double[0];
	2215	double[] e = new double[0];
	2216	double[,] z = new double[0,0];
	2217	double polim1 = 0;
	2218	double poli = 0;
	2219	double t = 0;
	2220
	2221	alpha = (double[])alpha.Clone();
	2222	beta = (double[])beta.Clone();
	2223	info = 0;
	2224	x = new double[0];
	2225	w = new double[0];
	2226
	2227	if( n<2 )
	2228	{
	2229	info = -1;
	2230	return;
	2231	}
	2232	info = 1;
	2233
	2234	//
	2235	// Initialize, D[1:N], E[1:N]
	2236	//
	2237	n = n-1;
	2238	d = new double[n+1];
	2239	e = new double[n];
	2240	for(i=1; i<=n; i++)
	2241	{
	2242	d[i-1] = alpha[i-1];
	2243	if( (double)(beta[i])<=(double)(0) )
	2244	{
	2245	info = -2;
	2246	return;
	2247	}
	2248	e[i-1] = Math.Sqrt(beta[i]);
	2249	}
	2250
	2251	//
	2252	// Caclulate Pn(a), Pn-1(a), and D[N+1]
	2253	//
	2254	beta[0] = 0;
	2255	polim1 = 0;
	2256	poli = 1;
	2257	for(i=1; i<=n; i++)
	2258	{
	2259	t = (a-alpha[i-1])poli-beta[i-1]polim1;
	2260	polim1 = poli;
	2261	poli = t;
	2262	}
	2263	d[n] = a-beta[n]*polim1/poli;
	2264
	2265	//
	2266	// EVD
	2267	//
	2268	if( !evd.smatrixtdevd(ref d, e, n+1, 3, ref z) )
	2269	{
	2270	info = -3;
	2271	return;
	2272	}
	2273
	2274	//
	2275	// Generate
	2276	//
	2277	x = new double[n+1];
	2278	w = new double[n+1];
	2279	for(i=1; i<=n+1; i++)
	2280	{
	2281	x[i-1] = d[i-1];
	2282	w[i-1] = mu0*math.sqr(z[0,i-1]);
	2283	}
	2284	}
	2285
	2286
	2287	/*************************************************************************
	2288	Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
	2289	nodes.
	2290
	2291	INPUT PARAMETERS:
	2292	N - number of nodes, >=1
	2293
	2294	OUTPUT PARAMETERS:
	2295	Info - error code:
	2296	* -4 an error was detected when calculating
	2297	weights/nodes. N is too large to obtain
	2298	weights/nodes with high enough accuracy.
	2299	Try to use multiple precision version.
	2300	* -3 internal eigenproblem solver hasn't converged
	2301	* -1 incorrect N was passed
	2302	* +1 OK
	2303	X - array[0..N-1] - array of quadrature nodes,
	2304	in ascending order.
	2305	W - array[0..N-1] - array of quadrature weights.
	2306
	2307
	2308	-- ALGLIB --
	2309	Copyright 12.05.2009 by Bochkanov Sergey
	2310	*************************************************************************/
	2311	public static void gqgenerategausslegendre(int n,
	2312	ref int info,
	2313	ref double[] x,
	2314	ref double[] w)
	2315	{
	2316	double[] alpha = new double[0];
	2317	double[] beta = new double[0];
	2318	int i = 0;
	2319
	2320	info = 0;
	2321	x = new double[0];
	2322	w = new double[0];
	2323
	2324	if( n<1 )
	2325	{
	2326	info = -1;
	2327	return;
	2328	}
	2329	alpha = new double[n];
	2330	beta = new double[n];
	2331	for(i=0; i<=n-1; i++)
	2332	{
	2333	alpha[i] = 0;
	2334	}
	2335	beta[0] = 2;
	2336	for(i=1; i<=n-1; i++)
	2337	{
	2338	beta[i] = 1/(4-1/math.sqr(i));
	2339	}
	2340	gqgeneraterec(alpha, beta, beta[0], n, ref info, ref x, ref w);
	2341
	2342	//
	2343	// test basic properties to detect errors
	2344	//
	2345	if( info>0 )
	2346	{
	2347	if( (double)(x[0])<(double)(-1) \| (double)(x[n-1])>(double)(1) )
	2348	{
	2349	info = -4;
	2350	}
	2351	for(i=0; i<=n-2; i++)
	2352	{
	2353	if( (double)(x[i])>=(double)(x[i+1]) )
	2354	{
	2355	info = -4;
	2356	}
	2357	}
	2358	}
	2359	}
	2360
	2361
	2362	/*************************************************************************
	2363	Returns nodes/weights for Gauss-Jacobi quadrature on [-1,1] with weight
	2364	function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
	2365
	2366	INPUT PARAMETERS:
	2367	N - number of nodes, >=1
	2368	Alpha - power-law coefficient, Alpha>-1
	2369	Beta - power-law coefficient, Beta>-1
	2370
	2371	OUTPUT PARAMETERS:
	2372	Info - error code:
	2373	* -4 an error was detected when calculating
	2374	weights/nodes. Alpha or Beta are too close
	2375	to -1 to obtain weights/nodes with high enough
	2376	accuracy, or, may be, N is too large. Try to
	2377	use multiple precision version.
	2378	* -3 internal eigenproblem solver hasn't converged
	2379	* -1 incorrect N/Alpha/Beta was passed
	2380	* +1 OK
	2381	X - array[0..N-1] - array of quadrature nodes,
	2382	in ascending order.
	2383	W - array[0..N-1] - array of quadrature weights.
	2384
	2385
	2386	-- ALGLIB --
	2387	Copyright 12.05.2009 by Bochkanov Sergey
	2388	*************************************************************************/
	2389	public static void gqgenerategaussjacobi(int n,
	2390	double alpha,
	2391	double beta,
	2392	ref int info,
	2393	ref double[] x,
	2394	ref double[] w)
	2395	{
	2396	double[] a = new double[0];
	2397	double[] b = new double[0];
	2398	double alpha2 = 0;
	2399	double beta2 = 0;
	2400	double apb = 0;
	2401	double t = 0;
	2402	int i = 0;
	2403	double s = 0;
	2404
	2405	info = 0;
	2406	x = new double[0];
	2407	w = new double[0];
	2408
	2409	if( (n<1 \| (double)(alpha)<=(double)(-1)) \| (double)(beta)<=(double)(-1) )
	2410	{
	2411	info = -1;
	2412	return;
	2413	}
	2414	a = new double[n];
	2415	b = new double[n];
	2416	apb = alpha+beta;
	2417	a[0] = (beta-alpha)/(apb+2);
	2418	t = (apb+1)*Math.Log(2)+gammafunc.lngamma(alpha+1, ref s)+gammafunc.lngamma(beta+1, ref s)-gammafunc.lngamma(apb+2, ref s);
	2419	if( (double)(t)>(double)(Math.Log(math.maxrealnumber)) )
	2420	{
	2421	info = -4;
	2422	return;
	2423	}
	2424	b[0] = Math.Exp(t);
	2425	if( n>1 )
	2426	{
	2427	alpha2 = math.sqr(alpha);
	2428	beta2 = math.sqr(beta);
	2429	a[1] = (beta2-alpha2)/((apb+2)*(apb+4));
	2430	b[1] = 4(alpha+1)(beta+1)/((apb+3)*math.sqr(apb+2));
	2431	for(i=2; i<=n-1; i++)
	2432	{
	2433	a[i] = 0.25(beta2-alpha2)/(ii(1+0.5apb/i)(1+0.5(apb+2)/i));
	2434	b[i] = 0.25(1+alpha/i)(1+beta/i)(1+apb/i)/((1+0.5(apb+1)/i)(1+0.5(apb-1)/i)math.sqr(1+0.5apb/i));
	2435	}
	2436	}
	2437	gqgeneraterec(a, b, b[0], n, ref info, ref x, ref w);
	2438
	2439	//
	2440	// test basic properties to detect errors
	2441	//
	2442	if( info>0 )
	2443	{
	2444	if( (double)(x[0])<(double)(-1) \| (double)(x[n-1])>(double)(1) )
	2445	{
	2446	info = -4;
	2447	}
	2448	for(i=0; i<=n-2; i++)
	2449	{
	2450	if( (double)(x[i])>=(double)(x[i+1]) )
	2451	{
	2452	info = -4;
	2453	}
	2454	}
	2455	}
	2456	}
	2457
	2458
	2459	/*************************************************************************
	2460	Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf) with
	2461	weight function W(x)=Power(x,Alpha)*Exp(-x)
	2462
	2463	INPUT PARAMETERS:
	2464	N - number of nodes, >=1
	2465	Alpha - power-law coefficient, Alpha>-1
	2466
	2467	OUTPUT PARAMETERS:
	2468	Info - error code:
	2469	* -4 an error was detected when calculating
	2470	weights/nodes. Alpha is too close to -1 to
	2471	obtain weights/nodes with high enough accuracy
	2472	or, may be, N is too large. Try to use
	2473	multiple precision version.
	2474	* -3 internal eigenproblem solver hasn't converged
	2475	* -1 incorrect N/Alpha was passed
	2476	* +1 OK
	2477	X - array[0..N-1] - array of quadrature nodes,
	2478	in ascending order.
	2479	W - array[0..N-1] - array of quadrature weights.
	2480
	2481
	2482	-- ALGLIB --
	2483	Copyright 12.05.2009 by Bochkanov Sergey
	2484	*************************************************************************/
	2485	public static void gqgenerategausslaguerre(int n,
	2486	double alpha,
	2487	ref int info,
	2488	ref double[] x,
	2489	ref double[] w)
	2490	{
	2491	double[] a = new double[0];
	2492	double[] b = new double[0];
	2493	double t = 0;
	2494	int i = 0;
	2495	double s = 0;
	2496
	2497	info = 0;
	2498	x = new double[0];
	2499	w = new double[0];
	2500
	2501	if( n<1 \| (double)(alpha)<=(double)(-1) )
	2502	{
	2503	info = -1;
	2504	return;
	2505	}
	2506	a = new double[n];
	2507	b = new double[n];
	2508	a[0] = alpha+1;
	2509	t = gammafunc.lngamma(alpha+1, ref s);
	2510	if( (double)(t)>=(double)(Math.Log(math.maxrealnumber)) )
	2511	{
	2512	info = -4;
	2513	return;
	2514	}
	2515	b[0] = Math.Exp(t);
	2516	if( n>1 )
	2517	{
	2518	for(i=1; i<=n-1; i++)
	2519	{
	2520	a[i] = 2*i+alpha+1;
	2521	b[i] = i*(i+alpha);
	2522	}
	2523	}
	2524	gqgeneraterec(a, b, b[0], n, ref info, ref x, ref w);
	2525
	2526	//
	2527	// test basic properties to detect errors
	2528	//
	2529	if( info>0 )
	2530	{
	2531	if( (double)(x[0])<(double)(0) )
	2532	{
	2533	info = -4;
	2534	}
	2535	for(i=0; i<=n-2; i++)
	2536	{
	2537	if( (double)(x[i])>=(double)(x[i+1]) )
	2538	{
	2539	info = -4;
	2540	}
	2541	}
	2542	}
	2543	}
	2544
	2545
	2546	/*************************************************************************
	2547	Returns nodes/weights for Gauss-Hermite quadrature on (-inf,+inf) with
	2548	weight function W(x)=Exp(-x*x)
	2549
	2550	INPUT PARAMETERS:
	2551	N - number of nodes, >=1
	2552
	2553	OUTPUT PARAMETERS:
	2554	Info - error code:
	2555	* -4 an error was detected when calculating
	2556	weights/nodes. May be, N is too large. Try to
	2557	use multiple precision version.
	2558	* -3 internal eigenproblem solver hasn't converged
	2559	* -1 incorrect N/Alpha was passed
	2560	* +1 OK
	2561	X - array[0..N-1] - array of quadrature nodes,
	2562	in ascending order.
	2563	W - array[0..N-1] - array of quadrature weights.
	2564
	2565
	2566	-- ALGLIB --
	2567	Copyright 12.05.2009 by Bochkanov Sergey
	2568	*************************************************************************/
	2569	public static void gqgenerategausshermite(int n,
	2570	ref int info,
	2571	ref double[] x,
	2572	ref double[] w)
	2573	{
	2574	double[] a = new double[0];
	2575	double[] b = new double[0];
	2576	int i = 0;
	2577
	2578	info = 0;
	2579	x = new double[0];
	2580	w = new double[0];
	2581
	2582	if( n<1 )
	2583	{
	2584	info = -1;
	2585	return;
	2586	}
	2587	a = new double[n];
	2588	b = new double[n];
	2589	for(i=0; i<=n-1; i++)
	2590	{
	2591	a[i] = 0;
	2592	}
	2593	b[0] = Math.Sqrt(4*Math.Atan(1));
	2594	if( n>1 )
	2595	{
	2596	for(i=1; i<=n-1; i++)
	2597	{
	2598	b[i] = 0.5*i;
	2599	}
	2600	}
	2601	gqgeneraterec(a, b, b[0], n, ref info, ref x, ref w);
	2602
	2603	//
	2604	// test basic properties to detect errors
	2605	//
	2606	if( info>0 )
	2607	{
	2608	for(i=0; i<=n-2; i++)
	2609	{
	2610	if( (double)(x[i])>=(double)(x[i+1]) )
	2611	{
	2612	info = -4;
	2613	}
	2614	}
	2615	}
	2616	}
	2617
	2618
	2619	}
	2620	public class gkq
	2621	{
	2622	/*************************************************************************
	2623	Computation of nodes and weights of a Gauss-Kronrod quadrature formula
	2624
	2625	The algorithm generates the N-point Gauss-Kronrod quadrature formula with
	2626	weight function given by coefficients alpha and beta of a recurrence
	2627	relation which generates a system of orthogonal polynomials:
	2628
	2629	P-1(x) = 0
	2630	P0(x) = 1
	2631	Pn+1(x) = (x-alpha(n))Pn(x) - beta(n)Pn-1(x)
	2632
	2633	and zero moment Mu0
	2634
	2635	Mu0 = integral(W(x)dx,a,b)
	2636
	2637
	2638	INPUT PARAMETERS:
	2639	Alpha alpha coefficients, array[0..floor(3*K/2)].
	2640	Beta beta coefficients, array[0..ceil(3*K/2)].
	2641	Beta[0] is not used and may be arbitrary.
	2642	Beta[I]>0.
	2643	Mu0 zeroth moment of the weight function.
	2644	N number of nodes of the Gauss-Kronrod quadrature formula,
	2645	N >= 3,
	2646	N = 2*K+1.
	2647
	2648	OUTPUT PARAMETERS:
	2649	Info - error code:
	2650	* -5 no real and positive Gauss-Kronrod formula can
	2651	be created for such a weight function with a
	2652	given number of nodes.
	2653	* -4 N is too large, task may be ill conditioned -
	2654	x[i]=x[i+1] found.
	2655	* -3 internal eigenproblem solver hasn't converged
	2656	* -2 Beta[i]<=0
	2657	* -1 incorrect N was passed
	2658	* +1 OK
	2659	X - array[0..N-1] - array of quadrature nodes,
	2660	in ascending order.
	2661	WKronrod - array[0..N-1] - Kronrod weights
	2662	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	2663	corresponding to extended Kronrod nodes).
	2664
	2665	-- ALGLIB --
	2666	Copyright 08.05.2009 by Bochkanov Sergey
	2667	*************************************************************************/
	2668	public static void gkqgeneraterec(double[] alpha,
	2669	double[] beta,
	2670	double mu0,
	2671	int n,
	2672	ref int info,
	2673	ref double[] x,
	2674	ref double[] wkronrod,
	2675	ref double[] wgauss)
	2676	{
	2677	double[] ta = new double[0];
	2678	int i = 0;
	2679	int j = 0;
	2680	double[] t = new double[0];
	2681	double[] s = new double[0];
	2682	int wlen = 0;
	2683	int woffs = 0;
	2684	double u = 0;
	2685	int m = 0;
	2686	int l = 0;
	2687	int k = 0;
	2688	double[] xgtmp = new double[0];
	2689	double[] wgtmp = new double[0];
	2690	int i_ = 0;
	2691
	2692	alpha = (double[])alpha.Clone();
	2693	beta = (double[])beta.Clone();
	2694	info = 0;
	2695	x = new double[0];
	2696	wkronrod = new double[0];
	2697	wgauss = new double[0];
	2698
	2699	if( n%2!=1 \| n<3 )
	2700	{
	2701	info = -1;
	2702	return;
	2703	}
	2704	for(i=0; i<=(int)Math.Ceiling((double)(3*(n/2))/(double)2); i++)
	2705	{
	2706	if( (double)(beta[i])<=(double)(0) )
	2707	{
	2708	info = -2;
	2709	return;
	2710	}
	2711	}
	2712	info = 1;
	2713
	2714	//
	2715	// from external conventions about N/Beta/Mu0 to internal
	2716	//
	2717	n = n/2;
	2718	beta[0] = mu0;
	2719
	2720	//
	2721	// Calculate Gauss nodes/weights, save them for later processing
	2722	//
	2723	gq.gqgeneraterec(alpha, beta, mu0, n, ref info, ref xgtmp, ref wgtmp);
	2724	if( info<0 )
	2725	{
	2726	return;
	2727	}
	2728
	2729	//
	2730	// Resize:
	2731	// * A from 0..floor(3n/2) to 0..2n
	2732	// * B from 0..ceil(3n/2) to 0..2n
	2733	//
	2734	ta = new double[(int)Math.Floor((double)(3*n)/(double)2)+1];
	2735	for(i_=0; i_<=(int)Math.Floor((double)(3*n)/(double)2);i_++)
	2736	{
	2737	ta[i_] = alpha[i_];
	2738	}
	2739	alpha = new double[2*n+1];
	2740	for(i_=0; i_<=(int)Math.Floor((double)(3*n)/(double)2);i_++)
	2741	{
	2742	alpha[i_] = ta[i_];
	2743	}
	2744	for(i=(int)Math.Floor((double)(3n)/(double)2)+1; i<=2n; i++)
	2745	{
	2746	alpha[i] = 0;
	2747	}
	2748	ta = new double[(int)Math.Ceiling((double)(3*n)/(double)2)+1];
	2749	for(i_=0; i_<=(int)Math.Ceiling((double)(3*n)/(double)2);i_++)
	2750	{
	2751	ta[i_] = beta[i_];
	2752	}
	2753	beta = new double[2*n+1];
	2754	for(i_=0; i_<=(int)Math.Ceiling((double)(3*n)/(double)2);i_++)
	2755	{
	2756	beta[i_] = ta[i_];
	2757	}
	2758	for(i=(int)Math.Ceiling((double)(3n)/(double)2)+1; i<=2n; i++)
	2759	{
	2760	beta[i] = 0;
	2761	}
	2762
	2763	//
	2764	// Initialize T, S
	2765	//
	2766	wlen = 2+n/2;
	2767	t = new double[wlen];
	2768	s = new double[wlen];
	2769	ta = new double[wlen];
	2770	woffs = 1;
	2771	for(i=0; i<=wlen-1; i++)
	2772	{
	2773	t[i] = 0;
	2774	s[i] = 0;
	2775	}
	2776
	2777	//
	2778	// Algorithm from Dirk P. Laurie, "Calculation of Gauss-Kronrod quadrature rules", 1997.
	2779	//
	2780	t[woffs+0] = beta[n+1];
	2781	for(m=0; m<=n-2; m++)
	2782	{
	2783	u = 0;
	2784	for(k=(m+1)/2; k>=0; k--)
	2785	{
	2786	l = m-k;
	2787	u = u+(alpha[k+n+1]-alpha[l])t[woffs+k]+beta[k+n+1]s[woffs+k-1]-beta[l]*s[woffs+k];
	2788	s[woffs+k] = u;
	2789	}
	2790	for(i_=0; i_<=wlen-1;i_++)
	2791	{
	2792	ta[i_] = t[i_];
	2793	}
	2794	for(i_=0; i_<=wlen-1;i_++)
	2795	{
	2796	t[i_] = s[i_];
	2797	}
	2798	for(i_=0; i_<=wlen-1;i_++)
	2799	{
	2800	s[i_] = ta[i_];
	2801	}
	2802	}
	2803	for(j=n/2; j>=0; j--)
	2804	{
	2805	s[woffs+j] = s[woffs+j-1];
	2806	}
	2807	for(m=n-1; m<=2*n-3; m++)
	2808	{
	2809	u = 0;
	2810	for(k=m+1-n; k<=(m-1)/2; k++)
	2811	{
	2812	l = m-k;
	2813	j = n-1-l;
	2814	u = u-(alpha[k+n+1]-alpha[l])t[woffs+j]-beta[k+n+1]s[woffs+j]+beta[l]*s[woffs+j+1];
	2815	s[woffs+j] = u;
	2816	}
	2817	if( m%2==0 )
	2818	{
	2819	k = m/2;
	2820	alpha[k+n+1] = alpha[k]+(s[woffs+j]-beta[k+n+1]*s[woffs+j+1])/t[woffs+j+1];
	2821	}
	2822	else
	2823	{
	2824	k = (m+1)/2;
	2825	beta[k+n+1] = s[woffs+j]/s[woffs+j+1];
	2826	}
	2827	for(i_=0; i_<=wlen-1;i_++)
	2828	{
	2829	ta[i_] = t[i_];
	2830	}
	2831	for(i_=0; i_<=wlen-1;i_++)
	2832	{
	2833	t[i_] = s[i_];
	2834	}
	2835	for(i_=0; i_<=wlen-1;i_++)
	2836	{
	2837	s[i_] = ta[i_];
	2838	}
	2839	}
	2840	alpha[2n] = alpha[n-1]-beta[2n]*s[woffs+0]/t[woffs+0];
	2841
	2842	//
	2843	// calculation of Kronrod nodes and weights, unpacking of Gauss weights
	2844	//
	2845	gq.gqgeneraterec(alpha, beta, mu0, 2*n+1, ref info, ref x, ref wkronrod);
	2846	if( info==-2 )
	2847	{
	2848	info = -5;
	2849	}
	2850	if( info<0 )
	2851	{
	2852	return;
	2853	}
	2854	for(i=0; i<=2*n-1; i++)
	2855	{
	2856	if( (double)(x[i])>=(double)(x[i+1]) )
	2857	{
	2858	info = -4;
	2859	}
	2860	}
	2861	if( info<0 )
	2862	{
	2863	return;
	2864	}
	2865	wgauss = new double[2*n+1];
	2866	for(i=0; i<=2*n; i++)
	2867	{
	2868	wgauss[i] = 0;
	2869	}
	2870	for(i=0; i<=n-1; i++)
	2871	{
	2872	wgauss[2*i+1] = wgtmp[i];
	2873	}
	2874	}
	2875
	2876
	2877	/*************************************************************************
	2878	Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre
	2879	quadrature with N points.
	2880
	2881	GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed table) is
	2882	used depending on machine precision and number of nodes.
	2883
	2884	INPUT PARAMETERS:
	2885	N - number of Kronrod nodes, must be odd number, >=3.
	2886
	2887	OUTPUT PARAMETERS:
	2888	Info - error code:
	2889	* -4 an error was detected when calculating
	2890	weights/nodes. N is too large to obtain
	2891	weights/nodes with high enough accuracy.
	2892	Try to use multiple precision version.
	2893	* -3 internal eigenproblem solver hasn't converged
	2894	* -1 incorrect N was passed
	2895	* +1 OK
	2896	X - array[0..N-1] - array of quadrature nodes, ordered in
	2897	ascending order.
	2898	WKronrod - array[0..N-1] - Kronrod weights
	2899	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	2900	corresponding to extended Kronrod nodes).
	2901
	2902
	2903	-- ALGLIB --
	2904	Copyright 12.05.2009 by Bochkanov Sergey
	2905	*************************************************************************/
	2906	public static void gkqgenerategausslegendre(int n,
	2907	ref int info,
	2908	ref double[] x,
	2909	ref double[] wkronrod,
	2910	ref double[] wgauss)
	2911	{
	2912	double eps = 0;
	2913
	2914	info = 0;
	2915	x = new double[0];
	2916	wkronrod = new double[0];
	2917	wgauss = new double[0];
	2918
	2919	if( (double)(math.machineepsilon)>(double)(1.0E-32) & (((((n==15 \| n==21) \| n==31) \| n==41) \| n==51) \| n==61) )
	2920	{
	2921	info = 1;
	2922	gkqlegendretbl(n, ref x, ref wkronrod, ref wgauss, ref eps);
	2923	}
	2924	else
	2925	{
	2926	gkqlegendrecalc(n, ref info, ref x, ref wkronrod, ref wgauss);
	2927	}
	2928	}
	2929
	2930
	2931	/*************************************************************************
	2932	Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi
	2933	quadrature on [-1,1] with weight function
	2934
	2935	W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
	2936
	2937	INPUT PARAMETERS:
	2938	N - number of Kronrod nodes, must be odd number, >=3.
	2939	Alpha - power-law coefficient, Alpha>-1
	2940	Beta - power-law coefficient, Beta>-1
	2941
	2942	OUTPUT PARAMETERS:
	2943	Info - error code:
	2944	* -5 no real and positive Gauss-Kronrod formula can
	2945	be created for such a weight function with a
	2946	given number of nodes.
	2947	* -4 an error was detected when calculating
	2948	weights/nodes. Alpha or Beta are too close
	2949	to -1 to obtain weights/nodes with high enough
	2950	accuracy, or, may be, N is too large. Try to
	2951	use multiple precision version.
	2952	* -3 internal eigenproblem solver hasn't converged
	2953	* -1 incorrect N was passed
	2954	* +1 OK
	2955	* +2 OK, but quadrature rule have exterior nodes,
	2956	x[0]<-1 or x[n-1]>+1
	2957	X - array[0..N-1] - array of quadrature nodes, ordered in
	2958	ascending order.
	2959	WKronrod - array[0..N-1] - Kronrod weights
	2960	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	2961	corresponding to extended Kronrod nodes).
	2962
	2963
	2964	-- ALGLIB --
	2965	Copyright 12.05.2009 by Bochkanov Sergey
	2966	*************************************************************************/
	2967	public static void gkqgenerategaussjacobi(int n,
	2968	double alpha,
	2969	double beta,
	2970	ref int info,
	2971	ref double[] x,
	2972	ref double[] wkronrod,
	2973	ref double[] wgauss)
	2974	{
	2975	int clen = 0;
	2976	double[] a = new double[0];
	2977	double[] b = new double[0];
	2978	double alpha2 = 0;
	2979	double beta2 = 0;
	2980	double apb = 0;
	2981	double t = 0;
	2982	int i = 0;
	2983	double s = 0;
	2984
	2985	info = 0;
	2986	x = new double[0];
	2987	wkronrod = new double[0];
	2988	wgauss = new double[0];
	2989
	2990	if( n%2!=1 \| n<3 )
	2991	{
	2992	info = -1;
	2993	return;
	2994	}
	2995	if( (double)(alpha)<=(double)(-1) \| (double)(beta)<=(double)(-1) )
	2996	{
	2997	info = -1;
	2998	return;
	2999	}
	3000	clen = (int)Math.Ceiling((double)(3*(n/2))/(double)2)+1;
	3001	a = new double[clen];
	3002	b = new double[clen];
	3003	for(i=0; i<=clen-1; i++)
	3004	{
	3005	a[i] = 0;
	3006	}
	3007	apb = alpha+beta;
	3008	a[0] = (beta-alpha)/(apb+2);
	3009	t = (apb+1)*Math.Log(2)+gammafunc.lngamma(alpha+1, ref s)+gammafunc.lngamma(beta+1, ref s)-gammafunc.lngamma(apb+2, ref s);
	3010	if( (double)(t)>(double)(Math.Log(math.maxrealnumber)) )
	3011	{
	3012	info = -4;
	3013	return;
	3014	}
	3015	b[0] = Math.Exp(t);
	3016	if( clen>1 )
	3017	{
	3018	alpha2 = math.sqr(alpha);
	3019	beta2 = math.sqr(beta);
	3020	a[1] = (beta2-alpha2)/((apb+2)*(apb+4));
	3021	b[1] = 4(alpha+1)(beta+1)/((apb+3)*math.sqr(apb+2));
	3022	for(i=2; i<=clen-1; i++)
	3023	{
	3024	a[i] = 0.25(beta2-alpha2)/(ii(1+0.5apb/i)(1+0.5(apb+2)/i));
	3025	b[i] = 0.25(1+alpha/i)(1+beta/i)(1+apb/i)/((1+0.5(apb+1)/i)(1+0.5(apb-1)/i)math.sqr(1+0.5apb/i));
	3026	}
	3027	}
	3028	gkqgeneraterec(a, b, b[0], n, ref info, ref x, ref wkronrod, ref wgauss);
	3029
	3030	//
	3031	// test basic properties to detect errors
	3032	//
	3033	if( info>0 )
	3034	{
	3035	if( (double)(x[0])<(double)(-1) \| (double)(x[n-1])>(double)(1) )
	3036	{
	3037	info = 2;
	3038	}
	3039	for(i=0; i<=n-2; i++)
	3040	{
	3041	if( (double)(x[i])>=(double)(x[i+1]) )
	3042	{
	3043	info = -4;
	3044	}
	3045	}
	3046	}
	3047	}
	3048
	3049
	3050	/*************************************************************************
	3051	Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
	3052
	3053	Reduction to tridiagonal eigenproblem is used.
	3054
	3055	INPUT PARAMETERS:
	3056	N - number of Kronrod nodes, must be odd number, >=3.
	3057
	3058	OUTPUT PARAMETERS:
	3059	Info - error code:
	3060	* -4 an error was detected when calculating
	3061	weights/nodes. N is too large to obtain
	3062	weights/nodes with high enough accuracy.
	3063	Try to use multiple precision version.
	3064	* -3 internal eigenproblem solver hasn't converged
	3065	* -1 incorrect N was passed
	3066	* +1 OK
	3067	X - array[0..N-1] - array of quadrature nodes, ordered in
	3068	ascending order.
	3069	WKronrod - array[0..N-1] - Kronrod weights
	3070	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	3071	corresponding to extended Kronrod nodes).
	3072
	3073	-- ALGLIB --
	3074	Copyright 12.05.2009 by Bochkanov Sergey
	3075	*************************************************************************/
	3076	public static void gkqlegendrecalc(int n,
	3077	ref int info,
	3078	ref double[] x,
	3079	ref double[] wkronrod,
	3080	ref double[] wgauss)
	3081	{
	3082	double[] alpha = new double[0];
	3083	double[] beta = new double[0];
	3084	int alen = 0;
	3085	int blen = 0;
	3086	double mu0 = 0;
	3087	int k = 0;
	3088	int i = 0;
	3089
	3090	info = 0;
	3091	x = new double[0];
	3092	wkronrod = new double[0];
	3093	wgauss = new double[0];
	3094
	3095	if( n%2!=1 \| n<3 )
	3096	{
	3097	info = -1;
	3098	return;
	3099	}
	3100	mu0 = 2;
	3101	alen = (int)Math.Floor((double)(3*(n/2))/(double)2)+1;
	3102	blen = (int)Math.Ceiling((double)(3*(n/2))/(double)2)+1;
	3103	alpha = new double[alen];
	3104	beta = new double[blen];
	3105	for(k=0; k<=alen-1; k++)
	3106	{
	3107	alpha[k] = 0;
	3108	}
	3109	beta[0] = 2;
	3110	for(k=1; k<=blen-1; k++)
	3111	{
	3112	beta[k] = 1/(4-1/math.sqr(k));
	3113	}
	3114	gkqgeneraterec(alpha, beta, mu0, n, ref info, ref x, ref wkronrod, ref wgauss);
	3115
	3116	//
	3117	// test basic properties to detect errors
	3118	//
	3119	if( info>0 )
	3120	{
	3121	if( (double)(x[0])<(double)(-1) \| (double)(x[n-1])>(double)(1) )
	3122	{
	3123	info = -4;
	3124	}
	3125	for(i=0; i<=n-2; i++)
	3126	{
	3127	if( (double)(x[i])>=(double)(x[i+1]) )
	3128	{
	3129	info = -4;
	3130	}
	3131	}
	3132	}
	3133	}
	3134
	3135
	3136	/*************************************************************************
	3137	Returns Gauss and Gauss-Kronrod nodes for quadrature with N points using
	3138	pre-calculated table. Nodes/weights were computed with accuracy up to
	3139	1.0E-32 (if MPFR version of ALGLIB is used). In standard double precision
	3140	accuracy reduces to something about 2.0E-16 (depending on your compiler's
	3141	handling of long floating point constants).
	3142
	3143	INPUT PARAMETERS:
	3144	N - number of Kronrod nodes.
	3145	N can be 15, 21, 31, 41, 51, 61.
	3146
	3147	OUTPUT PARAMETERS:
	3148	X - array[0..N-1] - array of quadrature nodes, ordered in
	3149	ascending order.
	3150	WKronrod - array[0..N-1] - Kronrod weights
	3151	WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
	3152	corresponding to extended Kronrod nodes).
	3153
	3154
	3155	-- ALGLIB --
	3156	Copyright 12.05.2009 by Bochkanov Sergey
	3157	*************************************************************************/
	3158	public static void gkqlegendretbl(int n,
	3159	ref double[] x,
	3160	ref double[] wkronrod,
	3161	ref double[] wgauss,
	3162	ref double eps)
	3163	{
	3164	int i = 0;
	3165	int ng = 0;
	3166	int[] p1 = new int[0];
	3167	int[] p2 = new int[0];
	3168	double tmp = 0;
	3169
	3170	x = new double[0];
	3171	wkronrod = new double[0];
	3172	wgauss = new double[0];
	3173	eps = 0;
	3174
	3175
	3176	//
	3177	// these initializers are not really necessary,
	3178	// but without them compiler complains about uninitialized locals
	3179	//
	3180	ng = 0;
	3181
	3182	//
	3183	// Process
	3184	//
	3185	ap.assert(((((n==15 \| n==21) \| n==31) \| n==41) \| n==51) \| n==61, "GKQNodesTbl: incorrect N!");
	3186	x = new double[n];
	3187	wkronrod = new double[n];
	3188	wgauss = new double[n];
	3189	for(i=0; i<=n-1; i++)
	3190	{
	3191	x[i] = 0;
	3192	wkronrod[i] = 0;
	3193	wgauss[i] = 0;
	3194	}
	3195	eps = Math.Max(math.machineepsilon, 1.0E-32);
	3196	if( n==15 )
	3197	{
	3198	ng = 4;
	3199	wgauss[0] = 0.129484966168869693270611432679082;
	3200	wgauss[1] = 0.279705391489276667901467771423780;
	3201	wgauss[2] = 0.381830050505118944950369775488975;
	3202	wgauss[3] = 0.417959183673469387755102040816327;
	3203	x[0] = 0.991455371120812639206854697526329;
	3204	x[1] = 0.949107912342758524526189684047851;
	3205	x[2] = 0.864864423359769072789712788640926;
	3206	x[3] = 0.741531185599394439863864773280788;
	3207	x[4] = 0.586087235467691130294144838258730;
	3208	x[5] = 0.405845151377397166906606412076961;
	3209	x[6] = 0.207784955007898467600689403773245;
	3210	x[7] = 0.000000000000000000000000000000000;
	3211	wkronrod[0] = 0.022935322010529224963732008058970;
	3212	wkronrod[1] = 0.063092092629978553290700663189204;
	3213	wkronrod[2] = 0.104790010322250183839876322541518;
	3214	wkronrod[3] = 0.140653259715525918745189590510238;
	3215	wkronrod[4] = 0.169004726639267902826583426598550;
	3216	wkronrod[5] = 0.190350578064785409913256402421014;
	3217	wkronrod[6] = 0.204432940075298892414161999234649;
	3218	wkronrod[7] = 0.209482141084727828012999174891714;
	3219	}
	3220	if( n==21 )
	3221	{
	3222	ng = 5;
	3223	wgauss[0] = 0.066671344308688137593568809893332;
	3224	wgauss[1] = 0.149451349150580593145776339657697;
	3225	wgauss[2] = 0.219086362515982043995534934228163;
	3226	wgauss[3] = 0.269266719309996355091226921569469;
	3227	wgauss[4] = 0.295524224714752870173892994651338;
	3228	x[0] = 0.995657163025808080735527280689003;
	3229	x[1] = 0.973906528517171720077964012084452;
	3230	x[2] = 0.930157491355708226001207180059508;
	3231	x[3] = 0.865063366688984510732096688423493;
	3232	x[4] = 0.780817726586416897063717578345042;
	3233	x[5] = 0.679409568299024406234327365114874;
	3234	x[6] = 0.562757134668604683339000099272694;
	3235	x[7] = 0.433395394129247190799265943165784;
	3236	x[8] = 0.294392862701460198131126603103866;
	3237	x[9] = 0.148874338981631210884826001129720;
	3238	x[10] = 0.000000000000000000000000000000000;
	3239	wkronrod[0] = 0.011694638867371874278064396062192;
	3240	wkronrod[1] = 0.032558162307964727478818972459390;
	3241	wkronrod[2] = 0.054755896574351996031381300244580;
	3242	wkronrod[3] = 0.075039674810919952767043140916190;
	3243	wkronrod[4] = 0.093125454583697605535065465083366;
	3244	wkronrod[5] = 0.109387158802297641899210590325805;
	3245	wkronrod[6] = 0.123491976262065851077958109831074;
	3246	wkronrod[7] = 0.134709217311473325928054001771707;
	3247	wkronrod[8] = 0.142775938577060080797094273138717;
	3248	wkronrod[9] = 0.147739104901338491374841515972068;
	3249	wkronrod[10] = 0.149445554002916905664936468389821;
	3250	}
	3251	if( n==31 )
	3252	{
	3253	ng = 8;
	3254	wgauss[0] = 0.030753241996117268354628393577204;
	3255	wgauss[1] = 0.070366047488108124709267416450667;
	3256	wgauss[2] = 0.107159220467171935011869546685869;
	3257	wgauss[3] = 0.139570677926154314447804794511028;
	3258	wgauss[4] = 0.166269205816993933553200860481209;
	3259	wgauss[5] = 0.186161000015562211026800561866423;
	3260	wgauss[6] = 0.198431485327111576456118326443839;
	3261	wgauss[7] = 0.202578241925561272880620199967519;
	3262	x[0] = 0.998002298693397060285172840152271;
	3263	x[1] = 0.987992518020485428489565718586613;
	3264	x[2] = 0.967739075679139134257347978784337;
	3265	x[3] = 0.937273392400705904307758947710209;
	3266	x[4] = 0.897264532344081900882509656454496;
	3267	x[5] = 0.848206583410427216200648320774217;
	3268	x[6] = 0.790418501442465932967649294817947;
	3269	x[7] = 0.724417731360170047416186054613938;
	3270	x[8] = 0.650996741297416970533735895313275;
	3271	x[9] = 0.570972172608538847537226737253911;
	3272	x[10] = 0.485081863640239680693655740232351;
	3273	x[11] = 0.394151347077563369897207370981045;
	3274	x[12] = 0.299180007153168812166780024266389;
	3275	x[13] = 0.201194093997434522300628303394596;
	3276	x[14] = 0.101142066918717499027074231447392;
	3277	x[15] = 0.000000000000000000000000000000000;
	3278	wkronrod[0] = 0.005377479872923348987792051430128;
	3279	wkronrod[1] = 0.015007947329316122538374763075807;
	3280	wkronrod[2] = 0.025460847326715320186874001019653;
	3281	wkronrod[3] = 0.035346360791375846222037948478360;
	3282	wkronrod[4] = 0.044589751324764876608227299373280;
	3283	wkronrod[5] = 0.053481524690928087265343147239430;
	3284	wkronrod[6] = 0.062009567800670640285139230960803;
	3285	wkronrod[7] = 0.069854121318728258709520077099147;
	3286	wkronrod[8] = 0.076849680757720378894432777482659;
	3287	wkronrod[9] = 0.083080502823133021038289247286104;
	3288	wkronrod[10] = 0.088564443056211770647275443693774;
	3289	wkronrod[11] = 0.093126598170825321225486872747346;
	3290	wkronrod[12] = 0.096642726983623678505179907627589;
	3291	wkronrod[13] = 0.099173598721791959332393173484603;
	3292	wkronrod[14] = 0.100769845523875595044946662617570;
	3293	wkronrod[15] = 0.101330007014791549017374792767493;
	3294	}
	3295	if( n==41 )
	3296	{
	3297	ng = 10;
	3298	wgauss[0] = 0.017614007139152118311861962351853;
	3299	wgauss[1] = 0.040601429800386941331039952274932;
	3300	wgauss[2] = 0.062672048334109063569506535187042;
	3301	wgauss[3] = 0.083276741576704748724758143222046;
	3302	wgauss[4] = 0.101930119817240435036750135480350;
	3303	wgauss[5] = 0.118194531961518417312377377711382;
	3304	wgauss[6] = 0.131688638449176626898494499748163;
	3305	wgauss[7] = 0.142096109318382051329298325067165;
	3306	wgauss[8] = 0.149172986472603746787828737001969;
	3307	wgauss[9] = 0.152753387130725850698084331955098;
	3308	x[0] = 0.998859031588277663838315576545863;
	3309	x[1] = 0.993128599185094924786122388471320;
	3310	x[2] = 0.981507877450250259193342994720217;
	3311	x[3] = 0.963971927277913791267666131197277;
	3312	x[4] = 0.940822633831754753519982722212443;
	3313	x[5] = 0.912234428251325905867752441203298;
	3314	x[6] = 0.878276811252281976077442995113078;
	3315	x[7] = 0.839116971822218823394529061701521;
	3316	x[8] = 0.795041428837551198350638833272788;
	3317	x[9] = 0.746331906460150792614305070355642;
	3318	x[10] = 0.693237656334751384805490711845932;
	3319	x[11] = 0.636053680726515025452836696226286;
	3320	x[12] = 0.575140446819710315342946036586425;
	3321	x[13] = 0.510867001950827098004364050955251;
	3322	x[14] = 0.443593175238725103199992213492640;
	3323	x[15] = 0.373706088715419560672548177024927;
	3324	x[16] = 0.301627868114913004320555356858592;
	3325	x[17] = 0.227785851141645078080496195368575;
	3326	x[18] = 0.152605465240922675505220241022678;
	3327	x[19] = 0.076526521133497333754640409398838;
	3328	x[20] = 0.000000000000000000000000000000000;
	3329	wkronrod[0] = 0.003073583718520531501218293246031;
	3330	wkronrod[1] = 0.008600269855642942198661787950102;
	3331	wkronrod[2] = 0.014626169256971252983787960308868;
	3332	wkronrod[3] = 0.020388373461266523598010231432755;
	3333	wkronrod[4] = 0.025882133604951158834505067096153;
	3334	wkronrod[5] = 0.031287306777032798958543119323801;
	3335	wkronrod[6] = 0.036600169758200798030557240707211;
	3336	wkronrod[7] = 0.041668873327973686263788305936895;
	3337	wkronrod[8] = 0.046434821867497674720231880926108;
	3338	wkronrod[9] = 0.050944573923728691932707670050345;
	3339	wkronrod[10] = 0.055195105348285994744832372419777;
	3340	wkronrod[11] = 0.059111400880639572374967220648594;
	3341	wkronrod[12] = 0.062653237554781168025870122174255;
	3342	wkronrod[13] = 0.065834597133618422111563556969398;
	3343	wkronrod[14] = 0.068648672928521619345623411885368;
	3344	wkronrod[15] = 0.071054423553444068305790361723210;
	3345	wkronrod[16] = 0.073030690332786667495189417658913;
	3346	wkronrod[17] = 0.074582875400499188986581418362488;
	3347	wkronrod[18] = 0.075704497684556674659542775376617;
	3348	wkronrod[19] = 0.076377867672080736705502835038061;
	3349	wkronrod[20] = 0.076600711917999656445049901530102;
	3350	}
	3351	if( n==51 )
	3352	{
	3353	ng = 13;
	3354	wgauss[0] = 0.011393798501026287947902964113235;
	3355	wgauss[1] = 0.026354986615032137261901815295299;
	3356	wgauss[2] = 0.040939156701306312655623487711646;
	3357	wgauss[3] = 0.054904695975835191925936891540473;
	3358	wgauss[4] = 0.068038333812356917207187185656708;
	3359	wgauss[5] = 0.080140700335001018013234959669111;
	3360	wgauss[6] = 0.091028261982963649811497220702892;
	3361	wgauss[7] = 0.100535949067050644202206890392686;
	3362	wgauss[8] = 0.108519624474263653116093957050117;
	3363	wgauss[9] = 0.114858259145711648339325545869556;
	3364	wgauss[10] = 0.119455763535784772228178126512901;
	3365	wgauss[11] = 0.122242442990310041688959518945852;
	3366	wgauss[12] = 0.123176053726715451203902873079050;
	3367	x[0] = 0.999262104992609834193457486540341;
	3368	x[1] = 0.995556969790498097908784946893902;
	3369	x[2] = 0.988035794534077247637331014577406;
	3370	x[3] = 0.976663921459517511498315386479594;
	3371	x[4] = 0.961614986425842512418130033660167;
	3372	x[5] = 0.942974571228974339414011169658471;
	3373	x[6] = 0.920747115281701561746346084546331;
	3374	x[7] = 0.894991997878275368851042006782805;
	3375	x[8] = 0.865847065293275595448996969588340;
	3376	x[9] = 0.833442628760834001421021108693570;
	3377	x[10] = 0.797873797998500059410410904994307;
	3378	x[11] = 0.759259263037357630577282865204361;
	3379	x[12] = 0.717766406813084388186654079773298;
	3380	x[13] = 0.673566368473468364485120633247622;
	3381	x[14] = 0.626810099010317412788122681624518;
	3382	x[15] = 0.577662930241222967723689841612654;
	3383	x[16] = 0.526325284334719182599623778158010;
	3384	x[17] = 0.473002731445714960522182115009192;
	3385	x[18] = 0.417885382193037748851814394594572;
	3386	x[19] = 0.361172305809387837735821730127641;
	3387	x[20] = 0.303089538931107830167478909980339;
	3388	x[21] = 0.243866883720988432045190362797452;
	3389	x[22] = 0.183718939421048892015969888759528;
	3390	x[23] = 0.122864692610710396387359818808037;
	3391	x[24] = 0.061544483005685078886546392366797;
	3392	x[25] = 0.000000000000000000000000000000000;
	3393	wkronrod[0] = 0.001987383892330315926507851882843;
	3394	wkronrod[1] = 0.005561932135356713758040236901066;
	3395	wkronrod[2] = 0.009473973386174151607207710523655;
	3396	wkronrod[3] = 0.013236229195571674813656405846976;
	3397	wkronrod[4] = 0.016847817709128298231516667536336;
	3398	wkronrod[5] = 0.020435371145882835456568292235939;
	3399	wkronrod[6] = 0.024009945606953216220092489164881;
	3400	wkronrod[7] = 0.027475317587851737802948455517811;
	3401	wkronrod[8] = 0.030792300167387488891109020215229;
	3402	wkronrod[9] = 0.034002130274329337836748795229551;
	3403	wkronrod[10] = 0.037116271483415543560330625367620;
	3404	wkronrod[11] = 0.040083825504032382074839284467076;
	3405	wkronrod[12] = 0.042872845020170049476895792439495;
	3406	wkronrod[13] = 0.045502913049921788909870584752660;
	3407	wkronrod[14] = 0.047982537138836713906392255756915;
	3408	wkronrod[15] = 0.050277679080715671963325259433440;
	3409	wkronrod[16] = 0.052362885806407475864366712137873;
	3410	wkronrod[17] = 0.054251129888545490144543370459876;
	3411	wkronrod[18] = 0.055950811220412317308240686382747;
	3412	wkronrod[19] = 0.057437116361567832853582693939506;
	3413	wkronrod[20] = 0.058689680022394207961974175856788;
	3414	wkronrod[21] = 0.059720340324174059979099291932562;
	3415	wkronrod[22] = 0.060539455376045862945360267517565;
	3416	wkronrod[23] = 0.061128509717053048305859030416293;
	3417	wkronrod[24] = 0.061471189871425316661544131965264;
	3418	wkronrod[25] = 0.061580818067832935078759824240055;
	3419	}
	3420	if( n==61 )
	3421	{
	3422	ng = 15;
	3423	wgauss[0] = 0.007968192496166605615465883474674;
	3424	wgauss[1] = 0.018466468311090959142302131912047;
	3425	wgauss[2] = 0.028784707883323369349719179611292;
	3426	wgauss[3] = 0.038799192569627049596801936446348;
	3427	wgauss[4] = 0.048402672830594052902938140422808;
	3428	wgauss[5] = 0.057493156217619066481721689402056;
	3429	wgauss[6] = 0.065974229882180495128128515115962;
	3430	wgauss[7] = 0.073755974737705206268243850022191;
	3431	wgauss[8] = 0.080755895229420215354694938460530;
	3432	wgauss[9] = 0.086899787201082979802387530715126;
	3433	wgauss[10] = 0.092122522237786128717632707087619;
	3434	wgauss[11] = 0.096368737174644259639468626351810;
	3435	wgauss[12] = 0.099593420586795267062780282103569;
	3436	wgauss[13] = 0.101762389748405504596428952168554;
	3437	wgauss[14] = 0.102852652893558840341285636705415;
	3438	x[0] = 0.999484410050490637571325895705811;
	3439	x[1] = 0.996893484074649540271630050918695;
	3440	x[2] = 0.991630996870404594858628366109486;
	3441	x[3] = 0.983668123279747209970032581605663;
	3442	x[4] = 0.973116322501126268374693868423707;
	3443	x[5] = 0.960021864968307512216871025581798;
	3444	x[6] = 0.944374444748559979415831324037439;
	3445	x[7] = 0.926200047429274325879324277080474;
	3446	x[8] = 0.905573307699907798546522558925958;
	3447	x[9] = 0.882560535792052681543116462530226;
	3448	x[10] = 0.857205233546061098958658510658944;
	3449	x[11] = 0.829565762382768397442898119732502;
	3450	x[12] = 0.799727835821839083013668942322683;
	3451	x[13] = 0.767777432104826194917977340974503;
	3452	x[14] = 0.733790062453226804726171131369528;
	3453	x[15] = 0.697850494793315796932292388026640;
	3454	x[16] = 0.660061064126626961370053668149271;
	3455	x[17] = 0.620526182989242861140477556431189;
	3456	x[18] = 0.579345235826361691756024932172540;
	3457	x[19] = 0.536624148142019899264169793311073;
	3458	x[20] = 0.492480467861778574993693061207709;
	3459	x[21] = 0.447033769538089176780609900322854;
	3460	x[22] = 0.400401254830394392535476211542661;
	3461	x[23] = 0.352704725530878113471037207089374;
	3462	x[24] = 0.304073202273625077372677107199257;
	3463	x[25] = 0.254636926167889846439805129817805;
	3464	x[26] = 0.204525116682309891438957671002025;
	3465	x[27] = 0.153869913608583546963794672743256;
	3466	x[28] = 0.102806937966737030147096751318001;
	3467	x[29] = 0.051471842555317695833025213166723;
	3468	x[30] = 0.000000000000000000000000000000000;
	3469	wkronrod[0] = 0.001389013698677007624551591226760;
	3470	wkronrod[1] = 0.003890461127099884051267201844516;
	3471	wkronrod[2] = 0.006630703915931292173319826369750;
	3472	wkronrod[3] = 0.009273279659517763428441146892024;
	3473	wkronrod[4] = 0.011823015253496341742232898853251;
	3474	wkronrod[5] = 0.014369729507045804812451432443580;
	3475	wkronrod[6] = 0.016920889189053272627572289420322;
	3476	wkronrod[7] = 0.019414141193942381173408951050128;
	3477	wkronrod[8] = 0.021828035821609192297167485738339;
	3478	wkronrod[9] = 0.024191162078080601365686370725232;
	3479	wkronrod[10] = 0.026509954882333101610601709335075;
	3480	wkronrod[11] = 0.028754048765041292843978785354334;
	3481	wkronrod[12] = 0.030907257562387762472884252943092;
	3482	wkronrod[13] = 0.032981447057483726031814191016854;
	3483	wkronrod[14] = 0.034979338028060024137499670731468;
	3484	wkronrod[15] = 0.036882364651821229223911065617136;
	3485	wkronrod[16] = 0.038678945624727592950348651532281;
	3486	wkronrod[17] = 0.040374538951535959111995279752468;
	3487	wkronrod[18] = 0.041969810215164246147147541285970;
	3488	wkronrod[19] = 0.043452539701356069316831728117073;
	3489	wkronrod[20] = 0.044814800133162663192355551616723;
	3490	wkronrod[21] = 0.046059238271006988116271735559374;
	3491	wkronrod[22] = 0.047185546569299153945261478181099;
	3492	wkronrod[23] = 0.048185861757087129140779492298305;
	3493	wkronrod[24] = 0.049055434555029778887528165367238;
	3494	wkronrod[25] = 0.049795683427074206357811569379942;
	3495	wkronrod[26] = 0.050405921402782346840893085653585;
	3496	wkronrod[27] = 0.050881795898749606492297473049805;
	3497	wkronrod[28] = 0.051221547849258772170656282604944;
	3498	wkronrod[29] = 0.051426128537459025933862879215781;
	3499	wkronrod[30] = 0.051494729429451567558340433647099;
	3500	}
	3501
	3502	//
	3503	// copy nodes
	3504	//
	3505	for(i=n-1; i>=n/2; i--)
	3506	{
	3507	x[i] = -x[n-1-i];
	3508	}
	3509
	3510	//
	3511	// copy Kronrod weights
	3512	//
	3513	for(i=n-1; i>=n/2; i--)
	3514	{
	3515	wkronrod[i] = wkronrod[n-1-i];
	3516	}
	3517
	3518	//
	3519	// copy Gauss weights
	3520	//
	3521	for(i=ng-1; i>=0; i--)
	3522	{
	3523	wgauss[n-2-2*i] = wgauss[i];
	3524	wgauss[1+2*i] = wgauss[i];
	3525	}
	3526	for(i=0; i<=n/2; i++)
	3527	{
	3528	wgauss[2*i] = 0;
	3529	}
	3530
	3531	//
	3532	// reorder
	3533	//
	3534	tsort.tagsort(ref x, n, ref p1, ref p2);
	3535	for(i=0; i<=n-1; i++)
	3536	{
	3537	tmp = wkronrod[i];
	3538	wkronrod[i] = wkronrod[p2[i]];
	3539	wkronrod[p2[i]] = tmp;
	3540	tmp = wgauss[i];
	3541	wgauss[i] = wgauss[p2[i]];
	3542	wgauss[p2[i]] = tmp;
	3543	}
	3544	}
	3545
	3546
	3547	}
	3548	}
	3549

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences