数论

OI&ACM, 数学
  1. 1 【bzoj 2186】[Sdoi2008]沙拉公主的困惑
  2. 2 【COGS 405】[NOIP2009] Hankson的趣味题
  3. 3 【COGS 1265】[NOIP2012] 同余方程
  4. 4 【COGS 630】【bzoj 2749】[NOIP2011] 计算系数
  5. 5 【COGS 792】[HAOI2012] 外星人
  6. 6 【COGS 662】【bzoj 1041】[HAOI2008] 圆上的整点
  7. 7 [【BZOJ 2705】【COGS 931】 [SDOI2012]Longge的问题]
  8. 8 【COGS 2057】[ZLXOI2015]殉国

看了一天vfk的PPT

1 【bzoj 2186】[Sdoi2008]沙拉公主的困惑

1~N!中与M!互质的数

\(ans=\varphi(m!)\times \frac {n!}{m!} \mod R\)
因为所有小于m!且与m!互质的数加上m!的整数倍都与m!互质，而其他数都不与m!互质。
\(f(m)=\frac {\varphi(m!)}{m!}=\sum_{gcd(pi,m)=1}\frac {pi-1}{pi}\)
\(\therefore ans=f(m)\times n! \mod R\)
欧拉筛法预处理f()函数
tips：
1 . 如果m<R<=n的话ans=0，因为ans为R的倍数
2 . 求f()时不乘R，不然没法求R的逆元，因为gcd(R,R)!=-1
3 . 求i在mod R意义下的逆元：\(ni[i]=i^{R-2}\)，或者用exgcd

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
#include<iostream>
#include<cstdio>
using namespace std;
const int N=10000001;
int R,T,n[10001],m[10001],maxn,tot,prime[N],ni[N],fac[N],ans[N];//factorial 阶乘 
bool not_p[N];
int power(int a,int b)
{
	int re=1;
	for(;b;b>>=1,a=(long long)a*a%R)
	  if(b&1) re=(long long)re*a%R;
	return re;
}
void Pre()
{
	not_p[1]=1;fac[1]=1;ans[1]=1;
	for(int i=2;i<=maxn;i++)
	{
		fac[i]=(i!=R)?((long long)fac[i-1]*i%R):fac[i-1];
		if(!not_p[i])
		{
			prime[++tot]=i;ni[i]=power(i,R-2);
			ans[i]=(long long)ans[i-1]*(i-1)%R;
			if(i!=R) ans[i]=(long long)ans[i]*ni[i]%R;
		}
		else ans[i]=ans[i-1];
		for(int j=1,k;j<=tot&&prime[j]*i<=maxn;j++)
		{
			not_p[prime[j]*i]=1;
			if(!(i%prime[j])) break;
		}
	}
}
int main()
{
	scanf("%d%d",&T,&R);
	for(int i=1;i<=T;i++) scanf("%d%d",&n[i],&m[i]),maxn=max(maxn,n[i]);
	Pre();
	for(int i=1;i<=T;i++)
	{
		if(R>m[i]&&R<=n[i]) {printf("0\n");continue;}
		printf("%lld\n",(long long)fac[n[i]]*ans[m[i]]%R);
	}
}

2 【COGS 405】[NOIP2009] Hankson的趣味题

已知正整数a0,a1,b0,b1，求正整数x的个数，满足：
gcd(x,a0)=a1;lim(x,b0)=b1;

首先表示暴力枚举也能过。。。（COGS开O2就是快）

正解：分解质因数
设质因数p
x=p^a×...
a0=p^c1×...
a1=p^c2×...
b0=p^c3×...
b1=p^c4×...
则min(a,c1)=c2，max(a,c3)=c4；

若c1==c2,则a>=c1; --①
若c1!=c2,则a=c2; --②
若c3==c4,则a<=c3; --③
若c3!=c4,则a=c4; --④

