郑州大学2024年寒假训练 Day7：数论

感觉这一块讲的有点太少了，只有辗转相除法，拓展欧几里得定理，素数筛，快速幂和逆元五个内容。数论的内容远远不止这些。不过一个视频也讲不了太多东西，讲的还是数学，也是没有办法。一边看题一边说吧。

辗转相除法：

辗转相除法又叫欧几里得算法，用来求解两个数的最大公因数（GCD，Greatest Common Divisor）。内容是 $gcd(a,b)=gcd(b,a\%b)$。中国有一个类似的算法叫做更相减损术，内容是 $gcd(a,b)=gcd(a,a-b)=gcd(b,a-b)$ 其实内容上和欧几里得算法是一个东西。

证明 $gcd(a,b)=gcd(b,a\%b)$：
不妨假设 $a>b$。设 $d=gcd(a,b)$ ，那么 $a=k_1\ast d,b=k_2\ast d$ ， 则有 $a\%b=a-m\ast b=k_1\ast d-m\ast k_2\ast d=(k_1-m\ast k_2)\ast d$，余数仍然是最大公因数 $d$ 的倍数。
还有种证明方法是这样的：$a=q\ast b+r$，其中 $0\le r<b$，那么 $a\%b=r$。因为 $a,b$ 的最大公因数 $d$ 可以整除 $a,b$，即 $d|a$ 且 $d|b$ ，所以 $d|(q\ast b)$，所以 $d|(r=a-q\ast b)$，也就是说明余数 $r$ 是 $d$ 的倍数。

而 $(a,b)\to (b,a\%b)$ 是不断减小的，因此这样算下去的话，最终就会得到 $a\%b=0$，此时的 $a$ 就是最大公因数。这个辗转相除的过程最多不会超过 $log$ 次，因此求解 $gcd$ 的时间复杂度是 $O(logn)$

代码片段如下：
这里不需要判断ab谁大谁小并交换，因为再调用一次 gcd(b,a%b) 就实现了交换。

int gcd(int a,int b){return (b)?gcd(b,a%b):a;
}

更神秘的写法：
其实就是每次先对a模上b，然后交换ab，b^=a^=b^=实现的其实就是交换功能。

int gcd(int,a,int b){while(b)b^=a^=b^=a%=b;return a;
}

---

拓展欧几里得定理：

也就是exgcd（Extended gcd）。它主要是用来求解一次线性同余方程。也就是形如 $a\ast x\equiv c(\mathrm{mod}p)$ 求 $x$ 的通解问题。这一篇讲的比较好

对这种问题，我们可以稍加转化，得到 $a\ast x+p\ast y=c$ 的等式，换种写法就是 $a\ast x+b\ast y=c$，变成了二元一次方程求通解问题。在解决这个问题之前，我们可以算出 $a\ast x+b\ast y=gcd(a,b)$ 问题的通解，然后通过这个算式可以得到我们所求算式的通解。那么怎么求 $a\ast x+b\ast y=gcd(a,b)$ 的通解呢。

贝祖定理（$B\acute{e}zout$ 定理）

内容是 ：对任意整数 $a,b$，存在一对整数 $x,y$，满足方程 $a\ast x+b\ast y=gcd(a,b)$ 。

证明使用归纳法。当 $b=0$ 时，$x=1,y=0$ 显然是方程的一组解（$gcd(a,0)=a$）。

当 $b>0$，设我们知道方程 $b\ast x+(a\%b)\ast y=gcd(b,a\%b)=gcd(a,b)$ 的解，求 $a\ast x+b\ast y=gcd(a,b)$ 的解。$$b\ast x+(a\%b)\ast y=gcd(b,a\%b)$$$$b\ast x+(a-b\ast \lfloor \frac{a}{b}\rfloor )\ast y=gcd(a,b)$$$$a\ast y+b\ast x-b\ast \lfloor \frac{a}{b}\rfloor \ast y=gcd(a,b)$$$$a\ast y+b\ast (x-\lfloor \frac{a}{b}\rfloor \ast y)=gcd(a,b)$$

设 $b\ast x+(a\%b)\ast y=gcd(b,a\%b)$ 一组特解是 $(x_0,y_0)$ ，可以发现，$a\ast x+b\ast y=gcd(a,b)$ 一组可行解就是 $x=y_0,y=x_0-\lfloor \frac{a}{b}\rfloor \ast y_0$。我们用辗转相除法向下走到尽头，从最下面的 $a\ast 1+b\ast 0=gcd(a,0)$ 一步一步向上推理，就可以得到原本式子的一组特解。

