Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/gkq.cs @ 2806

Last change on this file since 2806 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)

File size: 40.5 KB
Line 
1/*************************************************************************
2Copyright (c) 2005-2009, 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
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class gkq
26    {
27        /*************************************************************************
28        Computation of nodes and weights of a Gauss-Kronrod quadrature formula
29
30        The algorithm generates the N-point Gauss-Kronrod quadrature formula  with
31        weight  function  given  by  coefficients  alpha  and beta of a recurrence
32        relation 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 zero moment Mu0
39
40            Mu0 = integral(W(x)dx,a,b)
41
42
43        INPUT PARAMETERS:
44            Alpha       –   alpha coefficients, array[0..floor(3*K/2)].
45            Beta        –   beta coefficients,  array[0..ceil(3*K/2)].
46                            Beta[0] is not used and may be arbitrary.
47                            Beta[I]>0.
48            Mu0         –   zeroth moment of the weight function.
49            N           –   number of nodes of the Gauss-Kronrod quadrature formula,
50                            N >= 3,
51                            N =  2*K+1.
52
53        OUTPUT PARAMETERS:
54            Info        -   error code:
55                            * -5    no real and positive Gauss-Kronrod formula can
56                                    be created for such a weight function  with  a
57                                    given number of nodes.
58                            * -4    N is too large, task may be ill  conditioned -
59                                    x[i]=x[i+1] found.
60                            * -3    internal eigenproblem solver hasn't converged
61                            * -2    Beta[i]<=0
62                            * -1    incorrect N was passed
63                            * +1    OK
64            X           -   array[0..N-1] - array of quadrature nodes,
65                            in ascending order.
66            WKronrod    -   array[0..N-1] - Kronrod weights
67            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
68                            corresponding to extended Kronrod nodes).
69
70          -- ALGLIB --
71             Copyright 08.05.2009 by Bochkanov Sergey
72        *************************************************************************/
73        public static void gkqgeneraterec(double[] alpha,
74            double[] beta,
75            double mu0,
76            int n,
77            ref int info,
78            ref double[] x,
79            ref double[] wkronrod,
80            ref double[] wgauss)
81        {
82            double[] ta = new double[0];
83            int i = 0;
84            int j = 0;
85            double[] t = new double[0];
86            double[] s = new double[0];
87            int wlen = 0;
88            int woffs = 0;
89            double u = 0;
90            int m = 0;
91            int l = 0;
92            int k = 0;
93            double[] xgtmp = new double[0];
94            double[] wgtmp = new double[0];
95            int i_ = 0;
96
97            alpha = (double[])alpha.Clone();
98            beta = (double[])beta.Clone();
99
100            if( n%2!=1 | n<3 )
101            {
102                info = -1;
103                return;
104            }
105            for(i=0; i<=(int)Math.Ceiling((double)(3*(n/2))/(double)(2)); i++)
106            {
107                if( (double)(beta[i])<=(double)(0) )
108                {
109                    info = -2;
110                    return;
111                }
112            }
113            info = 1;
114           
115            //
116            // from external conventions about N/Beta/Mu0 to internal
117            //
118            n = n/2;
119            beta[0] = mu0;
120           
121            //
122            // Calculate Gauss nodes/weights, save them for later processing
123            //
124            gq.gqgeneraterec(ref alpha, ref beta, mu0, n, ref info, ref xgtmp, ref wgtmp);
125            if( info<0 )
126            {
127                return;
128            }
129           
130            //
131            // Resize:
132            // * A from 0..floor(3*n/2) to 0..2*n
133            // * B from 0..ceil(3*n/2)  to 0..2*n
134            //
135            ta = new double[(int)Math.Floor((double)(3*n)/(double)(2))+1];
136            for(i_=0; i_<=(int)Math.Floor((double)(3*n)/(double)(2));i_++)
137            {
138                ta[i_] = alpha[i_];
139            }
140            alpha = new double[2*n+1];
141            for(i_=0; i_<=(int)Math.Floor((double)(3*n)/(double)(2));i_++)
142            {
143                alpha[i_] = ta[i_];
144            }
145            for(i=(int)Math.Floor((double)(3*n)/(double)(2))+1; i<=2*n; i++)
146            {
147                alpha[i] = 0;
148            }
149            ta = new double[(int)Math.Ceiling((double)(3*n)/(double)(2))+1];
150            for(i_=0; i_<=(int)Math.Ceiling((double)(3*n)/(double)(2));i_++)
151            {
152                ta[i_] = beta[i_];
153            }
154            beta = new double[2*n+1];
155            for(i_=0; i_<=(int)Math.Ceiling((double)(3*n)/(double)(2));i_++)
156            {
157                beta[i_] = ta[i_];
158            }
159            for(i=(int)Math.Ceiling((double)(3*n)/(double)(2))+1; i<=2*n; i++)
160            {
161                beta[i] = 0;
162            }
163           
164            //
165            // Initialize T, S
166            //
167            wlen = 2+n/2;
168            t = new double[wlen];
169            s = new double[wlen];
170            ta = new double[wlen];
171            woffs = 1;
172            for(i=0; i<=wlen-1; i++)
173            {
174                t[i] = 0;
175                s[i] = 0;
176            }
177           
178            //
179            // Algorithm from Dirk P. Laurie, "Calculation of Gauss-Kronrod quadrature rules", 1997.
180            //
181            t[woffs+0] = beta[n+1];
182            for(m=0; m<=n-2; m++)
183            {
184                u = 0;
185                for(k=(m+1)/2; k>=0; k--)
186                {
187                    l = m-k;
188                    u = u+(alpha[k+n+1]-alpha[l])*t[woffs+k]+beta[k+n+1]*s[woffs+k-1]-beta[l]*s[woffs+k];
189                    s[woffs+k] = u;
190                }
191                for(i_=0; i_<=wlen-1;i_++)
192                {
193                    ta[i_] = t[i_];
194                }
195                for(i_=0; i_<=wlen-1;i_++)
196                {
197                    t[i_] = s[i_];
198                }
199                for(i_=0; i_<=wlen-1;i_++)
200                {
201                    s[i_] = ta[i_];
202                }
203            }
204            for(j=n/2; j>=0; j--)
205            {
206                s[woffs+j] = s[woffs+j-1];
207            }
208            for(m=n-1; m<=2*n-3; m++)
209            {
210                u = 0;
211                for(k=m+1-n; k<=(m-1)/2; k++)
212                {
213                    l = m-k;
214                    j = n-1-l;
215                    u = u-(alpha[k+n+1]-alpha[l])*t[woffs+j]-beta[k+n+1]*s[woffs+j]+beta[l]*s[woffs+j+1];
216                    s[woffs+j] = u;
217                }
218                if( m%2==0 )
219                {
220                    k = m/2;
221                    alpha[k+n+1] = alpha[k]+(s[woffs+j]-beta[k+n+1]*s[woffs+j+1])/t[woffs+j+1];
222                }
223                else
224                {
225                    k = (m+1)/2;
226                    beta[k+n+1] = s[woffs+j]/s[woffs+j+1];
227                }
228                for(i_=0; i_<=wlen-1;i_++)
229                {
230                    ta[i_] = t[i_];
231                }
232                for(i_=0; i_<=wlen-1;i_++)
233                {
234                    t[i_] = s[i_];
235                }
236                for(i_=0; i_<=wlen-1;i_++)
237                {
238                    s[i_] = ta[i_];
239                }
240            }
241            alpha[2*n] = alpha[n-1]-beta[2*n]*s[woffs+0]/t[woffs+0];
242           
243            //
244            // calculation of Kronrod nodes and weights, unpacking of Gauss weights
245            //
246            gq.gqgeneraterec(ref alpha, ref beta, mu0, 2*n+1, ref info, ref x, ref wkronrod);
247            if( info==-2 )
248            {
249                info = -5;
250            }
251            if( info<0 )
252            {
253                return;
254            }
255            for(i=0; i<=2*n-1; i++)
256            {
257                if( (double)(x[i])>=(double)(x[i+1]) )
258                {
259                    info = -4;
260                }
261            }
262            if( info<0 )
263            {
264                return;
265            }
266            wgauss = new double[2*n+1];
267            for(i=0; i<=2*n; i++)
268            {
269                wgauss[i] = 0;
270            }
271            for(i=0; i<=n-1; i++)
272            {
273                wgauss[2*i+1] = wgtmp[i];
274            }
275        }
276
277
278        /*************************************************************************
279        Returns   Gauss   and   Gauss-Kronrod   nodes/weights  for  Gauss-Legendre
280        quadrature with N points.
281
282        GKQLegendreCalc (calculation) or  GKQLegendreTbl  (precomputed  table)  is
283        used depending on machine precision and number of nodes.
284
285        INPUT PARAMETERS:
286            N           -   number of Kronrod nodes, must be odd number, >=3.
287
288        OUTPUT PARAMETERS:
289            Info        -   error code:
290                            * -4    an  error   was   detected   when  calculating
291                                    weights/nodes.  N  is  too  large   to  obtain
292                                    weights/nodes  with  high   enough   accuracy.
293                                    Try  to   use   multiple   precision  version.
294                            * -3    internal eigenproblem solver hasn't converged
295                            * -1    incorrect N was passed
296                            * +1    OK
297            X           -   array[0..N-1] - array of quadrature nodes, ordered in
298                            ascending order.
299            WKronrod    -   array[0..N-1] - Kronrod weights
300            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
301                            corresponding to extended Kronrod nodes).
302
303
304          -- ALGLIB --
305             Copyright 12.05.2009 by Bochkanov Sergey
306        *************************************************************************/
307        public static void gkqgenerategausslegendre(int n,
308            ref int info,
309            ref double[] x,
310            ref double[] wkronrod,
311            ref double[] wgauss)
312        {
313            double eps = 0;
314
315            if( (double)(AP.Math.MachineEpsilon)>(double)(1.0E-32) & (n==15 | n==21 | n==31 | n==41 | n==51 | n==61) )
316            {
317                info = 1;
318                gkqlegendretbl(n, ref x, ref wkronrod, ref wgauss, ref eps);
319            }
320            else
321            {
322                gkqlegendrecalc(n, ref info, ref x, ref wkronrod, ref wgauss);
323            }
324        }
325
326
327        /*************************************************************************
328        Returns   Gauss   and   Gauss-Kronrod   nodes/weights   for   Gauss-Jacobi
329        quadrature on [-1,1] with weight function
330
331            W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
332
333        INPUT PARAMETERS:
334            N           -   number of Kronrod nodes, must be odd number, >=3.
335            Alpha       -   power-law coefficient, Alpha>-1
336            Beta        -   power-law coefficient, Beta>-1
337
338        OUTPUT PARAMETERS:
339            Info        -   error code:
340                            * -5    no real and positive Gauss-Kronrod formula can
341                                    be created for such a weight function  with  a
342                                    given number of nodes.
343                            * -4    an  error  was   detected   when   calculating
344                                    weights/nodes. Alpha or  Beta  are  too  close
345                                    to -1 to obtain weights/nodes with high enough
346                                    accuracy, or, may be, N is too large.  Try  to
347                                    use multiple precision version.
348                            * -3    internal eigenproblem solver hasn't converged
349                            * -1    incorrect N was passed
350                            * +1    OK
351                            * +2    OK, but quadrature rule have exterior  nodes,
352                                    x[0]<-1 or x[n-1]>+1
353            X           -   array[0..N-1] - array of quadrature nodes, ordered in
354                            ascending order.
355            WKronrod    -   array[0..N-1] - Kronrod weights
356            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
357                            corresponding to extended Kronrod nodes).
358
359
360          -- ALGLIB --
361             Copyright 12.05.2009 by Bochkanov Sergey
362        *************************************************************************/
363        public static void gkqgenerategaussjacobi(int n,
364            double alpha,
365            double beta,
366            ref int info,
367            ref double[] x,
368            ref double[] wkronrod,
369            ref double[] wgauss)
370        {
371            int clen = 0;
372            double[] a = new double[0];
373            double[] b = new double[0];
374            double alpha2 = 0;
375            double beta2 = 0;
376            double apb = 0;
377            double t = 0;
378            int i = 0;
379            double s = 0;
380
381            if( n%2!=1 | n<3 )
382            {
383                info = -1;
384                return;
385            }
386            if( (double)(alpha)<=(double)(-1) | (double)(beta)<=(double)(-1) )
387            {
388                info = -1;
389                return;
390            }
391            clen = (int)Math.Ceiling((double)(3*(n/2))/(double)(2))+1;
392            a = new double[clen];
393            b = new double[clen];
394            for(i=0; i<=clen-1; i++)
395            {
396                a[i] = 0;
397            }
398            apb = alpha+beta;
399            a[0] = (beta-alpha)/(apb+2);
400            t = (apb+1)*Math.Log(2)+gammafunc.lngamma(alpha+1, ref s)+gammafunc.lngamma(beta+1, ref s)-gammafunc.lngamma(apb+2, ref s);
401            if( (double)(t)>(double)(Math.Log(AP.Math.MaxRealNumber)) )
402            {
403                info = -4;
404                return;
405            }
406            b[0] = Math.Exp(t);
407            if( clen>1 )
408            {
409                alpha2 = AP.Math.Sqr(alpha);
410                beta2 = AP.Math.Sqr(beta);
411                a[1] = (beta2-alpha2)/((apb+2)*(apb+4));
412                b[1] = 4*(alpha+1)*(beta+1)/((apb+3)*AP.Math.Sqr(apb+2));
413                for(i=2; i<=clen-1; i++)
414                {
415                    a[i] = 0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
416                    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)*AP.Math.Sqr(1+0.5*apb/i));
417                }
418            }
419            gkqgeneraterec(a, b, b[0], n, ref info, ref x, ref wkronrod, ref wgauss);
420           
421            //
422            // test basic properties to detect errors
423            //
424            if( info>0 )
425            {
426                if( (double)(x[0])<(double)(-1) | (double)(x[n-1])>(double)(+1) )
427                {
428                    info = 2;
429                }
430                for(i=0; i<=n-2; i++)
431                {
432                    if( (double)(x[i])>=(double)(x[i+1]) )
433                    {
434                        info = -4;
435                    }
436                }
437            }
438        }
439
440
441        /*************************************************************************
442        Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
443
444        Reduction to tridiagonal eigenproblem is used.
445
446        INPUT PARAMETERS:
447            N           -   number of Kronrod nodes, must be odd number, >=3.
448
449        OUTPUT PARAMETERS:
450            Info        -   error code:
451                            * -4    an  error   was   detected   when  calculating
452                                    weights/nodes.  N  is  too  large   to  obtain
453                                    weights/nodes  with  high   enough   accuracy.
454                                    Try  to   use   multiple   precision  version.
455                            * -3    internal eigenproblem solver hasn't converged
456                            * -1    incorrect N was passed
457                            * +1    OK
458            X           -   array[0..N-1] - array of quadrature nodes, ordered in
459                            ascending order.
460            WKronrod    -   array[0..N-1] - Kronrod weights
461            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
462                            corresponding to extended Kronrod nodes).
463
464          -- ALGLIB --
465             Copyright 12.05.2009 by Bochkanov Sergey
466        *************************************************************************/
467        public static void gkqlegendrecalc(int n,
468            ref int info,
469            ref double[] x,
470            ref double[] wkronrod,
471            ref double[] wgauss)
472        {
473            double[] alpha = new double[0];
474            double[] beta = new double[0];
475            int alen = 0;
476            int blen = 0;
477            double mu0 = 0;
478            int k = 0;
479            int i = 0;
480
481            if( n%2!=1 | n<3 )
482            {
483                info = -1;
484                return;
485            }
486            mu0 = 2;
487            alen = (int)Math.Floor((double)(3*(n/2))/(double)(2))+1;
488            blen = (int)Math.Ceiling((double)(3*(n/2))/(double)(2))+1;
489            alpha = new double[alen];
490            beta = new double[blen];
491            for(k=0; k<=alen-1; k++)
492            {
493                alpha[k] = 0;
494            }
495            beta[0] = 2;
496            for(k=1; k<=blen-1; k++)
497            {
498                beta[k] = 1/(4-1/AP.Math.Sqr(k));
499            }
500            gkqgeneraterec(alpha, beta, mu0, n, ref info, ref x, ref wkronrod, ref wgauss);
501           
502            //
503            // test basic properties to detect errors
504            //
505            if( info>0 )
506            {
507                if( (double)(x[0])<(double)(-1) | (double)(x[n-1])>(double)(+1) )
508                {
509                    info = -4;
510                }
511                for(i=0; i<=n-2; i++)
512                {
513                    if( (double)(x[i])>=(double)(x[i+1]) )
514                    {
515                        info = -4;
516                    }
517                }
518            }
519        }
520
521
522        /*************************************************************************
523        Returns Gauss and Gauss-Kronrod nodes for quadrature with N  points  using
524        pre-calculated table. Nodes/weights were  computed  with  accuracy  up  to
525        1.0E-32 (if MPFR version of ALGLIB is used). In standard double  precision
526        accuracy reduces to something about 2.0E-16 (depending  on your compiler's
527        handling of long floating point constants).
528
529        INPUT PARAMETERS:
530            N           -   number of Kronrod nodes.
531                            N can be 15, 21, 31, 41, 51, 61.
532
533        OUTPUT PARAMETERS:
534            X           -   array[0..N-1] - array of quadrature nodes, ordered in
535                            ascending order.
536            WKronrod    -   array[0..N-1] - Kronrod weights
537            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
538                            corresponding to extended Kronrod nodes).
539
540
541          -- ALGLIB --
542             Copyright 12.05.2009 by Bochkanov Sergey
543        *************************************************************************/
544        public static void gkqlegendretbl(int n,
545            ref double[] x,
546            ref double[] wkronrod,
547            ref double[] wgauss,
548            ref double eps)
549        {
550            int i = 0;
551            int ng = 0;
552            int[] p1 = new int[0];
553            int[] p2 = new int[0];
554            double tmp = 0;
555
556            System.Diagnostics.Debug.Assert(n==15 | n==21 | n==31 | n==41 | n==51 | n==61, "GKQNodesTbl: incorrect N!");
557            x = new double[n-1+1];
558            wkronrod = new double[n-1+1];
559            wgauss = new double[n-1+1];
560            for(i=0; i<=n-1; i++)
561            {
562                x[i] = 0;
563                wkronrod[i] = 0;
564                wgauss[i] = 0;
565            }
566            eps = Math.Max(AP.Math.MachineEpsilon, 1.0E-32);
567            if( n==15 )
568            {
569                ng = 4;
570                wgauss[0] = 0.129484966168869693270611432679082;
571                wgauss[1] = 0.279705391489276667901467771423780;
572                wgauss[2] = 0.381830050505118944950369775488975;
573                wgauss[3] = 0.417959183673469387755102040816327;
574                x[0] = 0.991455371120812639206854697526329;
575                x[1] = 0.949107912342758524526189684047851;
576                x[2] = 0.864864423359769072789712788640926;
577                x[3] = 0.741531185599394439863864773280788;
578                x[4] = 0.586087235467691130294144838258730;
579                x[5] = 0.405845151377397166906606412076961;
580                x[6] = 0.207784955007898467600689403773245;
581                x[7] = 0.000000000000000000000000000000000;
582                wkronrod[0] = 0.022935322010529224963732008058970;
583                wkronrod[1] = 0.063092092629978553290700663189204;
584                wkronrod[2] = 0.104790010322250183839876322541518;
585                wkronrod[3] = 0.140653259715525918745189590510238;
586                wkronrod[4] = 0.169004726639267902826583426598550;
587                wkronrod[5] = 0.190350578064785409913256402421014;
588                wkronrod[6] = 0.204432940075298892414161999234649;
589                wkronrod[7] = 0.209482141084727828012999174891714;
590            }
591            if( n==21 )
592            {
593                ng = 5;
594                wgauss[0] = 0.066671344308688137593568809893332;
595                wgauss[1] = 0.149451349150580593145776339657697;
596                wgauss[2] = 0.219086362515982043995534934228163;
597                wgauss[3] = 0.269266719309996355091226921569469;
598                wgauss[4] = 0.295524224714752870173892994651338;
599                x[0] = 0.995657163025808080735527280689003;
600                x[1] = 0.973906528517171720077964012084452;
601                x[2] = 0.930157491355708226001207180059508;
602                x[3] = 0.865063366688984510732096688423493;
603                x[4] = 0.780817726586416897063717578345042;
604                x[5] = 0.679409568299024406234327365114874;
605                x[6] = 0.562757134668604683339000099272694;
606                x[7] = 0.433395394129247190799265943165784;
607                x[8] = 0.294392862701460198131126603103866;
608                x[9] = 0.148874338981631210884826001129720;
609                x[10] = 0.000000000000000000000000000000000;
610                wkronrod[0] = 0.011694638867371874278064396062192;
611                wkronrod[1] = 0.032558162307964727478818972459390;
612                wkronrod[2] = 0.054755896574351996031381300244580;
613                wkronrod[3] = 0.075039674810919952767043140916190;
614                wkronrod[4] = 0.093125454583697605535065465083366;
615                wkronrod[5] = 0.109387158802297641899210590325805;
616                wkronrod[6] = 0.123491976262065851077958109831074;
617                wkronrod[7] = 0.134709217311473325928054001771707;
618                wkronrod[8] = 0.142775938577060080797094273138717;
619                wkronrod[9] = 0.147739104901338491374841515972068;
620                wkronrod[10] = 0.149445554002916905664936468389821;
621            }
622            if( n==31 )
623            {
624                ng = 8;
625                wgauss[0] = 0.030753241996117268354628393577204;
626                wgauss[1] = 0.070366047488108124709267416450667;
627                wgauss[2] = 0.107159220467171935011869546685869;
628                wgauss[3] = 0.139570677926154314447804794511028;
629                wgauss[4] = 0.166269205816993933553200860481209;
630                wgauss[5] = 0.186161000015562211026800561866423;
631                wgauss[6] = 0.198431485327111576456118326443839;
632                wgauss[7] = 0.202578241925561272880620199967519;
633                x[0] = 0.998002298693397060285172840152271;
634                x[1] = 0.987992518020485428489565718586613;
635                x[2] = 0.967739075679139134257347978784337;
636                x[3] = 0.937273392400705904307758947710209;
637                x[4] = 0.897264532344081900882509656454496;
638                x[5] = 0.848206583410427216200648320774217;
639                x[6] = 0.790418501442465932967649294817947;
640                x[7] = 0.724417731360170047416186054613938;
641                x[8] = 0.650996741297416970533735895313275;
642                x[9] = 0.570972172608538847537226737253911;
643                x[10] = 0.485081863640239680693655740232351;
644                x[11] = 0.394151347077563369897207370981045;
645                x[12] = 0.299180007153168812166780024266389;
646                x[13] = 0.201194093997434522300628303394596;
647                x[14] = 0.101142066918717499027074231447392;
648                x[15] = 0.000000000000000000000000000000000;
649                wkronrod[0] = 0.005377479872923348987792051430128;
650                wkronrod[1] = 0.015007947329316122538374763075807;
651                wkronrod[2] = 0.025460847326715320186874001019653;
652                wkronrod[3] = 0.035346360791375846222037948478360;
653                wkronrod[4] = 0.044589751324764876608227299373280;
654                wkronrod[5] = 0.053481524690928087265343147239430;
655                wkronrod[6] = 0.062009567800670640285139230960803;
656                wkronrod[7] = 0.069854121318728258709520077099147;
657                wkronrod[8] = 0.076849680757720378894432777482659;
658                wkronrod[9] = 0.083080502823133021038289247286104;
659                wkronrod[10] = 0.088564443056211770647275443693774;
660                wkronrod[11] = 0.093126598170825321225486872747346;
661                wkronrod[12] = 0.096642726983623678505179907627589;
662                wkronrod[13] = 0.099173598721791959332393173484603;
663                wkronrod[14] = 0.100769845523875595044946662617570;
664                wkronrod[15] = 0.101330007014791549017374792767493;
665            }
666            if( n==41 )
667            {
668                ng = 10;
669                wgauss[0] = 0.017614007139152118311861962351853;
670                wgauss[1] = 0.040601429800386941331039952274932;
671                wgauss[2] = 0.062672048334109063569506535187042;
672                wgauss[3] = 0.083276741576704748724758143222046;
673                wgauss[4] = 0.101930119817240435036750135480350;
674                wgauss[5] = 0.118194531961518417312377377711382;
675                wgauss[6] = 0.131688638449176626898494499748163;
676                wgauss[7] = 0.142096109318382051329298325067165;
677                wgauss[8] = 0.149172986472603746787828737001969;
678                wgauss[9] = 0.152753387130725850698084331955098;
679                x[0] = 0.998859031588277663838315576545863;
680                x[1] = 0.993128599185094924786122388471320;
681                x[2] = 0.981507877450250259193342994720217;
682                x[3] = 0.963971927277913791267666131197277;
683                x[4] = 0.940822633831754753519982722212443;
684                x[5] = 0.912234428251325905867752441203298;
685                x[6] = 0.878276811252281976077442995113078;
686                x[7] = 0.839116971822218823394529061701521;
687                x[8] = 0.795041428837551198350638833272788;
688                x[9] = 0.746331906460150792614305070355642;
689                x[10] = 0.693237656334751384805490711845932;
690                x[11] = 0.636053680726515025452836696226286;
691                x[12] = 0.575140446819710315342946036586425;
692                x[13] = 0.510867001950827098004364050955251;
693                x[14] = 0.443593175238725103199992213492640;
694                x[15] = 0.373706088715419560672548177024927;
695                x[16] = 0.301627868114913004320555356858592;
696                x[17] = 0.227785851141645078080496195368575;
697                x[18] = 0.152605465240922675505220241022678;
698                x[19] = 0.076526521133497333754640409398838;
699                x[20] = 0.000000000000000000000000000000000;
700                wkronrod[0] = 0.003073583718520531501218293246031;
701                wkronrod[1] = 0.008600269855642942198661787950102;
702                wkronrod[2] = 0.014626169256971252983787960308868;
703                wkronrod[3] = 0.020388373461266523598010231432755;
704                wkronrod[4] = 0.025882133604951158834505067096153;
705                wkronrod[5] = 0.031287306777032798958543119323801;
706                wkronrod[6] = 0.036600169758200798030557240707211;
707                wkronrod[7] = 0.041668873327973686263788305936895;
708                wkronrod[8] = 0.046434821867497674720231880926108;
709                wkronrod[9] = 0.050944573923728691932707670050345;
710                wkronrod[10] = 0.055195105348285994744832372419777;
711                wkronrod[11] = 0.059111400880639572374967220648594;
712                wkronrod[12] = 0.062653237554781168025870122174255;
713                wkronrod[13] = 0.065834597133618422111563556969398;
714                wkronrod[14] = 0.068648672928521619345623411885368;
715                wkronrod[15] = 0.071054423553444068305790361723210;
716                wkronrod[16] = 0.073030690332786667495189417658913;
717                wkronrod[17] = 0.074582875400499188986581418362488;
718                wkronrod[18] = 0.075704497684556674659542775376617;
719                wkronrod[19] = 0.076377867672080736705502835038061;
720                wkronrod[20] = 0.076600711917999656445049901530102;
721            }
722            if( n==51 )
723            {
724                ng = 13;
725                wgauss[0] = 0.011393798501026287947902964113235;
726                wgauss[1] = 0.026354986615032137261901815295299;
727                wgauss[2] = 0.040939156701306312655623487711646;
728                wgauss[3] = 0.054904695975835191925936891540473;
729                wgauss[4] = 0.068038333812356917207187185656708;
730                wgauss[5] = 0.080140700335001018013234959669111;
731                wgauss[6] = 0.091028261982963649811497220702892;
732                wgauss[7] = 0.100535949067050644202206890392686;
733                wgauss[8] = 0.108519624474263653116093957050117;
734                wgauss[9] = 0.114858259145711648339325545869556;
735                wgauss[10] = 0.119455763535784772228178126512901;
736                wgauss[11] = 0.122242442990310041688959518945852;
737                wgauss[12] = 0.123176053726715451203902873079050;
738                x[0] = 0.999262104992609834193457486540341;
739                x[1] = 0.995556969790498097908784946893902;
740                x[2] = 0.988035794534077247637331014577406;
741                x[3] = 0.976663921459517511498315386479594;
742                x[4] = 0.961614986425842512418130033660167;
743                x[5] = 0.942974571228974339414011169658471;
744                x[6] = 0.920747115281701561746346084546331;
745                x[7] = 0.894991997878275368851042006782805;
746                x[8] = 0.865847065293275595448996969588340;
747                x[9] = 0.833442628760834001421021108693570;
748                x[10] = 0.797873797998500059410410904994307;
749                x[11] = 0.759259263037357630577282865204361;
750                x[12] = 0.717766406813084388186654079773298;
751                x[13] = 0.673566368473468364485120633247622;
752                x[14] = 0.626810099010317412788122681624518;
753                x[15] = 0.577662930241222967723689841612654;
754                x[16] = 0.526325284334719182599623778158010;
755                x[17] = 0.473002731445714960522182115009192;
756                x[18] = 0.417885382193037748851814394594572;
757                x[19] = 0.361172305809387837735821730127641;
758                x[20] = 0.303089538931107830167478909980339;
759                x[21] = 0.243866883720988432045190362797452;
760                x[22] = 0.183718939421048892015969888759528;
761                x[23] = 0.122864692610710396387359818808037;
762                x[24] = 0.061544483005685078886546392366797;
763                x[25] = 0.000000000000000000000000000000000;
764                wkronrod[0] = 0.001987383892330315926507851882843;
765                wkronrod[1] = 0.005561932135356713758040236901066;
766                wkronrod[2] = 0.009473973386174151607207710523655;
767                wkronrod[3] = 0.013236229195571674813656405846976;
768                wkronrod[4] = 0.016847817709128298231516667536336;
769                wkronrod[5] = 0.020435371145882835456568292235939;
770                wkronrod[6] = 0.024009945606953216220092489164881;
771                wkronrod[7] = 0.027475317587851737802948455517811;
772                wkronrod[8] = 0.030792300167387488891109020215229;
773                wkronrod[9] = 0.034002130274329337836748795229551;
774                wkronrod[10] = 0.037116271483415543560330625367620;
775                wkronrod[11] = 0.040083825504032382074839284467076;
776                wkronrod[12] = 0.042872845020170049476895792439495;
777                wkronrod[13] = 0.045502913049921788909870584752660;
778                wkronrod[14] = 0.047982537138836713906392255756915;
779                wkronrod[15] = 0.050277679080715671963325259433440;
780                wkronrod[16] = 0.052362885806407475864366712137873;
781                wkronrod[17] = 0.054251129888545490144543370459876;
782                wkronrod[18] = 0.055950811220412317308240686382747;
783                wkronrod[19] = 0.057437116361567832853582693939506;
784                wkronrod[20] = 0.058689680022394207961974175856788;
785                wkronrod[21] = 0.059720340324174059979099291932562;
786                wkronrod[22] = 0.060539455376045862945360267517565;
787                wkronrod[23] = 0.061128509717053048305859030416293;
788                wkronrod[24] = 0.061471189871425316661544131965264;
789                wkronrod[25] = 0.061580818067832935078759824240055;
790            }
791            if( n==61 )
792            {
793                ng = 15;
794                wgauss[0] = 0.007968192496166605615465883474674;
795                wgauss[1] = 0.018466468311090959142302131912047;
796                wgauss[2] = 0.028784707883323369349719179611292;
797                wgauss[3] = 0.038799192569627049596801936446348;
798                wgauss[4] = 0.048402672830594052902938140422808;
799                wgauss[5] = 0.057493156217619066481721689402056;
800                wgauss[6] = 0.065974229882180495128128515115962;
801                wgauss[7] = 0.073755974737705206268243850022191;
802                wgauss[8] = 0.080755895229420215354694938460530;
803                wgauss[9] = 0.086899787201082979802387530715126;
804                wgauss[10] = 0.092122522237786128717632707087619;
805                wgauss[11] = 0.096368737174644259639468626351810;
806                wgauss[12] = 0.099593420586795267062780282103569;
807                wgauss[13] = 0.101762389748405504596428952168554;
808                wgauss[14] = 0.102852652893558840341285636705415;
809                x[0] = 0.999484410050490637571325895705811;
810                x[1] = 0.996893484074649540271630050918695;
811                x[2] = 0.991630996870404594858628366109486;
812                x[3] = 0.983668123279747209970032581605663;
813                x[4] = 0.973116322501126268374693868423707;
814                x[5] = 0.960021864968307512216871025581798;
815                x[6] = 0.944374444748559979415831324037439;
816                x[7] = 0.926200047429274325879324277080474;
817                x[8] = 0.905573307699907798546522558925958;
818                x[9] = 0.882560535792052681543116462530226;
819                x[10] = 0.857205233546061098958658510658944;
820                x[11] = 0.829565762382768397442898119732502;
821                x[12] = 0.799727835821839083013668942322683;
822                x[13] = 0.767777432104826194917977340974503;
823                x[14] = 0.733790062453226804726171131369528;
824                x[15] = 0.697850494793315796932292388026640;
825                x[16] = 0.660061064126626961370053668149271;
826                x[17] = 0.620526182989242861140477556431189;
827                x[18] = 0.579345235826361691756024932172540;
828                x[19] = 0.536624148142019899264169793311073;
829                x[20] = 0.492480467861778574993693061207709;
830                x[21] = 0.447033769538089176780609900322854;
831                x[22] = 0.400401254830394392535476211542661;
832                x[23] = 0.352704725530878113471037207089374;
833                x[24] = 0.304073202273625077372677107199257;
834                x[25] = 0.254636926167889846439805129817805;
835                x[26] = 0.204525116682309891438957671002025;
836                x[27] = 0.153869913608583546963794672743256;
837                x[28] = 0.102806937966737030147096751318001;
838                x[29] = 0.051471842555317695833025213166723;
839                x[30] = 0.000000000000000000000000000000000;
840                wkronrod[0] = 0.001389013698677007624551591226760;
841                wkronrod[1] = 0.003890461127099884051267201844516;
842                wkronrod[2] = 0.006630703915931292173319826369750;
843                wkronrod[3] = 0.009273279659517763428441146892024;
844                wkronrod[4] = 0.011823015253496341742232898853251;
845                wkronrod[5] = 0.014369729507045804812451432443580;
846                wkronrod[6] = 0.016920889189053272627572289420322;
847                wkronrod[7] = 0.019414141193942381173408951050128;
848                wkronrod[8] = 0.021828035821609192297167485738339;
849                wkronrod[9] = 0.024191162078080601365686370725232;
850                wkronrod[10] = 0.026509954882333101610601709335075;
851                wkronrod[11] = 0.028754048765041292843978785354334;
852                wkronrod[12] = 0.030907257562387762472884252943092;
853                wkronrod[13] = 0.032981447057483726031814191016854;
854                wkronrod[14] = 0.034979338028060024137499670731468;
855                wkronrod[15] = 0.036882364651821229223911065617136;
856                wkronrod[16] = 0.038678945624727592950348651532281;
857                wkronrod[17] = 0.040374538951535959111995279752468;
858                wkronrod[18] = 0.041969810215164246147147541285970;
859                wkronrod[19] = 0.043452539701356069316831728117073;
860                wkronrod[20] = 0.044814800133162663192355551616723;
861                wkronrod[21] = 0.046059238271006988116271735559374;
862                wkronrod[22] = 0.047185546569299153945261478181099;
863                wkronrod[23] = 0.048185861757087129140779492298305;
864                wkronrod[24] = 0.049055434555029778887528165367238;
865                wkronrod[25] = 0.049795683427074206357811569379942;
866                wkronrod[26] = 0.050405921402782346840893085653585;
867                wkronrod[27] = 0.050881795898749606492297473049805;
868                wkronrod[28] = 0.051221547849258772170656282604944;
869                wkronrod[29] = 0.051426128537459025933862879215781;
870                wkronrod[30] = 0.051494729429451567558340433647099;
871            }
872           
873            //
874            // copy nodes
875            //
876            for(i=n-1; i>=n/2; i--)
877            {
878                x[i] = -x[n-1-i];
879            }
880           
881            //
882            // copy Kronrod weights
883            //
884            for(i=n-1; i>=n/2; i--)
885            {
886                wkronrod[i] = wkronrod[n-1-i];
887            }
888           
889            //
890            // copy Gauss weights
891            //
892            for(i=ng-1; i>=0; i--)
893            {
894                wgauss[n-2-2*i] = wgauss[i];
895                wgauss[1+2*i] = wgauss[i];
896            }
897            for(i=0; i<=n/2; i++)
898            {
899                wgauss[2*i] = 0;
900            }
901           
902            //
903            // reorder
904            //
905            tsort.tagsort(ref x, n, ref p1, ref p2);
906            for(i=0; i<=n-1; i++)
907            {
908                tmp = wkronrod[i];
909                wkronrod[i] = wkronrod[p2[i]];
910                wkronrod[p2[i]] = tmp;
911                tmp = wgauss[i];
912                wgauss[i] = wgauss[p2[i]];
913                wgauss[p2[i]] = tmp;
914            }
915        }
916    }
917}
Note: See TracBrowser for help on using the repository browser.