source: trunk/sources/ALGLIB/gkq.cs @ 2625

/*************************************************************************
Copyright (c) 2005-2009, Sergey Bochkanov (ALGLIB project).

>>> SOURCE LICENSE >>>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation (www.fsf.org); either version 2 of the 
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License is available at
http://www.fsf.org/licensing/licenses

>>> END OF LICENSE >>>
*************************************************************************/

using System;

namespace alglib
{
    public class gkq
    {
        /*************************************************************************
        Computation of nodes and weights of a Gauss-Kronrod quadrature formula

        The algorithm generates the N-point Gauss-Kronrod quadrature formula  with
        weight  function  given  by  coefficients  alpha  and beta of a recurrence
        relation which generates a system of orthogonal polynomials:

            P-1(x)   =  0
            P0(x)    =  1
            Pn+1(x)  =  (x-alpha(n))*Pn(x)  -  beta(n)*Pn-1(x)

        and zero moment Mu0

            Mu0 = integral(W(x)dx,a,b)


        INPUT PARAMETERS:
            Alpha       –   alpha coefficients, array[0..floor(3*K/2)].
            Beta        –   beta coefficients,  array[0..ceil(3*K/2)].
                            Beta[0] is not used and may be arbitrary.
                            Beta[I]>0.
            Mu0         –   zeroth moment of the weight function.
            N           –   number of nodes of the Gauss-Kronrod quadrature formula,
                            N >= 3,
                            N =  2*K+1.

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -5    no real and positive Gauss-Kronrod formula can
                                    be created for such a weight function  with  a
                                    given number of nodes.
                            * -4    N is too large, task may be ill  conditioned -
                                    x[i]=x[i+1] found.
                            * -3    internal eigenproblem solver hasn't converged
                            * -2    Beta[i]<=0
                            * -1    incorrect N was passed
                            * +1    OK
            X           -   array[0..N-1] - array of quadrature nodes,
                            in ascending order.
            WKronrod    -   array[0..N-1] - Kronrod weights
            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
                            corresponding to extended Kronrod nodes).

          -- ALGLIB --
             Copyright 08.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gkqgeneraterec(double[] alpha,
            double[] beta,
            double mu0,
            int n,
            ref int info,
            ref double[] x,
            ref double[] wkronrod,
            ref double[] wgauss)
        {
            double[] ta = new double[0];
            int i = 0;
            int j = 0;
            double[] t = new double[0];
            double[] s = new double[0];
            int wlen = 0;
            int woffs = 0;
            double u = 0;
            int m = 0;
            int l = 0;
            int k = 0;
            double[] xgtmp = new double[0];
            double[] wgtmp = new double[0];
            int i_ = 0;

            alpha = (double[])alpha.Clone();
            beta = (double[])beta.Clone();

            if( n%2!=1 | n<3 )
            {
                info = -1;
                return;
            }
            for(i=0; i<=(int)Math.Ceiling((double)(3*(n/2))/(double)(2)); i++)
            {
                if( (double)(beta[i])<=(double)(0) )
                {
                    info = -2;
                    return;
                }
            }
            info = 1;
            
            //
            // from external conventions about N/Beta/Mu0 to internal
            //
            n = n/2;
            beta[0] = mu0;
            
            //
            // Calculate Gauss nodes/weights, save them for later processing
            //
            gq.gqgeneraterec(ref alpha, ref beta, mu0, n, ref info, ref xgtmp, ref wgtmp);
            if( info<0 )
            {
                return;
            }
            
            //
            // Resize:
            // * A from 0..floor(3*n/2) to 0..2*n
            // * B from 0..ceil(3*n/2)  to 0..2*n
            //
            ta = new double[(int)Math.Floor((double)(3*n)/(double)(2))+1];
            for(i_=0; i_<=(int)Math.Floor((double)(3*n)/(double)(2));i_++)
            {
                ta[i_] = alpha[i_];
            }
            alpha = new double[2*n+1];
            for(i_=0; i_<=(int)Math.Floor((double)(3*n)/(double)(2));i_++)
            {
                alpha[i_] = ta[i_];
            }
            for(i=(int)Math.Floor((double)(3*n)/(double)(2))+1; i<=2*n; i++)
            {
                alpha[i] = 0;
            }
            ta = new double[(int)Math.Ceiling((double)(3*n)/(double)(2))+1];
            for(i_=0; i_<=(int)Math.Ceiling((double)(3*n)/(double)(2));i_++)
            {
                ta[i_] = beta[i_];
            }
            beta = new double[2*n+1];
            for(i_=0; i_<=(int)Math.Ceiling((double)(3*n)/(double)(2));i_++)
            {
                beta[i_] = ta[i_];
            }
            for(i=(int)Math.Ceiling((double)(3*n)/(double)(2))+1; i<=2*n; i++)
            {
                beta[i] = 0;
            }
            
            //
            // Initialize T, S
            //
            wlen = 2+n/2;
            t = new double[wlen];
            s = new double[wlen];
            ta = new double[wlen];
            woffs = 1;
            for(i=0; i<=wlen-1; i++)
            {
                t[i] = 0;
                s[i] = 0;
            }
            
            //
            // Algorithm from Dirk P. Laurie, "Calculation of Gauss-Kronrod quadrature rules", 1997.
            //
            t[woffs+0] = beta[n+1];
            for(m=0; m<=n-2; m++)
            {
                u = 0;
                for(k=(m+1)/2; k>=0; k--)
                {
                    l = m-k;
                    u = u+(alpha[k+n+1]-alpha[l])*t[woffs+k]+beta[k+n+1]*s[woffs+k-1]-beta[l]*s[woffs+k];
                    s[woffs+k] = u;
                }
                for(i_=0; i_<=wlen-1;i_++)
                {
                    ta[i_] = t[i_];
                }
                for(i_=0; i_<=wlen-1;i_++)
                {
                    t[i_] = s[i_];
                }
                for(i_=0; i_<=wlen-1;i_++)
                {
                    s[i_] = ta[i_];
                }
            }
            for(j=n/2; j>=0; j--)
            {
                s[woffs+j] = s[woffs+j-1];
            }
            for(m=n-1; m<=2*n-3; m++)
            {
                u = 0;
                for(k=m+1-n; k<=(m-1)/2; k++)
                {
                    l = m-k;
                    j = n-1-l;
                    u = u-(alpha[k+n+1]-alpha[l])*t[woffs+j]-beta[k+n+1]*s[woffs+j]+beta[l]*s[woffs+j+1];
                    s[woffs+j] = u;
                }
                if( m%2==0 )
                {
                    k = m/2;
                    alpha[k+n+1] = alpha[k]+(s[woffs+j]-beta[k+n+1]*s[woffs+j+1])/t[woffs+j+1];
                }
                else
                {
                    k = (m+1)/2;
                    beta[k+n+1] = s[woffs+j]/s[woffs+j+1];
                }
                for(i_=0; i_<=wlen-1;i_++)
                {
                    ta[i_] = t[i_];
                }
                for(i_=0; i_<=wlen-1;i_++)
                {
                    t[i_] = s[i_];
                }
                for(i_=0; i_<=wlen-1;i_++)
                {
                    s[i_] = ta[i_];
                }
            }
            alpha[2*n] = alpha[n-1]-beta[2*n]*s[woffs+0]/t[woffs+0];
            
            //
            // calculation of Kronrod nodes and weights, unpacking of Gauss weights
            //
            gq.gqgeneraterec(ref alpha, ref beta, mu0, 2*n+1, ref info, ref x, ref wkronrod);
            if( info==-2 )
            {
                info = -5;
            }
            if( info<0 )
            {
                return;
            }
            for(i=0; i<=2*n-1; i++)
            {
                if( (double)(x[i])>=(double)(x[i+1]) )
                {
                    info = -4;
                }
            }
            if( info<0 )
            {
                return;
            }
            wgauss = new double[2*n+1];
            for(i=0; i<=2*n; i++)
            {
                wgauss[i] = 0;
            }
            for(i=0; i<=n-1; i++)
            {
                wgauss[2*i+1] = wgtmp[i];
            }
        }


        /*************************************************************************
        Returns   Gauss   and   Gauss-Kronrod   nodes/weights  for  Gauss-Legendre
        quadrature with N points.

        GKQLegendreCalc (calculation) or  GKQLegendreTbl  (precomputed  table)  is
        used depending on machine precision and number of nodes.

        INPUT PARAMETERS:
            N           -   number of Kronrod nodes, must be odd number, >=3.

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -4    an  error   was   detected   when  calculating
                                    weights/nodes.  N  is  too  large   to  obtain
                                    weights/nodes  with  high   enough   accuracy.
                                    Try  to   use   multiple   precision  version.
                            * -3    internal eigenproblem solver hasn't converged
                            * -1    incorrect N was passed
                            * +1    OK
            X           -   array[0..N-1] - array of quadrature nodes, ordered in
                            ascending order.
            WKronrod    -   array[0..N-1] - Kronrod weights
            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
                            corresponding to extended Kronrod nodes).


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gkqgenerategausslegendre(int n,
            ref int info,
            ref double[] x,
            ref double[] wkronrod,
            ref double[] wgauss)
        {
            double eps = 0;

            if( (double)(AP.Math.MachineEpsilon)>(double)(1.0E-32) & (n==15 | n==21 | n==31 | n==41 | n==51 | n==61) )
            {
                info = 1;
                gkqlegendretbl(n, ref x, ref wkronrod, ref wgauss, ref eps);
            }
            else
            {
                gkqlegendrecalc(n, ref info, ref x, ref wkronrod, ref wgauss);
            }
        }


        /*************************************************************************
        Returns   Gauss   and   Gauss-Kronrod   nodes/weights   for   Gauss-Jacobi
        quadrature on [-1,1] with weight function

            W(x)=Power(1-x,Alpha)*Power(1+x,Beta).

        INPUT PARAMETERS:
            N           -   number of Kronrod nodes, must be odd number, >=3.
            Alpha       -   power-law coefficient, Alpha>-1
            Beta        -   power-law coefficient, Beta>-1

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -5    no real and positive Gauss-Kronrod formula can
                                    be created for such a weight function  with  a
                                    given number of nodes.
                            * -4    an  error  was   detected   when   calculating
                                    weights/nodes. Alpha or  Beta  are  too  close
                                    to -1 to obtain weights/nodes with high enough
                                    accuracy, or, may be, N is too large.  Try  to
                                    use multiple precision version.
                            * -3    internal eigenproblem solver hasn't converged
                            * -1    incorrect N was passed
                            * +1    OK
                            * +2    OK, but quadrature rule have exterior  nodes,
                                    x[0]<-1 or x[n-1]>+1
            X           -   array[0..N-1] - array of quadrature nodes, ordered in
                            ascending order.
            WKronrod    -   array[0..N-1] - Kronrod weights
            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
                            corresponding to extended Kronrod nodes).


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gkqgenerategaussjacobi(int n,
            double alpha,
            double beta,
            ref int info,
            ref double[] x,
            ref double[] wkronrod,
            ref double[] wgauss)
        {
            int clen = 0;
            double[] a = new double[0];
            double[] b = new double[0];
            double alpha2 = 0;
            double beta2 = 0;
            double apb = 0;
            double t = 0;
            int i = 0;
            double s = 0;

            if( n%2!=1 | n<3 )
            {
                info = -1;
                return;
            }
            if( (double)(alpha)<=(double)(-1) | (double)(beta)<=(double)(-1) )
            {
                info = -1;
                return;
            }
            clen = (int)Math.Ceiling((double)(3*(n/2))/(double)(2))+1;
            a = new double[clen];
            b = new double[clen];
            for(i=0; i<=clen-1; i++)
            {
                a[i] = 0;
            }
            apb = alpha+beta;
            a[0] = (beta-alpha)/(apb+2);
            t = (apb+1)*Math.Log(2)+gammafunc.lngamma(alpha+1, ref s)+gammafunc.lngamma(beta+1, ref s)-gammafunc.lngamma(apb+2, ref s);
            if( (double)(t)>(double)(Math.Log(AP.Math.MaxRealNumber)) )
            {
                info = -4;
                return;
            }
            b[0] = Math.Exp(t);
            if( clen>1 )
            {
                alpha2 = AP.Math.Sqr(alpha);
                beta2 = AP.Math.Sqr(beta);
                a[1] = (beta2-alpha2)/((apb+2)*(apb+4));
                b[1] = 4*(alpha+1)*(beta+1)/((apb+3)*AP.Math.Sqr(apb+2));
                for(i=2; i<=clen-1; i++)
                {
                    a[i] = 0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
                    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));
                }
            }
            gkqgeneraterec(a, b, b[0], n, ref info, ref x, ref wkronrod, ref wgauss);
            
            //
            // test basic properties to detect errors
            //
            if( info>0 )
            {
                if( (double)(x[0])<(double)(-1) | (double)(x[n-1])>(double)(+1) )
                {
                    info = 2;
                }
                for(i=0; i<=n-2; i++)
                {
                    if( (double)(x[i])>=(double)(x[i+1]) )
                    {
                        info = -4;
                    }
                }
            }
        }


        /*************************************************************************
        Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.

        Reduction to tridiagonal eigenproblem is used.

        INPUT PARAMETERS:
            N           -   number of Kronrod nodes, must be odd number, >=3.

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -4    an  error   was   detected   when  calculating
                                    weights/nodes.  N  is  too  large   to  obtain
                                    weights/nodes  with  high   enough   accuracy.
                                    Try  to   use   multiple   precision  version.
                            * -3    internal eigenproblem solver hasn't converged
                            * -1    incorrect N was passed
                            * +1    OK
            X           -   array[0..N-1] - array of quadrature nodes, ordered in
                            ascending order.
            WKronrod    -   array[0..N-1] - Kronrod weights
            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
                            corresponding to extended Kronrod nodes).

          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gkqlegendrecalc(int n,
            ref int info,
            ref double[] x,
            ref double[] wkronrod,
            ref double[] wgauss)
        {
            double[] alpha = new double[0];
            double[] beta = new double[0];
            int alen = 0;
            int blen = 0;
            double mu0 = 0;
            int k = 0;
            int i = 0;

            if( n%2!=1 | n<3 )
            {
                info = -1;
                return;
            }
            mu0 = 2;
            alen = (int)Math.Floor((double)(3*(n/2))/(double)(2))+1;
            blen = (int)Math.Ceiling((double)(3*(n/2))/(double)(2))+1;
            alpha = new double[alen];
            beta = new double[blen];
            for(k=0; k<=alen-1; k++)
            {
                alpha[k] = 0;
            }
            beta[0] = 2;
            for(k=1; k<=blen-1; k++)
            {
                beta[k] = 1/(4-1/AP.Math.Sqr(k));
            }
            gkqgeneraterec(alpha, beta, mu0, n, ref info, ref x, ref wkronrod, ref wgauss);
            
            //
            // test basic properties to detect errors
            //
            if( info>0 )
            {
                if( (double)(x[0])<(double)(-1) | (double)(x[n-1])>(double)(+1) )
                {
                    info = -4;
                }
                for(i=0; i<=n-2; i++)
                {
                    if( (double)(x[i])>=(double)(x[i+1]) )
                    {
                        info = -4;
                    }
                }
            }
        }


        /*************************************************************************
        Returns Gauss and Gauss-Kronrod nodes for quadrature with N  points  using
        pre-calculated table. Nodes/weights were  computed  with  accuracy  up  to
        1.0E-32 (if MPFR version of ALGLIB is used). In standard double  precision
        accuracy reduces to something about 2.0E-16 (depending  on your compiler's
        handling of long floating point constants).

        INPUT PARAMETERS:
            N           -   number of Kronrod nodes.
                            N can be 15, 21, 31, 41, 51, 61.

        OUTPUT PARAMETERS:
            X           -   array[0..N-1] - array of quadrature nodes, ordered in
                            ascending order.
            WKronrod    -   array[0..N-1] - Kronrod weights
            WGauss      -   array[0..N-1] - Gauss weights (interleaved with zeros
                            corresponding to extended Kronrod nodes).


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gkqlegendretbl(int n,
            ref double[] x,
            ref double[] wkronrod,
            ref double[] wgauss,
            ref double eps)
        {
            int i = 0;
            int ng = 0;
            int[] p1 = new int[0];
            int[] p2 = new int[0];
            double tmp = 0;

            System.Diagnostics.Debug.Assert(n==15 | n==21 | n==31 | n==41 | n==51 | n==61, "GKQNodesTbl: incorrect N!");
            x = new double[n-1+1];
            wkronrod = new double[n-1+1];
            wgauss = new double[n-1+1];
            for(i=0; i<=n-1; i++)
            {
                x[i] = 0;
                wkronrod[i] = 0;
                wgauss[i] = 0;
            }
            eps = Math.Max(AP.Math.MachineEpsilon, 1.0E-32);
            if( n==15 )
            {
                ng = 4;
                wgauss[0] = 0.129484966168869693270611432679082;
                wgauss[1] = 0.279705391489276667901467771423780;
                wgauss[2] = 0.381830050505118944950369775488975;
                wgauss[3] = 0.417959183673469387755102040816327;
                x[0] = 0.991455371120812639206854697526329;
                x[1] = 0.949107912342758524526189684047851;
                x[2] = 0.864864423359769072789712788640926;
                x[3] = 0.741531185599394439863864773280788;
                x[4] = 0.586087235467691130294144838258730;
                x[5] = 0.405845151377397166906606412076961;
                x[6] = 0.207784955007898467600689403773245;
                x[7] = 0.000000000000000000000000000000000;
                wkronrod[0] = 0.022935322010529224963732008058970;
                wkronrod[1] = 0.063092092629978553290700663189204;
                wkronrod[2] = 0.104790010322250183839876322541518;
                wkronrod[3] = 0.140653259715525918745189590510238;
                wkronrod[4] = 0.169004726639267902826583426598550;
                wkronrod[5] = 0.190350578064785409913256402421014;
                wkronrod[6] = 0.204432940075298892414161999234649;
                wkronrod[7] = 0.209482141084727828012999174891714;
            }
            if( n==21 )
            {
                ng = 5;
                wgauss[0] = 0.066671344308688137593568809893332;
                wgauss[1] = 0.149451349150580593145776339657697;
                wgauss[2] = 0.219086362515982043995534934228163;
                wgauss[3] = 0.269266719309996355091226921569469;
                wgauss[4] = 0.295524224714752870173892994651338;
                x[0] = 0.995657163025808080735527280689003;
                x[1] = 0.973906528517171720077964012084452;
                x[2] = 0.930157491355708226001207180059508;
                x[3] = 0.865063366688984510732096688423493;
                x[4] = 0.780817726586416897063717578345042;
                x[5] = 0.679409568299024406234327365114874;
                x[6] = 0.562757134668604683339000099272694;
                x[7] = 0.433395394129247190799265943165784;
                x[8] = 0.294392862701460198131126603103866;
                x[9] = 0.148874338981631210884826001129720;
                x[10] = 0.000000000000000000000000000000000;
                wkronrod[0] = 0.011694638867371874278064396062192;
                wkronrod[1] = 0.032558162307964727478818972459390;
                wkronrod[2] = 0.054755896574351996031381300244580;
                wkronrod[3] = 0.075039674810919952767043140916190;
                wkronrod[4] = 0.093125454583697605535065465083366;
                wkronrod[5] = 0.109387158802297641899210590325805;
                wkronrod[6] = 0.123491976262065851077958109831074;
                wkronrod[7] = 0.134709217311473325928054001771707;
                wkronrod[8] = 0.142775938577060080797094273138717;
                wkronrod[9] = 0.147739104901338491374841515972068;
                wkronrod[10] = 0.149445554002916905664936468389821;
            }
            if( n==31 )
            {
                ng = 8;
                wgauss[0] = 0.030753241996117268354628393577204;
                wgauss[1] = 0.070366047488108124709267416450667;
                wgauss[2] = 0.107159220467171935011869546685869;
                wgauss[3] = 0.139570677926154314447804794511028;
                wgauss[4] = 0.166269205816993933553200860481209;
                wgauss[5] = 0.186161000015562211026800561866423;
                wgauss[6] = 0.198431485327111576456118326443839;
                wgauss[7] = 0.202578241925561272880620199967519;
                x[0] = 0.998002298693397060285172840152271;
                x[1] = 0.987992518020485428489565718586613;
                x[2] = 0.967739075679139134257347978784337;
                x[3] = 0.937273392400705904307758947710209;
                x[4] = 0.897264532344081900882509656454496;
                x[5] = 0.848206583410427216200648320774217;
                x[6] = 0.790418501442465932967649294817947;
                x[7] = 0.724417731360170047416186054613938;
                x[8] = 0.650996741297416970533735895313275;
                x[9] = 0.570972172608538847537226737253911;
                x[10] = 0.485081863640239680693655740232351;
                x[11] = 0.394151347077563369897207370981045;
                x[12] = 0.299180007153168812166780024266389;
                x[13] = 0.201194093997434522300628303394596;
                x[14] = 0.101142066918717499027074231447392;
                x[15] = 0.000000000000000000000000000000000;
                wkronrod[0] = 0.005377479872923348987792051430128;
                wkronrod[1] = 0.015007947329316122538374763075807;
                wkronrod[2] = 0.025460847326715320186874001019653;
                wkronrod[3] = 0.035346360791375846222037948478360;
                wkronrod[4] = 0.044589751324764876608227299373280;
                wkronrod[5] = 0.053481524690928087265343147239430;
                wkronrod[6] = 0.062009567800670640285139230960803;
                wkronrod[7] = 0.069854121318728258709520077099147;
                wkronrod[8] = 0.076849680757720378894432777482659;
                wkronrod[9] = 0.083080502823133021038289247286104;
                wkronrod[10] = 0.088564443056211770647275443693774;
                wkronrod[11] = 0.093126598170825321225486872747346;
                wkronrod[12] = 0.096642726983623678505179907627589;
                wkronrod[13] = 0.099173598721791959332393173484603;
                wkronrod[14] = 0.100769845523875595044946662617570;
                wkronrod[15] = 0.101330007014791549017374792767493;
            }
            if( n==41 )
            {
                ng = 10;
                wgauss[0] = 0.017614007139152118311861962351853;
                wgauss[1] = 0.040601429800386941331039952274932;
                wgauss[2] = 0.062672048334109063569506535187042;
                wgauss[3] = 0.083276741576704748724758143222046;
                wgauss[4] = 0.101930119817240435036750135480350;
                wgauss[5] = 0.118194531961518417312377377711382;
                wgauss[6] = 0.131688638449176626898494499748163;
                wgauss[7] = 0.142096109318382051329298325067165;
                wgauss[8] = 0.149172986472603746787828737001969;
                wgauss[9] = 0.152753387130725850698084331955098;
                x[0] = 0.998859031588277663838315576545863;
                x[1] = 0.993128599185094924786122388471320;
                x[2] = 0.981507877450250259193342994720217;
                x[3] = 0.963971927277913791267666131197277;
                x[4] = 0.940822633831754753519982722212443;
                x[5] = 0.912234428251325905867752441203298;
                x[6] = 0.878276811252281976077442995113078;
                x[7] = 0.839116971822218823394529061701521;
                x[8] = 0.795041428837551198350638833272788;
                x[9] = 0.746331906460150792614305070355642;
                x[10] = 0.693237656334751384805490711845932;
                x[11] = 0.636053680726515025452836696226286;
                x[12] = 0.575140446819710315342946036586425;
                x[13] = 0.510867001950827098004364050955251;
                x[14] = 0.443593175238725103199992213492640;
                x[15] = 0.373706088715419560672548177024927;
                x[16] = 0.301627868114913004320555356858592;
                x[17] = 0.227785851141645078080496195368575;
                x[18] = 0.152605465240922675505220241022678;
                x[19] = 0.076526521133497333754640409398838;
                x[20] = 0.000000000000000000000000000000000;
                wkronrod[0] = 0.003073583718520531501218293246031;
                wkronrod[1] = 0.008600269855642942198661787950102;
                wkronrod[2] = 0.014626169256971252983787960308868;
                wkronrod[3] = 0.020388373461266523598010231432755;
                wkronrod[4] = 0.025882133604951158834505067096153;
                wkronrod[5] = 0.031287306777032798958543119323801;
                wkronrod[6] = 0.036600169758200798030557240707211;
                wkronrod[7] = 0.041668873327973686263788305936895;
                wkronrod[8] = 0.046434821867497674720231880926108;
                wkronrod[9] = 0.050944573923728691932707670050345;
                wkronrod[10] = 0.055195105348285994744832372419777;
                wkronrod[11] = 0.059111400880639572374967220648594;
                wkronrod[12] = 0.062653237554781168025870122174255;
                wkronrod[13] = 0.065834597133618422111563556969398;
                wkronrod[14] = 0.068648672928521619345623411885368;
                wkronrod[15] = 0.071054423553444068305790361723210;
                wkronrod[16] = 0.073030690332786667495189417658913;
                wkronrod[17] = 0.074582875400499188986581418362488;
                wkronrod[18] = 0.075704497684556674659542775376617;
                wkronrod[19] = 0.076377867672080736705502835038061;
                wkronrod[20] = 0.076600711917999656445049901530102;
            }
            if( n==51 )
            {
                ng = 13;
                wgauss[0] = 0.011393798501026287947902964113235;
                wgauss[1] = 0.026354986615032137261901815295299;
                wgauss[2] = 0.040939156701306312655623487711646;
                wgauss[3] = 0.054904695975835191925936891540473;
                wgauss[4] = 0.068038333812356917207187185656708;
                wgauss[5] = 0.080140700335001018013234959669111;
                wgauss[6] = 0.091028261982963649811497220702892;
                wgauss[7] = 0.100535949067050644202206890392686;
                wgauss[8] = 0.108519624474263653116093957050117;
                wgauss[9] = 0.114858259145711648339325545869556;
                wgauss[10] = 0.119455763535784772228178126512901;
                wgauss[11] = 0.122242442990310041688959518945852;
                wgauss[12] = 0.123176053726715451203902873079050;
                x[0] = 0.999262104992609834193457486540341;
                x[1] = 0.995556969790498097908784946893902;
                x[2] = 0.988035794534077247637331014577406;
                x[3] = 0.976663921459517511498315386479594;
                x[4] = 0.961614986425842512418130033660167;
                x[5] = 0.942974571228974339414011169658471;
                x[6] = 0.920747115281701561746346084546331;
                x[7] = 0.894991997878275368851042006782805;
                x[8] = 0.865847065293275595448996969588340;
                x[9] = 0.833442628760834001421021108693570;
                x[10] = 0.797873797998500059410410904994307;
                x[11] = 0.759259263037357630577282865204361;
                x[12] = 0.717766406813084388186654079773298;
                x[13] = 0.673566368473468364485120633247622;
                x[14] = 0.626810099010317412788122681624518;
                x[15] = 0.577662930241222967723689841612654;
                x[16] = 0.526325284334719182599623778158010;
                x[17] = 0.473002731445714960522182115009192;
                x[18] = 0.417885382193037748851814394594572;
                x[19] = 0.361172305809387837735821730127641;
                x[20] = 0.303089538931107830167478909980339;
                x[21] = 0.243866883720988432045190362797452;
                x[22] = 0.183718939421048892015969888759528;
                x[23] = 0.122864692610710396387359818808037;
                x[24] = 0.061544483005685078886546392366797;
                x[25] = 0.000000000000000000000000000000000;
                wkronrod[0] = 0.001987383892330315926507851882843;
                wkronrod[1] = 0.005561932135356713758040236901066;
                wkronrod[2] = 0.009473973386174151607207710523655;
                wkronrod[3] = 0.013236229195571674813656405846976;
                wkronrod[4] = 0.016847817709128298231516667536336;
                wkronrod[5] = 0.020435371145882835456568292235939;
                wkronrod[6] = 0.024009945606953216220092489164881;
                wkronrod[7] = 0.027475317587851737802948455517811;
                wkronrod[8] = 0.030792300167387488891109020215229;
                wkronrod[9] = 0.034002130274329337836748795229551;
                wkronrod[10] = 0.037116271483415543560330625367620;
                wkronrod[11] = 0.040083825504032382074839284467076;
                wkronrod[12] = 0.042872845020170049476895792439495;
                wkronrod[13] = 0.045502913049921788909870584752660;
                wkronrod[14] = 0.047982537138836713906392255756915;
                wkronrod[15] = 0.050277679080715671963325259433440;
                wkronrod[16] = 0.052362885806407475864366712137873;
                wkronrod[17] = 0.054251129888545490144543370459876;
                wkronrod[18] = 0.055950811220412317308240686382747;
                wkronrod[19] = 0.057437116361567832853582693939506;
                wkronrod[20] = 0.058689680022394207961974175856788;
                wkronrod[21] = 0.059720340324174059979099291932562;
                wkronrod[22] = 0.060539455376045862945360267517565;
                wkronrod[23] = 0.061128509717053048305859030416293;
                wkronrod[24] = 0.061471189871425316661544131965264;
                wkronrod[25] = 0.061580818067832935078759824240055;
            }
            if( n==61 )
            {
                ng = 15;
                wgauss[0] = 0.007968192496166605615465883474674;
                wgauss[1] = 0.018466468311090959142302131912047;
                wgauss[2] = 0.028784707883323369349719179611292;
                wgauss[3] = 0.038799192569627049596801936446348;
                wgauss[4] = 0.048402672830594052902938140422808;
                wgauss[5] = 0.057493156217619066481721689402056;
                wgauss[6] = 0.065974229882180495128128515115962;
                wgauss[7] = 0.073755974737705206268243850022191;
                wgauss[8] = 0.080755895229420215354694938460530;
                wgauss[9] = 0.086899787201082979802387530715126;
                wgauss[10] = 0.092122522237786128717632707087619;
                wgauss[11] = 0.096368737174644259639468626351810;
                wgauss[12] = 0.099593420586795267062780282103569;
                wgauss[13] = 0.101762389748405504596428952168554;
                wgauss[14] = 0.102852652893558840341285636705415;
                x[0] = 0.999484410050490637571325895705811;
                x[1] = 0.996893484074649540271630050918695;
                x[2] = 0.991630996870404594858628366109486;
                x[3] = 0.983668123279747209970032581605663;
                x[4] = 0.973116322501126268374693868423707;
                x[5] = 0.960021864968307512216871025581798;
                x[6] = 0.944374444748559979415831324037439;
                x[7] = 0.926200047429274325879324277080474;
                x[8] = 0.905573307699907798546522558925958;
                x[9] = 0.882560535792052681543116462530226;
                x[10] = 0.857205233546061098958658510658944;
                x[11] = 0.829565762382768397442898119732502;
                x[12] = 0.799727835821839083013668942322683;
                x[13] = 0.767777432104826194917977340974503;
                x[14] = 0.733790062453226804726171131369528;
                x[15] = 0.697850494793315796932292388026640;
                x[16] = 0.660061064126626961370053668149271;
                x[17] = 0.620526182989242861140477556431189;
                x[18] = 0.579345235826361691756024932172540;
                x[19] = 0.536624148142019899264169793311073;
                x[20] = 0.492480467861778574993693061207709;
                x[21] = 0.447033769538089176780609900322854;
                x[22] = 0.400401254830394392535476211542661;
                x[23] = 0.352704725530878113471037207089374;
                x[24] = 0.304073202273625077372677107199257;
                x[25] = 0.254636926167889846439805129817805;
                x[26] = 0.204525116682309891438957671002025;
                x[27] = 0.153869913608583546963794672743256;
                x[28] = 0.102806937966737030147096751318001;
                x[29] = 0.051471842555317695833025213166723;
                x[30] = 0.000000000000000000000000000000000;
                wkronrod[0] = 0.001389013698677007624551591226760;
                wkronrod[1] = 0.003890461127099884051267201844516;
                wkronrod[2] = 0.006630703915931292173319826369750;
                wkronrod[3] = 0.009273279659517763428441146892024;
                wkronrod[4] = 0.011823015253496341742232898853251;
                wkronrod[5] = 0.014369729507045804812451432443580;
                wkronrod[6] = 0.016920889189053272627572289420322;
                wkronrod[7] = 0.019414141193942381173408951050128;
                wkronrod[8] = 0.021828035821609192297167485738339;
                wkronrod[9] = 0.024191162078080601365686370725232;
                wkronrod[10] = 0.026509954882333101610601709335075;
                wkronrod[11] = 0.028754048765041292843978785354334;
                wkronrod[12] = 0.030907257562387762472884252943092;
                wkronrod[13] = 0.032981447057483726031814191016854;
                wkronrod[14] = 0.034979338028060024137499670731468;
                wkronrod[15] = 0.036882364651821229223911065617136;
                wkronrod[16] = 0.038678945624727592950348651532281;
                wkronrod[17] = 0.040374538951535959111995279752468;
                wkronrod[18] = 0.041969810215164246147147541285970;
                wkronrod[19] = 0.043452539701356069316831728117073;
                wkronrod[20] = 0.044814800133162663192355551616723;
                wkronrod[21] = 0.046059238271006988116271735559374;
                wkronrod[22] = 0.047185546569299153945261478181099;
                wkronrod[23] = 0.048185861757087129140779492298305;
                wkronrod[24] = 0.049055434555029778887528165367238;
                wkronrod[25] = 0.049795683427074206357811569379942;
                wkronrod[26] = 0.050405921402782346840893085653585;
                wkronrod[27] = 0.050881795898749606492297473049805;
                wkronrod[28] = 0.051221547849258772170656282604944;
                wkronrod[29] = 0.051426128537459025933862879215781;
                wkronrod[30] = 0.051494729429451567558340433647099;
            }
            
            //
            // copy nodes
            //
            for(i=n-1; i>=n/2; i--)
            {
                x[i] = -x[n-1-i];
            }
            
            //
            // copy Kronrod weights
            //
            for(i=n-1; i>=n/2; i--)
            {
                wkronrod[i] = wkronrod[n-1-i];
            }
            
            //
            // copy Gauss weights
            //
            for(i=ng-1; i>=0; i--)
            {
                wgauss[n-2-2*i] = wgauss[i];
                wgauss[1+2*i] = wgauss[i];
            }
            for(i=0; i<=n/2; i++)
            {
                wgauss[2*i] = 0;
            }
            
            //
            // reorder
            //
            tsort.tagsort(ref x, n, ref p1, ref p2);
            for(i=0; i<=n-1; i++)
            {
                tmp = wkronrod[i];
                wkronrod[i] = wkronrod[p2[i]];
                wkronrod[p2[i]] = tmp;
                tmp = wgauss[i];
                wgauss[i] = wgauss[p2[i]];
                wgauss[p2[i]] = tmp;
            }
        }
    }
}