Context Navigation

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

Visit:

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

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

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

Download in other formats:

Update cookies preferences