BZOJ 2219 数论之神 BSGS+CRT

题意:链接

方法: BSGS+CRT

解析:

这道题有什么区别呢?

就是他取模的值不是一个质数了

这怎么办呢?

我们来把这个数分解质因数。

P=piti

然后对于每一个 piti 我们单独计算。

最后用中国剩余定理合并

这其实跟想象中的中国剩余定理没有什么关系,其实是其中的一个性质,就是 mod p 的解数总数量等于每一个mod piti 的乘积。

这是一个性质,别问我证明。

原方程 xA=B(mod C)

转化为 xA=B(mod Ci)>Ci 是质数或是质数的幂次。

我们可以单独计算。

这跟之前的题类似。

不过有一点需要注意的地方。

就是如果直接跑最后求乘积的话是错的。

为什么呢,因为这个解的范围显然改变了。原来为C,不过现在变成了每个Ci的范围,显然不对。

况且,对于B来说,其可能比Ci要大,显然是错误的。

那么对于一个Ci,我们需要知道什么呢,知道它的质因子Pi,幂次Ti

首先先将这个B取余一下Ci,让这个方程正常一点。

然后这里就有问题了。

如果B整除Ci怎么办呢?

整除Ci的话,那么我们不妨设所有答案的gcd为 Pik ,明确一下,如果要是解最多的话显然我们让这个k越小越好。

回归下方程 xA=0(mod PiTi)

PiAk=0(mod PiTi)

Ak>=Ti

所以这个k最小是多少呢? Ti1A+1

这应该是比较好想的?

然后这种情况我们就讨论完了。

现在来讨论B不为零的时候。

这时候我们可以把B换一个形式。

比如 PipW

然后呢,此时等式变成了什么

xA=PipW(mod PiTi)

既然这样的话,右边的 Pip 一定可以均分到左边对吧?

如果不均分的话,显然二者是不会同余的。

所以不均分的话显然该方程无解。

然后我们来考虑均分的情况

对于该方程,我们可以把左右及模数共有的部分给约去。

所以变成了 yA=W(mod PiTik)gcd(y,Pi)=1

显然所有的子方程都可以化成这个东西。

然后来考虑这个东西怎么求解。

mod=PiTik

然后我们可以求出模数的原根g;

然后呢我们可以求出来W在模mod下对g的指标。

不妨设为ind_W;

我们再设y在模mod下对g的指标为 indy

则原方程转化为

Aindy=indW(mod φ(mod))

为什么取原根的优越性就体现出来吧?

因为该方程的取余后的值是可以覆盖[1,φ(mod))的所有值,也就是可以将这个方程进一步转化。

至于这个方程怎么解?

很简单吧,这就是 Ax=B(mod C) 嘛,exgcd秒解。

但是这里有一个问题,也就是说我们如此解出来的方案的数量是否与原方程相同呢?

xA=PipW(mod PiTi)

yA=W(mod PiTik)gcd(y,Pi)=1

就这两个方程,我们看到,差异性就体现在模的数部分,这会有什么影响呢?

别忘了,如果我们求出来所有的解的话,正常的话,应该是在 (mod PiTi) 意义下的,但是现在变成了 (mod PiTik) 意义下的解..所以我们显然会少计算解,少计算了多少呢?

先看我们把x缩小了多少

显然是y显然是 xPipA 这么多,所以我们解出来的东西实际应该在 (mod PiTipA) 这个范围内但是我们解的范围是多少呢? PiTik 这个范围内的,并且我们知道,k=p即我们求出了 PiTip 这个范围内的

所以我们要将该方程最终解中乘以个 PiTipAPiTip PippA

然后我们解完所有的方程后合并就是最终解。

这题解我写的好累啊。

popoqqq爷有个例子,可以去看看,但是有点bug,细心的人自己发现吧。

