Context Navigation

source: branches/1888_OaaS/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.6.0/ALGLIB-3.6.0/integration.cs @ 16811

Visit:

Last change on this file since 16811 was 8418, checked in by abeham, 12 years ago
#1905: Added ALGLIB 3.6.0 to ExtLibs
File size: 131.4 KB

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

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

Download in other formats:

Update cookies preferences