BSGS/exBSGS 学习笔记【同余】

作者: wjyyy 分类: BSGS,分块,同余,学习笔记 发布时间: 2019-06-06 16:43

点击量:48

〇、前言

BSGS 是用来求形如 $a^x\equiv b\pmod p$ 的最小非负整数解的算法。

当 $\gcd(a,p)=1$ 时可以直接使用 BSGS,当 $\gcd(a,p)\ne 1$ 时需要用到 exBSGS 进行转化。

参考资料:

PoPoQQQ 《原根与指标》

ButterflyDew 数论ex

一、分块

对于 $a^x$,在 $\text{mod}$ 质数 $p$ 的意义下有 $\{1,2,\ldots,p-1\}$ 共 $p-1$ 种取值。

考虑把 $a^0,a^1,\ldots,a^{p-2}$ 分成 $\left\lceil\frac p{\lfloor\sqrt p\rfloor}\right\rceil$ 段,每一段的长度为 $\lfloor\sqrt p\rfloor$。预处理出 $a^1,a^2,\ldots,a^{\lfloor\sqrt p\rfloor}$ 和 $a^{\sqrt p},a^{2\sqrt p},\ldots,a^{\left\lceil\frac p{\lfloor\sqrt p\rfloor}\right\rceil \sqrt p}$,就可以用这些元素组成 $a^0,a^1,\ldots,a^{p-2}$ 中的任意一个数了。

即把 $x$ 表示成 $x=k\lfloor\sqrt p\rfloor-r$ 的形式,减法可以防止后面求逆元的麻烦。

二、meet-in-middle

在 $a^x\equiv b\pmod p$ 中,我们已知 $a,b,p$,因此 $a^x$ 的值是已知的。

在上面一段中,把 $x$ 表示成 $x=k\lfloor\sqrt p\rfloor-r$,则 $a^{k\lfloor\sqrt p\rfloor-r}\equiv b\pmod p$。

移项得 $a^{k\lfloor\sqrt p\rfloor}\equiv a^rb\pmod p$。

如果我们把 $a^1b,a^2b,\ldots,a^{\lfloor\sqrt p\rfloor}b$ 放在哈希表 $A$ 中,然后遍历 $a^{\lfloor\sqrt p\rfloor},a^{2\lfloor\sqrt p\rfloor},\ldots,a^{\left\lceil\frac p{\lfloor\sqrt p\rfloor}\right\rceil \sqrt p}$。

可以看出,等式右边已经被存在表中了,现在考虑等式左边。如果存在 $a^x\equiv b\pmod p$,那么就一定存在一个 $k$,可以在哈希表中找到和 $a^{k\lfloor\sqrt p\rfloor}$ 相等的元素。

现在要解决的问题是找到最小非负整数解。

对于大块来说,$a^{\lfloor\sqrt p\rfloor},a^{2\lfloor\sqrt p\rfloor},\ldots,a^{\left\lceil\frac p{\lfloor\sqrt p\rfloor}\right\rceil \sqrt p}$ 从小到大枚举就可以了,在哈希表中查到的顺序就不一定了。

事实上,我们用链表或前向星作为哈希表,从小到大插入 $a^1,a^2,\ldots,a^{\lfloor\sqrt p\rfloor}$。在同一链表中如果有多个符合条件的,那么就选择最后面一个;在同一前向星表中如果有多个符合条件的,那么取最前面一个。

三、BSGS

当 $\gcd(a,p)=1$ 时,$a^0,a^1,\ldots,a^{p-2}$ 刚好是 $p-1$ 个不同的数,此时直接用上面的 meet-in-middle 就可以解决了。

注意一点,我的代码风格习惯把在小段把 $a$ 滚动到 $a^{\lfloor\sqrt p\rfloor}$ 然后直接用。此时如果滚动的是 $a^ib$,那么后面还要除以 $b$。这时 $b$ 可能没有模 $p$ 的逆元(如 $b=0$),因此只滚动 $a^i$ 比较好,每次手动乘 $b$ 不会很慢。

