IntegerMath.h

Go to the documentation of this file.

/* This file is part of gbLAB.
 *
 * gbLAB is distributed without any warranty under the MIT License.
 */
 
 
#ifndef gbLAB_IntegerMath_h_
#define gbLAB_IntegerMath_h_
 
#include <numeric>
#include <Eigen/Dense>
#include <vector>
#include <iostream>
#include <algorithm>
#include <deque>
 
 
 
namespace gbLAB
{
    template <typename IntScalarType>
    struct IntegerMath
    {
        inline static IntScalarType positive_modulo(IntScalarType i, IntScalarType n) {
            return (i % n + n) % n;
        }
 
        static IntScalarType sgn(const IntScalarType& a)
        {
            return a<0? -1 : 1;
        }
        
        static IntScalarType gcd(const IntScalarType &a, const IntScalarType &b)
        {
            const IntScalarType absA(abs(a));
            const IntScalarType absB(abs(b));
            return absB > 0 ? gcd(absB, absA % absB) : (absA > 0 ? absA : 1);
        }
        
        template<typename T>
        static IntScalarType gcd(const Eigen::MatrixBase<T>& a)
        {
            switch (a.size())
            {
                case 0:
                {
                    throw std::runtime_error("gcd: array size is zero\n");
                    return 0;
                    break;
                }
                    
                case 1:
                {
                    return a(0);
                    break;
                }
                    
                case 2:
                {
                    return gcd(a(0),a(1));
                    break;
                }
                    
                default:
                {
                    IntScalarType temp(a(0));
                    for(long k=1;k<a.size();++k)
                    {
                        if (temp==0 && a(k)==0 && k!=a.size()-1)
                            temp= 0;
                        else
                            temp=gcd(temp,a(k));
                    }
                    return temp;
                    break;
                }
            }
        }
        
        static IntScalarType lcm(const IntScalarType &a, const IntScalarType &b)
        {
            return a * b / gcd(a, b);
        }
 
        template<typename T>
        static IntScalarType lcm(const Eigen::MatrixBase<T>& a)
        {
            switch (a.size())
            {
                case 0:
                {
                    throw std::runtime_error("lcm: array size is zero\n");
                    return 0;
                    break;
                }
 
                case 1:
                {
                    return a(0);
                    break;
                }
 
                case 2:
                {
                    return lcm(a(0),a(1));
                    break;
                }
 
                default:
                {
                    IntScalarType temp(a(0));
                    for(long k=1;k<a.size();++k)
                    {
                        if (temp==0 && a(k)==0 && k!=a.size()-1)
                            temp= 0;
                        else
                            temp=lcm(temp,a(k));
                    }
                    return temp;
                    break;
                }
            }
 
        }
 
        static Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic>
        integerGramSchmidt(const Eigen::Vector<IntScalarType,Eigen::Dynamic>& a)
        {
            int dim= a.size();
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> matrixA;
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,1> tempx;
            matrixA.resize(1,1);
            matrixA << a(0);
 
            for(int i=0;i<dim-1;i++) {
                matrixA.conservativeResize(i + 1, i + 2);
                // zero the new column
                matrixA.col(i + 1) = Eigen::Vector<IntScalarType, Eigen::Dynamic>::Zero(i + 1);
 
                matrixA(0, i + 1) = a(i + 1);
                tempx.conservativeResize(i + 2, 1);
 
                // fill the new row with tempx
                if (i != 0)
                {
                    tempx(i + 1) = 0;
                    matrixA.row(i) = tempx;
                }
 
                Eigen::Vector<IntScalarType,Eigen::Dynamic> det=
                        Eigen::Vector<IntScalarType,Eigen::Dynamic>::Zero(i+2);
                for(int j=0;j<i+2; j++)
                {
                    Eigen::MatrixXd minor= Eigen::MatrixXd::Zero(i+1,i+1);
                    std::vector<IntScalarType> indicesToKeep(i+2);
                    std::iota(std::begin(indicesToKeep),std::end(indicesToKeep),0);
                    indicesToKeep.erase(std::remove(indicesToKeep.begin(),indicesToKeep.end(),j),indicesToKeep.end());
 
                    Eigen::Vector<IntScalarType,Eigen::Dynamic> indicesToKeepVector;
                    indicesToKeepVector= Eigen::Map<Eigen::Vector<IntScalarType,Eigen::Dynamic>>(indicesToKeep.data(),i+1);
                    minor= matrixA(Eigen::indexing::all, indicesToKeepVector).template cast<double>();
                    det(j)= round(minor.determinant());
                }
 
                IntScalarType e= gcd(det);
 
 
                if (det == Eigen::Vector<IntScalarType,Eigen::Dynamic>::Zero(i+2))
                {
                    tempx= Eigen::Vector<IntScalarType,Eigen::Dynamic>::Zero(i+2);
                    tempx(i)= 1;
                    continue;
                }
 
                for (int j = 0; j < i + 2; j++)
                    tempx(j) = pow(-1, j + 1) * det(j) / e;
            }
            matrixA.conservativeResize(dim,dim);
            matrixA.row(dim-1)= tempx;
 
            return matrixA.block(1,0,dim-1,dim);
        }
 