代码:

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define N 140345
#define INF 0x7fffffffffffffffll
using namespace std;
typedef long long ll;
int T,cnt,tot,totc;
int head[N+10];
ll a,b,c;
ll prime[N];
ll the_prime_c[N];
ll cnt_num_c[N];
struct node
{
    ll from,to,val;
    int next;
}edge[N+10];
void cutit(ll x)
{
    totc=0;
    for(ll i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            the_prime_c[++totc]=i;
            cnt_num_c[totc]=1;
            x/=i;
            while(x%i==0){x/=i;cnt_num_c[totc]++;}
        }
    }
    if(x!=1)the_prime_c[++totc]=x,cnt_num_c[totc]=1;
}
void init()
{
    memset(head,-1,sizeof(head));
    cnt=1;
}
void edgeadd(ll from,ll to,ll val)
{
    edge[cnt].from=from,edge[cnt].to=to,edge[cnt].val=val;
    edge[cnt].next=head[from];
    head[from]=cnt++;
}
void exgcd(ll a,ll b,ll &x,ll &y,ll &gcd)
{
    if(!b)
    {
        x=1,y=0,gcd=a;
        return;
    }
    exgcd(b,a%b,y,x,gcd);
    y=y-a/b*x;
}
void get_factor(ll x)
{
    tot=0;
    for(int i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            prime[++tot]=i;
            while(x%i==0)x/=i;
        }
    }
    if(x!=1)prime[++tot]=x; 
}
ll quick_my(ll x,ll y,ll mod)
{
    ll ret=1;
    while(y)
    {
        if(y&1)ret=(ret*x)%mod;
        x=(x*x)%mod;
        y>>=1;
    }
    return ret;
}
ll get_inv(ll x,ll mod)
{
    ll X,Y,GCD;
    exgcd(x,mod,X,Y,GCD);
    return (X%mod+mod)%mod;
}
ll get_phi(ll x)
{
    ll ret=x,tmp=x;
    for(ll i=2;i*i<=tmp;i++)
    {
        if(tmp%i==0)
        {
            ret=ret*(i-1)/i;
            while(tmp%i==0){tmp/=i;}
        }
    }
    if(tmp!=1)ret=ret*(tmp-1)/tmp;
    return ret;
}
int check(ll x,ll MOD,ll PHI)
{
    for(int i=1;i<=tot;i++)
    {
        if(quick_my(x,PHI/prime[i],MOD)==1)return 0;
    }
    return 1;
}
ll find_primitive_root(ll x)
{
    ll tmp=get_phi(x);
    get_factor(tmp);
    for(int i=2;;i++)
    {
        if(check(i,x,tmp))return i;
    }
}
ll Baby_Step_Giant_Step(ll A,ll B,ll C)
{
    init();
    ll m=(int)ceil(sqrt(C));
    ll k=1;
    for(int i=0;i<m;i++)
    {
        int flag=1;
        for(int j=head[k%N];j!=-1;j=edge[j].next)
        {
            if(edge[j].val==k){flag=0;break;}
        }
        if(flag)edgeadd(k%N,i,k);
        k=k*A%C;
    }
    ll invk=get_inv(k,C);
    ll invD=1;
    for(int i=0;i<=m;i++)
    {
        ll tmpB=B*invD%C;
        for(int j=head[tmpB%N];j!=-1;j=edge[j].next)
        {
            if(edge[j].val==tmpB)return edge[j].to+i*m;
        }
        invD=invD*invk%C;
    }
    return -1;
}
ll work(ll a,ll b,ll pri,ll num)
{
    ll B=quick_my(pri,num,INF);b%=B;
    if(b==0)return quick_my(pri,num-((num-1)/a+1),INF);
    ll cnt_b=0;
    while(b%pri==0)
    {
        cnt_b++,b/=pri,num--;
    }
    if(cnt_b%a!=0)return 0;
    ll g=find_primitive_root(B);
    ll indB=Baby_Step_Giant_Step(g,b,B);
    ll X,Y,GCD;
    exgcd(a,get_phi(B),X,Y,GCD);
    // A*ind_g(A)=ind_g(B)(mod φ(C1));
    if(indB%GCD!=0)return 0;
    return GCD*quick_my(pri,cnt_b-cnt_b/a,INF);
}
int main()
{
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%lld%lld",&a,&b,&c);
        c=c*2+1;
        cutit(c);
        ll ans=1;
        for(int i=1;i<=totc;i++)
        {
            ans*=work(a,b,the_prime_c[i],cnt_num_c[i]);
        }
        printf("%lld\n",ans);
    }
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值