得到特解后会想得到 $a\ast x+b\ast y=gcd(a,b)$ 的通解，这个好说， $a\ast x$ 多一个 $\frac{a\ast b}{gcd(a,b)}$，那么 $b\ast y$ 就少一个 $\frac{a\ast b}{gcd(a,b)}$ 就行了。也就是说，$x$ 每加一个 $\frac{b}{gcd(a,b)}$，$y$ 就减少一个 $\frac{a}{gcd(a,b)}$，可以维持等式不变。写成式子就是，通解 $\overline{x}=x+k\ast \frac{b}{gcd(a,b)},\overline{y}=y-k\ast \frac{a}{gcd(a,b)}$。令 $t1=\frac{b}{gcd(a,b)}$，那么就是 $\overline{x}=x+k\ast t1$，这相当于通解模 $t1$ 同余特解，即 $\overline{x}\equiv x(\mathrm{mod}t1)$，因此求 $x$ 的最小非负整数解就是 $(x\%t1+t1)\%t1$，同理 $y$ 的通解。

代码片段如下：

int exgcd(int a,int b,int &x,int &y){if(!b){x=1;y=0;return a;}int d=exgcd(b,a%b,x,y);int z=x;x=y;y=z-(a/b)*y;return d;
}

压行：

int exgcd(int a,int b,int &x,int &y){if(!b){x=1;y=0;return a;}int d=exgcd(b,a%b,y,x);//这里相当于交换了xy的值y-=a/b*x;return d;
}

另外提一嘴，对任意一个方程 $a\ast x+b\ast y=d=gcd(a,b)$ ，等式两边同时除以最大公因数 $d$，就得到了 $\frac{a}{d}\ast x+\frac{b}{d}\ast y=1$ ，解和上面是一样的。也就是说当 $a,b$ 互质时 $a\ast x+b\ast y=1$ 一定有解。反之也成立（不然你 $a,b$ 都是 $d$ 的倍数，$x,y$ 是整数，它们乘积的和也一定是 $d$ 的倍数，但是 $1$ 不包含任何大于1的因数，因此 $a,b$ 不可能同时是任何一个 $d$ 的倍数）。

拓展欧几里得值域分析：

因为在这里面无法取模，否则会导致出错，所以会担心运算过程中出现超大的数导致溢出。不过安心，所有x和y的中间值的绝对值都小于 ∣ m a x { a , b } ∣ |max\{a,b\}| ∣max{a,b}∣。证明还是数学归纳法。详细的证明步骤在上面给出的博客中详细证明了。不多赘述。

裴蜀定理

内容：设 $a,b$ 为正整数，则关于 $a,b$ 的方程 $a\ast x+b\ast y=c$ 有整数解当且仅当 $c$ 是 $gcd(a,b)$ 的倍数。

也就是说如果 $c$ 不能被 $gcd(a,b)$ 整除，方程无整数解，否则有整数解。不证明了。

我们通过上面的方式得到了 $a\ast x+b\ast y=gcd(a,b)$ 一组特解，现在要得到 $a\ast x+b\ast y=c$ 通解。先算特解，再找到通解的之间的变化量即可。这个好说，假设 $c=k\ast d$ ，那么在等式 $a\ast x+b\ast y=d=gcd(a,b)$ 两边同时乘以 $k=\frac{c}{d}$ 就完事了，式子就变成了 $a\ast (\frac{c}{d}\cdot x)+b\ast (\frac{c}{d}\cdot y)=c$。新的特解就是 $x=\frac{c}{d}\cdot x,y=\frac{c}{d}\cdot y$。不过由于系数没变，所以变化量是不变的，还是 $t1=\frac{b}{d},t2=\frac{a}{d}$。

总结，通过拓展欧几里得算法可以算出等式 $a\ast x+b\ast y=d=gcd(a,b)$ 一组特解 $x_0,y_0$，那么等式 $a\ast x+b\ast y=c$ 一组特解是 $x=c/d\ast x_0,y_0=c/d\ast y_0$，通解 $\overline{x}$ 的变化量是 $t1=b/d$，通解 $\overline{y}$ 的变化量是 $t2=a/d$。等式 $a\ast x+b\ast y=c$ 的通解为 $\overline{x}=x+k\ast t1,\overline{y}=y-k\ast t2$

