初等数论笔记-算法竞赛相关-Ⅱ

Miller Rabin

$O(logn)$的大质数检验，用快速乘mul的话是$O(log^2n)$

引理1：设$p$是一个素数，$a$是一个正整数且$p\nmid a$，则$a^{p-1}\equiv 1(mod \;p)$。

引理2：如果$p$是一个素数，且$0<x<p$，则方程$x^2\equiv 1(mod\;p)$​的解为$x=1,p-1$

算法流程：

1.若为2，则true

2.若为非2的偶数或1，则false

3.若x不满足$a^{x-1}\equiv 1(mod \;x)$，则​不为质数

4.75%概率判断一个数x是否为质数：//ps为什么是75%，不知道，玄学~~

令$x-1=u·2^t$​​，其中$u$​​​为奇数

则$a^{x-1}\equiv 1(mod \;x)\implies (a^u)^{2^t}\equiv 1(mod \;x)$​

我们令$b=(a^u)mod \;x$​​​，然后进行t次循环，每次循环都让$b=(b·b)mod\;x$​​，这样t次循环后$b=a^{u·2^t}=a^{x-1}$​了

此时假如x是质数，由引理2可知b=1或b=x-1，否则x不是质数

t次循环结束后，由引理1得，若b不为1，则x不是素数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

inline ll mul(ll x,ll y, ll p){
	return (x*y-(ll)(x/(long double)p*y+1e-3)*p+p)%p;
}
bool Miller(ll x){
	//const static ll a[8]={2,3,7,11,233,331,24251,9875321};
	//const static int a[10]={2,3,7,11,233,331,39,61,24251,9875321} -1e14
	//const static int a[12]={2,3,5,7,11,233,331,19,23,29,31,37}; -1e9
	if(x==2) return true;
	if(!(x&1)||x<2) return false;
	ll u=x-1;
	while(!(u&1)){ //如果为偶数，就一直除以2，直到u为奇数
		u>>=1;
	}
	for(int i=0;i<8;++i){
		if(x==a[i]) return true;
		//	ll a=rand()%(x-1)+1;//也可以rand 
		ll t=u,m=qpow(a[i],u,x);
		while(t!=x-1&&m!=1&&m!=x-1){
			m=mul(m,m,x);//long long的平方 
			m=m*m%x;
			t<<=1;
		}
		if(m!=x-1&&!(t&1)) return false;
	}
	return true;
}

Pollard’s Rho

$O(n^{\frac{1}{4}})$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

ll PollardRho(ll n){
	if(n==4) return 2;
	if(Miller(n))
	return n;
	while (1){
		ll c=rand()%(n-1)+1;
		auto f=[=](ll x) { return mul(x,x,n)+c%n; };
		ll t=0,r=0,p=1,q;
		do{
			for(int i=0;i<128;++i){
				t=f(t),r=f(f(r));
				q=mul(p,abs(t-r),n);
				if(!q)
					break;
				p=q;
			}
			ll d=__gcd(p,n);
			if(d>1)
			return d;
		} while(t!=r);
	}
}

Euclideanoid

对式子

$f(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor$

在低于线性的复杂度下【$O(logn)$​​】进行计算。由于其复杂度与欧几里得算法相当，故得名类欧几里得算法。由于该算法名字在目前暂时没有英文翻译（可能有，但我没查到），在请教过一位大佬后冒昧将其翻译成Euclidean-oid。下面是算法的推导过程，如有笔误请指正。

当$a≥c$或者$b≥c$时，

$\begin{aligned} f(a,\;b,\;c,\;n)&=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor\\ &=\sum_{i=0}^n\lfloor\frac{(\lfloor \frac{a}{c}\rfloor c+a\%c) i+ (\lfloor \frac{b}{c}\rfloor c+b \% c) }{c} \rfloor\\ &=\frac{n(n+1)}{2} \lfloor\frac{a}{c}\rfloor+(n+1)\lfloor\frac{b}{c}\rfloor+f(a\%c,\;b\%c,\;c,\;n) \end{aligned}$

对于$a<c$且$b<c$的情况，

$\begin{aligned} f(a,\;b,\;c,\;n)&=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor\\ &=\sum_{i=0}^n \sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}1\\ &=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\sum_{i=0}^n \;[j<\lfloor\frac{ai+b}{c}\rfloor]\\ \end{aligned}$

考虑将$i$​提出来，对于 $j<\lfloor\frac{ai+b}{c}\rfloor\iff j+1≤\frac{ai+b}{c} \iff jc+c-b≤ai \iff i>\lfloor \frac{jc+c-b-1}{a}\rfloor $​

