gsSturmSolver.h

#include <gsPolynomial/gsMonomialPoly.h>


namespace gismo{

template<class T>
class gsSturmSolver
{
    

};

        template<typename T>
        gsMatrix<T> FindRootIntervals(gsMonomialPoly<T> & poly, T eps)
    {
                gsMatrix<T> root(1, 3);
                root.setZero();

                root = CalcIntervals(poly, eps, root);
                root.conservativeResize(root.rows() - 1, 3);
                root = SortRoots(root);

                return root;
        }
        

        template<typename T>
        gsMatrix<T> CalcIntervals(const gsMonomialPoly<T> & poly, T eps, gsMatrix<T> rootIntervals)
        {
                gsMonomialPoly<T> copyPoly = poly;
                gsVector<T> coeff = copyPoly.coefs();
                int degree = copyPoly.basis().size() - 1;

                // checks if the input polynomial has 0 as a root
                int divByX = 0;

                while (!copyPoly.isNull() && copyPoly.coefs()(0) == 0)
        {
                        rootIntervals(rootIntervals.rows() - 1, 0) = 0;
                        rootIntervals(rootIntervals.rows() - 1, 1) = 0;
                        rootIntervals(rootIntervals.rows() - 1, 2) += 1;

                        gsMonomialBasis<T> monomialBasis(1);
                        gsMatrix<T> coefs(2, 1);
                        coefs << 0, 1;
                        gsMonomialPoly<T> x(monomialBasis, coefs);
                        copyPoly = PolyDivision(copyPoly, x, 1);
                        divByX = 1;
                }

                if (divByX == 1)
        {
                        rootIntervals.conservativeResize(rootIntervals.rows() + 1, 3);
                }

                //if the polynomial is not constant, the sturm sequence is calculated
                if (! copyPoly.isConstant() )
        {
                        gsMatrix<T> sturmSeq(degree + 1, degree + 1);
                        sturmSeq = SturmSequence(copyPoly);

                        //checks if the polynomial has multiple roots
                        if (sturmSeq.bottomRows(1).isZero())
                        {
                                gsVector< T >vec = sturmSeq.row(sturmSeq.rows() - 2);
                                gsMonomialPoly<T >poly(copyPoly.basis(), vec);
                                rootIntervals = CalcIntervals(poly, eps, rootIntervals);
                                rootIntervals = CalcIntervals(PolyDivision(copyPoly, poly, 1), eps, rootIntervals);
                        }
                        // the polynomial doesn't have multiple roots
                        else{
                                RootIsolation(T(-CauchyBound(copyPoly)), CauchyBound(copyPoly),
                              eps, sturmSeq, rootIntervals);

                        }
                }
                return rootIntervals;
        }

        
        template <typename T>
        void RootIsolation(T leftBound, T rightBound, T eps, gsMatrix<T> & sturmSeq, gsMatrix<T> & rootIntervals){

                int roots = SignChanges(HornerEval(sturmSeq, leftBound)) - SignChanges(HornerEval(sturmSeq, rightBound));

                //no roots in this interval
                if (roots == 0) { return; }

                //one root and the length of the intervals is <= 0.01
                else if (roots == 1 && rightBound - leftBound <= 0.01)
        {
                        rootIntervals(rootIntervals.rows() - 1, 0) = leftBound;
                        rootIntervals(rootIntervals.rows() - 1, 1) = rightBound;
                        rootIntervals(rootIntervals.rows() - 1, 2) = 1;
                        rootIntervals.conservativeResize(rootIntervals.rows() + 1, 3);
                }
                //more roots and the length of the interval is smaller than the border eps
                else if (roots > 1 && rightBound - leftBound < eps)
        {
                        rootIntervals(rootIntervals.rows() - 1, 0) = leftBound;
                        rootIntervals(rootIntervals.rows() - 1, 1) = rightBound;
                        rootIntervals(rootIntervals.rows() - 1, 2) = roots;
                        rootIntervals.conservativeResize(rootIntervals.rows() + 1, 3);
                }

                //more roots and the length of the interval is greater than eps
                else
        {
                        T bound = (rightBound + leftBound) / 2;
                        gsVector<T> f = sturmSeq.row(0);
                        T root_check = HornerEvalPoly(f, bound);
            
                        // the middle of the interval is a root
                        if (root_check == 0)
                        {
                                roots = SignChanges(HornerEval(sturmSeq, T(bound - eps/2)))
                      - SignChanges(HornerEval(sturmSeq, T(bound + eps/2)));
                                if (roots > 1)
                                {
                                        rootIntervals(rootIntervals.rows() - 1, 0) = bound - eps / 2;
                                        rootIntervals(rootIntervals.rows() - 1, 1) = bound + eps / 2;
                                        rootIntervals(rootIntervals.rows() - 1, 2) = roots;
                                        rootIntervals.conservativeResize(rootIntervals.rows() + 1, 3);
                                }
                                else
                                {
                                        rootIntervals(rootIntervals.rows() - 1, 0) = bound;
                                        rootIntervals(rootIntervals.rows() - 1, 1) = bound;
                                        rootIntervals(rootIntervals.rows() - 1, 2) = 1;
                                        rootIntervals.conservativeResize(rootIntervals.rows() + 1, 3);
                                }
                                RootIsolation(leftBound, T(bound - eps/2), eps, sturmSeq, rootIntervals);
                                RootIsolation(T(bound + eps/2), rightBound, eps, sturmSeq, rootIntervals);
                        }
                        // the middle of the interval is not a root, so a bisection is done
                        else
                        {
                                RootIsolation(leftBound, bound, eps, sturmSeq, rootIntervals);
                                RootIsolation(bound, rightBound, eps, sturmSeq, rootIntervals);

                        }
                }
        }