---

素数筛

试除法（开方暴力）

我们想要知道一个数 $x$ 是不是素数，朴素的想法是枚举 $2\sim x-1$ ，暴力地尝试 $x$ 能不能被其中一个整除。这样是 $O(n)$ 的。一个显而易见的性质是，如果一个数 $x$ 是合数，那么它一定存在一个质因数 $p$，而且 $2\le p\le \sqrt{x}$。

反证法证明：假设合数 $x$ 最小质因数为 $p$，如果上述性质不满足条件，说明 $p>\sqrt{x}$，因为 $x$ 为合数，因此 $x$ 存在至少一个另一个质因数 $q$，根据题设 $p$ 是最小的质因数，因此 $q>=q$，那么有 $p\ast q\ge p^2>(\sqrt{x}{)}^2=x$，两个因数的乘积大于了这个数，矛盾，因此 $2\le p\le \sqrt{x}$。

所以我们枚举 $2\sim \lfloor \sqrt{x}\rfloor$ 就行了，时间复杂度是 $O(\sqrt{n})$。代码如下：

bool isp(int x){if(x==1)return false;for(int i=2;i*i<=x;i++)if(x%i==0)return false;return true;
}

埃式筛

我们从 $2$ 往后看每个数，如果这个数是质数，直接去标记一下它所有的倍数，这样就可以把有这个质因数的合数标记掉。如果一个数是合数，那么它一定存在小于它的质因子。如果我们对每个遇到的质因数都这样做，那么我们在看到一个合数之前，它一定会被标记掉。我们从前往后看到的所有没有标记的数都是素数。所以我们从 $1$ 到 $n$ 走一遍这种操作，就可以得到 $1\sim n$ 中所有的素数。

这种像筛子一样一次一次把合数筛出来的过程就是埃式筛，它筛数的次数约有：$n/2+n/3+n/5+\cdots +n/p+\dots (p是质数)\le n/1+n/2+n/3+\cdots +n/n$ ，当 $n\to \infty$ 时，$n/1+n/2+n/3+\cdots +n/n=nlogn$（这东西学名叫调和级数，想看证明可以搜），因此埃式筛的筛的次数也就是时间复杂度为 $O(nlogn)$，而且是常数很小的 $nlogn$。代码如下：

bool vis[maxn];
void eratosthenes(int n){vis[1]=true;for(int i=2;i<=n;i++){if(!vis[i])for(int j=2*i;j<=n;j+=i)vis[j]=true;}
}

欧拉筛

上面发现一个合数可能会被筛去多次，有没有办法使得每个数只被筛去一次？有。

我们想要一个合数只被它最小的质因数与和它对应的次大的因数（最大因数就是它本身）筛去。当我们看到一个数的时候，枚举目前已经找到的质数，质数与我们看到的这个数相乘得到的数就是以枚举的质数为最小质因子的数。当我们枚举的质数已经可以整除我们现在看到的这个数的时候，就说明后面的质数不再小于等于我们看到的这个数的最小质因数了，就可以退出了。

通过这种想法，任何一个合数都一定会被它最小的质因数与和它对应的次大的因数筛去一次。所以时间复杂度是 $O(n)$ 的。代码如下。

int prime[maxn],cnt;
bool vis[maxn];
void Eular(int n){for(int i=2;i<=n;i++){if(!vis[i])prime[++cnt]=i;for(int j=1,p=prime[1];j<=cnt && i*p<=n;p=prime[++j]){vis[i*p]=true;if(i%p==0)break;}}
}

这个筛的过程一定要融会贯通，因为它每个筛到的数都一定是这个数最小的质因数来筛，这个性质很强，所以对这个东西略微修改就可以同时处理很多有用的信息。比如每个数的最小质因数，每个数所有的质因数，计算欧拉函数，莫比乌斯函数等等。

---

快速幂

普通算 $a^b$ 的方式就是把 $b$ 个 $a$ 乘起来，时间复杂度是 $O(b)$ 的。

发现一个东西，就是 $a\ast a=a^2,a^2\ast a^2=a^4,a^4\ast a^4=a^8,...$ 通过这种方式可以快速算出指数很大的 $a$ 的幂，然后我们再用预先算出来的这些幂指数去凑出答案的幂指数就行了。

