给出一组正整数$x,n,r$，使得$r^2\equiv x(mod\: n)$，求出所有满足该等式的$r$。

假设有另一个解$r'$满足条件，则有$r^2-r'^2=kn$

因式分解，得$(r+r')(r-r')=kn$

将$n$分解成$a*b$，则有$\left\{\begin{matrix}r+r'=xa\\ r-r'=yb\end{matrix}\right.$

两式相加得$2r=xa+yb$，这是一个二元线性不定方程，可用扩欧求出x的通解。

假设已经求出了$x$的通解$x=x_{0}+k\Delta x$，

由于$r+r'=xa$，所以$r'=xa-r=(x_{0}+k\Delta x)a-r=x_{0}a-r+k(a\Delta x)$,

 设$\Delta t=a\Delta x$，则$r'_{0}=((x_{0}a-r)\%\Delta t+\Delta t)\%\Delta t为$r'$的第一个非负整数解

 因此$r'$的通解为$r'=r'_{0}+k\Delta t$

枚举所有的$a,b$，将所有$r'$的可行解插入一个集合里就行了。

 1 #include<bits/stdc++.h>
 2 
 3 using namespace std;
 4 typedef long long ll;
 5 ll x,n,r,ka;
 6 set<ll> st;
 7 void exgcd(ll a,ll b,ll& x,ll& y,ll& g) {
 8     if(!b)x=1,y=0,g=a;
 9     else exgcd(b,a%b,y,x,g),y-=x*(a/b);
10 }
11 
12 void solve(ll a,ll b) {
13     ll c=2*r,x,y,g;
14     exgcd(a,b,x,y,g);
15     if(c%g)return;
16     x*=c/g;
17     ll dx=abs(b/g);
18     ll dt=dx*a;
19     ll t=((a*x-r)%dt+dt)%dt;
20     for(; t<n; t+=dt)st.insert(t);
21 }
22 
23 int main() {
24     while(scanf("%lld%lld%lld",&x,&n,&r)&&x) {
25         st.clear();
26         for(ll i=1; i*i<=n; ++i)if(n%i==0)solve(i,n/i),solve(n/i,i);
27         printf("Case %lld:",++ka);
28         for(ll i:st)printf(" %lld",i);
29         printf("\n");
30     }
31     return 0;
32 }