        // Finds such x and y, that a  x + b  y = gcd(a, b)
        //https://github.com/ADJA/algos/blob/master/NumberTheory/DiophantineEquation.cpp
        static IntScalarType extended_gcd(IntScalarType a, IntScalarType b, IntScalarType &x, IntScalarType &y)
        {
            if (b == 0)
            {
                x = 1;
                y = 0;
                return a;
            }
            IntScalarType x1, y1;
            IntScalarType g = extended_gcd(b, a % b, x1, y1);
            x = y1;
            y = x1 - (a / b) * y1;
            return g;
        }
 
        // Solves equation a  x + b  y = c, writes answer to x and y
        static void solveDiophantine2vars(IntScalarType a, IntScalarType b, IntScalarType c, IntScalarType &x, IntScalarType &y)
        {
 
            IntScalarType g = extended_gcd(a, b, x, y);
 
            if (c % g != 0)
            {
                std::cout << a << "  " << b << "  " << c << "   " << g << std::endl;
                puts("Impossible");
                exit(0);
            }
 
            c /= g;
 
            x = x * c ;
            y = y * c ;
        }
 
        // Find u such that a_1 u_1 + a_2 u_2 + .... + a_n u_n = 1
        // Method 0:
        // "A fast algorithm to find reduced hyperplane unit cells and solve N-dimensional Bézout's identities."
        //  Acta Crystallographica Section A: Foundations and Advances 77.5 (2021).
        // template<typename ArrayType>
        // static ArrayType solveBezout(const ArrayType& a)
        template<typename T>
        static Eigen::Vector<IntScalarType,Eigen::Dynamic> solveBezout(const Eigen::MatrixBase<T>& a)
        {
            int n= a.size();
            if (n<2) throw std::runtime_error("the size of arrays should be at least two");
            if (abs(IntegerMath<IntScalarType>::gcd(a)) != 1 || a.isZero())
                throw std::runtime_error("No solution since the gcd is not 1 or -1.");
 
            Eigen::Vector<IntScalarType,Eigen::Dynamic> u(n);
 
            IntScalarType p,q;
            IntScalarType g12= extended_gcd(a(0),a(1),p,q); // g12 - gcd of a(0) and a(1)
 
            // form the n-1 sized vector na= {g,a2,...,a_{n-1}}
            Eigen::Vector<IntScalarType,Eigen::Dynamic> na(n-1);
            na(0)= g12; na(Eigen::seq(1,n-2))= a(Eigen::seq(2,n-1));
 
            Eigen::Vector<IntScalarType,Eigen::Dynamic> k(n-1);
            if (n==2)
            {
                IntScalarType g= extended_gcd(a(0),a(1),u(0),u(1));
                if (g<0) u= -u;
                return u;
            }
            else
            {
                k = solveBezout(na);
                // {p*k0,q*k0,k1,..,k_{n-2}}
                u(0) = p * k(0);
                u(1) = q * k(0);
                u(Eigen::seq(2, n - 1)) = k(Eigen::seq(1, n - 2));
                return u;
            }
        }
 