举个例子，如果我们要算 $a^{10}$ ，那么我们可以先算出来 $a^1,a^2,a^4,a^8$，然后 $a^{10}=a^2\ast a^8$。当指数很大的时候，优势就更大了，比如 $a^{1025}=a^{1024}\ast a^1$，而处理到 $a^{1024}$ 只需要十次。

其实这东西有点像二进制，我们处理出来的是 $a^{1_{(2)}},a^{10_{(2)}},a^{100_{(2)}},a^{1000_{(2)}}$ ，相乘其实就是幂指数相加，通过幂指数相加得到任何一个我们想要的幂指数。比如 $a^{10}=a^{1010_{(2)}}=a^{1000_{(2)}+10_{(2)}}=a^{1000_{(2)}}\ast a^{10_{(2)}}$。想要凑出指数为 $b$，二进制位最多 $lo{g}_{2}b$ 位，所以快速幂的时间复杂度为 $O(lo{g}_{2}b)$

代码如下：

ll qpow(ll a,ll b,ll mod){ll base=a%mod,ans=1;while(b){if(b&1){ans*=base;ans%=mod;}base*=base;base%=mod;b>>=1;}return ans;
}

---

逆元

这个东西主要是为了解决模意义下不能相除的问题。那么我们怎么才能实现相除呢？

假设 $b$ 在取模意义下相当于 $\frac{1}{a}$ 或者说 $a^{-1}$，那么有 $a\ast b\equiv 1(\mathrm{mod}p)$。

费马小定理：

内容：若 $p$ 是质数，且 $a$ 不是 $p$ 的倍数，则有：$a^{p-1}\equiv 1(\mathrm{mod}p)$

它是欧拉定理的一个特例，不过在大部分情况下已经够用。这里要特别注意限定条件必须满足。因为 $a^{p-1}\equiv 1(\mathrm{mod}p)$，所以 $a^{p-2}\ast a\equiv 1(\mathrm{mod}p)$，这里 $a^{p-2}$ 就实现了上面的 $b$ 的功能，因此 $a^{p-2}\equiv \frac{1}{a}(\mathrm{mod}p)$。

简单来说，如果模数是质数 $p$，要除以 $a$，那么就直接乘以 $a^{p-2}$ 就可以实现相同的作用了。

---

A - 最大公约数和最小公倍数问题

思路：

如果把 $x_0,y_0$ 质因数分解，就得到了 $x_0=p_1^{c_1}\ast p_2^{c_2}\ast p_3^{c_3}\ast \cdots \ast p_i^{c_i},y_0=p_1^{t_1}\ast p_2^{t_2}\ast p_3^{t_3}\ast \cdots \ast p_j^{t_j}$，而且 $t_i\ge c_i$。对 $P,Q$ 来说，它们肯定都有 $x_0$ 的质因数部分，而从 $x_0$ 到 $y_0$ 的这一块质因数部分完全不重合。也就是说从 $\frac{y_0}{x_0}=p_1^{m_1}\ast p_2^{m_2}\ast p_3^{m_3}\ast \cdots \ast p_j^{m_j}$ ，取出若干种质因数$p_k^{m_k}$ 分给 $P$，剩下的分给 $Q$ 就是一个可行的解。（每种质因数只能分给 $P,Q$ 其中一个，否则分给 $P,Q$ 的部分会出现重合）

因为一种质因数要么分给 $P$，要么分给 $Q$，每种质因数有两种选择，如果 $\frac{y_0}{x_0}$ 有 $x$ 种质因数，那么分法也就是答案就有 $2^x$ 种。

code：

#include <iostream>
#include <cstdio>
using namespace std;int x,y;int main(){cin>>x>>y;if(y%x){cout<<0;return 0;}x=y/x;int cnt=0;for(int i=2;i*i<=x;i++)if(x%i==0){cnt++;while(x%i==0)x/=i;}if(x>1)cnt++;cout<<(1ll<<cnt);return 0;
}

---

B - 二元一次不定方程 (exgcd)

思路：

求解形如 $a\ast x+b\ast y=c$ 方程的通解已经在上面详细描述过了，这里尝试算出正整数解的个数和 $x,y$ 的最大和最小正整数解。

