二次同余方程模合数的一般解法

0.不讨论复杂情况的解释
对于一般的二次同余方程形如

可以通过配方化为下式

可通过换元，得到

解出X的取值，然后用2ax+b回带，用扩展欧几里得解线性同余方程就可以得到方程本来的解
注意上式得到了一般二次同余方程的形式，即

下文讨论上式的一般解法，并给出代码进行解释

1.朴素的二次同余方程的解法
因为，因此很明显可以直接枚举。时间复杂度，如果m稍微大一些效率就会很低

2.对于1的第一个优化
由数论知识可得，对于0中的方程，我们可以分解m的质因数，即

对于每一个方程，用1中的朴素解法来暴力枚举，时间复杂度优化到了
然后对于同余方程组的每一组解都做一个中国剩余定理就可以求得x

3.对于1的第二个优化
由数论知识可得，方程 的解数为（假设方程有解）：
当mi为奇质数时，方程有2个解，假设一个是x0，则另一个是mi^pi-x0
当mi等于2且pi大于2时，方程有4个解，假设一个是x0，则另外三个是mi^pi-x0，x0+mi^(pi-1),mi^(pi-1)-x0
因此，如果原方程有解，那么只要枚举到一个解，那么就可以返回了。

4.对于2的优化
由2得，即使分解了质因数，但要枚举质因数的幂还是太困难，还是容易超时
对于2中同余方程组的一个方程 ，考虑这个方程 ，假设我们求出了这个方程的一个解x0，那么就可以解出这个方程的一个解。
具体方法是在求出了x0之后，有
这是一个线性同余方程了，解出k的值，可得，x1是方程 的一个解。
这样一直推下去，可以得到原方程的一个解。
利用3中结论，可得另外的解。
而解一个二次同余方程模质数，还是采用1中的枚举法，这样时间复杂度是

5.对于4的优化
实际上，解一个二次同余方程模奇质数的方法以及很普遍了，不懂的可以看看大神acdreamers的博客http://blog.csdn.net/acdreamers/article/details/10182281。
对于模2的情况，这个直接枚举就行了。
有一个小优化，详见代码。

6.代码
因为现在我没有在网上找到哪个题可以交二次同余方程模合数的。。因此现在代码可能有bug。

#include<cstdio>
#include<cstdlib>
#include<set>
#include<algorithm>
#define ll long long
#define pll pair<ll,ll>
using namespace std;
ll fac[70][2],eqans[70][5],tempa[70],m[70],ans[1000100],n,a,cnt,w;
void frac(ll n)//分解模
{
    for(int i=2;i*i<=n;++i)
        if(n%i==0)
        {
            while(n%i==0)
            {
                fac[cnt][0]=i;
                ++fac[cnt][1];
                n/=i;
            }
            ++cnt;
        }
    if(n>1)
    {
        fac[cnt][0]=n;
        ++fac[cnt++][1];
    }
}
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 mul(ll a,ll b,ll p)
{
    ll res=0;
    while(b)
    {
        if(b&1)res=(res+a)%p;
        a=(a+a)%p;
        b>>=1;
    }
    return res;
}
ll pow(ll a,ll b,ll p)
{
    ll res=1%p;
    while(b)
    {
        if(b&1)res=mul(res,a,p);
        a=mul(a,a,p);
        b>>=1;
    }
    return res;
}
pll mul(pll a,pll b,ll p)//二次域上乘法
{
    pll res;
    res.first=(a.first*b.first%p+a.second*b.second%p*w%p)%p;
    res.second=(a.first*b.second%p+a.second*b.first%p)%p;
    return res;
}
pll pow(pll a,ll b,ll p)//二次域上快速幂
{
    pll res=pll(1,0);
    while(b)
    {
        if(b&1)res=mul(res,a,p);
        a=mul(a,a,p);
        b>>=1;
    }
    return res;
}
ll work(ll a,ll p)//解二次同余方程模奇质数
{
    if(a==0)return 0;
    if(pow(a,(p-1)/2,p)==p-1)return -1;
    ll t;
    while(1)
    {
        t=rand()%p;
        w=((t*t-a)%p+p)%p;
        if(pow(w,(p-1)/2,p)==p-1)break;
    }
    pll tmp=pll(t,1);
    return pow(tmp,(p+1)/2,p).first;
}
ll work2(ll a)//枚举二次同余方程模8同余的解
{
    for(int i=0;i<8;++i)
        if(i*i%8==a%8)
            return i;
    return -1;
}
ll solve(ll a,ll p,ll k)//求二次同余方程模质数的幂的一个解
{
    ll prep=1,nowp=p,w;
    if(p==2)//质数为2的时候，a=0则0是一个解，a=1则1是一个解，a=2 or a=3无解
    {
        if(a==0)return 0;
        if(a==1)return 1;
        if(a==2||a==3)return -1;
        w=work2(a);prep=2;nowp=8;k-=2;//此时只用讨论模8(模4的已经讨论完了)
    }
    else w=work(a%p,p);
    if(w==-1)return -1;
    for(int i=1;i<k;++i)//模质数升幂
    {
        prep*=p;
        nowp*=p;
        ll d,x,y,c;
        c=((a-w*w)%nowp+nowp)%nowp;
        exgcd(2*w*prep,nowp,d,x,y);
        if(c%d==0)w=(w+(c/d*x%nowp+nowp)%(nowp/d)*prep)%nowp;
        else return -1;
    }
    return w;
}
ll CRT()//合并方程组
{
    ll res=0;
    for(int i=0;i<cnt;++i)
    {
        ll x,y,d,Mi=n/m[i];
        exgcd(Mi,m[i],d,x,y);
        res=(res+Mi*x*tempa[i])%n;
    }
    if(res<0)res+=n;
    return res;
}
void dfs(int pos)//决定选取哪些方程做CRT
{
    if(pos==cnt)
    {
        ans[++ans[0]]=CRT();
        return;
    }
    m[pos]=pow(fac[pos][0],fac[pos][1],n+1);
    for(int i=0;i<eqans[pos][19];++i)
    {
        tempa[pos]=eqans[pos][i];
        dfs(pos+1);
    }
}
void quadratic_congruence()//二次同余主函数
{
    for(int i=0;i<cnt;++i)//对于每个质数的幂求二次同余
    {
        ll p=pow(fac[i][0],fac[i][1],n+1),t=solve(a%p,fac[i][0],fac[i][1]);
        if(t==-1)//如果有一个无解就返回
        {
            puts("No Solution");
            return;
        }
        set<ll>s;
        set<ll>::iterator iter;
        if(t==0)//以前的程序有bug,如果整除的话那么不只有2或4个解
        {
            ll pp=pow(fac[i][0],(fac[i][1]+1)/2,n+1);
            for(int i=0;pp*i<p;++i)
                s.insert(pp*i);
        }
        else
        {
            s.insert(t);
            s.insert(p-t);
            if(fac[i][0]==2&&fac[i][1]>1)s.insert((t+p/2)%p);
            if(fac[i][0]==2&&fac[i][1]>1)s.insert(((p/2-t)+p)%p);
        }
        for(iter=s.begin();iter!=s.end();iter++)
            eqans[i][eqans[i][19]++]=*iter;
    }
    dfs(0);
}
int main()
{
    scanf("%I64d%I64d",&a,&n);
    frac(n);
    quadratic_congruence();
    sort(ans+1,ans+ans[0]+1);
    for(int i=1;i<=ans[0];++i)
        printf("%I64d%c",ans[i],i==ans[0]?'\n':' ');
}