<commit_before>#ifndef EIGEN_LAURIEGAUTSCHI_H
#define EIGEN_LAURIEGAUTSCHI_H

namespace Eigen {

    /**
     * \ingroup NumericalIntegration_Module
     *
     * \class LaurieGautschi
     *
     * \tparam _Scalar floating point type
     *
     * \brief This class computes Kronrod abscissae & weights for arbitrary precision.
     *
     * Based on work of <NAME> and <NAME>.
     * <NAME> (1997). Calculation of Gauss-Kronrod Quadrature Rules.
     * Mathematics of Computation, 66(219), 1133-1145
     * Created by <NAME>, November 7, 2011.
     * Ported to C++/Eigen by <NAME>, <NAME>, and <NAME> September 2014
     *
     * \TODO Ensure only appropriates types are used for Scalar, e.g. prohibit integers.
     */
    template<typename _Scalar>
    class LaurieGautschi
    {
    public:
        typedef _Scalar Scalar;
        typedef typename Eigen::Matrix<Scalar,Eigen::Dynamic,1>  VectorType;
        typedef typename Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> MatrixType;
        typedef typename VectorType::Index IndexType;
        typedef typename Eigen::SelfAdjointEigenSolver<MatrixType> SelfAdjointEigenSolverType;

        /**
         * \brief Recurrence coefficients for monic Jacobi polynomials.
         *
         * This method generates the first N recurrence
         * coefficients for monic Jacobi polynomials with parameters
         * alpha and beta. These are orthogonal on [-1,1] relative to the
         * weight function w(t)=(1-t)^a(1+t)^b. The N alpha-coefficients
         * are stored in \a alphaOut, the n beta-coefficients in \a betaOut.
         * http://en.wikipedia.org/wiki/Jacobi_polynomials
         *
         * Created by <NAME>, 6-22-1998; edited by <NAME>, 4-4-2002.
         * Ported to C++/Eigen by <NAME>, <NAME> and  <NAME>, September 2014
         *
         * \param[in] N Number of recurrence coefficients
         * \param[in] alpha Alpha parameter of the Jacobi-polynomials
         * \param[in] beta Beta parameter of the Jacobi-polynomials
         * \param[in,out] alphaOut N alpha-coefficients
         * \param[in,out] betaOut N beta-coefficients
         */
        static void r_jacobi(const IndexType N,const Scalar alpha,const Scalar beta,
            VectorType& alphaOut, VectorType& betaOut)
        {
            //TODO : make use the eigen assert facilities
            assert(alpha > Scalar(-1));
            assert(beta > Scalar(-1));
            assert(alphaOut.rows() == betaOut.rows());
            assert(alphaOut.rows() > 0);
            assert(N <= alphaOut.rows());

            using std::pow;
            using std::tgamma;
            alphaOut(0) = (beta-alpha)/(alpha+beta+Scalar(2.));
            betaOut(0) = pow(Scalar(2.),(alpha+beta+Scalar(1.))) * tgamma(alpha+Scalar(1.)) * tgamma(beta+Scalar(1.)) / tgamma(alpha+beta+Scalar(2.));

            for(IndexType n=1;n<N;++n)
            {
                Scalar nAlphaBeta = Scalar(2.)*n+alpha+beta;
                alphaOut(n) = (beta*beta - alpha*alpha) / (nAlphaBeta*(nAlphaBeta+Scalar(2.)));
                betaOut(n) =  Scalar(4.)*(n+alpha)*(n+beta)*n*(n+alpha+beta) / (nAlphaBeta*nAlphaBeta*(nAlphaBeta+Scalar(1.))*(nAlphaBeta-Scalar(1.)));
            }
        }