假设我们已经算到了 $\overline{x}=x+k\ast t1,\overline{y}=y-k\ast t2$。$x,y$ 的最小正整数解好算，先看 $x$ ，因为是正整数解，不是非负整数解，不包括0，解的值域应该是 $[1,t1]$，所以我们的计算式应该是 $x_{min}=((x-1)\%t1+t1)\%t1+1$，同理 $y_{min}=((y-1)\%t2+t2)\%t2+1$。

而 $x$ 取到最大正整数解时，$y$ 此时应该正好取到最小正整数解，因此我们把 $y_{min}$ 代入，算出此时的 $x$ 即可。同理 $y$。

正整数解的个数其实就看一个变量就可以了。比如看 $x$ ，其实就是 $x_{min}$ 和 $x_{max}$ 中间差了几个 $t1$ 。

code：

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;ll T,a,b,c;ll exgcd(ll a,ll b,ll &x,ll &y){if(!b){x=1;y=0;return a;}ll d=exgcd(b,a%b,x,y);ll z=x;x=y;y=z-(a/b)*y;return d;
}int main(){cin>>T;while(T--){scanf("%lld%lld%lld",&a,&b,&c);ll x,y,d=exgcd(a,b,x,y);if(c%d){printf("-1\n");continue;}x=c/d*x;y=c/d*y;ll t1=b/d,t2=a/d;ll xmi=((x-1)%t1+t1)%t1+1,ymi=((y-1)%t2+t2)%t2+1,xmx,ymx;xmx=(c-ymi*b)/a;ymx=(c-xmi*a)/b;if(xmx<=0 || ymx<=0){printf("%lld %lld\n",xmi,ymi);}else {printf("%lld %lld %lld %lld %lld\n",(xmx-xmi)/t1+1,xmi,ymi,xmx,ymx);}}return 0;
}

---

C - 质因数分解

思路：

枚举所有可能的因数，得到的第一个因数对应的另一个因数就是答案。

根据上面提到的性质，如果一个数 $x$ 是合数，那么它一定存在一个质因数 $p$，而且 $2\le p\le \sqrt{x}$。我们最多找到 $\sqrt{x}$ 就能找到答案了。

code：

#include <iostream>
#include <cstdio>
using namespace std;int n;int main(){cin>>n;for(int i=2;i*i<=n;i++)if(n%i==0){cout<<n/i;return 0;}return 0;
}

---

D - 质数筛

思路：

写作素数筛，读作开方暴力，直接暴力验证每个数是不是质数就行了。

code：

#include <iostream>
#include <cstdio>
using namespace std;
const int maxn=105;int n,a[maxn];bool isp(int x){if(x==1)return false;for(int i=2;i*i<=x;i++)if(x%i==0)return false;return true;
}int main(){cin>>n;for(int i=1;i<=n;i++)cin>>a[i];for(int i=1;i<=n;i++)if(isp(a[i]))cout<<a[i]<<" ";return 0;
}

---

E - Prime Distance

题意：

给你两个数 $L,U$ ，$1\le L<U\le 2147483647$ 且 $U-L\le 1000000$ 。问 $L\sim U$ 中相邻的两个质数中差值最小和最大的一对质数是什么。

思路：

直接跑筛法是跑不动的，但是根据一个数 $x$ 如果是合数，那么一定存在一个小于等于 $\sqrt{x}$ 的质因数的性质可以知道， $L\sim U$ 中的所有合数最大的最小质因数不会超过 $\sqrt{U}$ ，我们直接预处理处理 $\sqrt{U}$ 以内的所有素数，然后对每个素数筛掉 $L\sim U$ 中它的倍数，$L\sim U$ 中剩下的没有被筛掉的数就是素数。

需要注意当 $L=1$ 时，会有一个 $1$ 不会被筛掉，而它既不是素数也不是合数，因此需要特判去除。

code：

