Context Navigation

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

Visit:

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

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

Download in other formats:

Update cookies preferences