        /**
         * \brief Recurrence coefficients for monic Jacobi polynomials on [0,1].
         *
         * This method generates the first N recurrence
         * coefficients for monic Jacobi polynomials on [0,1] with
         * parameters alpha and beta. These are orthogonal on [0,1] relative
         * to the weight function w(t)=(1-t)^alpha t^beta. The N alpha-
         * coefficients are stored in \a alphaOut, the N beta-
         * coefficients in \a betaOut.
         * http://en.wikipedia.org/wiki/Jacobi_polynomials
         *
         * Created by <NAME>, 6-22-1998; edited by <NAME>, 4-4-2002.
         * Ported to C++/Eigen by <NAME>, <NAME> and  <NAME>, September 2014
         *
         * \param[in] N Number of recurrence coefficients
         * \param[in] alpha Alpha parameter of the Jacobi-polynomials
         * \param[in] beta Beta parameter of the Jacobi-polynomials
         * \param[in,out] alphaOut N alpha-coefficients
         * \param[in,out] betaOut N beta-coefficients
         */
        static void r_jacobi_01(const IndexType N,const Scalar alpha,const Scalar beta,
            VectorType& alphaOut, VectorType& betaOut)
        {
            //TODO : make use the eigen assert facilities
            assert(alpha > Scalar(-1));
            assert(beta > Scalar(-1));
            assert(alphaOut.rows() == betaOut.rows());
            assert(alphaOut.rows() > 0);
            assert(N <= alphaOut.rows());

            r_jacobi(N, alpha, beta, alphaOut, betaOut);

            for(IndexType n=0; n<N; ++n)
            {
                alphaOut(n) = (Scalar(1.)+alphaOut(n))/Scalar(2.);
            }

            using std::pow;
            betaOut(0) = betaOut(0)/pow(Scalar(2),alpha+beta+Scalar(1.));

            for(IndexType n=1; n<N; ++n)
            {
                betaOut(n) = betaOut(n)/Scalar(4.);
            }

        }

        /**
         * \brief Jacobi-Kronrod matrix.
         *
         * This method produces the alpha- and beta-elements in
         * the Jacobi-Kronrod matrix of order 2N+1 for the weight
         * function (or measure) w. The input data for the weight
         * function w are the recurrence coefficients of the associated
         * orthogonal polynomials, which are stored in \a alphaIn and \a betaIn .
         * At least ceil(3*N/2)+1 coefficients should be provided.
         * The 2N+1 alpha- and beta-elements are returned in \a alpha and \a beta
         * respectively.
         *
         * Created by <NAME>, 6-22.1998
         * Edited by <NAME>, November 7, 2011
         * Ported to C++/Eigen by <NAME>, <NAME> and  <NAME>, September 2014
         *
         * \param[in] N Number of nodes
         * \param[in,out] alphaIn The recurrence coefficients of the associated orthogonal polynomials
         * \param[in,out] betaIn The recurrence coefficients of the associated orthogonal polynomials
         * \param[in,out] alpha Alpha-elements in the Jacobi-Kronrod matrix of order 2N+1
         * \param[in,out] beta Beta-elements in the Jacobi-Kronrod matrix of order 2N+1
         *
         */
        static void r_kronrod(const IndexType N,VectorType const& alphaIn, VectorType const& betaIn,
            VectorType& alpha, VectorType& beta)
        {
            //TODO : make use the eigen assert facilities
            assert(alphaIn.rows() == betaIn.rows());
            assert(alphaIn.rows() > 0);
            assert(alphaIn.rows() >= ceil(3*N/2)+1 );
            assert(alpha.rows() == 2*N+1);
            assert(beta.rows() == 2*N+1);

            alpha = VectorType::Zero(2*N+1);
            beta = VectorType::Zero(2*N+1);

            for(IndexType k=0; k<=floor(3*N/2)+1; ++k)
            {
                alpha(k) = alphaIn(k);
            }

            for(IndexType k=0; k<=ceil(3*N/2)+1; ++k)
            {
                beta(k) = betaIn(k);
            }

            VectorType sigma = VectorType::Zero(floor(N/2)+2);
            VectorType tempVector = VectorType::Zero(floor(N/2)+2);

            tempVector(1) = beta(N+1);

            for(IndexType m=0; m<N-2+1; ++m)
            {
                Scalar u = 0;
                for(IndexType k=floor((m+1)/2); k>=0;--k)
                {
                    IndexType l = m-k;
                    u = u + ( alpha(k+N+1)-alpha(l) )*tempVector(k+1) + beta(k+N+1)*sigma(k) - beta(l)*sigma(k+1);
                    sigma(k+1) = u;
                }

                VectorType swap = sigma;
                sigma = tempVector;
                tempVector = swap;
            }

            for(IndexType j=floor(N/2); j>=0; --j)
            {
                sigma(j+1) = sigma(j);
            }

            for(IndexType m = N-1; m<2*N-3+1; ++m)
            {
                IndexType k = m+1-N;
                IndexType j = 0;
                Scalar u = 0;
                for(k=m+1-N; k<floor((m-1)/2)+1; ++k)
                {
                    IndexType l = m-k;
                    j = N-1-l;
                    u = u - (alpha(k+N+1)-alpha(l))*tempVector(j+1) - beta(k+N+1)*sigma(j+1) + beta(l)*sigma(j+2);
                    sigma(j+1) = u;
                }

                k=floor((m+1)/2);

                if(m % 2 == 0)
                {
                    alpha(k+N+1) = alpha(k) + (sigma(j+1)-beta(k+N+1)*sigma(j+2)) / tempVector(j+2);
                }
                else
                {
                    beta(k+N+1) = sigma(j+1) / sigma(j+2);
                }

                VectorType swap = sigma;
                sigma = tempVector;
                tempVector = swap;
            }

            alpha(2*N) = alpha(N-1)-beta(2*N)*sigma(1)/tempVector(1);
        }