#include <iostream>
#include <cstdio>
#include <vector>
using namespace std;
const int maxn=1e5+5;
typedef long long ll;
const ll inf=1e18;ll L,U;int prime[maxn],cnt;
bool tvis[maxn*50];
void Eular(int n){tvis[1]=true;for(int i=2;cnt<=n;i++){if(!tvis[i])prime[++cnt]=i;for(int j=1,p=prime[1];j<=cnt && i*p<maxn*50;p=prime[++j]){tvis[i*p]=true;if(i%p==0)break;}}
}int main(){Eular(1e5);while(cin>>L>>U){vector<bool> vis(U-L+5);if(L==1)vis[0]=true;for(ll i=1,p=prime[1];i<=cnt && p*p<=U;p=prime[++i])for(ll j=max((L+p-1)/p,2ll);j*p<=U;j++)vis[j*p-L]=true;ll dmi=inf,dmx=0,l1,r1,l2,r2;for(ll i=0,lst=-1;i<=U-L;i++){if(!vis[i]){if(~lst){if(i-lst<dmi){dmi=i-lst;l1=lst;r1=i;}if(i-lst>dmx){dmx=i-lst;l2=lst;r2=i;}lst=i;}else lst=i;}}if(dmx)printf("%lld,%lld are closest, %lld,%lld are most distant.\n",L+l1,L+r1,L+l2,L+r2);else printf("There are no adjacent primes.\n");}return 0;
}

---

F - 线性筛素数

思路：

素数筛板子，$prime$ 数组中记录的就是从 $1\sim n$ 的所有素数，$prime[k]$ 就是第 $k$ 大的素数。

code：

#include <iostream>
#include <cstdio>
using namespace std;
const int maxn=1e7+5;int n,q;int prime[maxn],cnt;
bool vis[maxn*10];
void Eular(int n){for(int i=2;i<=n;i++){if(!vis[i])prime[++cnt]=i;for(int j=1,p=prime[1];j<=cnt && i*p<=n;p=prime[++j]){vis[i*p]=true;if(i%p==0)break;}}
}int main(){cin>>n>>q;Eular(n);for(int i=1,k;i<=q;i++){cin>>k;cout<<prime[k]<<endl;}return 0;
}

---

G - 又是毕业季II

思路：

好题。洛谷 P1414 又是毕业季II看不懂可以去题解区翻题解看

取出 $k$ 个数的最大公因数，我们算一遍 $看$ 个人的最大公因数是 $k\ast log$ 次数的，再枚举选择数…不用想了，复杂度爆炸。

那怎么得到 $k$ 个人的最大公因数。如果我们知道前 $k-1$ 个人的最大公因数，如果第 $k$ 个人有因数是这个最大公因数，那么最大公因数不变，否则是更小的需要所有人都有的因数。

嗯？如果我们一开始就对每个数标记一下以它为因数的人有几个不就好了，假设记录数组叫 $vis$。这样处理一下的话，如果我们找 $k$ 个人的最大公因数，就直接在 $1\sim inf=10^6$ 中找最大的 $i$，满足 $vis[i]\ge k$ 即可。

code：

#include <iostream>
#include <cstdio>
using namespace std;
const int maxn=1e4+5;
const int maxf=1e6+5;int n,vis[maxf];//是i的倍数的有多少个 int ans[maxn];int main(){cin>>n;for(int i=1,t;i<=n;i++){cin>>t;for(int p=1;p*p<=t;p++)if(t%p==0){vis[p]++;if(p*p!=t)vis[t/p]++;}}for(int i=1e6;i>=1;i--)if(!ans[vis[i]])ans[vis[i]]=i;for(int i=n,mx;i>=1;i--){mx=max(mx,ans[i]);ans[i]=mx;}for(int i=1;i<=n;i++)cout<<ans[i]<<endl;return 0;
}

---

H - 快速幂

思路：

快速幂板题

code：

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;ll a,b,p;ll qpow(ll a,ll b,ll mod){ll base=a%mod,ans=1;while(b){if(b&1){ans*=base;ans%=mod;}base*=base;base%=mod;b>>=1;}return ans;
}int main(){cin>>a>>b>>p;printf("%lld^%lld mod %lld=%lld\n",a,b,p,qpow(a,b,p));return 0;
}

---

I - Carmichael Numbers

题意：

Alvaro 认为，与传统海鲜饭相比，加密海鲜饭的主要优势在于，除了你自己，没人知道你吃的是什么。

一个合数如果能够通过费马素性检验：$$a^{n}modn=a$$就认为他是个一个 $Carmichaelnumber$。

现在给你一个数，问你它是不是 $Carmichaelnumber$

思路：

可以先把合数筛出来。如果给的是合数，再进行一下费马素性检验，如果都能通过就输出 $Carmichaelnumber$。