当①&&③时，需满足c1<=a<=c3，此时ans×=c3-c1+1；（乘法原理）
当①&&④时，需满足a=c4>=c1；
当②&&③时，需满足a=c2<=c3；
当②&&④时，需满足a=c2=c4；
如果不满足，则ans=0。 

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
#include<cstdio>
using namespace std;
int T,ans,a0,a1,b0,b1;
int gcd(int x,int y)
{
	return y?gcd(y,x%y):x;
}
int main()
{
	freopen("son.in","r",stdin);
	freopen("son.out","w",stdout);
	scanf("%d",&T);
	while(T--)
	{
		ans=0;
		scanf("%d%d%d%d",&a0,&a1,&b0,&b1);
		for(int x=1;x*x<=b1;x++)
		  if(!(b1%x))
		  {
		  	if(gcd(x,a0)==a1&&(long long)gcd(x,b0)*b1==(long long)x*b0) ans++;
		  	if(x*x!=b1&&gcd(b1/x,a0)==a1&&gcd(b1/x,b0)==b0/x) ans++;
		  }  
		printf("%d\n",ans);
	}
}

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
#include<cstdio>
using namespace std;
const int N=50001;
int T,ans,a0,a1,b0,b1,tot,prime[N];
bool not_p[N];
void pre()
{
	not_p[1]=1;
	for(int i=2;i<N;i++)
	{
		if(!not_p[i]) prime[++tot]=i;
		for(int j=1;j<=tot&&prime[j]*i<N;j++)
		{
			not_p[prime[j]*i]=1;
			if(!(i%prime[j])) break;
		}
	}
}
bool calc(int x)
{
	int c1=0,c2=0,c3=0,c4=0;
	while(!(a0%x)) a0/=x,c1++;
	while(!(a1%x)) a1/=x,c2++;
	while(!(b0%x)) b0/=x,c3++;
	while(!(b1%x)) b1/=x,c4++;
	if(c1==c2&&c3==c4&&c1<=c3) ans*=c3-c1+1;
	else if((c1==c2&&((c3==c4&&c1>c3)||(c3!=c4&&c4<c1)))||
	        (c1!=c2&&((c3==c4&&c2>c3)||(c3!=c4&&c2!=c4))))  ans=0;
	return ans;
}
int main()
{
	freopen("son.in","r",stdin);
	freopen("son.out","w",stdout);
	scanf("%d",&T);pre();
	while(T--)
	{
		ans=1;
		scanf("%d%d%d%d",&a0,&a1,&b0,&b1);
		for(int i=1,x=prime[i];x*x<=b1;i++,x=prime[i])
		if(!(b1%x))
		  if(!calc(x)) 
		    break;
		if(b1!=1) calc(b1);
		printf("%d\n",ans);
	}
}

3 【COGS 1265】[NOIP2012] 同余方程

求关于 x 的同余方程 ax ≡ 1 (mod b)的最小正整数解。

ax+by=1;
解exgcd，求最小的x

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include<cstdio>
#include<cmath>
using namespace std;
int a,b,x,y;
void exgcd(int a,int b,int &x,int &y)
{
	if(!b) { x=1; y=0; }
	else { exgcd(b,a%b,y,x); y-=x*(a/b); }
}
int main()
{
	freopen("mod.in","r",stdin);
	freopen("mod.out","w",stdout);
	scanf("%d%d",&a,&b);
	exgcd(a,b,x,y);
	x%=b; while(x<=0) x+=b;
	printf("%d",x);
}

4 【COGS 630】【bzoj 2749】[NOIP2011] 计算系数

给定一个多项式 \((ax+by)^k\) ，请求出多项式展开后 \(x^ny^m\) 项的系数mod 10007的结果。

