Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.9.0/ALGLIB-3.9.0/integration.cs @ 12790

Visit:

Last change on this file since 12790 was 12790, checked in by gkronber, 9 years ago
#2435: updated alglib to version 3.9.0
File size: 134.4 KB

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

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

Download in other formats:

Update cookies preferences