NumberTheory

gistfile1.cpp

namespace NT {
    typedef long long LL;
    const int SIEVE_SIZE=(int)1e7;
    const int MOD=(int)1e9+7;

    int minp[SIEVE_SIZE+1],prime[SIEVE_SIZE+1]; //the least prime factor, the prime list    
    //O(n) get all prime numbers <= SIEVE_SIZE, return the number of primes
    int Sieve(int n=SIEVE_SIZE) {
        int c=0;
        for(int i=2;i<=n;++i) {
            if(!minp[i])    prime[c++]=i,minp[i]=i;
            for(int j=0;j<c && i*prime[j]<=n;++j) {
                minp[i*prime[j]]=prime[j];
                if(i%prime[j]==0)   break;
            }
        }
        return c;
    }

    vector<int> pl,lo,eu,rv,mu;     //[0,SIZE]
    //O(n) get pl[] (prime list), lo(minimum prime factor), eu[](euler function), rv[](the inverse of MOD)
    void FastSieve(int SIZE,int mod=MOD){
        lo=eu=rv=mu=vector<int>(SIZE+1);
        lo[1]=eu[1]=rv[1]=mu[1]=1;
        for(int x=2;x<=SIZE;x++){
            rv[x]=rv[mod%x]*LL(mod-mod/x)%mod;
            if(!lo[x]) pl.push_back(lo[x]=x),eu[x]=x-1,mu[x]=-1;
            for(size_t i=0;i<pl.size() && x*pl[i]<=SIZE;i++){
                int t=x*pl[i];
                lo[t]=pl[i];
                eu[t]=eu[x]*(pl[i]-(x%pl[i]!=0));
                if(x%pl[i]==0) {
                    mu[t]=0;    break;
                }else
                    mu[t]=-mu[x];
            }
        }
    }
    void FastSievePl(int SIZE) {
        lo=vector<int>(SIZE+1);
        lo[1]=1;
        for(int x=2;x<=SIZE;x++){            
            if(!lo[x]) pl.push_back(lo[x]=x);
            for(size_t i=0;i<pl.size() && x*pl[i]<=SIZE;i++){
                lo[x*pl[i]]=pl[i];                
                if(x%pl[i]==0)  break;                    
            }
        }    
    }
    //O(n) get rv[](inverse of mod)
    void FastSieveRv(int SIZE,int mod) {
        rv=vector<int>(SIZE+1); 
        rv[1]=1;
        for(int x=2;x<=SIZE;x++)
            rv[x]=rv[mod%x]*LL(mod-mod/x)%mod;
    }
    //O(n) get eu[](euler function)
    void FastSieveEu(int SIZE) {
        lo=eu=vector<int>(SIZE+1);
        lo[1]=eu[1]=1;
        for(int x=2;x<=SIZE;x++){            
            if(!lo[x]) pl.push_back(lo[x]=x),eu[x]=x-1;
            for(size_t i=0;i<pl.size() && x*pl[i]<=SIZE;i++){
                int t=x*pl[i];
                lo[t]=pl[i];
                eu[t]=eu[x]*(pl[i]-(x%pl[i]!=0));
                if(x%pl[i]==0)  break;                    
            }
        }  
    }
    //O(n) get mu[](mobius function)
    void FastSieveMu(int SIZE) {
        lo=mu=vector<int>(SIZE+1);
        lo[1]=mu[1]=1;
        for(int x=2;x<=SIZE;x++){        
            if(!lo[x]) pl.push_back(lo[x]=x),mu[x]=-1;
            for(size_t i=0;i<pl.size() && x*pl[i]<=SIZE;i++){
                int t=x*pl[i];
                lo[t]=pl[i];                
                if(x%pl[i]==0) {
                    mu[t]=0;    break;
                }else
                    mu[t]=-mu[x];
            }
        }
    }
    //O(sqrt(n)) get all prime factors of n
    vector<LL> Factorize(LL n){ //Make sure to call FastSievePl() before
        vector<LL> u;
        int i,t=sqrt(n+1);
        for(i=0;pl[i]<=t;i++) if(n%pl[i]==0){
            do{
                n/=pl[i];
                u.push_back(pl[i]);
            }while(n%pl[i]==0);
            t=sqrt(n+1);
        }
        if(n>1) u.push_back(n);
        return u;
    }
    //O(sqrt(n)) test if n is a prime number
    bool IsPrime(LL n) {    //Make sure to call FastSievePl() before
        if(n<(int)pl.size()) return n>=2 && lo[n]==n;
        int i,t=sqrt(n+1);
        for(i=0;pl[i]<=t;i++) if(n%pl[i]==0) return false;
        return true;
    }

    /*  Number Theory Functions */

