高次同余笔记(三):离散对数和原根

我们来看这个方程:
这里写图片描述
a,b,c在int内,c是质数。
求x在[0,c-1]内所有的解。
这个怎么搞?
那我们换个方程:
这里写图片描述
这个方程 的解很明显是
这里写图片描述
但是我们换个角度,因为开根号这个操作数论里面不好搞。
这样有
这里写图片描述
通过这个式子可以算出
这里写图片描述
这样x就是一个指数式子,不含根号了。
但是e这个东西数论里面没有啊。
注意到这个底数的选取是任意的,那么随便选一个底数记作G。
这样有这里写图片描述
但是这个底数的选取也是有要求的,起码这个底数使得x属于[1,c-1]的时候G^x%c要取遍[1,c-1]的所有数。不然换成原来的指数式子是有问题的。
比如下面一个表,当G=3,c=7的时候:

  • x 1 2 3 4 5 6
  • G^x 3 2 6 4 5 1

然后,符合条件的G被称为模c的原根。
这样取对数的操作就是离散对数了。
然后我们再看logGb这个怎么求。
x=logGb->G^x=b(mod c)
这个是什么东西?
不就是之前讲过的baby-step-giant-step吗?
这样方程中的a已知,logGb可求,c已知,这样原方程被转化成一个线性同余式子,用exgcd可求logGx。
知道了logGx后,x=G^logGx,用快速幂可以求出x。
如果上面那个线性同余有多解,那么最后的x也是有多解的。
思路就是这样。
下面来说原根怎么求。
G是模c的原根的充要条件是G^(c-1) = 1 (mod c)在G < c的时候当且仅当指数为c-1的时候成立。
这样容易想到枚举一个G,对于每一个G,枚举i从0到c-1,计算G^(c-1),如果前面就出现1了那么G不是模c的原根。
有一个优化,就是i的枚举不需要从0到c-1,只需要枚举(c-1)/c的一个质因子,对于每个质因子都这么枚举就快了。

扩展是当c为合数的时候,这个我还没弄懂怎么CRT合并,哪位大神来讲解一下。
贴一个HDU3930

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<map>
#define ll long long
using namespace std;
const int N=1000010;
ll fac[70],ans[N],a,b,c;
bool isprime[N];
int prime[N/12],pcnt,fcnt,cas;
inline void getprime(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!isprime[i])prime[pcnt++]=i;
        for(int j=0;j<pcnt&&prime[j]*i<=n;j++)
        {
            isprime[prime[j]*i]=true;
            if(i%prime[j]==0)break;
        }
    }
}
inline void divide(ll n)
{
    fcnt=0;
    ll t=sqrt(n);
    for(int i=0;prime[i]<=t&&i<pcnt;++i)
        if(n%prime[i]==0)
        {
            fac[fcnt++]=prime[i];
            while(n%prime[i]==0)n/=prime[i];
        }
    if(n>1)fac[fcnt++]=n;
}
inline ll mul(ll n,ll m,ll p)
{
    ll ans=0;
    while(m)
    {
        if(m%2)ans=(ans+n)%p;
        n=2*n%p;
        m/=2;
    }
    return ans;
}
inline ll pow(ll n,ll m,ll p)
{
    ll ans=1%p;
    while(m)
    {
        if(m%2)ans=mul(ans,n,p);
        n=mul(n,n,p);
        m/=2;
    }
    return ans;
}
ll primitive_root(ll p)
{
    divide(p-1);
    for(int r=2;;++r)
    {
        bool flag=1;
        for(int i=0;i<fcnt;++i)
        {
            ll t=(p-1)/fac[i];
            if(pow(r,t,p)==1)
            {
                flag=0;
                break;
            }
        }
        if(flag)return r;
    }
}
void exgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
    if(!b){x=1;d=a;y=0;}
    else
    {
        exgcd(b,a%b,d,y,x);
        y-=a/b*x;
    }
}
ll BSGS(ll a,ll b,ll p)
{
    map<ll,int>mp;
    ll m=ceil(sqrt(p)),d=1,val=1,g,x,y;
    for(int i=0;i<m;++i)
    {
        mp.insert(make_pair(val,i));
        val=val*a%p;
    }
    for(int i=0;i<m;++i)
    {
        exgcd(d,p,g,x,y);
        x=(b/g*x%p+p)%(p/g);
        if(mp.count(x))return i*m+mp[x];
        d=d*val%p;
    }
    return -1;
}
void work(ll a,ll b,ll c)
{
    ll rt=primitive_root(c);
    ll logrtb=BSGS(rt,b,c);
    ll d,x,y;
    exgcd(a,c-1,d,x,y);
    if(logrtb%d){puts("-1");return;}
    ll t1=logrtb/d,t2=(c-1)/d;
    ans[0]=(x*t1%t2+t2)%t2;
    for(int i=1;i<d;++i)
        ans[i]=ans[i-1]+t2;
    for(int i=0;i<d;++i)
        ans[i]=pow(rt,ans[i],c);
    sort(ans,ans+d);
    for(int i=0;i<d;++i)
        printf("%I64d\n",ans[i]);
}
int main()
{
    getprime(N);
    while(~scanf("%I64d%I64d%I64d",&a,&c,&b))
    {
        printf("case%d:\n",++cas);
        work(a,b,c);
    }
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值