        /**
         * \brief Gauss-Kronrod quadrature formula.<commit_msg>Add std scope resolution to floor() and ceil() function calls to eliminate a compiler error thrown by a new install of Ubuntu 15.04 with all packages updated and upgraded.<commit_after>#ifndef EIGEN_LAURIEGAUTSCHI_H
#define EIGEN_LAURIEGAUTSCHI_H

namespace Eigen {

    /**
     * \ingroup NumericalIntegration_Module
     *
     * \class LaurieGautschi
     *
     * \tparam _Scalar floating point type
     *
     * \brief This class computes Kronrod abscissae & weights for arbitrary precision.
     *
     * Based on work of <NAME> and <NAME>.
     * <NAME> (1997). Calculation of Gauss-Kronrod Quadrature Rules.
     * Mathematics of Computation, 66(219), 1133-1145
     * Created by <NAME>, November 7, 2011.
     * Ported to C++/Eigen by <NAME>, <NAME>, and <NAME> September 2014
     *
     * \TODO Ensure only appropriates types are used for Scalar, e.g. prohibit integers.
     */
    template<typename _Scalar>
    class LaurieGautschi
    {
    public:
        typedef _Scalar Scalar;
        typedef typename Eigen::Matrix<Scalar,Eigen::Dynamic,1>  VectorType;
        typedef typename Eigen::Matrix<Scalar,Eigen::Dynamic, Eigen::Dynamic> MatrixType;
        typedef typename VectorType::Index IndexType;
        typedef typename Eigen::SelfAdjointEigenSolver<MatrixType> SelfAdjointEigenSolverType;