    //O(log) get gcd(a,b)
    LL gcd(LL a,LL b) {
        return !b?a:gcd(b,a%b);
    }
    //O(log) find (x,y) st. ax+by=gcd(a,b)
    template<class T>
    T extgcd(T a,T b,T &x,T &y) {
        if(!b) return x=1,y=0,a;
        T u=extgcd(b,a%b,y,x);
        return y-=a/b*x,u;
    }
    //O(log) get x's multiplicative inverse mod m
    // if and only if gcd(x,m)==1, otherwise the inverse not exists,return 0    
    LL Inv(LL x,LL m) {
        LL r,v;
        return extgcd(x,m,r,v)==1?(r+m)%m:0;
    }
    //O(log(e)) get a^e (mod m)
    LL Pow(LL a,LL e,int m=MOD) {
        if(m==1)    return 0;
        LL s=1%m;
        for(a=a%m;e;a=a*a%m,e>>=1)
            if(e&1) s=s*a%m;
        return s;
    }
    //O(sqrt(n)) get eu(n)(euler function)
    int Eu(int n){
        int i,m=n;
        for(i=2;i*i<=n;i++) if(n%i==0)
            for(m=m/i*(i-1);n%i==0;n/=i);
        if(n>1) m=m/n*(n-1);
        return m;
    }

    //Discrete logarithm        return a^x %p = b the minimum solution x, if no solution, return -1
    //O(sqrt(p) * log(sqrt(p)))     p can be a composite number    
    int log(int a, int r, int m){
        if(r>=m) return -1;
        int i,g,x,c=0,at=int(2+sqrt(m));
        for(i=0,x=1%m;i<50;i++,x=LL(x)*a%m) if(x==r) return i;
        for(g=x=1;__gcd(int(LL(x)*a%m),m)!=g;c++) g=__gcd(x=LL(x)*a%m,m);
        if(r%g) return -1;
        if(x==r) return c;
        static unordered_map<int,int> u;    u.clear();
        g=Eu(m/g),u[x]=0;
        g=Pow(a,g-at%g,m);
        for(i=1;i<at;i++){ // Baby Step
            u.insert(MP(x=LL(x)*a%m,i));
            if(x==r) return c+i;
        }
        for(i=1;i<at;i++){ // Giant Step
            unordered_map<int,int>::iterator t=u.find(r=LL(r)*g%m);
            if(t!=u.end()) return c+i*at+t->second;
        }   
        return -1;
    }

    //O(log)  Linear Diophantine equations
    //find a solution of the equation : ax+by=c, the all solution is (x-k*(b/gcd(a,b)), y+k*(a/gcd(a,b)))
    //if no solution ,return false
    bool DioQuation(LL a,LL b,LL c,LL &x,LL &y) {
        LL r=extgcd(a,b,x,y);
        if(c%r) return false;
        x*=c/r; y*=c/r;
        return true;
    }

    //O(n*gcd)  Linear congruence equations     x%m[i]=r[i]  return the minimum x, return -1 if no solution
    //Can handle the case when m[i] and m[j] are not coprime
    //if r[i]>=m[i], return -1
    template<class T>
    T Equations(int n,const vector<T> &m,const vector<T> &r){
        T ans=0,lcm=1,x,y;
        for(int i=0;i<n;i++){
            T g=extgcd(lcm,m[i],x,y);
            if((r[i]-ans)%g || r[i]>=m[i]) return -1;
            x*=(r[i]-ans)/g; y=m[i]/g;
            ans+=(x%y+y)%y*lcm;
            ans%=lcm*=y;
        }
        return ans;
    }
    //Linear congruence equations     a[i]*x%m[i]=r[i]  return the minimum x, return -1 if no solution
    template<class T>
    T Equations(int n,vector<T> a,vector<T> m,vector<T> r) {        
        for(int i=0;i<n;++i) {
            T g=gcd(a[i],m[i]);
            if(r[i]%g || r[i]>=m[i])     return -1;
            a[i]/=g,m[i]/=g,r[i]/=g,r[i]=(LL)r[i]*Inv(a[i],m[i])%m[i];                        
        }        
        return Equations(n,m,r);
    }

    //O(log)    solve linear module equation ax=r(mod m)    return all solutions in the range [0,m)
    template<class T>
    vector<T> LinearEquation(T a,T r,T m) {
        vector<T> sol;
        T d,e,x,y;
        d=extgcd(a,m,x,y);
        if(r%d) return 0;        
        e=(x*(r/d)%m+m)%m;
        for(int i=0;i<d;++i)
            sol.push_back((e+i*(m/d))%m);
        return d;
    }

    //O(log(n)) get the last non-zero digit of n!
    const int GF[]={6, 6, 2, 6, 4, 4, 4, 8, 4, 6};
    const int GT[]={6, 2, 4, 8};
    const int G[] ={1, 1, 2, 6, 4};
    int LastDigit(int n) {
        if(n<5) return G[n];
        return GF[n%10]*GT[n/5*3%4]%10*LastDigit(n/5)%10;
    }
    int LastDigit(char *buf) {
        const int mod[20]={1,1,2,6,4,2,2,4,2,8,4,4,8,4,6,8,8,6,8,2},MAXN=10000;
        int len=strlen(buf),i,c,ret=1;    static int a[MAXN];
        if (len==1)         return mod[buf[0] - '0'];
        for (i=0;i<len;i++) a[i] = buf[len - 1 - i] - '0';        
        for (;len;len-=!a[len - 1]) {
            ret=ret*mod[a[1]%2*10+a[0]]%5;            
            for(c=0,i=len-1;i>=0;--i)
                c=c*10+a[i],a[i]=c/5,c%=5;
        }
        return ret+ret%2*5;        
    }
}