Context Navigation

source: branches/DataAnalysis Refactoring/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.1.0/ALGLIB-3.1.0/integration.cs @ 5779

Visit:

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

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

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

Download in other formats:

Update cookies preferences