NOI Section 2 数论
目录
- 错拍问题
- 同余方程
- 容斥原理
- 康托展开
- 生成函数
- 组合数&CRT
- 乘法逆元
- min25
- Min-Max容斥
- 不等式证明
- 质因数分解
- 斐波那契数列
- 多项式
- 常系数齐次线性递推
- 高斯消元
- 牛顿迭代
- 插值
错排问题
用\(f(n)\)表示\(n\)个元素的错排问题，有公式：
\[f(n)=(n-1)[f(n-1)+f(n-2)],n\ge 3
\]
特别地，\(f(1)=0,f(2)=1,f(3)=2\)
P1595 信封问题
problem
求\(n\)个信不装回原来信封的情况数。
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <queue>
#include <vector>
using namespace std;
const int maxn=1e6+10;
const int maxm=5e5+10;
const int mod=1e9+7;
#define int long long
#define inf 0x3f3f3f3f3f
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int f[maxn];
int n;
signed main(){
n=read();
f[1]=0,f[2]=1;
for(int i=3;i<=n;i++) f[i]=(i-1)*(f[i-1]+f[i-2]);
printf("%lld",f[n]);
return 0;
}
同余方程
P1082[LOJ#2605][NOIp2012] 同余方程
problem
给定\(a,b\)，求最小的正整数\(x\)使得\(ax\equiv 1(\operatorname{mod} b)\)
constraint
\(N\le 2\times 10^9\)
solution
拓展欧几里得。
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=5e6+10;
const int maxm=5e5+10;
const int mod=998244353;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int exgcd(int a,int b,int &x,int &y){
if(!b) return x=1,y=0,a;
int d=exgcd(b,a%b,x,y);
int temp=x;x=y,y=temp-a/b*y;
return d;
}
int n,m,x,y;
signed main(){
n=read(),m=read();
exgcd(n,m,x,y);
printf("%lld",(x+m)%m);
return 0;
}
P1403 [AHOI2005]约数研究
problem
定义\(f(N)=x\)，\(x\)是\(N\)的约数个数。如\(12\)的约数有\(1,2,3,4,6,12\)，则\(f(12)=6\).
现给定\(n\)，求
\[\sum_{i=1}^nf(i)
\]
constraint
\(1\le n\le 10^6\)
solution
\([1,n]\)内有约数\(i\)的个数是\(\left\lfloor\dfrac{n}{i}\right\rfloor\)， 那么答案就是\(\sum\limits_{i=1}^n\left\lfloor\dfrac{n}{i}\right\rfloor\)
然而我们发现有很多\(\left\lfloor\dfrac{n}{i}\right\rfloor\)都是一样的，对于这些一样的数我们每次算一次，似乎很浪费时间
所以我们每次\(i\)跳到\(\lfloor\frac nj\rfloor=\lfloor\frac ni\rfloor+1\)这样的\(j\)上,对于中间一样的数一次性算掉
复杂度\(O(2\sqrt n)\)
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
using namespace std;
const int maxn=1e6+10;
const int maxm=5e5+10;
const int mod=1e9+7;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int n,ans;
signed main(){
n=read();
for(int i=1,j;i<=n;i=j+1) j=n/(n/i),ans+=(n/i)*(j-i+1);
printf("%lld",ans);
return 0;
}
P2158 [SDOI2008]仪仗队
problem
在一个\(n\times n\)的方阵中，原点能看到多少个点？
constraint
\(1\le n \le 4\times 10^5\)
solution
一个点\((x,y)\)能够被看到，当且仅当\(\gcd(x,y)=1\),所以也就是求
\[\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=1]
\]
由于整个图形沿着\(y=x\)对称，且\(y=x\)上只能看到\((1,1)\)，我们就可以只考虑一个三角形内即可
所以答案变成
\[2\sum_{i=1}^n\sum_{j=i}^n[\gcd(i,j)=1]+1,n\ge 2
\]
由于欧拉函数\(\phi (n)\)的定义是\(\le n\)的数中与\(n\)互质的个数，即\(\phi(n)=\sum_{i=1}^n[\gcd(i,n)=1]\)
所以原式子改为
\[2\sum_{i=1}^n\phi(i)+1,n\ge 2
\]
求欧拉函数前缀和即可
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
using namespace std;
const int maxn=1e6+10;
const int maxm=5e5+10;
const int mod=1e9+7;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int n,ans,p[maxn];
signed main(){
n=read();
for(int i=1;i<=n;i++) p[i]=i;
for(int i=2;i<=n;i++) if(p[i]==i) for(int j=i;j<=n;j+=i) p[j]=p[j]*(i-1)/i;
for(int i=1;i<n;i++) ans+=p[i];
printf("%lld\n",(n==1)?0:ans<<1|1);
return 0;
}
P2455 [SDOI2006]线性方程组
problem
求方程组
已知 \(n\) 元线性一次方程组。
\[\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases}
\]
请根据输入的数据，编程输出方程组的解的情况。
constraint
\(n\le 50\)
solution
高斯消元。
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstring>
using namespace std;
const int maxn=1e6+10;
const int maxm=5e5+10;
const int mod=1e9+7;
#define int long long
#define random(a,b) ((a)+rand()%((b)-(a)+1))
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
double fread(){
double x;scanf("%lf",&x);return x;
}
int n,b[maxn],fg;
double a[120][120],f[120];
signed main(){
srand(time(NULL));
n=read();
for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) a[i][j]=fread();
for(int i=1;i<=n;i++) for(int j=i+1;j<=n;j++) {
int x=random(1,1000);
if(x&1) for(int k=1;k<=n+1;k++) swap(a[i][k],a[j][k]);
}
for(int i=1;i<=n;i++){ //高斯消元
int mx=i;
for(int j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[mx][i])) mx=j;
if(!a[mx][i]) continue;
for(int j=1;j<=n+1;j++) swap(a[i][j],a[mx][j]);
for(int j=1;j<=n;j++) if(j!=i){
double t=a[j][i]/a[i][i];
for(int k=1;k<=n+1;k++) a[j][k]-=a[i][k]*t;
}
}
for(int i=1;i<=n;i++){
if(a[i][i]==0){
if(a[i][n+1]) return !printf("-1");
else fg=1;
continue;
}
f[i]=a[i][n+1]/a[i][i];
if(!f[i]) f[i]=0.00;
}
if(fg) return !printf("0");
for(int i=1;i<=n;i++) printf("x%lld=%.2lf\n",i,f[i]);
return 0;
}
P4296 [AHOI2007]密码箱
problem
给定\(n\)，求所有\(x\in[1,n]\)使得\(x^2\equiv n(\operatorname{mod} 1)\)
constraint
\(n\le 2\times 10^9\)
solution
我们有\(x^2-1=kn,(x-1)(x+1)=kn\)，可以构造\(a,b\)使得\(a\mid (x+1),b\mid (x-1)\)或者\(a\mid(x-1),b\mid(x+1)\)，我们可以枚举\(a,b\)找到所有的\(x\)
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=1e6+10;
const int maxm=5e5+10;
const int mod=1e9+7;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int b,n;
set<int> st;
signed main(){
n=read();
if(n==1) return !printf("None");
st.insert(1);
for(int a=1;a<=sqrt(n);a++){
if(!(n%a)) {
b=n/a;
for(int i=b+1;i<=n;i+=b) if(!((i+1)%a)) st.insert(i);
for(int i=b-1;i<=n;i+=b) if(!((i-1)%a)) st.insert(i);
}
}
for(std::set<int>::iterator i=st.begin();i!=st.end();i++) printf("%lld\n",*i);
return 0;
}
P1641 [SCOI2010]生成字符串
problem
把\(n\)个\(1\)和\(m\)个\(0\)排成一列，使得任意前\(k\)个字符的\(1\)的个数不能少于\(0\)的个数，求方案数模\(20100403\)
constraint
\(n\le 10^6\)
solution
在二维平面上，把\((x,y)\)中的\(x\)看做1和0数量的和，\(y\)看做1和0数量的差
则：
- 向右上走，表示这一位选\(1\)
- 向左下走，表示这一位选\(0\)
那么走到\((n+m,n-m)\)是所有的方案数，相当于从\(n+m\)步中选\(m\)步是向右上走，即\(C(n+m,m)\)
其中不符合条件的是行走过程中走到了\(y=-1\)这条线，由于对称性可知相当于从\((0,2)\)走到了\((n+m,n-m)\)，方案数是\(C(n+m,m-1)\)
相减就是答案
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=5e6+10;
const int maxm=5e5+10;
const int mod=20100403;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int n,m;
int qpow(int a,int b){
int ret=1;
while(b){
if(b&1) ret=ret*a%mod;
a=a*a%mod;
b>>=1;
}
return ret;
}
int fac[maxn],inv[maxn];
int c(int x,int y){
int z=fac[x]*inv[y]%mod;
return z*inv[x-y]%mod;
}
signed main(){
n=read(),m=read();
fac[0]=1;
for(int i=1;i<=n+m;i++) fac[i]=fac[i-1]*i%mod;
inv[n+m]=qpow(fac[n+m],mod-2);
for(int i=n+m-1;i>=0;i--) inv[i]=inv[i+1]*(i+1)%mod;
printf("%lld\n",(c(n+m,m)-c(n+m,m-1)+mod)%mod);
return 0;
}
康托展开
problem
给定\(N\)，求给定的\(1\sim N\)的排列是\(N\)全排列的第几个。
constraint
\(N\le 10^6\)
solution
\[f(1\sim N)=1+\sum_{i=1}^Na_i\times(n-i)!
\]
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=5e6+10;
const int maxm=5e5+10;
const int mod=998244353;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int n,a[maxn],c[maxn];
int fac[maxn],ans;
int lowbit(int x){return x&(-x);}
void add(int x,int k){
while(x<=n) c[x]+=k,x+=lowbit(x);
}
int sum(int x){
int t=0;
while(x) t+=c[x],x-=lowbit(x);
return t;
}
signed main(){
n=read();
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=fac[i-1]*i%mod;
for(int i=1;i<=n;i++) add(i,1);
for(int i=1;i<=n;i++){
a[i]=read();
ans=(ans+((sum(a[i])-1)*fac[n-i])%mod)%mod;
add(a[i],-1);
}
printf("%lld\n",ans+1);
return 0;
}
组合数学
P1128 [HNOI2001]求正整数
problem
给定\(n\)，求含有\(n\)个不同因子的最小正整数\(m\)。
constraint
\(n\le 5\times 10^4\)
solution
- 错误的做法：贪心
反例：输入为8时，答案应是\(24=2^3\times 2\)，而非\(30=2\times 3\times 5\) - 正确的做法：dp（dfs+剪枝也可以）
考虑\(f(i,j)\)表示包含前\(j\)个质因数，因子个数为\(i\)的最小数。因为：
\(\tau(n)=\prod_i(k_i+1), \rm{where}\;n=\prod_i p_i^{k_i}\)
考虑到枚举最后一个质数的指数，有转移方程：
\[f_{i, j} = \min_{k|i}\left(f_{\frac ik, j-1}p_j^{k-1}\right)
\]
考虑取对数，dp变成了
\[f_{i, j} = \min_{k|i}\left(f_{\frac ik, j-1}+(k-1)\log p_j\right)
\]
最后要求结果的时候只要找到转移的方向把对应的\(p_j^{k-1}\)乘上去就好了。需要使用高精乘低精
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=5e6+10;
const int maxm=5e5+10;
const int mod=998244353;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
const int p[20]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71};
double logp[20];
double f[600][20];
int d[maxn],cnt,A[maxn],len;
int n,m;
void mul(int x){
int v=0;
for(int i=0;i<len;i++) v=(A[i]=A[i]*x+v)/10,A[i]%=10;
while(v) A[len++]=v%10,v/=10;
}
signed main(){
m=0,n=read();
for(int i=1;i<=n;i++) if(n%i==0) d[m++]=i;
for(int i=0;i<20;i++) f[0][i]=0.0;
for(int i=0;i<20;i++) logp[i]=log(p[i]);
for(int i=1;i<m;i++){
for(int k=0;k<20;k++) f[i][k]=1e9;
for(int j=0;j<i;j++) if(d[i]%d[j]==0){
int t=d[i]/d[j];
for(int k=1;k<20;k++) f[i][k]=min(f[i][k],f[j][k-1]+logp[k-1]*(t-1));
}
}
A[0]=len=1;
int j=0;
for(int i=0;i<20;i++) if(f[m-1][i]<f[m-1][j]) j=i;
for(int i=m-1,nxt;i;i=nxt,--j){
for(nxt=0;d[i]%d[nxt]||f[i][j]<f[nxt][j-1]
+logp[j-1]*(d[i]/d[nxt]-1)-1e-5;++nxt);
for(int k=1;k<d[i]/d[nxt];++k) mul(p[j-1]);
}
while(len--) printf("%lld",A[len]);
return 0;
}
P2303 [SDOI2012] Longge 的问题
problem
给定\(n\)，求\(\sum\limits_{i=1}^n\gcd(i,n)\)
constraint
\(n\le 2^{32}\)
solution
拆式子：
\[\sum_{i=1}^n\gcd(i,n)=\sum_{j\mid n} j\sum_{i=1}^n[\gcd(i,n)=j]=\sum_{j\mid n}j\times\sum_{i=1}^{\frac{n}{d}}[\gcd(i,\dfrac{n}{j})=1] =\sum_{j\mid n}j\phi(\dfrac{n}{j})
\]
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=5e6+10;
const int maxm=5e5+10;
const int mod=998244353;
const int N=(1<<16)+5;
#define int long long
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int n,tot,p[N];
bool flag[N];
void sieve(int n){
for(int i=2;i<=n;i++) {
if(!flag[i]) p[++tot]=i;
for(int j=1;j<=tot&&i*p[j]<=n;j++){
flag[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
}
int phi(int x){
int ans=x;
for(int i=1;i<=tot&&p[i]*p[i]<=x;i++){
if(x%p[i]) continue;
ans=ans/p[i]*(p[i]-1);
while(x%p[i]==0) x/=p[i];
}
if(x>1) ans=ans/x*(x-1);
return ans;
}
signed main(){
int n,ans=0;
n=read();
sieve((int)sqrt(n));
for(int i=1;i*i<=n;i++){
if(n%i==0){
ans+=i*phi(n/i);
if(i*i!=n) ans+=n/i*phi(i);
}
}
printf("%lld\n",ans);
return 0;
}
P2522 [HAOI2011]Problem b
problem
给定\((a,b,c,d,k)\)，求：
\[\sum_{i=a}^b\sum_{j=c}^d[\gcd(i,j)=k]
\]
\(n\)组数据。
constraint
\(a,b,c,d,n,k\le 5\times 10^4\)
solution
推式子+整除分块。
首先是一个容斥：\(Ans(a,b,c,d)=Ans(1,b,1,d)-Ans(1,a-1,1,d)-Ans(1,b,1,c-1)+Ans(1,a-1,1,c-1)\)
所以我们只需要求\(Ans(1,x,1,y)\)即可，也就是\(\sum\limits_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=k]\)
接下来开始推式子,设
\[F(n)=\sum_{k\mid n}f(k)=\left\lfloor \dfrac{a}{n}\right\rfloor \left\lfloor \dfrac{b}{n}\right\rfloor
\]
(这个式子是根据莫比乌斯反演得来的)
根据莫比乌斯反演，我们得到：
\(f(n)=\sum_{n\mid k}\mu(\left\lfloor \dfrac{k}{n}\right\rfloor )F(k)\)
答案\(\rm{Ans}=f(d)\)
于是
\[\texttt{Ans} =\sum_{d\mid k}\mu(\left\lfloor \dfrac{k}{d}\right\rfloor )F(k)
\]
枚举\(\left\lfloor \dfrac{k}{d}\right\rfloor\)设为\(t\)，则
\[\texttt{Ans}=\sum_{t=1}^{\min(a,b)}\mu(t)\left\lfloor \dfrac{a}{td}\right\rfloor\left\lfloor \dfrac{b}{td}\right\rfloor
\]
运用整除分块。
- \(\texttt{long long}\)就\(\rm{TLE}\)，改成\(\rm{int}\)就\(\rm{AC}\)?!离谱...
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=5e6+10;
const int mod=998244353;
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
int a,b,c,d,k,x,y,n,ans,t;
int gcd(int x,int y){
return !y?x:gcd(y,x%y);
}
bool vis[maxn];
int prim[maxn],mu[maxn],sum[maxn],cnt;
void get_mu(int n){
mu[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]) mu[i]=-1,prim[++cnt]=i;
for(int j=1;j<=cnt&&i*prim[j]<=n;j++){
vis[i*prim[j]]=1;
if(i%prim[j]==0) break;
else mu[i*prim[j]]=-mu[i];
}
}
for(int i=1;i<=n;i++) sum[i]=sum[i-1]+mu[i];
}
int rep;
int calc(int a,int b){
rep=min(a,b),ans=0;
for(int l=1,r;l<=rep;l=r+1){
r=min(a/(a/l),b/(b/l));
ans+=(1ll*a/(1ll*l*k))*(1ll*b/(1ll*l*k))*(sum[r]-sum[l-1]);
}
return ans;
}
signed main(){
t=read();
get_mu(50000);
while(t--){
ans=0;
a=read(),b=read(),c=read(),d=read(),k=read();
printf("%lld\n",calc(b,d)-calc(b,c-1)-calc(a-1,d)+calc(a-1,c-1));
}
return 0;
}
P2257 YY的GCD
problem
给定\((a,b)\)，求\(\sum_{i=1}^a\sum_{j=1}^b[\gcd(i,j)\in\texttt P]\)，\(\texttt{P}\)是质数集。
constraint
\(T=10^4,N,M\le 10^7\)
solution
莫比乌斯反演+整除分块。题解咕咕了。
- 也许是第一道卡常题，头一次把
#define int long long删掉，改成局部long long才过，见识到了见识到了
code
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <set>
#include <cstring>
using namespace std;
const int maxn=2e7+10;
const int mod=998244353;
int read(){
int a=0,op=1;char c=getchar();
while(c<'0'||c>'9') {if(c=='-') op=-1;c=getchar();}
while(c>='0'&&c<='9') a*=10,a+=c^48,c=getchar();
return a*op;
}
bool vis[maxn];
int prim[maxn],mu[maxn],g[maxn],cnt;
long long sum[maxn];
void get_mu(int n){
mu[1]=1;
for(int i=2;i<=n;i++){
if(!vis[i]) mu[i]=-1,prim[++cnt]=i;
for(int j=1;j<=cnt&&prim[j]*i<=n;j++){
vis[i*prim[j]]=1;
if(i%prim[j]==0) break;
else mu[prim[j]*i]=-mu[i];
}
}
for(int j=1;j<=cnt;j++) for(int i=1;i*prim[j]<=n;i++) g[i*prim[j]]+=mu[i];
for(int i=1;i<=n;i++) sum[i]=sum[i-1]+(int)g[i];
}
int T,n,m;
long long ans;
signed main(){
T=read();
get_mu(10000000);
while(T--){
n=read(),m=read();
if(n>m) swap(n,m);
ans=0;
for(int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
}
printf("%lld\n",ans);
}
return 0;
}
problem
constraint
solution
code
浙公网安备 33010602011771号