二元一次不定方程
定义:
a
,
b
,
c
a,b,c
a,b,c是整数,
a
b
≠
0
ab≠0
ab=0,那么形如
a
x
+
b
y
=
c
ax+by=c
ax+by=c的方程称为二元一次不定方程。
定理:设
a
,
b
a,b
a,b是整数,且
d
=
(
a
,
b
)
d=(a,b)
d=(a,b),如果
d
∣
c
d|c
d∣c,那么方程存在无数多个整数解,否则方程不存在整数解。
(
(
a
,
b
)
=
g
c
d
(
a
,
b
)
)
((a,b)=gcd(a,b))
((a,b)=gcd(a,b))
二元一次不定方程和同余方程之间可以相互转换的,例如在
a
>
0
,
b
>
0
a>0,b>0
a>0,b>0的条件下,求解二元一次不定方程
a
x
+
b
y
=
c
ax+by=c
ax+by=c与求解同余方程
a
x
≡
c
(
m
o
d
b
)
ax≡c(mod\ b)
ax≡c(mod b)是相同的。
求解同余方程使用扩展欧几里得算法,二元一次不定方程同样可以使用扩展欧几里得算法。
由扩展欧几里得可知:设
a
a
a和
b
b
b不全为0,则存在整数
x
,
y
x,y
x,y,使得:
g
c
d
(
a
,
b
)
=
a
x
+
b
y
gcd(a,b)=ax+by
gcd(a,b)=ax+by
所以当
d
∣
c
d|c
d∣c时,
a
x
+
b
y
=
c
ax+by=c
ax+by=c存在整数解。
通过扩展欧几里得算法求出特解
x
0
,
y
0
x_0,y_0
x0,y0,那么通解为
x
0
+
n
×
b
d
,
y
0
−
n
×
a
d
(
n
∈
Z
)
x_0+n\times\frac{b}{d},y_0-n\times\frac{a}{d}\ (n∈Z)
x0+n×db,y0−n×da (n∈Z).
void exgcd(int a,int b,int&d,int&x,int&y)//扩展欧几里得算法
{
if(b==0)
{
x=1,y=0;
d=a;
}
else
{
exgcd(b,a%b,d,x,y);
int t=x;
x=y;
y=t-a/b*y ;
}
}
n n n元一次不定方程
定理:
n
n
n元一次不定方程
a
1
x
1
+
a
2
x
2
+
.
.
.
+
a
n
x
n
=
c
(
a
1
,
a
2
,
.
.
.
,
a
n
,
c
∈
N
)
a_1x_1+a_2x_2+...+a_nx_n=c(a_1,a_2,...,a_n,c∈N)
a1x1+a2x2+...+anxn=c(a1,a2,...,an,c∈N),有解的充分必要条件是
(
a
1
,
a
2
,
.
.
.
,
a
n
)
∣
c
(a_1,a_2,...,a_n)|c
(a1,a2,...,an)∣c.
解
n
n
n元一次不定方程可以转化为求二元一次不定方程,
n
n
n元一次不定方程
a
1
x
1
+
a
2
x
2
+
.
.
.
+
a
n
x
n
=
c
a_1x_1+a_2x_2+...+a_nx_n=c
a1x1+a2x2+...+anxn=c,可以转为为下面形式:
S
{
a
1
x
1
+
a
2
x
2
=
d
1
t
1
d
1
t
1
+
a
3
x
3
=
d
2
t
2
d
2
t
2
+
a
4
x
4
=
d
3
t
3
.
.
.
d
n
−
1
t
n
−
1
+
a
n
x
n
=
c
S\left\{ \begin{aligned} &a_1x_1+a_2x_2=d_1t_1&\\ &d_1t_1+a_3x_3=d_2t_2&\\ &d_2t_2+a_4x_4=d_3t_3&\\ &\ \ \ \ .&\\ &\ \ \ \ .&\\ &\ \ \ \ .&\\ &d_{n-1}t_{n-1}+a_nx_n=c&\\ \end{aligned} \right.
S⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1x1+a2x2=d1t1d1t1+a3x3=d2t2d2t2+a4x4=d3t3 . . .dn−1tn−1+anxn=c
(
d
n
−
1
,
a
n
)
=
(
a
1
,
a
2
,
.
.
.
,
a
n
)
(d_{n-1},an)=(a_1,a_2,...,a_n)
(dn−1,an)=(a1,a2,...,an),所以上述定理成立。
求出最后一个方程的解,然后代入倒数第二个方程,依此类推可以求出方程的所有解。
对于
m
m
m个
n
n
n元一次不定方程组成的方程组,其中
m
<
n
m<n
m<n,可以消去
m
−
1
m-1
m−1个未知数,从而消去
m
−
1
m-1
m−1个不定方程,将方程组转化为一个
n
−
m
+
1
n-m+1
n−m+1元的一次不定方程。然后依照
n
n
n元一次不定方程求解。
解不定方程简单应用
The Balance(poj 2142)
有一天平,两个质量为
a
a
a和
b
b
b的砝码,数量不限,且两边都可以放,要求称出质量为
c
c
c的物品,即两边平衡。尽量少的使用砝码,当砝码数量相同时总重量最小。
输入为
a
,
b
,
c
a,b,c
a,b,c.
思路:如果a小于b交换a,b。显然可以构造不定方程
a
x
+
b
y
=
c
ax+by=c
ax+by=c,此时
x
x
x和
y
y
y不一定为正整数,可以先求出一组特解
x
0
,
y
0
x_0,y_0
x0,y0,通解就为
x
0
+
n
×
b
d
,
y
0
+
n
×
a
d
(
n
∈
Z
)
x_0+n\times\frac{b}{d},y_0+n\times\frac{a}{d}\ (n∈Z)
x0+n×db,y0+n×da (n∈Z).因为要求砝码数尽量少,即
m
i
n
(
∣
x
∣
+
∣
y
∣
)
min(|x|+|y|)
min(∣x∣+∣y∣),代入通解为
m
i
n
(
∣
x
0
+
n
×
b
d
∣
+
∣
y
0
−
n
×
a
d
∣
)
min(|x_0+n\times\frac{b}{d}|+|y_0-n\times\frac{a}{d}|)
min(∣x0+n×db∣+∣y0−n×da∣),因为
a
>
b
a>b
a>b,显然最小值应在
y
0
d
a
\frac{y_0d}{a}
ay0d附近,可以for循环枚举前后的几个值,然后取
m
i
n
(
∣
x
∣
+
∣
y
∣
)
min(|x|+|y|)
min(∣x∣+∣y∣),当
m
i
n
(
∣
x
∣
+
∣
y
∣
)
min(|x|+|y|)
min(∣x∣+∣y∣)相等时
m
i
n
(
a
∣
x
∣
+
b
∣
y
∣
)
min(a|x|+b|y|)
min(a∣x∣+b∣y∣)
代码如下:
#include <iostream>
#include <algorithm>
#include <math.h>
using namespace std;
typedef long long ll;
const int maxx=1e6+7;
void exgcd(ll a,ll b,ll&d,ll&x,ll&y)
{
if(!b)
{
x=1,y=0;
d=a;
}
else
{
exgcd(b,a%b,d,x,y);
ll temp=x;
x=y;
y=temp-a/b*y;
}
}
int main()
{
ll a,b,c,d,x0,y0;
ios::sync_with_stdio(false);
while(cin>>a>>b>>c&&(a||b||c))
{
ll x=0x3f3f3f3f,y=0x3f3f3f3f,flag=0;
if(a<b) {swap(a,b);flag=1;}
exgcd(a,b,d,x0,y0);
a=a/d;
b=b/d;
x0=x0*c/d;
y0=y0*c/d;
for(ll i=y0/a-10;i<=y0/a+10;i++)
{
ll q=x0+i*b,p=y0-i*a;
q=q<0?(-q):q; p=p<0?(-p):p;
if(q+p<x+y)
x=q,y=p;
else if(q+p==x+y&&(q*a+p*b)<(x*a+y*b))
x=q,y=p;
}
if(!flag)
cout<<x<<" "<<y<<endl;
else cout<<y<<" "<<x<<endl;
}
}
跳骚(poj 1091)
Z城市居住着很多只跳蚤。在Z城市周六生活频道有一个娱乐节目。一只跳蚤将被请上一个高空钢丝的正中央。钢丝很长,可以看作是无限长。节目主持人会给该跳蚤发一张卡片。卡片上写有N+1个自然数。其中最后一个是M,而前N个数都不超过M,卡片上允许有相同的数字。跳蚤每次可以从卡片上任意选择一个自然数S,然后向左,或向右跳S个单位长度。而他最终的任务是跳到距离他左边一个单位长度的地方,并捡起位于那里的礼物。
比如当N=2,M=18时,持有卡片(10, 15, 18)的跳蚤,就可以完成任务:他可以先向左跳10个单位长度,然后再连向左跳3次,每次15个单位长度,最后再向右连跳3次,每次18个单位长度。而持有卡片(12, 15, 18)的跳蚤,则怎么也不可能跳到距他左边一个单位长度的地方。
当确定N和M后,显然一共有M^N张不同的卡片。现在的问题是,在这所有的卡片中,有多少张可以完成任务。
输入两个整数N和M(N <= 15 , M <= 100000000)。
思路:假设卡片上面数字为
a
1
,
a
2
,
.
.
.
a
n
,
M
a_1,a_2,...a_n,M
a1,a2,...an,M,跳骚选择卡片的次数为
x
1
,
x
2
,
.
.
.
x
n
,
x
n
+
1
x_1,x_2,...x_n,x_{n+1}
x1,x2,...xn,xn+1,那么有不定方程
a
1
x
1
+
a
2
x
2
+
.
.
.
+
a
n
x
n
+
M
x
n
+
1
=
1
a_1x_1+a_2x_2+...+a_nx_n+Mx_{n+1}=1
a1x1+a2x2+...+anxn+Mxn+1=1,即
(
a
1
,
a
2
,
.
.
.
,
a
n
,
M
)
=
1
(a_1,a_2,...,a_n,M)=1
(a1,a2,...,an,M)=1,可以计算出
(
a
1
,
a
2
,
.
.
.
,
a
n
,
M
)
≠
1
(a_1,a_2,...,a_n,M)≠1
(a1,a2,...,an,M)=1的所有情况,然后用
M
n
M^n
Mn减去即可。计算
(
a
1
,
a
2
,
.
.
.
,
a
n
,
M
)
≠
1
(a_1,a_2,...,a_n,M)≠1
(a1,a2,...,an,M)=1的所有情况可以使用容斥定理计算,如果gcd不等于1的话,那么这些数肯定都至少有一个M的质因子,假设d是M的一个质因子,那么小于M的可以整除d的数个数量就是
M
d
\frac{M}{d}
dM个,从其中任选n个,因为可以重复选择,所以数量就是
M
d
n
\frac{M}{d}^n
dMn。
所以
a
n
s
=
⋃
i
=
1
r
M
d
i
n
(
r
为
M
的
质
因
子
数
量
)
ans=\bigcup\limits_{i=1}^{r}\frac{M}{d_i}^n\ (r为M的质因子数量)
ans=i=1⋃rdiMn (r为M的质因子数量)
由容斥原理可知:
∣
⋃
i
=
1
n
S
i
∣
=
∑
C
⊆
M
n
(
−
1
)
s
i
z
e
(
C
)
−
1
∣
⋂
T
⊆
C
T
∣
\left|\bigcup\limits_{i=1}^{n}S_i\right|=\sum\limits_{C\subseteq M}^{n}(-1)^{size(C)-1}\left|\bigcap\limits_{T\subseteq C}T\right|
∣∣∣∣i=1⋃nSi∣∣∣∣=C⊆M∑n(−1)size(C)−1∣∣∣∣∣T⊆C⋂T∣∣∣∣∣
#include <iostream>
#include <algorithm>
#include <math.h>
using namespace std;
typedef long long ll;
const int maxx=1e6+7;
ll n,m,p[35],cnt,sum,ans;
ll qpow(ll a,ll b)
{
ll ans=1;
while(b)
{
if(b&1)
ans=ans*a;
a=a*a;
b>>=1;
}
return ans;
}
void dfs(int now,int id,int over,ll product)
{
if(now==over)
{
sum+=qpow(m/product,n);
return;
}
for(int i=id;i<=cnt;i++)
{
dfs(now+1,i+1,over,product*p[i]);
}
}
int main()
{
ios::sync_with_stdio(0);
cin>>n>>m;
ll a=m;
for(int i=2;i*i<=a;i++)
if(a%i==0)
{
p[++cnt]=i;
while(a%i==0)
a/=i;
}
if(a>1) p[++cnt]=a;
for(int i=1;i<=cnt;i++)
{
sum=0;
dfs(0,1,i,1);
if(i&1)
ans+=sum;
else ans-=sum;
}
cout<<qpow(m,n)-ans<<endl;
}
特殊的不定方程
-
毕达哥拉斯定理:平面上的一个直角三角形中,两个直角边边长的平方加起来等于斜边长的平方。如果设直角三角形的两条直角边长度分别是 a a a和 b b b,斜边长度是 c c c,那么可以用数学语言表达: a 2 + b 2 = c 2 a^2+b^2=c^2 a2+b2=c2.
满足上面这个方程的正整数三元组被称为毕达哥拉斯三元组。
如果一个毕达哥拉斯三元组 x , y , z x,y,z x,y,z,满足 ( x , y , z ) = 1 (x,y,z)=1 (x,y,z)=1,那么这个毕达哥拉斯三元组被称为本原的毕达哥拉斯三元组。
定理:正整数 a , b , c a,b,c a,b,c构成一个本原毕达哥拉斯三元组且y为偶数,当且仅当存在互质的正整数 m , n ( m > n ) m,n(m>n) m,n(m>n)其中 m m m为奇数, n n n为偶数;或者 n n n为奇数, m m m为偶数,并且满足:- a = m 2 + n 2 a=m^2+n^2 a=m2+n2
- b = 2 m n b=2mn b=2mn
- c = m 2 − n 2 c=m^2-n^2 c=m2−n2
-
费马大定理:方程 x n + y n = z n x^n+y^n=z^n xn+yn=zn无非0整数解,其中 n n n为整数且 n ≥ 3 n≥3 n≥3.
-
佩尔方程:形如 x 2 − d y 2 = 1 x^2-dy^2=1 x2−dy2=1( d > 1 d>1 d>1且不为完全平方数)的不定方程称为佩尔方程。
佩尔方程当 d d d为完全平方数时, ( x + d y ) ( x − d y ) = 1 (x+\sqrt{d}y)(x-\sqrt{d}y)=1 (x+dy)(x−dy)=1,显然方程无解。
解佩尔方程:
若佩尔方程最小特解是(x_1,y_1),那么可有迭代公式:
x n = x n − 1 x 1 + d y n − 1 y 1 x_n=x_{n-1}x_1+dy_{n-1}y_1 xn=xn−1x1+dyn−1y1
y n = x n − 1 y 1 + y n − 1 x 1 y_n=x_{n-1}y_1+y_{n-1}x_1 yn=xn−1y1+yn−1x1
即:
[ x k y k ] = [ x 1 d y 1 y 1 x 1 ] k − 1 [ x 1 y 1 ] \begin{bmatrix} x_k\\ y_k\\ \end{bmatrix}={\begin{bmatrix} x_1&dy_1\\ y_1&x_1\\ \end{bmatrix}}^{k-1}\begin{bmatrix} x_1\\ y_1\\ \end{bmatrix} [xkyk]=[x1y1dy1x1]k−1[x1y1]
我们可以通过矩阵快速幂求出第k个解。
对于特解可以暴力求解:
由 x 2 − d y 2 = 1 x^2-dy^2=1 x2−dy2=1,得 x 2 = d y 2 + 1 x^2=dy^2+1 x2=dy2+1,即 x = d y 2 + 1 x=\sqrt{dy^2+1} x=dy2+1,那么令 y = 1 y=1 y=1,代入求解 x x x,然后判断是否满足,否则就 y = y + 1 y=y+1 y=y+1,直到求出满足的特解。
void find()
{
int x,y=1;
while(1)
{
x=(long long)sqrt(d*y*y+1);
if(x*x-d*y*y==1)
break;
y++;
}
}
特殊的不定方程的应用
Street Numbers(poj 1320)
题意:求前十组两个整数不相等
n
,
m
n,m
n,m,满足
(
n
<
m
)
(n<m)
(n<m),使得:
1
+
2
+
.
.
.
+
(
n
−
1
)
=
(
n
+
1
)
+
.
.
.
+
m
1+2+...+(n-1) =(n+1)+...+m
1+2+...+(n−1)=(n+1)+...+m
思路:
n
(
n
−
1
)
2
=
(
m
−
n
)
(
m
+
n
+
1
)
2
\frac{n(n-1)}{2}=\frac{(m-n)(m+n+1)}{2}
2n(n−1)=2(m−n)(m+n+1),化简可得
(
2
m
+
1
)
2
−
8
n
2
=
1
(2m+1)^2-8n^2=1
(2m+1)2−8n2=1
佩尔方程递推求解即可。
#include <iostream>
#include <algorithm>
using namespace std;
const int maxx=1e6+7;
int main()
{
ios::sync_with_stdio(0);
int x0=3,y0=1,px=3,py=1,d=8,x,y;
for(int i=1;i<=10;i++)
{
x=px*x0+d*py*y0;
y=px*y0+py*x0;
cout<<y<<" "<<(x-1)/2<<endl;
px=x,py=y;
}
}
Right-angled Triangle(fzu 1669)
求满足以a、b为直角边,c为斜边,而且满足a + b + c <= L的直角三角形的个数。
思路:根据本原毕达哥拉斯三元组满足的条件枚举n,m,求出本原毕达哥拉斯三元组后这个三元组的倍数也是毕达哥拉斯三元组。枚举所有即可。
#include <iostream>
#include <algorithm>
using namespace std;
const int maxx=1e6+7;
int main()
{
ios::sync_with_stdio(0);
int l;
while(cin>>l)
{
int ans=0;
for(int i=1;i*i<=l;i++)
{
for(int j=i+1;j*j<=l;j++)
{
if(i%2!=j%2)
{
int x=j*j-i*i,y=i*j*2,z=i*i+j*j;
if(__gcd(x,__gcd(y,z))==1)
{
for(int k=1;k*(z+y+x)<=l;k++)
{
ans++;
}
}
}
}
}
cout<<ans<<endl;
}
}
No more tricks, Mr Nanguo (hdu 3292)
给出d和k,求满足
x
2
−
d
y
2
=
1
x^2-dy^2=1
x2−dy2=1的第k小的x的值。
思路:佩尔方程先暴力求出最小解,然后矩阵快速幂求解即可。
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;
typedef long long ll;
const int maxx=1e6+7,mod=8191;
struct matrix
{
ll a[5][5];
}m;
matrix times(matrix x,matrix y,int n)
{
matrix temp;
memset(temp.a,0,sizeof temp);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
temp.a[i][j]=(temp.a[i][j]+x.a[i][k]*y.a[k][j])%mod;
return temp;
}
matrix qpow(matrix a,ll b,int n)
{
matrix ans;
memset(ans.a,0,sizeof ans.a);
for(int i=1;i<=n;i++)
ans.a[i][i]=1;
while(b)
{
if(b&1) ans=times(ans,a,n);
a=times(a,a,n);
b>>=1;
}
return ans;
}
int main()
{
ios::sync_with_stdio(0);
int n,k;
while(cin>>n>>k)
{
int temp=sqrt(n+0.0);
if(temp*temp==n)
cout<<"No answers can meet such conditions"<<endl;
else
{
ll x,y=1;
while(true)
{
x=(ll)sqrt(n*y*y+1);
if(x*x-y*y*n==1)
break;
y++;
}
m.a[1][1]=x%mod;
m.a[1][2]=n*y%mod;
m.a[2][1]=y%mod;
m.a[2][2]=x%mod;
matrix ans=qpow(m,k-1,2);
x=x*ans.a[1][1]%mod+y*ans.a[1][2]%mod;
cout<<x%mod<<endl;
}
}
}