Context Navigation

source: branches/2435-alglib_3_15/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.15.0/ALGLIB-3.15.0/integration.cs @ 16684

Visit:

Last change on this file since 16684 was 16684, checked in by gkronber, 5 years ago
#2435: added 3.15.0 version of alglib as external library (build from source)
File size: 140.9 KB

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

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

Download in other formats:

Update cookies preferences