source: trunk/sources/ALGLIB/gq.cs @ 2636

/*************************************************************************
Copyright (c) 2005-2007, 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 gq
    {
        /*************************************************************************
        Computation of nodes and weights for a Gauss quadrature formula

        The algorithm generates the N-point Gauss 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 zeroth moment Mu0

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

        INPUT PARAMETERS:
            Alpha   –   array[0..N-1], alpha coefficients
            Beta    –   array[0..N-1], beta coefficients
                        Zero-indexed element is not used and may be arbitrary.
                        Beta[I]>0.
            Mu0     –   zeroth moment of the weight function.
            N       –   number of nodes of the quadrature formula, N>=1

        OUTPUT PARAMETERS:
            Info    -   error code:
                        * -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.
            W       -   array[0..N-1] - array of quadrature weights.

          -- ALGLIB --
             Copyright 2005-2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgeneraterec(ref double[] alpha,
            ref double[] beta,
            double mu0,
            int n,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            int i = 0;
            double[] d = new double[0];
            double[] e = new double[0];
            double[,] z = new double[0,0];

            if( n<1 )
            {
                info = -1;
                return;
            }
            info = 1;
            
            //
            // Initialize
            //
            d = new double[n+1];
            e = new double[n+1];
            for(i=1; i<=n-1; i++)
            {
                d[i] = alpha[i-1];
                if( (double)(beta[i])<=(double)(0) )
                {
                    info = -2;
                    return;
                }
                e[i] = Math.Sqrt(beta[i]);
            }
            d[n] = alpha[n-1];
            
            //
            // EVD
            //
            if( !tdevd.tridiagonalevd(ref d, e, n, 3, ref z) )
            {
                info = -3;
                return;
            }
            
            //
            // Generate
            //
            x = new double[n];
            w = new double[n];
            for(i=1; i<=n; i++)
            {
                x[i-1] = d[i];
                w[i-1] = mu0*AP.Math.Sqr(z[1,i]);
            }
        }


        /*************************************************************************
        Computation of nodes and weights for a Gauss-Lobatto quadrature formula

        The algorithm generates the N-point Gauss-Lobatto quadrature formula  with
        weight function given by coefficients alpha and beta of a recurrence 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 zeroth moment Mu0

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

        INPUT PARAMETERS:
            Alpha   –   array[0..N-2], alpha coefficients
            Beta    –   array[0..N-2], beta coefficients.
                        Zero-indexed element is not used, may be arbitrary.
                        Beta[I]>0
            Mu0     –   zeroth moment of the weighting function.
            A       –   left boundary of the integration interval.
            B       –   right boundary of the integration interval.
            N       –   number of nodes of the quadrature formula, N>=3
                        (including the left and right boundary nodes).

        OUTPUT PARAMETERS:
            Info    -   error code:
                        * -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.
            W       -   array[0..N-1] - array of quadrature weights.

          -- ALGLIB --
             Copyright 2005-2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgenerategausslobattorec(double[] alpha,
            double[] beta,
            double mu0,
            double a,
            double b,
            int n,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            int i = 0;
            double[] d = new double[0];
            double[] e = new double[0];
            double[,] z = new double[0,0];
            double pim1a = 0;
            double pia = 0;
            double pim1b = 0;
            double pib = 0;
            double t = 0;
            double a11 = 0;
            double a12 = 0;
            double a21 = 0;
            double a22 = 0;
            double b1 = 0;
            double b2 = 0;
            double alph = 0;
            double bet = 0;

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

            if( n<=2 )
            {
                info = -1;
                return;
            }
            info = 1;
            
            //
            // Initialize, D[1:N+1], E[1:N]
            //
            n = n-2;
            d = new double[n+3];
            e = new double[n+2];
            for(i=1; i<=n+1; i++)
            {
                d[i] = alpha[i-1];
            }
            for(i=1; i<=n; i++)
            {
                if( (double)(beta[i])<=(double)(0) )
                {
                    info = -2;
                    return;
                }
                e[i] = Math.Sqrt(beta[i]);
            }
            
            //
            // Caclulate Pn(a), Pn+1(a), Pn(b), Pn+1(b)
            //
            beta[0] = 0;
            pim1a = 0;
            pia = 1;
            pim1b = 0;
            pib = 1;
            for(i=1; i<=n+1; i++)
            {
                
                //
                // Pi(a)
                //
                t = (a-alpha[i-1])*pia-beta[i-1]*pim1a;
                pim1a = pia;
                pia = t;
                
                //
                // Pi(b)
                //
                t = (b-alpha[i-1])*pib-beta[i-1]*pim1b;
                pim1b = pib;
                pib = t;
            }
            
            //
            // Calculate alpha'(n+1), beta'(n+1)
            //
            a11 = pia;
            a12 = pim1a;
            a21 = pib;
            a22 = pim1b;
            b1 = a*pia;
            b2 = b*pib;
            if( (double)(Math.Abs(a11))>(double)(Math.Abs(a21)) )
            {
                a22 = a22-a12*a21/a11;
                b2 = b2-b1*a21/a11;
                bet = b2/a22;
                alph = (b1-bet*a12)/a11;
            }
            else
            {
                a12 = a12-a22*a11/a21;
                b1 = b1-b2*a11/a21;
                bet = b1/a12;
                alph = (b2-bet*a22)/a21;
            }
            if( (double)(bet)<(double)(0) )
            {
                info = -3;
                return;
            }
            d[n+2] = alph;
            e[n+1] = Math.Sqrt(bet);
            
            //
            // EVD
            //
            if( !tdevd.tridiagonalevd(ref d, e, n+2, 3, ref z) )
            {
                info = -3;
                return;
            }
            
            //
            // Generate
            //
            x = new double[n+2];
            w = new double[n+2];
            for(i=1; i<=n+2; i++)
            {
                x[i-1] = d[i];
                w[i-1] = mu0*AP.Math.Sqr(z[1,i]);
            }
        }


        /*************************************************************************
        Computation of nodes and weights for a Gauss-Radau quadrature formula

        The algorithm generates the N-point Gauss-Radau  quadrature  formula  with
        weight function given by the coefficients alpha and  beta  of a recurrence
        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 zeroth moment Mu0

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

        INPUT PARAMETERS:
            Alpha   –   array[0..N-2], alpha coefficients.
            Beta    –   array[0..N-1], beta coefficients
                        Zero-indexed element is not used.
                        Beta[I]>0
            Mu0     –   zeroth moment of the weighting function.
            A       –   left boundary of the integration interval.
            N       –   number of nodes of the quadrature formula, N>=2
                        (including the left boundary node).

        OUTPUT PARAMETERS:
            Info    -   error code:
                        * -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.
            W       -   array[0..N-1] - array of quadrature weights.


          -- ALGLIB --
             Copyright 2005-2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgenerategaussradaurec(double[] alpha,
            double[] beta,
            double mu0,
            double a,
            int n,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            int i = 0;
            double[] d = new double[0];
            double[] e = new double[0];
            double[,] z = new double[0,0];
            double polim1 = 0;
            double poli = 0;
            double t = 0;

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

            if( n<2 )
            {
                info = -1;
                return;
            }
            info = 1;
            
            //
            // Initialize, D[1:N], E[1:N]
            //
            n = n-1;
            d = new double[n+1+1];
            e = new double[n+1];
            for(i=1; i<=n; i++)
            {
                d[i] = alpha[i-1];
                if( (double)(beta[i])<=(double)(0) )
                {
                    info = -2;
                    return;
                }
                e[i] = Math.Sqrt(beta[i]);
            }
            
            //
            // Caclulate Pn(a), Pn-1(a), and D[N+1]
            //
            beta[0] = 0;
            polim1 = 0;
            poli = 1;
            for(i=1; i<=n; i++)
            {
                t = (a-alpha[i-1])*poli-beta[i-1]*polim1;
                polim1 = poli;
                poli = t;
            }
            d[n+1] = a-beta[n]*polim1/poli;
            
            //
            // EVD
            //
            if( !tdevd.tridiagonalevd(ref d, e, n+1, 3, ref z) )
            {
                info = -3;
                return;
            }
            
            //
            // Generate
            //
            x = new double[n+1];
            w = new double[n+1];
            for(i=1; i<=n+1; i++)
            {
                x[i-1] = d[i];
                w[i-1] = mu0*AP.Math.Sqr(z[1,i]);
            }
        }


        /*************************************************************************
        Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
        nodes.

        INPUT PARAMETERS:
            N           -   number of nodes, >=1

        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,
                            in ascending order.
            W           -   array[0..N-1] - array of quadrature weights.


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgenerategausslegendre(int n,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            double[] alpha = new double[0];
            double[] beta = new double[0];
            int i = 0;

            if( n<1 )
            {
                info = -1;
                return;
            }
            alpha = new double[n];
            beta = new double[n];
            for(i=0; i<=n-1; i++)
            {
                alpha[i] = 0;
            }
            beta[0] = 2;
            for(i=1; i<=n-1; i++)
            {
                beta[i] = 1/(4-1/AP.Math.Sqr(i));
            }
            gqgeneraterec(ref alpha, ref beta, beta[0], n, ref info, ref x, ref w);
            
            //
            // 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  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 nodes, >=1
            Alpha       -   power-law coefficient, Alpha>-1
            Beta        -   power-law coefficient, Beta>-1

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -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/Alpha/Beta was passed
                            * +1    OK
            X           -   array[0..N-1] - array of quadrature nodes,
                            in ascending order.
            W           -   array[0..N-1] - array of quadrature weights.


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgenerategaussjacobi(int n,
            double alpha,
            double beta,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            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<1 | (double)(alpha)<=(double)(-1) | (double)(beta)<=(double)(-1) )
            {
                info = -1;
                return;
            }
            a = new double[n];
            b = new double[n];
            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( n>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<=n-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));
                }
            }
            gqgeneraterec(ref a, ref b, b[0], n, ref info, ref x, ref w);
            
            //
            // 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  nodes/weights  for  Gauss-Laguerre  quadrature  on  [0,+inf) with
        weight function W(x)=Power(x,Alpha)*Exp(-x)

        INPUT PARAMETERS:
            N           -   number of nodes, >=1
            Alpha       -   power-law coefficient, Alpha>-1

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -4    an  error  was   detected   when   calculating
                                    weights/nodes. Alpha is 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/Alpha was passed
                            * +1    OK
            X           -   array[0..N-1] - array of quadrature nodes,
                            in ascending order.
            W           -   array[0..N-1] - array of quadrature weights.


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgenerategausslaguerre(int n,
            double alpha,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            double[] a = new double[0];
            double[] b = new double[0];
            double t = 0;
            int i = 0;
            double s = 0;

            if( n<1 | (double)(alpha)<=(double)(-1) )
            {
                info = -1;
                return;
            }
            a = new double[n];
            b = new double[n];
            a[0] = alpha+1;
            t = gammafunc.lngamma(alpha+1, ref s);
            if( (double)(t)>=(double)(Math.Log(AP.Math.MaxRealNumber)) )
            {
                info = -4;
                return;
            }
            b[0] = Math.Exp(t);
            if( n>1 )
            {
                for(i=1; i<=n-1; i++)
                {
                    a[i] = 2*i+alpha+1;
                    b[i] = i*(i+alpha);
                }
            }
            gqgeneraterec(ref a, ref b, b[0], n, ref info, ref x, ref w);
            
            //
            // test basic properties to detect errors
            //
            if( info>0 )
            {
                if( (double)(x[0])<(double)(0) )
                {
                    info = -4;
                }
                for(i=0; i<=n-2; i++)
                {
                    if( (double)(x[i])>=(double)(x[i+1]) )
                    {
                        info = -4;
                    }
                }
            }
        }


        /*************************************************************************
        Returns  nodes/weights  for  Gauss-Hermite  quadrature on (-inf,+inf) with
        weight function W(x)=Exp(-x*x)

        INPUT PARAMETERS:
            N           -   number of nodes, >=1

        OUTPUT PARAMETERS:
            Info        -   error code:
                            * -4    an  error  was   detected   when   calculating
                                    weights/nodes.  May be, N is too large. Try to
                                    use multiple precision version.
                            * -3    internal eigenproblem solver hasn't converged
                            * -1    incorrect N/Alpha was passed
                            * +1    OK
            X           -   array[0..N-1] - array of quadrature nodes,
                            in ascending order.
            W           -   array[0..N-1] - array of quadrature weights.


          -- ALGLIB --
             Copyright 12.05.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void gqgenerategausshermite(int n,
            ref int info,
            ref double[] x,
            ref double[] w)
        {
            double[] a = new double[0];
            double[] b = new double[0];
            double t = 0;
            int i = 0;
            double s = 0;

            if( n<1 )
            {
                info = -1;
                return;
            }
            a = new double[n];
            b = new double[n];
            for(i=0; i<=n-1; i++)
            {
                a[i] = 0;
            }
            b[0] = Math.Sqrt(4*Math.Atan(1));
            if( n>1 )
            {
                for(i=1; i<=n-1; i++)
                {
                    b[i] = 0.5*i;
                }
            }
            gqgeneraterec(ref a, ref b, b[0], n, ref info, ref x, ref w);
            
            //
            // test basic properties to detect errors
            //
            if( info>0 )
            {
                for(i=0; i<=n-2; i++)
                {
                    if( (double)(x[i])>=(double)(x[i+1]) )
                    {
                        info = -4;
                    }
                }
            }
        }
    }
}