注意一点,我的代码风格习惯把在小段把 $a$ 滚动到 $a^{\lfloor\sqrt p\rfloor}$ 然后直接用。此时如果滚动的是 $a^ib$,那么后面还要除以 $b$。这时 $b$ 可能没有模 $p$ 的逆元(如 $b=0$),因此只滚动 $a^i$ 比较好,每次手动乘 $b$ 不会很慢。

代码 of POJ 2417

#include<cstdio>
#include<cstring>
#include<map>
#include<cstdlib>
#include<cmath>
#define ll long long
#define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
char buf[2100000],*p1=buf,*p2=buf;
int read()
{
    int x=0;
    char ch=gc();
    while(ch<'0'||ch>'9')
    {
        if(ch==EOF)
            exit(0);
        ch=gc();
    }
    while(ch>='0'&&ch<='9')
    {
        x=x*10+(ch&15);
        ch=gc();
    }
    return x;
}
ll qpow(ll x,ll y,ll p)
{
    ll ans=1;
    while(y)
    {
        if(y&1)
            ans=ans*x%p;
        x=x*x%p;
        y>>=1;
    }
    return ans;
}
struct edge//前向星
{
    int n,nxt,x,v;
}e[100000];
int head[10007],ecnt=-1,cnt=0;
void add(int bs,int x,int v)
{
    e[++ecnt]=(edge){++cnt,head[bs],x,v};
    head[bs]=ecnt;
}
int Find(int n,int x)//找到符合条件的第一个就可以返回答案了
{
    for(int i=head[n];~i;i=e[i].nxt)
        if(e[i].x==x)
            return e[i].v;
    return 0;
}
int main()
{
#ifdef wjyyy
    freopen("a.in","r",stdin);
#endif
    int p,b,n;
    while(1)
    {
        for(int i=0;i<10007;++i)
            head[i]=-1;
        ecnt=-1;
        p=read();
        b=read();
        n=read();
        int t=sqrt(p),flag=0;
        ll f=n;//这里最好不要这样写 详见下一份代码
        for(int i=1;i<=t;++i)//插入 a^1,a^2,…
        {
            f=f*b%p;
            add(f%10007,f,i);
        }
        ll g=f=f*qpow(n,p-2,p)%p,k;
        for(ll q=t;q-t<p;q+=t,g=g*f%p)//在 a^{√p},a^{2√p},… 中找符合条件的答案
            if(k=Find(g%10007,g))
            {
                printf("%lld\n",q-k);
                flag=1;
                break;
            }
        if(!flag)
            puts("no solution");
    }
    return 0;
}

四、exBSGS

如果 $\gcd(a,p)\ne 1$,那么 $a^{k\lfloor\sqrt p\rfloor-r}\equiv b\pmod p$ 就不能把 $a^{-r}$ 移到右边了。

需要想办法把原式转化成底数模数互质的情况。

根据同余,$a^x\equiv b\pmod p\Leftrightarrow\exists k\in \mathbb Z,a^x+kp=b$。

由于 $\gcd(a,p)\ne 1$,令 $g=\gcd(a,p)$,则 $g|a^x+kp$,当后面一个式子成立时,必然有 $g|b$。

如果此时 $b\nmid a^x+kp$,那么同余方程无解,原方程无解。

否则等式两边同除以 $g$,得到 $\frac{a^x}g+\frac{kp}g=\frac bg$。

由于产生影响的是 $a$,不是 $a^x$,所以此时单独提出来一个 $\frac ag$ 作为“系数”,继续判断是否互质;同理,在 $\frac{kp}g$ 中产生影响的也只有 $p$,所以式子是这样的:$\frac aga^{x-1}+k\frac pg=\frac bg$。

转化回来就是 $\frac aga^{x-1}\equiv\frac bg\pmod{\frac pg}$。又回到了 $a^x\equiv b\pmod p$ 的形式。

一直这样转化,会出现上述无解的情况,或在某一步满足了 $\gcd(a,p)=1$。