根据多项式定理，\((ax+by)^n=\sum_{r=0}^n C_n^r (ax)^{n-r}(by)^r\)
\(\therefore ans=C_k^n\times a^n\times b^m \mod 10007\)
而\(C_k^n=\frac{k!}{n!(k-n)!}=\frac{\prod_{i=1}^{i<=n} k-n+i}{n!}\)

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
#include<cstdio>
using namespace std;
const int mod=10007;
int a,b,k,n,m;
int power(int x,int y)
{
	int re=1;
	for(;y;y>>=1,x=x*x%mod)
	  if(y&1) re=re*x%mod;
	return re;
}
int C(int x,int y)
{
	int re1=1,re2=1;
	if(y-x<x) x=y-x;
	for(int i=1;i<=x;i++)
	{
		re1=re1*(y-x+i)%mod;
		re2=re2*i%mod;
	}
	return re1*power(re2,mod-2)%mod;
}
int main()
{
	freopen("factor.in","r",stdin);
	freopen("factor.out","w",stdout);
	scanf("%d%d%d%d%d",&a,&b,&k,&n,&m);
	printf("%d\n",C(n,k)*power(a%mod,n)%mod*power(b%mod,m)%mod);
}

5 【COGS 792】[HAOI2012] 外星人

求最小的x，使\(\varphi^x(N)=1\)，即对N做\(\varphi(N)\)操作多少次后为1

由题目给的公式，对于素数x>2，phi(x)会变出好多2，而phi(2) = 1；所以求出每一个质数经过分解会变出多少个2来，2的个数就是ans。
令f[i]表示i分解出了几个2，用筛法：
i为质数，f[i]=f[i-1]；
否则，f[i×prime[j]=f[i]+f[prime[j]]；
tip：如果一开始n为奇数，则ans需+1，因为变出2需要第一步。

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
#include<cstdio>
using namespace std;
const int N=100001;
int prime[N],tot,f[N],T,m,x,y;
long long ans;
bool notp[N],yes;
void pre()
{
	notp[1]=1;f[1]=1;
	for(int i=2;i<N;i++)
	{
		if(!notp[i]) prime[++tot]=i,f[i]=f[i-1];
		for(int j=1;j<=tot&&prime[j]*i<N;j++)
		{
			notp[prime[j]*i]=1;
			f[prime[j]*i]=f[prime[j]]+f[i];
			if(!(i%prime[j])) break;
		}
	}
}
int main()
{
	freopen("alien.in","r",stdin);
	freopen("alien.out","w",stdout);
    pre();
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d",&m);yes=1;ans=0;
		while(m--)
		{
			scanf("%d%d",&x,&y);ans+=(long long)f[x]*y;
			if(x==2) yes=0;
		}  
		printf("%lld\n",ans+yes);
	}
}

6 【COGS 662】【bzoj 1041】[HAOI2008] 圆上的整点

平面上有一个圆,圆心坐标为(0,0)，半径为n.问圆周上有多少个x,y坐标均为整数的点。

法一：见hzwer 

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
#include<cstdio>
#include<cmath>
using namespace std;
long long r,ans;
long long gcd(long long x,long long y)
{
	return y?gcd(y,x%y):x;
}
int main()
{
	freopen("cir.in","r",stdin);
	freopen("cir.out","w",stdout);
	scanf("%lld",&r);
	for(long long d=1;d*d<=2*r;d++)
	  if(!((2*r)%d))
	  {
	  	for(long long a=1;a*a<=r/d;a++)
	  	{
	  		long long b=sqrt(2*r/d-a*a);
			if(a!=b&&a*a+b*b==2*r/d&&gcd(a*a,b*b)==1) 
			  ans++; 
	  	}
	  	if(d*d!=2*r)
	  	for(long long a=1;a*a<=d/2;a++)
	  	{
	  		long long b=sqrt(d-a*a);
	  		if(a!=b&&a*a+b*b==d&&gcd(a*a,b*b)==1) 
			  ans++;
	  	}
	  }
	printf("%lld\n",ans*4+4);
}

