BZOJ-4891 [Tjoi2017]龙舟(Pollard-Rho算法+exgcd)

题目描述

  给出一个长为 \(m\) 的序列 \(b=[b_1,b_2,\cdots,b_m]\)，和 \(n\) 个长为 \(m\) 的序列 \(a_1=[a_{11},a_{12},\cdots,a_{1m}],a_{2}=[a_{21},a_{22},\cdots,a_{2m}],\cdots,a_{n}=[a_{n1},a_{n2},\cdots,a_{nm}]\)。

  有 \(k\) 次询问，每次给出 \(x,M\)，求：

\[\Bigg(\frac{\displaystyle\prod_{i=1}^{m}b_i}{\displaystyle\prod_{i=1}^{m}a_{xi}}\Bigg)^{-1}\pmod{M} \]

  数据范围：\(n\leq 20,m\leq 10000,k\leq 50,M,a,b\leq 2\times 10^{18}\)。

分析

  如果上下分解质因子，约分后求值，时间复杂度为 \(O(kmM^{\frac{1}{4}})\)，显然不行。

  考虑直接把模数 \(M\) 分解，这样可以把分子分母提前约掉 \(M\) 的质因子，剩余的部分和 \(M\) 互质，扩展欧几里得算法求逆元即可；如果某一质因子有负数次幂，则它与 \(M\) 并不互质，输出 \(-1\)。

  时间复杂度 \(O(k(M^\frac{1}{4}+mcnt))\)。

代码

#include<bits/stdc++.h>
using namespace std;
#define int128 __int128
long long n,m,k;
long long b[100010],a[30][10010];
int cnt;
long long factor[100010],num[100010];
long long x,y,M;
long long ans,inv;
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return ;
    }
    exgcd(b,a%b,y,x);
    y=y-x*(a/b);
}
long long quick_mul(long long a,long long b,long long mod)
{
    long long ans=0;
    while(b)
    {
        if(b&1)
            ans=(ans+a)%mod;
        a=(a+a)%mod;
        b>>=1;
    }
    return ans;
}
long long quick_pow(long long a,long long b,long long mod)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
            ans=quick_mul(ans,a,mod)%mod;
        a=quick_mul(a,a,mod)%mod;
        b>>=1;
    }
    return ans;
}
long long max_factor;
bool Miller_Rabin(long long p)
{
    if(p<2)
        return 0;
    if(p==2)
        return 1;
    if(p==3)
        return 1;
    long long d=p-1,r=0;
    while(d%2==0)
    {
        r++;
        d>>=1;
    }
    for(long long k=0;k<10;k++)
    {
        long long a=rand()%(p-2)+2;
        long long x=quick_pow(a,d,p);
        if(x==1||x==p-1)
            continue;
        for(int i=0;i<r-1;i++)
        {
            x=(int128)x*x%p;
            if(x==p-1)
                break;
        }
        if(x!=p-1)
            return 0;
    }
    return 1;
}
long long f(long long x,long long c,long long n)
{
    return ((int128)x*x+c)%n;
}
long long Pollard_Rho(long long n)
{
    long long s=0,t=0;
    long long c=rand()%(n-1)+1;
    int step=0,goal=1;
    long long val=1;
    for(goal=1; ;goal<<=1,s=t,val=1)
    {
        for(step=1;step<=goal;step++)
        {
            t=f(t,c,n);
            val=(int128)val*abs(t-s)%n;
            if(step%127==0)
            {
                long long d=__gcd(val,n);
                if(d>1)
                    return d;
            }
        }
        long long d=__gcd(val,n);
        if(d>1)
            return d;
    }
}
void fac(long long x)
{
    if(x<=max_factor||x<2)
        return ;
    if(Miller_Rabin(x))
    {
        factor[++cnt]=x;
        //max_factor=max(max_factor,x);
        return ;
    }
    long long p=x;
    while(p>=x)
        p=Pollard_Rho(x);
    while(x%p==0)
        x=x/p;
    fac(x);fac(p);
}
long long solve(long long x,long long M)
{
    ans=1,inv=1;
    cnt=0;
    memset(num,0,sizeof(num));
    fac(M);
    sort(factor+1,factor+1+cnt);
    cnt=unique(factor+1,factor+1+cnt)-(factor+1);
    for(int i=1;i<=m;i++)
    {
        long long temp=b[i];
        for(int j=1;j<=cnt;j++)
        {
            while(temp%factor[j]==0)
            {
                num[j]++;
                temp=temp/factor[j];
            }
        }
        ans=quick_mul(ans,temp,M)%M;
    }
    for(int i=1;i<=m;i++)
    {
        long long temp=a[x][i];
        for(int j=1;j<=cnt;j++)
        {
            while(temp%factor[j]==0)
            {
                num[j]--;
                temp=temp/factor[j];
            }
        }
        inv=quick_mul(inv,temp,M)%M;
    }
    for(int i=1;i<=cnt;i++)
        if(num[i]<0)
            return -1;
    for(int i=1;i<=cnt;i++)
        if(num[i]!=0)
            ans=quick_mul(ans,quick_pow(factor[i],num[i],M),M);
    exgcd(inv,M,x,y);//inv*x+M*y=1
    inv=(x%M+M)%M;
    return quick_mul(ans,inv,M);
}
int main()
{
    srand((unsigned)time(NULL));
    cin>>n>>m>>k;
    for(int i=1;i<=m;i++)
        scanf("%lld",&b[i]);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++)
            scanf("%lld",&a[i][j]);
    while(k--)
    {
        scanf("%lld %lld",&x,&M);
        printf("%lld\n",solve(x,M));
    }
    return 0;
}