Miller-Rabin + Pollard-rho

Miller-Rabin + Pollard-Rho算法模板题2019 ACM-ICPC Asia Xuzhou E

题意：给一个大小为\(n\)的序列a,定义$ Z = \Pi_{i=1}^{n} a_i$。另给出 \(X ，Y\),若\(b_i = Z * X^{i}\)，求最大的\(i\)使得\(b_i | Y!\)。

题解：板子题，题解略。

#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=(a);i<(b);i++)
#define per(i,a,b) for(int i=(b)-1;i>=(a);i--)
#define pb push_back
#define db double
#define VI vector<int>
#define ll long long
#define ld long double
#define ull unsigned long long
#define PII pair<int,int>
#define fi first
#define se second
mt19937_64 mrand(random_device{}());
ll rnd(ll x){return mrand()%x;}
#define all(x) x.begin(),x.end()
ll gcd(ll a,ll b){return !b?a:gcd(b,a%b);}
const ll mod = 1e9+7;
ll qpow(ll a,ll b){ll t=1;assert(b>=0);for(;b;b>>=1){if(b&1)t=t*a%mod;a=a*a%mod;}return t;}

//时间复杂度7*log^3(n)
ll qmul(ll a,ll b,ll mod)//快速乘
{
    ll c=(ld)a/mod*b;
    ll res=(ull)a*b-(ull)c*mod;
    return (res+mod)%mod;
}
ll qpow(ll a,ll n,ll mod)//快速幂
{
    ll res=1;
    while(n)
    {
        if(n&1) res=qmul(res,a,mod);
        a=qmul(a,a,mod);
        n>>=1;
    }
    return res;
}
bool MRtest(ll n)//Miller Rabin Test
{
    if(n<3||n%2==0) return n==2;//特判
    ll u=n-1,t=0;
    while(u%2==0) u/=2,++t;
    //ll ud[] = {2,7,61};
    ll ud[]={2,325,9375,28178,450775,9780504,1795265022};
    for(ll a:ud)
    {
        ll v=qpow(a,u,n);
        if(v==1||v==n-1||v==0) continue;
        for(int j=1;j<=t;j++)
        {
            v=qmul(v,v,n);
            if(v==n-1&&j!=t){v=1;break;}//出现一个n-1，后面都是1，直接跳出
            if(v==1) return 0;//这里代表前面没有出现n-1这个解，二次检验失败
        }
        if(v!=1) return 0;//Fermat检验
    }
    return 1;
}

//找到n的一个因子，时间复杂度是O(n^(1/4))
//结合Miller-Rabin可用来求出最大素因子
//模板题：https://www.luogu.com.cn/problem/P4718
//进一步优化可以考虑倍增和固定上限结合
ll Pollard_Rho(ll N){
  if (N == 4)return 2;
  if (MRtest(N))return N;
  while(1){
    ll c = rnd(N-1) + 1;
    auto f = [=](ll x) { return ((__int128)x * x + c) % N; };
    ll t = 0, r = 0, p = 1, q;
    do{
      for (int i = 0; i < 128; ++i){ // 令固定距离C=128
        t = f(t), r = f(f(r));
        if (t == r || (q = (__int128)p * abs(t - r) % N) == 0) // 如果发现环，或者积即将为0，退出
          break;
        p = q;
      }
      ll d = gcd(p, N);
      if (d > 1)return d;
    }while (t != r);
  }
}

ll max_prime_factor(ll x){
  ll fac = Pollard_Rho(x);
  if (fac == x)return x;
  else return max(max_prime_factor(fac), max_prime_factor(x / fac));
}

const int N = 1e5+5;
ll a[N],x,y;
int n;
ll fac[30],idx;
int fac_cnt[30];

void chai(ll x){
  while(x!=1){
    ll temp = max_prime_factor(x);
    fac[++idx] = temp;
    while(x%temp==0){
      fac_cnt[idx]++;
      x/=temp;
    }
  }
}

ll cal(ll n,ll p){
  ll res = 0;
  while(n){
    res += n/p;
    n/=p;
  }
  return res;
}

ll solve(){
  ll ans = 4e18;
  for(int i=1;i<=idx;i++){
    ll& p = fac[i];
    ll ty = cal(y,p);
    for(int j=1;j<=n;j++){
      ty -= cal(a[j],p);
    }
    ans = min(ans,ty/fac_cnt[i]);
  }
  return ans;
}

int main(){
  int t;cin >> t;
  while(t--){
    idx = 0;
    memset(fac_cnt,0,sizeof fac_cnt);
    scanf("%d%lld%lld",&n,&x,&y);
    for(int i=1;i<=n;i++)scanf("%lld",a+i); 
    chai(x);
    ll ans = solve();
    cout << ans << endl;
  }
}