$\begin{aligned} f(a,\;b,\;c,\;n)&=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\sum_{i=0}^n \;[j<\lfloor\frac{ai+b}{c}\rfloor]\\ &=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor-1}\sum_{i=0}^n \;[i>\lfloor \frac{jc+c-b-1}{a}\rfloor]\\ &=\sum_{j=0}^{\lfloor\frac{an+b}{c}\rfloor-1} n-\lfloor \frac{jc+c-b-1}{a}\rfloor\\ &=nm-f(c,c-b-1,a,m-1) \end{aligned}$

其中，$m=\lfloor\frac{an+b}{c}\rfloor$

exEuclideanoid

在类欧几里得算法的基础上进行扩展，对$g$和$h$进行计算。

$\begin{aligned} g(a,b,c,n)&=\sum_{i=0}^ni\lfloor\frac{ai+b}{c}\rfloor\\ when \;a≥c||b≥c\;&=g(a\%c,\;b\%c,\;c,\;n)+\lfloor\frac{a}{c}\rfloor\frac{n(n+1)(2n+1)}{6}+\lfloor\frac{b}{c}\rfloor \frac{n(n+1)}{2}\\ else \;&=\frac{1}{2}\;[nm(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)] \end{aligned}$

其中，$m=\lfloor\frac{an+b}{c}\rfloor$

$\begin{aligned} h(a,b,c,n)&=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^2\\ when \;a≥c||b≥c,\;&=h(a\%c,b\%c,c,n)+2\lfloor\frac{b}{c}\rfloor f(a\%c,b\%c,c,n)+2\lfloor\frac{a}{c}\rfloor g(a\%c,b\%c,c,n)+\\ &\lfloor\frac{a}{c}\rfloor^2\frac{n(n+1)(2n+1)}{6}+ \lfloor\frac{b}{c}\rfloor^2 (n+1)+\lfloor\frac{a}{c}\rfloor\lfloor\frac{b}{c}\rfloor n(n+1)\\ else,\;&=nm(m+1)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n) \end{aligned}$

其中，$m=\lfloor\frac{an+b}{c}\rfloor$

P5170 【模板】类欧几里得算法 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

const int mod=998244353;
void pre(int n=1e7){
	return ;
}
ll exgcd(ll a,ll b,ll &x,ll &y){
	ll d=a;
	if(b){
		d=exgcd(b,a%b,y,x);
		y-=(a/b)*x;
	} else{
		x=1;y=0;
	}
	return d;
}
ll inv(ll a,ll mod){
	ll x,y;
	if(exgcd(a,mod,x,y)!=1)
		return -1;
	return (x%mod+mod)%mod;
}
const ll inv6=inv(6,mod),inv2=inv(2,mod);
struct node{
	int f,h,g;
	node(){f=g=h=0;}
};
node calc(ll a,ll b,ll c,ll n){
	node d,e;
	ll n1=n*(n+1)%mod*inv2%mod,n2=n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
	ll m=(a*n+b)/c;
	ll k1=a/c,k2=b/c;
	if(a==0){
		d.f=(n+1)*k2%mod;
		d.g=n1*k2%mod;
		d.h=(n+1)*k2%mod*k2%mod;
		return d;
	}
	if(a>=c||b>=c){
		e=calc(a%c,b%c,c,n);
		d.f=(e.f+(n+1)*k2%mod+n1*k1%mod)%mod;
		d.g=(e.g+k1*n2%mod+k2*n1%mod)%mod;
		d.h=(e.h+2*k2*e.f%mod+2*k1*e.g%mod+k1*k1%mod*n2%mod+k2*k2%mod*(n+1)%mod+k1*k2%mod*n%mod*(n+1)%mod)%mod;
		return d;
	}
	e=calc(c,c-b-1,a,m-1);
	d.f=((n*m%mod-e.f)%mod+mod)%mod;
	d.g=((m*n1%mod-inv2*e.h-inv2*e.f)%mod+mod)%mod;
	d.h=((n*m%mod*(m+1)%mod-2*e.g-2*e.f-d.f)%mod+mod)%mod;
	return d;
}
void solves(){
	ll n,a,b,c;cin>>n>>a>>b>>c;
	node ans=calc(a,b,c,n);
	cout<<ans.f<<" "<<ans.h<<" "<<ans.g<<endl;
}
int main(){
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
	int _OwO_=1;
	pre();
	cin>>_OwO_;
	for(int _=1;_<=_OwO_;++_){
		solves();
	}
}