        template<typename T>
        T CauchyBound(gsMonomialPoly<T> &poly)
        {
                gsVector<T> vec = poly.coefs();

                int lead = vec.size() - 1; // right-most non-zero
                for (int i = vec.size() - 1; i >= 0; i--)
        {
                        if (vec(i) != 0)
            {
                                lead = i;
                                break;
                        }
                }
        return 1 + (vec/vec[lead]).array().abs().maxCoeff();
        }

        template<typename T>
        gsMatrix<T> SturmSequence(gsMonomialPoly<T> & poly)
        {
                int deg =  poly.basis().size() - 1;
                gsMatrix<T> sturm(deg + 1, deg + 1);
                sturm.setZero();

                //first row = f
                sturm.row(0) = poly.coefs().transpose();

                //second row = f'
                sturm.row(1) = PolyDerivWithSameDegree(poly).coefs().transpose();

                if ( PolyDerivWithSameDegree(poly).isConstant() )
        {
                        sturm.conservativeResize(2, sturm.cols());
                        return sturm;
                }

                // i-th row = -rem(f_(i-2), f_(i-1)) 
                // until it is constant
                gsVector<T> dividend;
                gsVector<T> divisor;

                for (int i = 2; i <= deg; i++)
                {
                        dividend = sturm.row(i - 2);
                        divisor = sturm.row(i - 1);
                        gsMonomialPoly<T> poly1(poly.basis(), dividend);
                        gsMonomialPoly<T> poly2(poly.basis(), divisor);
                        gsMonomialPoly<T> poly3 = PolyDivision(poly1, poly2, 2);

                        if ( poly3.isNull() )
            {
                                sturm.conservativeResize(i + 1, sturm.cols());
                                return sturm;
                        }


                        for (int j = 0; j < poly3.coefs().size(); j++)
            {
                                sturm(i, j) -= poly3.coefs()(j);
                        }

                        if ( poly3.isConstant() )
            {
                                sturm.conservativeResize(i + 1, sturm.cols());
                                return sturm;
                        }

                }

                return sturm;
        }
        
        template<typename T>
        gsMonomialPoly<T> PolyDerivWithSameDegree(const gsMonomialPoly<T> & poly)
        {
                gsVector<T> deriv(poly.deg()+1);        
        //AM:
        //deriv.topRows(poly.deg()) = gsVector<int>::LinSpaced(poly.deg(),1,poly.deg()).array()
        //    * poly.coefs().bottomRows(poly.deg()).array();
        
                for (int i = 0; i < deriv.size() - 1; i++)
                {
                        deriv(i) = poly.coefs()(i + 1)*(i + 1);
                }
                deriv(deriv.size() - 1) = 0;
        return gsMonomialPoly<T>(poly.basis(), deriv);
        }