        template<typename T>
        static Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> ccum(const Eigen::MatrixBase<T>& qin)
        //static Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic>
        //ccum(const Eigen::Vector<IntScalarType,Eigen::Dynamic>& qin)
        {
            int n= qin.size();
            Eigen::Vector<IntScalarType,Eigen::Dynamic> q(n);
            q= qin;
 
            try
            {
                if (n==1)
                    throw std::runtime_error("Size of the input matrix for CCUM should be >1.");
                if (abs(gcd(qin)) != 1)
                    throw std::runtime_error("The absolute gcd of the input array is not 1.");
            }
            catch(std::runtime_error& e)
            {
                std::cerr << e.what();
            }
 
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> Q=
                    Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic>::Identity(n,n);
 
            // first ensure that q(0) is non-zero; otherwise swap
            bool swapFirst= false;
            int swapFirstElementWith;
            if (q(0)==0)
            {
                int elementIndex= -1;
                for (auto element : q)
                {
                    elementIndex++;
                    if (element !=0) {
                        // swap
                        IntScalarType temp= q(0);
                        q(0)= element;
                        q(elementIndex)= temp;
                        swapFirst = true;
                        swapFirstElementWith= elementIndex;
                        break;
                    }
                }
            }
            Q.col(0)= q;
 
            int elementIndex= -1;
            IntScalarType y;
            // ensure q(0) and q(1) are co-prime; otherwise swap
            for (const IntScalarType& element : q)
            {
                elementIndex++;
                if (elementIndex == 0) continue;
                if (!(q(0)==0 && element==0))
                    y= gcd(q(0),element);
                else
                    continue;
                if (abs(y)==1)
                {
                    // swap rows
                    int temp= Q(1,0);
                    Q(1,0)= Q(elementIndex,0);
                    Q(elementIndex,0)= temp;
 
                    IntScalarType u,v;
                    IntegerMath<IntScalarType>::solveDiophantine2vars(q(0),element,1,u,v);
                    //Q(0,elementIndex)=-v; Q(elementIndex,elementIndex)= u;
                    Q(0,1)=-v; Q(1,1)= u;
                    // unswap
                    Q.row(1).swap(Q.row(elementIndex));
 
                    if(swapFirst) Q.row(0).swap(Q.row(swapFirstElementWith));
                    return Q;
                }
            }
 
            // At this point, q(0) is not co-prime with any other elements of q
            if (!(q(0) ==0 && q(1) == 0))
                y= gcd(q(0),q(1));
            else
                y= 0;
            IntScalarType u,v;
            IntegerMath<IntScalarType>::solveDiophantine2vars(q(0),q(1),y,u,v);
 
            IntScalarType alpha,beta;
            IntegerMath<IntScalarType>::solveDiophantine2vars(u,v,1,alpha,beta);
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> M=
                    Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic>::Identity(n,n);
            M.block(0,0,2,2) << u,     v,
                               -beta, alpha;
 
            Eigen::Vector<IntScalarType,Eigen::Dynamic> t;
            t= M*q;
            elementIndex= -1;
            int swappedElementIndex;
            bool swapped= false;
            for (IntScalarType& element : t)
            {
                elementIndex++;
                if (elementIndex<2) continue;
                // Find a element in {t(2),...,} that y:=gcd(q(0), q(1)) does not divide
                // We are guaranteed to find such an element since gcd(q)=1
                if (element%y != 0)
                {
                    // swap the found element with t(1)
                    int temp= t(elementIndex);
                    t(elementIndex)= t(1);
                    t(1)= temp;
                    swappedElementIndex= elementIndex;
                    swapped= true;
                    break;
                }
                assert(elementIndex!=t.size());
            }
 
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> tempMatrix(n,n) ;
            tempMatrix= ccum(t);
 
            // inv(P)*ccum(t)
            // swap rows 0 and swappedElementIndex
            if (swapped) tempMatrix.row(1).swap(tempMatrix.row(swappedElementIndex));
 
            // inv(M)*inv(P)*ccum(t)
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> invM=
                    Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic>::Identity(n,n);
            invM.block(0,0,2,2) << alpha, -v,
                                   beta,   u;
            Eigen::Matrix<IntScalarType,Eigen::Dynamic,Eigen::Dynamic> output(n,n);
            output= invM*tempMatrix;
            if(swapFirst) output.row(0).swap(output.row(swapFirstElementWith));
            return output;
        }
 
 
        static std::deque<std::vector<IntScalarType>> comb(const int& n, const int& k)
        {
            std::deque<std::vector<IntScalarType>> output;
            std::string bitmask(k, 1); // K leading 1's
            bitmask.resize(n, 0); // N-K trailing 0's
 
            // print integers and permute bitmask
            do {
                std::vector<IntScalarType> kCombination(k);
                int arrIndex= 0;
 
                for (int i = 0; i < n; ++i) // [0..N-1] integers
                {
                    if (bitmask[i])
                    {
                        kCombination[arrIndex]= (IntScalarType) i;
                        arrIndex++;
                    }
                }
                assert(arrIndex==k);
                output.push_back(kCombination);
            } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
            return output;
        }
 
    };
 
    template <typename IntScalarType, int dim>
    class MatrixDimIExt : public Eigen::Matrix<IntScalarType,dim,dim>
    {
    public:
        static Eigen::Matrix<IntScalarType,dim,dim> adjoint(const Eigen::Matrix<IntScalarType,dim,dim>& input)
        {
            Eigen::Matrix<IntScalarType,dim,dim> output;
            for (int i=0; i<dim; i++)
            {
                // remove i-th row
                Eigen::Matrix<IntScalarType,dim-1,dim> temp1 ;
                temp1 << input(Eigen::seq(0,i-1),Eigen::all),
                         input(Eigen::seq(i+1,dim-1),Eigen::all);
                for (int j=0; j<dim; j++)
                {
                    // remove jth column
                    Eigen::Matrix<IntScalarType,dim-1,dim-1> temp2 ;
                    temp2 << temp1(Eigen::all,Eigen::seq(0,j-1)), temp1(Eigen::all,Eigen::seq(j+1,dim-1));
 
                    output(i,j)=(IntScalarType) std::pow(-1,i+j+2)*temp2.template cast<double>().determinant();
                }
            }
 
            return output.transpose();
        }
    };
 
    template<typename IntScalarType>
    class MatrixDimIExt<IntScalarType,1> : public Eigen::Matrix<IntScalarType,1,1>
    {
    public:
        static Eigen::Matrix<IntScalarType,1,1> adjoint(const Eigen::Matrix<IntScalarType,1,1>& input)
        {
            Eigen::Matrix<IntScalarType,1,1> output;
            output << 1;
            return output;
        }
 
    };
}
#endif