        /**
         * \brief Recurrence coefficients for monic Jacobi polynomials.
         *
         * This method generates the first N recurrence
         * coefficients for monic Jacobi polynomials with parameters
         * alpha and beta. These are orthogonal on [-1,1] relative to the
         * weight function w(t)=(1-t)^a(1+t)^b. The N alpha-coefficients
         * are stored in \a alphaOut, the n beta-coefficients in \a betaOut.
         * http://en.wikipedia.org/wiki/Jacobi_polynomials
         *
         * Created by <NAME>, 6-22-1998; edited by <NAME>, 4-4-2002.
         * Ported to C++/Eigen by <NAME>, <NAME> and  <NAME>, September 2014
         *
         * \param[in] N Number of recurrence coefficients
         * \param[in] alpha Alpha parameter of the Jacobi-polynomials
         * \param[in] beta Beta parameter of the Jacobi-polynomials
         * \param[in,out] alphaOut N alpha-coefficients
         * \param[in,out] betaOut N beta-coefficients
         */
        static void r_jacobi(const IndexType N,const Scalar alpha,const Scalar beta,
            VectorType& alphaOut, VectorType& betaOut)
        {
            //TODO : make use the eigen assert facilities
            assert(alpha > Scalar(-1));
            assert(beta > Scalar(-1));
            assert(alphaOut.rows() == betaOut.rows());
            assert(alphaOut.rows() > 0);
            assert(N <= alphaOut.rows());

            using std::pow;
            using std::tgamma;
            alphaOut(0) = (beta-alpha)/(alpha+beta+Scalar(2.));
            betaOut(0) = pow(Scalar(2.),(alpha+beta+Scalar(1.))) * tgamma(alpha+Scalar(1.)) * tgamma(beta+Scalar(1.)) / tgamma(alpha+beta+Scalar(2.));

            for(IndexType n=1;n<N;++n)
            {
                Scalar nAlphaBeta = Scalar(2.)*n+alpha+beta;
                alphaOut(n) = (beta*beta - alpha*alpha) / (nAlphaBeta*(nAlphaBeta+Scalar(2.)));
                betaOut(n) =  Scalar(4.)*(n+alpha)*(n+beta)*n*(n+alpha+beta) / (nAlphaBeta*nAlphaBeta*(nAlphaBeta+Scalar(1.))*(nAlphaBeta-Scalar(1.)));
            }
        }

        /**
         * \brief Recurrence coefficients for monic Jacobi polynomials on [0,1].
         *
         * This method generates the first N recurrence
         * coefficients for monic Jacobi polynomials on [0,1] with
         * parameters alpha and beta. These are orthogonal on [0,1] relative
         * to the weight function w(t)=(1-t)^alpha t^beta. The N alpha-
         * coefficients are stored in \a alphaOut, the N beta-
         * coefficients in \a betaOut.
         * http://en.wikipedia.org/wiki/Jacobi_polynomials
         *
         * Created by <NAME>, 6-22-1998; edited by <NAME>, 4-4-2002.
         * Ported to C++/Eigen by <NAME>, <NAME> and  <NAME>, September 2014
         *
         * \param[in] N Number of recurrence coefficients
         * \param[in] alpha Alpha parameter of the Jacobi-polynomials
         * \param[in] beta Beta parameter of the Jacobi-polynomials
         * \param[in,out] alphaOut N alpha-coefficients
         * \param[in,out] betaOut N beta-coefficients
         */
        static void r_jacobi_01(const IndexType N,const Scalar alpha,const Scalar beta,
            VectorType& alphaOut, VectorType& betaOut)
        {
            //TODO : make use the eigen assert facilities
            assert(alpha > Scalar(-1));
            assert(beta > Scalar(-1));
            assert(alphaOut.rows() == betaOut.rows());
            assert(alphaOut.rows() > 0);
            assert(N <= alphaOut.rows());

            r_jacobi(N, alpha, beta, alphaOut, betaOut);

            for(IndexType n=0; n<N; ++n)
            {
                alphaOut(n) = (Scalar(1.)+alphaOut(n))/Scalar(2.);
            }

            using std::pow;
            betaOut(0) = betaOut(0)/pow(Scalar(2),alpha+beta+Scalar(1.));

            for(IndexType n=1; n<N; ++n)
            {
                betaOut(n) = betaOut(n)/Scalar(4.);
            }

        }

