数论及其应用——解不定方程
在初等代数中，我们熟悉二元一次方程组的求法，但是很多时候，我们有n个变量，却没有n个彼此独立的方程，因此我们是无法给出方程组的唯一解的，即其解情况是不确定的。
定义1：a,b,c是整数，ab != 0，那么形为ax + by = c的方程成为二元一次不定方程。
定理1：设d = gcd(a,b)，如果d|c，那么方程存在无穷多个整数解，否则方程不存在整数解。
其实这条定理是和扩展欧几里得算法是有关系的。在先前的文章中，我们给出了扩展欧几里得定理，即存在整数对(x,y),使得gcd(a,b) = xa + yb。既然存在了一组整数解，我们再从几何的角度讲这个二元一次不等方程看成直线的解析式，就很好理解它为什么有无数组整数解了。
基于对二元不定方程的理解，我们可以马上把我们得到的结论推广到n元不定方程。
定理2：n元不定方程a1x1 + a2x2 +……+anxn = c(a1、a2、a3……an均为整数)有解的充要条件是gcd(a1,a2,……,an)|c。
基于这个定理，我们通过一个题目来实现它的应用。(Problem source : pku 1091)
Description
比如当N=2，M=18时，持有卡片(10, 15, 18)的跳蚤，就可以完成任务：他可以先向左跳10个单位长度，然后再连向左跳3次，每次15个单位长度，最后再向右连跳3次，每次18个单位长度。而持有卡片(12, 15, 18)的跳蚤，则怎么也不可能跳到距他左边一个单位长度的地方。
当确定N和M后，显然一共有M^N张不同的卡片。现在的问题是，在这所有的卡片中，有多少张可以完成任务。
Input
Output
#include<stdio.h> #define N 64 using namespace std; int bo[N],t; void divide(int m) { int i; t = 0; for(i = 2;i*i <= m;i++) { if(m%i == 0) { t++; bo[t] = i; while(m%i == 0) m /= i; } } if(m != 1) { t++; bo[t] = m; } } long long quick_multi(long long a , long long b) //快速计算a^b { long long ans = 1; while(b) { ans *= a; b--; } return ans; } long long ans , temp; int a[N],m,n; void dfs(int b,int ct ,int c) { int i , x; if(ct == c) { x = m; for(i = 1;i <= c;i++) x/= a[i]; temp += quick_multi(x , n); return; } for(i = b + 1;i<=t;i++) { a[ct + 1] = bo[i]; dfs(i , ct + 1 , c); } } int main() { int i; while(scanf("%d%d",&n,&m) != EOF) { ans = t = 0; divide(m); for(i = 1;i<=t;i++) { temp = 0; dfs(0,0,i); if(i&1) ans += temp; else ans -= temp; } ans = quick_multi(m,n) - ans; printf("%lld\n",ans); } return 0; }
介绍了判断n元一次不定方程是否有解的算法，这里我们介绍如何求解不定方程。首先我们来探讨求解形如ax + by = c的二元一次不定方程。
在学习扩展欧几里得算法的时候我们曾提及它在求逆元、解同余方程、不定方程中有着重要的作用。其实很好理解，扩展欧几里得是求解满足ax + by = gcd(a , b)的整数对(x , y)，这是二元不定方程的形式。我们稍作修改，可将不定方程转化成同余方程——ax = gcd(a , b) (mod y)。类似的，我们求解ax = 1(mod m)，就是求整数a在模m情况下的逆元。这就很解释扩展欧几里得为何在解决这三个问题上发挥着如此大的作用了。
下面我们来探讨利用扩展欧几里得求解二元一次不定方程的最小解问题。(Problem source : pku 2142)
Description
Input
Output
- You can measure dmg using x many amg weights and y many bmg weights.
- The total number of weights (x + y) is the smallest among those pairs of nonnegative integers satisfying the previous condition.
- The total mass of weights (ax + by) is the smallest among those pairs of nonnegative integers satisfying the previous two conditions.
题目大意：给出两种规格的砝码a、b，和一个整数n，现在请你求解出一种在天平上放置砝码的方案称出质量为n的物体，并且使得所用砝码数量最少，如果数量相同，则使得砝码的总质量最少。
数理分析：这是一道需要我们进行数学抽象的实际问题，我们不难想象，这道题目本质上是要求我们求解不定方程ax + by = n，并使得|x| + |y|最小，如果两种不同方案中使用的砝码数相同，则使得a|x| + b|y|最小。
这里x、y之所以添加绝对值，是考虑在实际问题中，利用正负号区别两种砝码放在在一个盘中还是放在不同的盘中，这是很好理解的。
那么下面我们应该关注的应该是如何求解这个不定方程。首先我们知道，利用扩展欧几里得算法我们可以得到ax + by = gcd(a , b)的一个特解，因此我们也不难得到ax + by = n的一个特解，因为我们知道这个方程有根的条件是gcd(a , b) | n，因此在ax + by = gcd(a ,b)的特解上乘以一个整数即可。
那么现在对于方程ax + by = n，我们手中有一组特解，然而我们需要找到符合题设要求的特解，如何实现呢？
显然，我们筛选符合要求的特解是需要在含有所有的解的解空间中找的，因此我们需要解决的问题是如何表示出这个二元不定方程的通解形式。 设(x0,y0)是二元不定方程ax + by = c的一组特解，并且满足gcd(a , b) = 0,那么x = x0 + bt ， y = y0 - at是该方程的通解形式。
证明如下：
1.将x ， y带入方程，等式成立，表明x，y是方程的解。
2.注意到gcd(a , b) = 1，这使得解空间中各解的间隔最小，因此一定会覆盖所有解，表明x，y是方程的通解。
而对于任意一个题目中给出的有解的二元不定方程ax + by = n，必然能够化成 a'x + b'y = n' 的形式，其中满足gcd(a',b') =1。那么通解我们就很容易表示出来了。
x = x0 + bt
y = y0 - at
不难发现，x与y都和整数t呈线性相关，t增大，|x|、|y|同时增大，那么我们要求(|x| + |y|)min，就需要让t尽量小，也就是求x的最小正整数解，并根据方程求出此时的y；求出y的最小正整数解，并根据方程求出此时的x，那么满足（|x|+|y|）min的方案一定会这两个方案中出现。如果相等，则继续比较ax + by即可。
而基于x = x0 + bt的通解形式，我们不难求出其最小整数解——xmin = (x0%b + b)%b，也许读者看到这个形式会感觉很熟悉，这其实是与前面求解同余方程的最小整数解呼应起来的，不过这里是通过通解的形式的角度来看的，实际上更好理解。
基于以上的数理分析，通过调用扩展欧几里得，不难完成编程实现。
ps：理论上比较两方案的砝码数然后比较砝码质量，但是实际编程中发现添加比较砝码总质量的过程会返回WA，下面给出AC代码(其中比较砝码总质量的过程被注释掉)。
#include<cstdio> #include<iostream> using namespace std; int gcd(int a , int b) { if(b==0) return a; return gcd(b,a%b); } void extend_Eulid(int a , int b,int &x,int &y) //解ax + by = 1 { if(b==0) { x = 1; y = 0; return; } extend_Eulid(b,a%b,x,y); int temp = x; x = y; y = temp - a/b*y; } int main() { int a , b , n; while(scanf("%d%d%d",&a,&b,&n) != EOF) { if(a + b + n == 0) break; int x , y; int gcdab = gcd(a,b); a /= gcdab; b /= gcdab; n /= gcdab; if(b > a) { int temp; temp = a; a = b; b = temp; } extend_Eulid(a , b , x ,y); int vx = x * n; vx = (vx % b + b)%b; int vy = (n - a * vx)/b; if(vy < 0) vy = -vy; y = y * n; y = (y % a + a) % a; x = (n - b * y)/a; if(x < 0) x = -x; if(x + y > vx + vy) //比较砝码数量 x = vx , y = vy; /*if(x + y == vx + vy) 比较砝码总质量 { if(a*x + b*y > a*vx + b*vy) x = vx , b = vy; }*/ printf("%d %d\n",x,y); } }
上文介绍了解n元一次不定方程的方法和如何判读n元一次方程是否有解之后，接下来我们将讨论几个特殊的不定方程。
毕达哥拉斯三元组：满足x^2 + y^2 = z^2的三整数x、y、z被称为毕达哥拉斯三元组。如果gcd(x , y , z) = 1，则称这个毕达哥拉斯三元组是本原的。
基于本原毕达哥拉斯三元组的定义，数学家经过推导论证，发现了本原毕达哥拉斯三元组有如下的性质(此处折叠证明)： x = m^2 - n^2 y = 2mn z = m^2 + n^2 其中m和n的奇偶性是相反的。即如果m是偶数则n为奇数，m为奇数则n为偶数。
基于对这个不定方程的引入，我们通过一个题目来编程实现它的应用。(Problem souce : pku 1305)
Description
Input
题目大意：给出整数n，表示a,b,c<= n，然后请你求解满足这一范围的本原毕达哥拉斯三元组的组数，并求出[1,n]中不涉及本原毕达哥拉斯三元组的整数个数。 数理分析：我们容易想到设置三层穷举并根据定义来判断a、b、c是否是本原毕达哥拉斯，但是难免超时。这里如果我们考虑应用本原毕达哥拉斯三元组的性质，只需枚举[1,n]范围内的m、n即可，在得到一组后，在将三个数同时扩大i倍，直到三元组中最大的数不满足所设范围。 编程实现：有了以上的分析，编程上只需设置两层穷举，至于求解[1,n]范围内不涉及本原毕达哥拉斯三元组的数字，只需开一个标记数组flag[]，在生成本原毕达哥拉斯三元组顺便标记最后筛选一遍便可以得到结果。
参考代码如下。
#include<iostream> #include<cstdio> #include<cstdlib> #include<cmath> #include<cstring> using namespace std; const int N=1000001; bool flag[N]; int gcd(int a,int b) { return b==0?a:gcd(b,a%b); } void solve(int t) { int tmp,m,n,i,ans1,ans2,x,y,z; ans1=ans2=0; memset(flag,0,sizeof(flag)); tmp=(int)sqrt(t+0.0); for(n=1;n<=tmp;n++) { for(m=n+1;m<=tmp;m++) { if(m*m+n*n>t) break; if(n%2!=m%2) { if(gcd(m,n)==1) { x=m*m-n*n; y=2*m*n; z=m*m+n*n; ans1++; for(i=1;;i++) { if(i*z>t) break; flag[i*x]=1; flag[i*y]=1; flag[i*z]=1; } } } } } for(i=1;i<=t;i++) if(!flag[i]) ans2++; printf("%d %d\n",ans1,ans2);} int main() { int n; while(~scanf("%d",&n)) solve(n); return 0; }
其实根据毕达哥拉斯三元组的定义我们不难发现，它是基于毕达哥拉斯定理的，即三元组中的整数是可以组成直角三角形的，那么下面我们就来看一道关于直角三角形的问题。(Problem source : fzu 1669)
Problem Description
A triangle is one of the basic shapes of geometry: a polygon with three corners or vertices and three sides or edges which are line segments. A triangle with vertices A, B, and C is denoted △ABC.Triangles can also be classified according to their internal angles, described below using degrees of arc:
- A right triangle (or right-angled triangle, formerly called a rectangled triangle) has one 90° internal angle (a right angle). The side opposite to the right angle is the hypotenuse; it is the longest side in the right triangle. The other two sides are the legs or catheti (singular: cathetus) of the triangle. Right triangles conform to the Pythagorean Theorem, wherein the sum of the squares of the two legs is equal to the square of the hypotenuse, i.e., a^2 + b^2 = c^2, where a and b are the legs and c is the hypotenuse.
- An oblique triangle has no internal angle equal to 90°.
- An obtuse triangle is an oblique triangle with one internal angle larger than 90° (an obtuse angle).
- An acute triangle is an oblique triangle with internal angles all smaller than 90° (three acute angles). An equilateral triangle is an acute triangle, but not all acute triangles are equilateral triangles.
What we consider here is very simple. Give you the length of L, you should calculate there are how many right-angled triangles such that a + b + c ≤ L where a and b are the legs and c is the hypotenuse. You should note that the three numbers a, b and c are all integers.
Input
There are multiply test cases. For each test case, the first line is an integer L(12≤L≤2000000), indicating the length of L.
Output
For each test case, output the number of right-angled triangles such that a + b + c ≤ L where a and b are the legs and c is the hypotenuse.
题目大意：给出整数L，求解周长小于L的所有直角三角形的种类数。
数理分析：从数论的角度去看这个问题，其实就是求解在满足x + y + z <= L的情况下，满足x^2 + y^2 = z^2 的整数组合(x , y , z)，也就是毕达哥拉斯三元组。
基于我们上文对本原毕达哥拉斯三元组的算法实现，我们想方设法将这个问题往上面的模型去靠拢。我们发现，类似于数论中的基本算术定理，其实任意一个毕达哥拉斯三元组都是可以有本原毕达哥拉斯三元组得到的，因此我们可以先找出本原毕达哥拉斯三元组的情况，让后将x、y、z同时乘以整数i，直到(x + y + z)*i>L停止记数。
有了上述的数理分析，在编程实现上，只需对记录本原毕达哥拉斯的模板上稍作改动即可。
参考代码如下。
#include<iostream> #include<cstdio> #include<cstdlib> #include<cmath> #include<cstring> using namespace std; int gcd(int a,int b) { return b==0?a:gcd(b,a%b); } void solve(int t) { int tmp,m,n,i,ans1,x,y,z; ans1 = 0; tmp=(int)sqrt(t+0.0); for(n=1;n<=tmp;n++) { for(m=n+1;m<=tmp;m++) { if(2*m*m+2*m*n>t) break; if(n%2!=m%2) { if(gcd(m,n)==1) { x=m*m-n*n; y=2*m*n; z=m*m+n*n; for(i = 1;;i++) { if(i*(x + y + z) <= t) ans1++; else break; } } } } } printf("%d\n",ans1); } int main() { int n; while(~scanf("%d",&n)) solve(n); return 0; }
今天将介绍一个数论中经典的不定方程——佩尔方程。
定义：形如x^2 - Dy^2 = 1(D > 1且D不是完全平方数)的方程叫做佩尔方程。
由定义我们很容易看到，如果D是完全平方数，则原方程可以写成(x + D^0.5 * y)*(x - D^0.5*y) = 1，由题目的限制条件，显然这个方程是无解的。
这里直接引入数论中的定理，折叠证明。(时间实在有限)
对于佩尔方程，设其最小整数解是(x1,x2)，则存在如下的线性递推关系，来表征第n组解。(注意，这里默认解空间是升序排列)
xn = x(n-1)x1 + Dy(n-1)y1
yn = x(n-1)y1 + y(n-1)x1
基于佩尔方程的解存在这样的递推关系，我们根据矩阵的相关定义，容易构造出矩阵的表达形式，随后我们便可以通过矩阵快速幂，来实现快速访问某个佩尔方程的解空间。
而很显然，为了实现上面的做法，我们需要非常关键的一组数据——该佩尔方程的最小解，如何求呢？
很容易想到的就是穷举法，在求得一组最小解后，即可通过上面矩阵表示的线性递推式，来访问佩尔方程的解空间了。
我们通过一个具体的题目来实践一下。(Problem source : hdu 3292)
题目大意：给出N、K，求解满足X^2 - Ny^2 = 1的解空间中，第K大的解X[k]。
数理分析：通过读题不难将这个实际问题转化到求解佩尔方程的数学模型上来。
编程实现：有了以上分析，再编程实现上只需调用矩阵快速幂算法便可轻松实现求解。
参考代码如下。
#include<iostream> #include<stdio.h> #include<cmath> #define MAXN 4 #define M 8191; using namespace std; typedef struct Matrax { int m[MAXN][MAXN]; }Matrax; Matrax per , d; int n; int x , y , D; void search() { y = 1; while(1) { x = (long long)sqrt(D * y * y + 1); if(x*x - D*y*y == 1) break; y++; } } void init() { d.m[0][0] = x%M; d.m[0][1] = D*y % M; d.m[1][0] = y%M; d.m[1][1] = x%M; for(int i = 0;i < 2;i++) for(int j = 0;j < 2;j++) per.m[i][j] = (i == j); } Matrax multi(Matrax a,Matrax b) { Matrax c; int k , i , j; for(i = 0;i < 2;i++) for(j = 0;j < 2;j++) { c.m[i][j] = 0; for(k = 0;k < 2;k++) { c.m[i][j] += a.m[i][k]*b.m[k][j]; } c.m[i][j]%=M; } return c; } Matrax power(int k) { Matrax p , ans = per; p = d; while(k) { if(k&1) { ans = multi(ans , p); k--; } k /=2; p = multi(p , p); } return ans; } int main() { int K; while(cin >> D >> K) { int ad = sqrt(D + 0.0); if(ad * ad == D) { printf("No answers can meet such conditions\n"); continue; } search(); n = 2; init(); d = power(K-1); int res = d.m[0][0]*x%M; res += d.m[0][1]*y%M; res %=M; printf("%d\n",res); } return 0; }