        template<typename T>
        gsMonomialPoly<T> PolyDivision(const gsMonomialPoly<T> & dividendPoly,
                                   const gsMonomialPoly<T> & divisorPoly, int num)
        {
                if (num != 1 && num != 2)
        {
                        std::cout << "ERROR: The number is not in the corresponding range. The dividend will be returned" << std::endl;
                        return dividendPoly;
                }
                gsMonomialPoly<T> dividendCopy = dividendPoly;
                gsMonomialPoly<T> divisorCopy = divisorPoly;

                // storing the coefs of dividend and divisor seperately in a gsVector/gsMatrix
                gsMatrix<T> dividend = dividendCopy.coefs().transpose();

                gsVector<T> divisor(dividend.size());
                divisor.setZero();
                for (int i = 0; i < divisorCopy.coefs().size(); i++)
        {
                        divisor(i) = divisorCopy.coefs()(i);
                }


                // getting the real degrees of dividend and divisor
                int deg_dividend = dividend.size() - 1;

                while (dividend(deg_dividend) == 0)
                {
                        deg_dividend--;
                }


                int deg_divisor = divisor.size() - 1;

                while (divisor(deg_divisor) == 0)
                {
                        deg_divisor--;
                }

                if (deg_divisor > deg_dividend)
        {
                        std::cout << "Error, the division is not possible. The degree of the divisor is greater than the degree of the dividend. The dividend will be returned";
                        return dividendPoly;
                }


                // Algorithm for polynomial long division
                gsVector<T> quot(dividend.size());
                gsVector<T> rem(dividend.size());

                quot.setZero();
                rem = dividend.transpose();
                int deg_rem = deg_dividend;
                
                int lead = 0;
                T lead_coeff = 0;

        gsVector<T> dummy;
                while (!(rem.maxCoeff() == 0 && rem.minCoeff() == 0) && deg_rem >= deg_divisor)
        {
                        lead       = deg_rem - deg_divisor;
                        lead_coeff = rem(deg_rem) / divisor(deg_divisor);
                        quot(lead) = quot(lead) + lead_coeff;
                        dummy      = ShiftRight(divisor, lead);
                        rem        = rem - dummy * lead_coeff; 
                        deg_rem    = deg_dividend;

                        while (rem(deg_rem) == 0 && deg_rem > 0)
                        {
                                deg_rem--;
                        }
                }

                // returning the right output
                if (num == 1){

                        gsMonomialBasis<T> monomial_basis(deg_dividend - deg_divisor);
                        gsMonomialPoly<T> polyQuot(monomial_basis, quot.head(deg_dividend - deg_divisor + 1));

                        return polyQuot;
                }
                else{

                        gsMonomialBasis<T> monomial_basis(deg_rem);
                        gsMonomialPoly<T> polyRem(monomial_basis, rem.head(deg_rem + 1));

                        return polyRem;
                }


        }

        template <typename T>
        gsVector<T> ShiftRight(const gsVector<T> & vec, int num)
        {        
                gsVector<T> shift_vec(vec.size());
        shift_vec.topRows(num).setZero();
        const index_t b = vec.size()-num;
        shift_vec.bottomRows(b) = vec.topRows(b);
                return shift_vec;
        }


        //returns the number of sign changes of a vector
        template <typename T>
        int SignChanges(const gsVector<T> & vec)
        {
        // AM:
        //int result = 0;
        //for (int i = 1; i < vec.size(); ++i)
        //    if (vec[i-1]*vec[i] .. 0) ++result;
            
                gsVector<T> vec_copy = vec;
                for (int i = 0; i < vec.size(); i++)
                {
                        vec_copy(i) = vec(i) / math::abs(vec(i));// 0, +1 or -1
                }

                int result = 0;
                for (int i = 0; i < vec_copy.size() - 1; i++)
                {
                        if (math::abs(vec_copy(i) + vec_copy(i + 1)) < 2) { result++; }

                }
                return result;
        }


        template <typename T>
        gsVector<T> HornerEval(gsMatrix<T> & mat, T num)
        {
                int col = mat.cols();
                gsVector<T> eval(mat.rows());
                eval = mat.col(col - 1);

                for (int i = 2; i <= mat.cols(); i++)
                {
                        eval = eval*num + mat.col(col - i);

                }
                return eval;
        }


        template <typename T>
        T HornerEvalPoly(const gsVector<T> & vec, T num)
        {
                T eval = vec(vec.size() - 1);

                for (int i = 2; i <= vec.size(); i++)
                {
                        eval = eval*num + vec(vec.size() - i);
                }
                return eval;
        }


        template<typename T>
        gsMatrix<T> SortRoots(gsMatrix<T> & unsort)
    {
                gsMatrix<T> sort(unsort.rows(), 3);
                sort.setZero();
                T min = -1;
                int position = -1;

                for (int i = 0; i < unsort.rows(); i++)
                {

                        min = unsort.col(0).minCoeff();
                        for (int j = 0; j < unsort.rows(); j++)
            {
                                if (min == unsort(j, 0)){
                                        position = j;
                                        break;
                                }
                        }

                        sort.row(i) = unsort.row(position);
                        unsort(position, 0) = INT_MAX;
                }
                
                for (int i = 0; i < sort.rows()-1; i++)
                {

                        if ((sort(i, 0) == sort(i + 1, 0) && sort(i, 1) <= sort(i + 1, 1)) ||
                                (sort(i, 0) <= sort(i + 1, 0) && sort(i, 1) == sort(i + 1, 1)))
            {
                                sort(i, 2) += 1;
                                for (int j = i+1; j < sort.rows()-1; j++){
                                        sort.row(j) = sort.row(j + 1);
                                }
                                
                                sort.conservativeResize(sort.rows()- 1, 3);
                                --i;
                        }       
                }
                return sort;
        }

}