Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.7.0/ALGLIB-3.7.0/integration.cs @ 14778

Visit:

Last change on this file since 14778 was 9268, checked in by mkommend, 12 years ago
#2007: Added AlgLib 3.7.0 to the external libraries.
File size: 134.0 KB

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

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

Download in other formats:

Update cookies preferences