code：

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=65005;int n;int prime[maxn],cnt;
bool vis[maxn];
void Eular(int n){vis[1]=true;for(int i=2;i<=n;i++){if(!vis[i])prime[++cnt]=i;for(int j=1,p=prime[1];j<=cnt && i*p<=n;p=prime[++j]){vis[i*p]=true;if(i%p==0)break;}}
}ll qpow(ll a,ll b,ll mod){ll base=a%mod,ans=1;while(b){if(b&1){ans*=base;ans%=mod;}base*=base;base%=mod;b>>=1;}return ans;
}int main(){Eular(65000);while(cin>>n,n){if(!vis[n]){printf("%d is normal.\n",n);continue;}bool f=true;for(int i=2;i<n;i++){if(qpow(i,n,n)!=i){f=false;//不满足素性检验 break;}}if(f)printf("The number %d is a Carmichael number.\n",n);else printf("%d is normal.\n",n);}return 0;
}

---

J - 因子和

思路：

洛谷 P1593 因子和

首先需要知道一个数 $a$ 的所有因数和是多少。把 $a$ 质因数分解，得到 $$a=p_1^{c_1}\ast p_2^{c_2}\ast p_3^{c_3}\ast \cdots \ast p_m^{c_m}$$ 那么它的所有因数的和就是 $$\begin{gathered} sum=(1+p_1+p_1^2+\cdots +p_1^{c_1})\cdot (1+p_2+p_2^2+\cdots +p_2^{c_2})\cdot \\ \cdots \cdot (1+p_m+p_m^2+\cdots +p_m^{c_m}) \end{gathered}$$

这个式子可以这样理解：我们从第一个括号里随便取一项，和第二个括号里随便取一项，以此类推，一直取到最后一个括号，这样就可以得到一个因数。因数一共有 ${\prod }_{i=1}^m(c_i+1)$ 项，也正好是这些选取方案，取法和因数一一对应。

那么 $a^b$ 的因数和就是：$$\begin{gathered} sum=(1+p_1+p_1^2+\cdots +p_1^{c_1}+\cdots +p_1^{b\ast c_1})\cdot (1+p_2+p_2^2+\cdots +p_2^{b\ast c_2})\cdot \\ \cdots \cdot (1+p_m+p_m^2+\cdots +p_m^{b\ast c_m}) \end{gathered}$$

第一个括号内是一个等比数列，首项为1，公比为 $p_1$。因此第一个括号内的内容的总和 $=\frac{1\ast (1-p_1^{b\ast c_1+1})}{1-p_1}=\frac{p_1^{b\ast c_1+1}-1}{p_1-1}$。对 $p_1^{b\ast c_1+1}$ 求快速幂，对 $p_1-1$ 算逆元就能快速求解了，对一个数，处理出它所有的质因数和对应个数，每个质因数用这个过程算出答案，结果累乘起来即可。

不过 $p_1-1$ 可以计算逆元的前提条件是 $p_1-1\not\equiv 0(\mathrm{mod}9901)$，因此当 $p_1\equiv 1(\mathrm{mod}9901)$ 时的情况单独处理，发现原式子中这个数就等于 $b\ast c_1+1$ 。

注意 $a=0$ 的情况需要特判。

code：

#include <iostream>
#include <cstdio>
#include <vector>
using namespace std;
typedef long long ll;
const ll mod=9901;
#define pii pair<int,int>ll a,b;
vector<pii> t;ll qpow(ll a,ll b){b%=mod-1;ll base=a%mod,ans=1;while(b){if(b&1){ans*=base;ans%=mod;}base*=base;base%=mod;b>>=1;}return ans;
}
ll inv(ll x){return qpow(x,mod-2);}int main(){cin>>a>>b;if(a==0){cout<<0<<endl;return 0;}for(int i=2,cnt;i*i<=a;i++){if(a%i==0){cnt=0;while(a%i==0){a/=i;cnt++;}t.push_back(pii(i,cnt));}}if(a>1)t.push_back(pii(a,1));ll ans=1;for(auto x:t){ll p=x.first,c=x.second;if(p%mod!=1)ans=ans*(qpow(p%mod,b*c+1)-1)%mod*inv(p-1)%mod;else ans=ans*(b*c+1)%mod;}cout<<(ans%mod+mod)%mod;return 0;
}