        /**
         * \brief Jacobi-Kronrod matrix.
         *
         * This method produces the alpha- and beta-elements in
         * the Jacobi-Kronrod matrix of order 2N+1 for the weight
         * function (or measure) w. The input data for the weight
         * function w are the recurrence coefficients of the associated
         * orthogonal polynomials, which are stored in \a alphaIn and \a betaIn .
         * At least ceil(3*N/2)+1 coefficients should be provided.
         * The 2N+1 alpha- and beta-elements are returned in \a alpha and \a beta
         * respectively.
         *
         * Created by <NAME>, 6-22.1998
         * Edited by <NAME>, November 7, 2011
         * Ported to C++/Eigen by <NAME>, <NAME> and  <NAME>, September 2014
         *
         * \param[in] N Number of nodes
         * \param[in,out] alphaIn The recurrence coefficients of the associated orthogonal polynomials
         * \param[in,out] betaIn The recurrence coefficients of the associated orthogonal polynomials
         * \param[in,out] alpha Alpha-elements in the Jacobi-Kronrod matrix of order 2N+1
         * \param[in,out] beta Beta-elements in the Jacobi-Kronrod matrix of order 2N+1
         *
         */
        static void r_kronrod(const IndexType N,VectorType const& alphaIn, VectorType const& betaIn,
            VectorType& alpha, VectorType& beta)
        {
            //TODO : make use the eigen assert facilities
            assert(alphaIn.rows() == betaIn.rows());
            assert(alphaIn.rows() > 0);
            assert(alphaIn.rows() >= std::ceil(3*N/2)+1 );
            assert(alpha.rows() == 2*N+1);
            assert(beta.rows() == 2*N+1);

            alpha = VectorType::Zero(2*N+1);
            beta = VectorType::Zero(2*N+1);

            for(IndexType k=0; k<=std::floor(3*N/2)+1; ++k)
            {
                alpha(k) = alphaIn(k);
            }

            for(IndexType k=0; k<=std::ceil(3*N/2)+1; ++k)
            {
                beta(k) = betaIn(k);
            }

            VectorType sigma = VectorType::Zero(std::floor(N/2)+2);
            VectorType tempVector = VectorType::Zero(std::floor(N/2)+2);

            tempVector(1) = beta(N+1);

            for(IndexType m=0; m<N-2+1; ++m)
            {
                Scalar u = 0;
                for(IndexType k=std::floor((m+1)/2); k>=0;--k)
                {
                    IndexType l = m-k;
                    u = u + ( alpha(k+N+1)-alpha(l) )*tempVector(k+1) + beta(k+N+1)*sigma(k) - beta(l)*sigma(k+1);
                    sigma(k+1) = u;
                }

                VectorType swap = sigma;
                sigma = tempVector;
                tempVector = swap;
            }

            for(IndexType j=std::floor(N/2); j>=0; --j)
            {
                sigma(j+1) = sigma(j);
            }

            for(IndexType m = N-1; m<2*N-3+1; ++m)
            {
                IndexType k = m+1-N;
                IndexType j = 0;
                Scalar u = 0;
                for(k=m+1-N; k<std::floor((m-1)/2)+1; ++k)
                {
                    IndexType l = m-k;
                    j = N-1-l;
                    u = u - (alpha(k+N+1)-alpha(l))*tempVector(j+1) - beta(k+N+1)*sigma(j+1) + beta(l)*sigma(j+2);
                    sigma(j+1) = u;
                }

                k=std::floor((m+1)/2);

                if(m % 2 == 0)
                {
                    alpha(k+N+1) = alpha(k) + (sigma(j+1)-beta(k+N+1)*sigma(j+2)) / tempVector(j+2);
                }
                else
                {
                    beta(k+N+1) = sigma(j+1) / sigma(j+2);
                }

                VectorType swap = sigma;
                sigma = tempVector;
                tempVector = swap;
            }

            alpha(2*N) = alpha(N-1)-beta(2*N)*sigma(1)/tempVector(1);
        }

        /**
         * \brief Gauss-Kronrod quadrature formula.
{'fraction_non_alphanumeric': 0.12652094918240378, 'fraction_numerical': 0.022903705335021746, 'mean_word_length': 1.775248281130634, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 75, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '866394', 'n_tokens_mistral': 5889, 'n_tokens_neox': 5423, 'n_words': 1385}