中高级数论 [欧拉函数线性筛，二次剩余]

欧拉函数线性筛

对于素数$p$,
$$\phi (p\ast i)=\{\begin{array}{ll}p-1& i=1\\ p\ast \phi (i)& p\mid i\\ (p-1)\ast \phi (i)& p\nmid i\end{array}$$

证明：

$i=1$ 显然

$i\ne 1$根据

$\phi (i)=i\prod \frac{p_i-1}{p_i}$

$p\mid i$时除掉$i$中的$p$

$p\nmid i$ 时除掉$p$

魔改欧拉筛素数即可

欧拉定理

当$(a,p)=1$

$a^{\phi (p)}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

即

$a^b\equiv a^{b\%\phi (p)}(\mathrm{m}\mathrm{o}\mathrm{d}p)$

扩展欧拉定理

当$b\ge \phi (p)$

$a^b\equiv a^{b\%\phi (p)+\phi (p)}(\mathrm{m}\mathrm{o}\mathrm{d}p)$

注意不要求互质

Miller Rabbin 算法

可以$O(lo{g}_x)$判断素数（有极小概率出错）

费马小定理：$p$为素数，$a\ne p$,$a^{p-1}\equiv 1(modp)$

我有一个对这个命题的十分美妙的证明,可惜这里空白太小,我写不下

二次探测定理：$p$为素数，$x^2\equiv 1(modp)$的解为$x=1$或$x=p-1$

证明见二次剩余

算法流程：设要判断的数为$x$

1. 特判$0,1,2$,偶数
2. 取一个较小素数$p$,设$a,b$满足$2^a\times b=x-1$
3. 设$x^{\prime }\equiv p^b(modx)$,不断取平方并对$x$取模，若结果为$1$，用二次探测判断。
4. 重复$a$次后$x^{\prime }\equiv p^{x-1}(modx)$，用费马小定理判断

例题：hdu2138 How many prime numbers

题意：给$N$个数，数数有多少个素数

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
typedef long long ll;
using namespace std;
int qpow(int a,int p,int m)
{
	int ans=1;
	while (p)
	{
		if (p&1) ans=((ll)ans*a)%m;
		a=((ll)a*a)%m,p>>=1;
	}
	return ans;
}
int pri[]={2,3,5,7,11,13,17,19,23,29};
bool check(int x)
{
	if (x==2) return true;
	if (!(x&1)) return false;
	int a=0,b=x-1;
	while (!(b&1)) ++a,b>>=1;
	for (int i=0;pri[i]<x&&pri[i]<29;i++)
	{
		int now=qpow(pri[i],b,x);
		for (int j=1;j<=a;j++)
		{
			int t=((ll)now*now)%x;
			if (t==1&&now!=1&&now!=x-1) return false;
			now=t;
		}
		if (now!=1) return false;
	}
	return true;
}
int main()
{
	int n;
	while (~scanf("%d",&n))
	{
		int ans=0,x;
		while (n--)
		{
			scanf("%d",&x);
			ans+=check(x);
		} 
		printf("%d\n",ans);
	}
	return 0;
}

关于底数选择：

以$2,3,7,61,24251$为底，$1e16$内唯一无法判断的合数为$46856248255981$

特判一下就可以了

代码不想改了

Pollard-Rho算法

这个算法就更玄学了，可以极快分解质因数。

有多快呢？$longlong$能跑几百个

该算法核心是找到这个数的一个非平凡因子（即不是$1$和本身，以下简称因子），然后递归就行了

如何找到一个因子呢？ 随机。

当然直接随机肯定不现实，需要一定的改进

①优秀的随机函数

定义随机函数$f_i=f_{i-1}^2+c$，其中$c$为随机常数。

至于为什么要选这个函数……鬼知道啊

②$gcd$

设给定的数为$x$，我们要找的是$d\mid x$

概率还是低了点

容易想到一个办法：取$gcd$，大于$1$就返回$gcd$的值

这样多了因数的倍数，很大的提高了概率

③倍增路径

每次取一段$2^k$个数，在模$x$乘起来，最后取$gcd$,$k$为当前段数

易知这样是等效的，避免了大量的$gcd$计算

使用了倍增，将$gcd$次数从$O(n)$降至$O(logn)$

另外不再使用随机值，而是用随机出来的值和上一段最后一个数作差，这样得到数的个数也大幅增加。

这样除非环特别短，否则环上两两的数都会统计入答案。

根据生日悖论，当环长为45时找不到答案的概率小于$1e-18$

而你的函数是根据$C$改变的，所以根本卡不了

所以如果出锅了，自觉去交一发Hash Killer III

④玄学优化

每一段跑$127$个数更新一次

亲测快了$3$倍

来自$luogu$题解区

原理就不得而知了

结合$Miller-Rabbin$食用更佳

例题：luogu模板

题意：给个数，是素数输出"Prime"，不是素数输出最大质因数

递归，记个全局变量。还可以剪枝。

运用你丰富的卡常技巧就可以通过此题

// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cstdlib>
#ifdef WIN32
#define lld "%I64d"
#else
#define lld "%lld"
#endif
typedef long long ll;
using namespace std;
inline ll read()
{
	ll ans=0;
	char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
	return ans;
}
ll base[]={2,3,7,61,24251};
inline ll qmul(const ll& a,const ll& b,const ll& m)
{
    ll ans=a*b-(ll)((long double)a*b/m+0.5)*m;
    return ans<0? ans+m:ans;
}
inline ll qpow(ll a,ll p,const ll& m)
{
    ll ans=1;
    while (p)
    {
        if (p&1) ans=qmul(ans,a,m);
        a=qmul(a,a,m),p>>=1;
    }
    return ans;
}
inline bool Miller_Rabbin(const ll& x)
{
    if (x==1) return false;
    if (x==2) return true;
    if (!(x&1)) return false;
    if (x==46856248255981ll) return false;
    int a=0;ll b=x-1;
    while (!(b&1)) ++a,b>>=1;
    for (register int i=0;i<5&&base[i]<x;i++)
    {
        ll now=qpow(base[i],b,x);
        for (register int j=1;j<=a;j++)
        {
            ll t=qmul(now,now,x);
            if (t==1&&now!=1&&now!=x-1) return false;
            now=t;
        }
        if (now!=1) return false;
    }
    return true;
}
ll ans;
inline ll f(const ll& x,const ll& c,const ll&m){return (qmul(x,x,m)+c)%m;}
inline ll gcd(ll x,ll y)
{
	static ll t;
	while (y>0) t=y,y=x%y,x=t;
	return x;
}
inline ll abso(const ll& x){return x>0? x:-x;}
inline ll Pollard_rho(const ll& x)
{
    ll s=0,t=0,c=rand()%(x-1)+1,d,len=1,m=1,i;
    for (;;len<<=1,s=t,m=1)
    {
        for (i=1;i<=len;i++)
        {
            t=f(t,c,x);
            m=qmul(m,abso(t-s),x);
            if (i%127==0&&(d=gcd(m,x))>1) return d;
        }
        if ((d=gcd(m,x))>1) return d;
    }
}
inline void fact(const ll& x)
{
    if (x<ans) return;
    if (Miller_Rabbin(x)) return (void)(ans=max(ans,x));
    ll d=Pollard_rho(x);
    if (Miller_Rabbin(d)) ans=max(ans,d);
    else fact(d);
    d=x/d;
    if (Miller_Rabbin(d)) ans=max(ans,d);
    else fact(d);	
}
int main()
{
    register int T;
    ll n,d;
    scanf("%d",&T);
    while (T--)
    {
        n=read();
        ans=0,fact(n),d=ans;
        if (n==d) printf("Prime\n");
        else printf(lld"\n",ans);
    }
    return 0;
}

原根

占个坑以后补

二次剩余

以下的$p$均为奇素数

定义：方程$x^2\equiv n(modp)$有解,称$n$是$modp$意义下的二次剩余，否则称非二次剩余

说人话：在$modp$下可以开根号

勒让德符号

$$(\frac{n}{p})=\{\begin{array}{ll}1& \text{n是p的二次剩余}\\ -1& \text{n是p的非二次剩余}\\ 0& p\mid n\end{array}$$

如无特别说明，下文分数+括号均为勒让德符号

引理

$(\frac{n}{p})\equiv n^{\frac{p-1}{2}}(modp)$

(强调：$p$是奇素数)

证明：

1. 当$n$是$p$的二次剩余，设$x^2\equiv n(modp)$,则$n^{\frac{p-1}{2}}\equiv x^{p-1}(modp)$。由费马小定理，$x^{p-1}\equiv 1(modp)$
2. 当$n$是$p$的非二次剩余，即不存在$x^2\equiv n(modp)$。对于任意的$a$,由于$a^2\ne n(modp)$,所以$a\ne a^{-1}n$。而逆元是互逆的，这样可以把$p-1$个数配成$\frac{p-1}{2}$对，每对数的积为$n$。将所有数乘起来，$(p-1)!\equiv n^{\frac{p-1}{2}}(modp)$。由威尔逊定理，$(p-1)!\equiv -1(modp)$
3. $p\mid n$,显然成立

引理

共有$\frac{p-1}{2}$个$n$是$p$下的二次剩余

证明：设两个数$u,v$满足$u\ne v,u^2\equiv v^2(modp)$,即$p\mid u^2-v^2=(u+v)(u-v)$

故$p\mid u+v$,反之同理

说人话：有且仅有模意义下的相反数平方相同。

所以$x^2$共有$\frac{p-1}{2}$个不同取值

所以随机出一个（非）二次剩余是$O(1)$的

引理

$(a+b{)}^p\equiv a^p+b^p(modp)$

二项式定理展开即可，不再证明。

重头戏开始。坐稳了，前方高能！

随机一数$a$，若$(\frac{a^2-n}{p})=-1$,令$w=\sqrt{a^2-n}$,则$(a+w{)}^{\frac{p+1}{2}}$是方程$x^2\equiv n(modp)$的根

什么？$a^2-n$不是非二次剩余吗？怎么开根？

没办法啊……那就类比虚数，定义$z=a+bw$吧（绝对不是人想出来的）

引理

$w^p\equiv -w(modp)$

证明：$w^p\equiv w^{p-1}w\equiv (a^2-n{)}^{\frac{p-1}{2}}w\equiv (\frac{a^2-n}{p})w\equiv -w（modp）$

最后证明：

$x^2\equiv (a+w{)}^{p+1}\equiv (a+w{)}^p(a+w)\equiv (a^p+w^p)(a+w)=(a-w)(a+w)=a^2-w^2\equiv n(modp)$

证毕……

哎等等？$x$是虚数怎么办？

假设$x$是虚数，设$x=a+bw$(不是上面的$a$) 由$x$定义，$a\ne 0$

$x^2=a^2+b^2{w}^2+2abw=a^2+a^2{b}^2-n{b}^2+2abw$是虚数，而$n$是实数。

假设不成立，所以$x$是实数(并不严谨)

例题：Timus OJ1132 Square Root

数据很小，暴力能过。但只有这题，将就了吧。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#include <cstdlib>
#define int long long
using namespace std;
struct complex{int x,y;};
int n,w2,p;
inline complex operator *(const complex& a,const complex& b){return (complex){(a.x*b.x%p+a.y*b.y%p*w2%p)%p,(a.x*b.y%p+a.y*b.x%p)%p};}
inline complex qpow(complex a,int b)
{
	complex ans;
	ans.x=1,ans.y=0;
	while (b)
	{
		if (b&1) ans=ans*a;
		a=a*a,b>>=1;
	}
	return ans;
}
inline int qpow(int a,int b)
{
	int ans=1;
	while (b)
	{
		if (b&1) ans=ans*a%p;
		a=a*a%p,b>>=1;
	}
	return ans;
}
inline int lerender(const int& x){return qpow(x,p>>1);}
int solve()
{
	if (lerender(n)==p-1) return -1;
	int a;
	while (true)
	{
		a=rand()%p;
		w2=(a*a-n)%p;
		if (w2<0) w2+=p;
		if (lerender(w2)==p-1) break;
	}
	complex res;
	res.x=a,res.y=1;
	res=qpow(res,(p+1)>>1);
	return res.x;
}
inline int read()
{
	int ans=0;
	char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=(ans<<3)+(ans<<1)+(c^48),c=getchar();
	return ans;
}
void write(const int& x)
{
	if (x<10) return(void)putchar(x^48);
	int t=x/10;
	write(t),putchar((x-(t<<1)-(t<<3))^48);
}
main()
{
	int T=read();
	while (T--)
	{
		n=read(),p=read();
		if (p==2) 
		{
			printf("1\n");
			continue;
		}
		int ans=solve();
		if (ans==-1) printf("No root\n");
		else
		{
			int ans2=p-ans;
			if (ans>ans2) swap(ans,ans2);
			write(ans),putchar(' '),write(ans2),putchar('\n');
		}
	}
	return 0;
}