数论及其应用——整除性问题
数论，顾名思义是将数作为研究对象的一个学科分支。而数论包括初等数论和高等数论，初等数论则主要研究整数的性质和方程的整数解(不要问我为什么都学高数了还在这搞初等数论，因为小学奥数经历完全空白啊……)。初等数论和其他离散数学分支一样，在计算机科学、信息安全(密码学)、离散控制系统和代数编码等许多领域有着广泛的应用。
这篇文章，将介绍关于数在整除性方面的一些基础的问题。
有关整除的一个简单问题。(Problem source : nefu 115)
数理分析：既然是关于斐波那契数列数的整除问题，那么我们就将斐波那契数列给写出来。
1，1，2，3，5，8，13，21，34 , 55，89 ， 144 ……
我们找3的倍数：n = 4 , 8，12……
我们找4的倍数：n = 6 ， 12……
于是我们做出猜想：当n是4的倍数，f[n]是3的倍数；当n是6的倍数，f[n]是4的倍数。
那么基于以上的猜想，当n既是4的倍数，又是6的倍数(即n是12的倍数，12为4、6的最小公倍数)，此时f[n]就既是3的倍数，也是4的倍数(即是12的倍数，12为3、4的最小公倍数)。
之所以要基于这样的猜想，是因为数据量太大，一一计算是不太可能的，于是我们就猜想对于整个斐波那契数列，3和4的倍数是有规律的分布的，即形成了循环节。
有了以上猜想，我们就很好来编程实现了。
简单的参考代码如下。
#include<stdio.h>
#include<iostream>
using namespace std;
int main()
{
int n;
while(cin >> n)
{
if(n%12 == 0) cout<<"YES"<<endl;
else if(n%4 ==0) cout<<'3'<<endl;
else if(n%6 ==0) cout<<'4'<<endl;
else cout<<"NO"<<endl;
}
}
通过这个简单的问题，我们其实就可以初步窥探数论这个领域的特点，它既简单又复杂。拿上面的例子，我们如果打出很多个斐波那契数，会发现都是符合我们做出的猜想的，但是归根结底，我们的猜想还是一个不完全归纳，如果想进行严谨的证明，问题的难度就大大提升了。这也许就是高等数论中那个著名的“黎曼猜想”经久以来还没有被数学家攻克的原因吧？也许不是？罢了，这也是猜想，哈哈。
通过这个问题我们可以看到，在解决数论问题以及在其他科学领域，有猜想的能力是非常重要的。猜想本身是对规律的一种猜测和概括，从而会打开一扇解决问题的大门。
我们再来看一道类似的问题。(Problem source : codeforces 248B)
题目大意：让你找到长度为n的最小整数，这个整数能够同时整除2,3,5,7。
数理分析：如果我们先忽视数据的大小来看这个问题，是很简单的。既然需要同时整除2,3,5,7，我们就计算它们的最小公倍数lcm(2,3,5,7) = 210。能整除210的最小的n位整数即为所求。这里编程上也很容易实现。
但是我们注意到n的取值范围是[1,10^5]，这么大的数据用整型保存显然不太现实，我们就需要考虑用字符了。我们先通过长整型的数据类型进行打表，观察到结果呈现出一定的规律分布。结果如下。
此时我们可以看到，除了n = 3，后面的结果的形式都可看做“1” + “一些0” + “周期变化的数” ， 而“一些零”我们会发现是呈递增分布的(这里将第二部分和第三部分交界处的0放入那一部分能够得到规律，需要一定的洞察力)，其正比于n-4。我们也能够找到这些周期分布的数给找出来，分别是05、08、17、05、20、11。基于此，我们就可以轻松得编程实现了。
参考代码如下。
#include<stdio.h> #include<cmath> using namespace std; int main() { char a[6][3] = {"05","08","17","02","20","11"}; int n; while(scanf("%d",&n) != EOF) { if(n <= 2) printf("-1\n"); else if (n == 3) printf("210\n"); else { int temp = (n - 4)%6; printf("1"); for(int i = 0;i < n - 4;i++) printf("0"); printf("%s0\n",a[temp]); } } }
通过上面两个问题，我们不难总结出来，在数论的一些问题中，虽然穷举往往能够解决一些看似简单的问题，但是随着数据的增大这个方法将不好用，所以我们必须通过归纳观察找到规律、线性的方法以降低程序的时间复杂度，以此来结局大数据的问题。
在整除性问题中，有两个很重要的概念——最小公倍数(用lcm(a,b)表示)和最大公约数(用gcd(a,b)表示)。关于这两者的定义不再累述，根据其名字也可以推断出其含义。我们今天就来简单的探讨一下这两者的计算方法和两者之间的关系。
最大公约数和最小公倍数之间的关系：lcm(a,b) = ab/gcd(a,b)。这个关系作为基本性质给出，证明也不难。如果将a、b分别表示成含有a、b最大公约数的乘积形式。a = m*gcd(a,b) , b = n*gcd(a,b) ，那么m和n之间显然是互质的，两个互质数的最小公倍数肯定是两个数的乘积，因此m、n的最小公倍数是m*n，因此，a和b的最小公倍数是m*n*gcd(a,b)，即ab/gcd(a,b)。
最大公约数和最小公倍数的计算：基于二者的关系式，我们只需要知道一个量，就能计算出另一个量。这里我们再探讨最大公约数的计算方法。
这里我们先给出这样一个证明：证a = qb + r中，a、b的公因子和b、r的公因子相同。
我们假设d|a 且 d|b。我们看到r = a - qb ，假如c = xa + yb，且d|a，d|b，那么d|c。这里是同样道理，所以d|r，加上之前的d|b，证明完毕。
所以我们求gcd(a,b)，可以转化成求gcd(b,a%b)，然后再转化成gcd(a%b,b%(a%b))，这样一直递推下去，最终会出现gcd(a,0)的情况，显然此时返回a即可。
基于上文的理论，我们通过简单的题目来实践一下。(Problem source : hdu 1108)
很直接的求最小公倍数，参考代码如下。
#include<cstdlib> #include<iostream> using namespace std; long long gcd(int a , int b) { if(b == 0) return a; else return gcd(b , a%b); } int main() { int m , n; while(cin >> m >> n) { cout << (m/gcd(m,n))*n<<endl; } return 0; }
我们再来看一道关于最小公倍数的问题。(Problem source : hdu 1019)
题目大意：求出n个整数的最小公因子。 数理分析：基于上面关于两个数最大公约数和最小公倍数的关系，我们很容易得到两个数的最小公因子，那么三个数呢？显然我们可以先求出两个的最小公因子lcm1，然后在求lcm1和第三个数的最小公因子lcm2，则lcm2就是这三个数的最小公因子。这种做法是显然正确的，无需证明。那么推广起来，n个数我们也采取两两计算来得到n个数的最小公因子。
参考代码如下。
#include<stdio.h> #include<iostream> using namespace std; long long gcd(int a , int b) { if(b == 0) return a; else return gcd(b , a%b); } int main() { int ncases; int n,i; long long lcm , temp; scanf("%d",&ncases); while(ncases--) { scanf("%d",&n); lcm = 1; for(i = 0;i < n;i++) { scanf("%I64d",&temp); lcm = lcm/gcd(lcm,temp)*temp; } printf("%I64d\n",lcm); } }
基于对n个整数求最小公倍数lcm算法的编程实现，我们再来看一道类似的问题。(Problem source : hdu 1788)
但是我们仔细来看这道题目，会发现其实中国剩余定理只是一个混淆视听的存在。
题目要求是满足一下等式的最小整数N
N%M1 = M1 - a
N%M2 = M2 - a
……
N%MI = MI - a
我们通过求余式转化，即可得到下面的式子。
(N + a)%M1 = 0
(N + a)%M2 = 0
……
(N + a)%MI = 0，显然这里的N+a其实就是M1、M2、M3……MI的最小公倍数，求N也就不是难事了。
这道题目告诉我们，在遇到解一元线性求余方程组的时候，虽然中国剩余定理是一种方法，但是在某种特殊情况下(比如这道题目)，我们是可以将问题简化处理的。
参考代码如下。
#include<stdio.h> #include<iostream> using namespace std; long long gcd(int a , int b) { if(b == 0) return a; else return gcd(b , a%b); } int main() { int m,i; int I,a; long long lcm , temp; while(scanf("%d%d",&I,&a) != EOF) { if(I == 0 && a == 0) break; lcm = 1; for(i = 0;i < I;i++) { scanf("%I64d",&temp); lcm = lcm/gcd(lcm,temp)*temp; } printf("%I64d\n",lcm - a);
我们再来看一道有关最大公约数和最小公倍数的题目。(Problem source : nefu116)
这里不过是结合的具体的情景，我们不难发现，两种剑法同时变化周期结束的时候，其实就是周期M，N的公倍数。而结合具体的情景，一出现两仪剑法就表示妖人死亡，所以我们要找的肯定是最小公倍数。
有了这层对具体情境的抽象，编程上就不难实现了，代码也与上题给出的一致。
#include<cstdlib>
#include<iostream>
using namespace std;
long long gcd(int a , int b)
{
if(b == 0)
return a;
else return gcd(b , a%b);
}
int main()
{
int m , n;
while(cin >> m >> n)
{
cout << (m/gcd(m,n))*n<<endl;
}
return 0;
}
基于在欧几里得算法，我们这里将简单介绍拓展欧几里得算法以及利用它能够解决的求解不定方程、求解同余方程、求解模的逆元等问题。
先给出拓展欧几里得定理：对于整数a，b，存在整数对x，y，是的gcd(a,b) = xa + yb。
①显然，当b = 0的时候，gcd(a,0) = 1*a + 0*b，x = 1 ， y = 0；
②而当ab ！= 0的时候，我们先假设gcd(a,b) = x1a + y1b。
同时有gcd(b,a % b) = x2b + y2(a%b)。
通过欧几里得算法的证明我们可以知道，gcd(a,b) = gcd(b,a % b)，由此不难推出——x1a + y1b = x2b + y2(a - (a/b)*b)。
进行展开化简，我们得到：x1a + y1b = y2a +y2(x2 - (a/b))b。
所以，x1 = y2 , y1 = (x2 - (a/b))y2。
那么我们可以看到，作为我们想要的结果，x1和y1的值是基于x2，y2的。如果重复上述过程，会发现x2，y2的值又是基于x3，y3的，这样便形成了递推，而最终这个过程会遇到①的情况，此时xi = 1 ，yi = 0，然后开始递归，最终能够得到x1和y1的值。
我们来结合一个问题来具体分析一下。(Problem source : pku 1061)
Description
Input
Output
数理分析：我们设步数为t，那么A青蛙t步后的坐标为x+mt，B青蛙为y+nt，由于是在地球上运动，所以我们迁移到数轴上，只要A、B的坐标之差是L的整数倍，这就表明它们在球体的同一位置，即相遇。
即，求解x + mt - y - nt = l*p(p为整数)。这个方程的t可取的最小值。
我们先对方程进行一些处理——(n - m)t + l*p = (x - y)。此时我们方程两边同mod l ，等式依然成立，方程变为(n - m)t = (x - y) (mod l)，这就是所谓的同余方程。
我们暂且回到(n - m)t + l*p = (x - y)这个方程，与我们的拓展欧几里得算法给出的式子gcd(a , b) = xa + yb进行联系思考，我们就会有一些启发。前者(n - m) 和 l已知，而t、p未知，我们想要求最小的t；后者a、b已知，我们通过拓展欧几里得来计算x、y。
再看等式的另一边，前者是(x - y),后者是gcd(a,b)这个差异显得就比较大。容易看出，如果(x-y) % gcd(a,b) = 0，那么t和p依然是有解，无非是在gcd(n - m, l)的基础上乘以(x-y)/gcd(n-m,l)，但是如果(x - y)%gcd(a,b)不为0呢？那么此时t、p是无解的，也就是对应着"Impossible"的情况了。
下面我么您详细分析一下(x - y)%gcd(n - m ,l) = 0的情况。这里为了表述方便，我们令M = gcd(num , l)。那么我们通过拓展欧几里得算法将会得到M = x(n - m) + yl中的x和y。
M = x(n - m) + yl ①
(x - y) = t(n - m) + pl ②
M = gcd(n-m , l) ③
通过这等式的联立观察，我们发现，t = x * (x - y)/M。但是还没有完，我们需要的是t的最小值，这里p又是不确定的，可能导致我们最终求出的t是扩大后的结果。
此时我们有了方程 (x - y) = t(n - m) + pl 一个解，但是我们需要得到该方程解集中的最小解，这该如何处理呢？
我们将上段的方程转化成同余方程：(x-y) = t(n-m) (mod l)，暂且放在这里，我们去研究另外一个简单的模型来找到同余方程解空间的分布规律。
给出一个同余方程 5x=1(mod3)， 解得x = 2,5,8,11,14......我们发现解空间有着相同的间隔——3，基于此，我们再更抽象的看问题。给出同余方程ax = b (mod n) ④
设解之间的间隔为dx，那么有a(x + dx) = b (mod n)⑤
④-⑤，得a*dx = 0(mod n)，那么显然a*dx是a和n的公倍数，这里我们当然要求最小的dx，所以计算a，n的最小公倍数——a*n/gcd(a,n)，由此我们看到，dx = n / gcd(a,n)。那么我们只要知道了解空间中的任意解x，只需计算 (x%dx + dx)%dx即可得到最正小解。
有了这个模型做基础，我们再回到我们的问题上来。对于同余方程(x-y) = t(n-m) (mod l)，dx = l/gcd(n-m,l)，那么(t%dx + dx )%dx就是我们所求的最小解。
编程实现：基于前文对于拓展欧几里得算法原理的引入，提到了递归的构造方法。这里需要注意的是gcd(a,b) = xa + yb中x与a、y与b在函数中要对应起来，否则可能会出现错误。
参考代码如下。
#include<cstdlib>
#include<iostream>
typedef long long int64;
using namespace std;
int64 exgcd(int64 a,int64 &x,int64 b,int64 &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
int64 r = exgcd(b,x,a%b,y);
int64 t = x;
x = y;
y = t-a/b*y;
return r;
}
int main()
{
int64 x,y,m,n,l;
int64 ar , br ,cr;
while(cin >> x >> y >> m >> n >> l)
{
int64 M = exgcd(n-m,ar,l,br);
if((x-y)%M || m == n)
cout << "Impossible" << endl;
else
{
int64 s = l/M;
ar = ar*((x-y)/M);
ar = (ar%s + s)%s;
cout << ar << endl;
}
}
}
通过上文介绍利用拓展欧几里得解决同余问题，我们在这里再来看一道有关同余方程的问题。(pku 2115)
Description
for (variable = A; variable != B; variable += C) statement;I.e., a loop which starts by setting variable to value A and while variable is not equal to B, repeats statement followed by increasing the variable by C. We want to know how many times does the statement get executed for particular values of A, B and C, assuming that all arithmetics is calculated in a k-bit unsigned integer type (with values 0 <= x < 2k) modulo 2k.
Input
The input is finished by a line containing four zeros.
Output
题目大意：我们熟知程序语言中的for循环语句，那么现在给出A，B，C，k的值，表示k进制数下，循环控制变量i从A开始，每次+C，到B结束，求你求解该循环执行的次数。
数理分析:值得注意的是，这里的k代表k位系统，即系统中的无符号整型的范围是[0,2^k]，如果某次循环超过2^k,则会从0开始，这与后面给出的A，B，C的范围也是呼应起来的，也将为问题的抽象做铺垫。
也正是基于这个限制条件，和储存系统的特性(超过最大值从0开始的循环性的机制)，我们很容易将其与取余相联系，我们能够将这个实际问题转化成抽象化的数论方程，即A + BX = C (mod M) ,其中2^k = M ①
列出了同余方程，我们就容易利用扩展欧几里得算法来对齐进行求解了。
对方程①进行等价转换，BX + MY = C - A , 首先我们求解BX + MY = gcd(B,M)的解，如果gcd(B,M)|(C - B),那么方程有解，否则无解。
类似pku1061，基于循环语句自身的特点，我们显然要求出x的最小值，而对于如何求出最小值，笔者在《数论及其应用——整除性问题》已经给出论述，这里不再累述。
有了以上数理分析，在编程中只需调用扩展欧几里得算法即可轻松实现。
参考代码如下。
#include<cstdio> #include<iostream> using namespace std; void exgcd(long long a,long long b , long long &d , long long &x,long long &y) //gcd(a,b) = ax +by , d = gcd(a,b) { if(b==0) { x = 1;y =0; d = a ; return ; } else { exgcd(b,a%b,d,x,y); long long temp = x; x = y; y = temp - (a/b)*y; } } int main() { long long x ,y , m , n , l; long long a , b , c , d, k; while(cin >> a >> b >> c >> k &&(a + b + c + k)) { long long temp = c; c = b - a; a = temp; b = (long long)1 << k; exgcd(a,b,d,x,y); if(c%d !=0) { cout << "FOREVER\n"; } else { long long ans = x * c/d; //最小解 long long temp = b/d; ans = ans % temp +temp; cout << ans % temp << endl; } } return 0; }
我们再来看一道有关整除的问题。(hdu 2197)
数理分析：对于长度为n的本原串，我们实在是不好找，而长度为n能够组成的所有情况我们是知道的——2^n，那么我们能够通过间接的找长度为n的非本原串来解决问题呢？设f(n)表示长度为n的本原串，g(n)表示长度为n的非本原串，则有f(n) = 2^n - g(n)。
那么我们下面来讨论g(n)的计算，基于题设关于非本原串的定义，我们知道每个非本原串都有一个“基元单位”，而这个非本原串其实就是由某个“基元单位”循环k次形成的，而k显然是整数，也就是说这个基元单位的长度ai(ai<n)是n的因子，我们尝试将用动态规划的思维将这个全局问题给子问题化，也就是说，g(n) = f(a1) + f(a2) + ……，其中a1、a2、a3……是n的非1和n的因子序列。
那么我们可以得到一个递归程式，如下。
f(n) = 2^n - f(a1) - f(a2) ……-f(ai)。
而该方程的初始情况很显然——f[1] = 2。
基于对整个计算过程中递归程式的给出，在具体计算中我们只需再简单的调用快速进行乘方运算的快速幂和试除法求解因子即可。
参考代码如下。
#include<iostream> #include<cstdio> #include<cstring> using namespace std; const int mod = 2008; long long quick_mod(long long a,long long b,long long m) { long long ans = 1; while(b) { if(b&1) { ans = (ans * a) % m; b--; } b /= 2; a = a * a % m; } return ans; } int get(int n) { if(n == 1) return 2; int sum = 0; for(int i = 2;i*i <= n;++i) { if(n%i == 0) { sum = (sum +get(i))%mod; if(n/i != i) sum = (sum + get(n/i))%mod; } } return (quick_mod(2,n,mod) - 2 - sum + mod)%mod; } int main() { int n; while(cin >> n) cout<<get(n)<<endl; return 0; }
可以看到，这道描述简单的问题，其实融合了数论中快速幂、试除法和基本的动态规划或者说递归思想，是一道很综合的好题目。
——<未完>