这时的式子是
$$
\frac{a^k}{g_1\times g_2\times \cdots\times g_k}a^{x-k}\equiv\frac b{g_1\times g_2\times \cdots\times g_k}\pmod{\frac c{g_1\times g_2\times \cdots\times g_k}}
$$
是 $A\cdot a^{x-k}\equiv B\pmod C$ 的形式。

这样的话直接用 BSGS 的方法就可以解决了。

此外,在这个式子中出现了 $a^{x-k}$,我们无法保证 $x-k\ge 0$。但是如果 $g_i>1$,则 $k\le \lfloor\log_2 C\rfloor$。

这样的话我们可以在程序开始的时候对 $i\in[0,\lfloor\log_2 C\rfloor]$检验 $a^i\equiv b\pmod c$ 是否成立,如果有成立可以直接输出。

这样做还有一个好处。当 $x=0$ 且 $b=1$ 时不受 $\gcd(a,p)$ 是否 $=1$ 的限制。可以直接在检验的时候判掉,不用特判了。注意有的题如果没有定义 $0^0=1$,那么认为 $0^x\equiv 1$ 无解;$0^x\equiv 0$ 的最小非负整数解是 $x=1$。

前面处理到互质的复杂度是 $O(\log n)$,后面是 BSGS,因此总复杂度为 $O(\sqrt n)$。

代码 of POJ 3243

#include<cmath>
#include<cstdio>
#include<cstring>
#define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#define ll long long
#define con x=read();z=read();k=read();continue
char buf[2100000],*p1=buf,*p2=buf;
int read()
{
    int x=0;
    char ch=gc();
    while(ch<'0'||ch>'9')
        ch=gc();
    while(ch>='0'&&ch<='9')
    {
        x=x*10+(ch&15);
        ch=gc();
    }
    return x;
}
struct edge
{
    int x,nxt,v;
}e[100000];
int head[10007],ecnt=-1;
void add(int x,int v)
{
    e[++ecnt]=(edge){x,head[x%10007],v};
    head[x%10007]=ecnt;
}
ll qpow(ll x,ll y,ll p)
{
    ll ans=1;
    while(y)
    {
        if(y&1)
            ans=ans*x%p;
        x=x*x%p;
        y>>=1;
    }
    return ans;
}
int Find(int x)
{
    for(int i=head[x%10007];~i;i=e[i].nxt)
        if(e[i].x==x)
            return e[i].v;
    return 0;
}
int bsgs(int A,int a,int b,int p)
{
    for(int i=0;i<10007;++i)
        head[i]=-1;
    ecnt=-1;
    int t=sqrt(p);
    ll f=1;
    for(int i=1;i<=t;++i)//maybe b doesn't have an inversion
    {
        f=f*a%p;
        add(f*b%p,i);
    }
    ll g=f*A%p,k;
    for(ll q=t;q-t<p;q+=t,g=g*f%p)
        if(k=Find(g))
            return q-k;
    return -1;
}
int gcd(int x,int y)
{
    while(y)
    {
        int t=y;
        y=x%y;
        x=t;
    }
    return x;
}
int main()
{
#ifdef wjyyy
    freopen("a.in","r",stdin);
#endif
    int x,z,k;
    x=read();
    z=read();
    k=read();
    while(x|z|k)
    {
        ll t=1%z;
        int flag=0;
        for(int i=0;i<=31;++i)
        {
            if(t==k)
            {
                printf("%d\n",i);
                flag=1;
                break;
            }
            t=t*x%z;
        }
        if(flag)
        {
            con;
        }
        int g=gcd(x,z),cnt=0;
        ll A=1;
        while(g>1)
        {
            if(k%g)
            {
                puts("No Solution");
                flag=1;
                break;
            }
            A=A*(x/g)%z;
            ++cnt;
            z/=g;
            k/=g;
            g=gcd(x,z);
        }
        if(flag)
        {
            con;
        }
        int ans=bsgs(A,x,k,z);
        if(ans==-1)
            puts("No Solution");
        else
            printf("%d\n",ans+cnt);
        con;
    }
    return 0;
}

说点什么

avatar
  Subscribe  
提醒
/* */