Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2434_crossvalidation/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.7.0/ALGLIB-3.7.0/integration.cs @ 16202

Last change on this file since 16202 was 9268, checked in by mkommend, 12 years ago

#2007: Added AlgLib 3.7.0 to the external libraries.

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