法二： 勾股数的性质（wiki）
勾股数：\(a^2+b^2=c^2\)
素勾股数：gcd(a,b,c)=1
设m>n，gcd(m,n)=1，m、n中有一个偶数
则\(a=m^2-n^2\)
\(b=2mn\)
\(c=m^2+n^2\)
所以枚举半径n的因数i，求c=i的素勾股数，即x<y且\(x^2+y^2=i^2\)且gcd(x,y)=1且x+y为奇数的个数
ans×8+4，因为枚举的x<y所以要×8，还有坐标轴上的点所以+4

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
#include<cstdio>
#include<cmath>
using namespace std;
int n;long long ans;
int gcd(int x,int y)
{
	return y?gcd(y,x%y):x;
}
void calc(int z)
{
	for(int x=1,y;2*x*x<z;x++)
	{
		y=sqrt(z-x*x);
		if((long long)x*x+y*y==z&&gcd(x,y)==1&&(x+y)&1) ans++;
	}
}
int main()
{
	freopen("cir.in","r",stdin);
	freopen("cir.out","w",stdout);
	scanf("%d",&n);
	for(int i=1;i*i<=n;i++)
	if(!(n%i))
	  {
		calc(i);
		if(i*i!=n) calc(n/i);
	  }
	printf("%lld\n",ans*8+4);
}

7 [【BZOJ 2705】【COGS 931】 [SDOI2012]Longge的问题]

求\(\sum_{1<=i<=n} gcd(i,n)\)

枚举n的约数k，令s(k)为满足gcd(m,n)=k,(1<=m<=n)的m的个数，则\(ans=\sum_{k\mid n}k\times s(k)\)
因为gcd(m,n)=k，所以\(gcd(\frac mk,\frac nk)=1\)，于是\(s(k)=\varphi(,\frac nk)\)
phi可以在根号的时间内求出

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
#include<cstdio>
using namespace std;
long long n,ans;
long long phi(long long x)
{
	long long y=x;
	for(long long i=2;i*i<=x;i++)
	  if(!(x%i))
	  {
	  	y=y/i*(i-1);
	  	while(!(x%i)) x/=i;
	  }
	if(x>1) y=y/x*(x-1);
	return y;
}
int main()
{
	freopen("gcdsum.in","r",stdin);
	freopen("gcdsum.out","w",stdout);
	scanf("%lld",&n);
	for(long long i=1;i*i<=n;i++)
	  if(!(n%i))
	  {
		ans+=i*phi(n/i);
		if(i*i!=n) ans+=n/i*phi(i);
	  }
	printf("%lld\n",ans);
}

8 【COGS 2057】[ZLXOI2015]殉国

给定a,b,c,求ax+by=c(x,y>=0)的最小最大的x+y，以及解的组数

求exgcd，由\(xx=x\times \frac cd+\frac bd \times t \ge 0,yy=y\times \frac cd-\frac ad \times t\ge 0\)
得到t的区间[l,r]
方法数=r-l+1;
因为方程线性，所以最值在l、r取到
要写long double!! 

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
#include<cstdio>
#include<cmath>
#define min(x,y) x<y?x:y
#define max(x,y) x>y?x:y
using namespace std;
long long a,b,c,d,x,y,xx,yy,ans1,ans2,ans;
void exgcd(long long a,long long b,long long &d,long long &x,long long &y)
{
	if(!b)  d=a,x=1,y=0;
	else exgcd(b,a%b,d,y,x),y-=x*(a/b);
}
int main()
{
	freopen("BlackHawk.in","r",stdin);
	freopen("BlackHawk.out","w",stdout);
	scanf("%lld%lld%lld",&a,&b,&c);
	exgcd(a,b,d,x,y);
	if(!(c%d))
	{
		xx=ceil((long double)-x/b*c);
		yy=floor((long double)y/a*c);
		ans=yy-xx+1;
		ans1=x*c/d+y*c/d+(b-a)/d*yy;
		ans2=x*c/d+y*c/d+(b-a)/d*xx;
	}
	if(ans<=0) printf("-1 -1\n0");
	else printf("%lld %lld\n%lld",min(ans1,ans2),max(ans1,ans2),ans);
}

数论