孙子定理也称为中国剩余定理。
《孙子算经》卷下第二十六题(“物不知数”问题):有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
孙子定理讲的是求解一元线性同余方程组的方法。
有人指出,孙子定理有以下5种解法:
1.枚举法
2.解不定方程法
3.逐级满足法
4.化为相同除数的同余式法
5.经典同余式方程组解法
据说最为简洁快速的是第4种方法。
这里介绍的程序为第5种解法。
要得到解x,首先需要计算模乘逆元。模乘逆元定义为:
满足 ab≡1(mod m),称b为a模乘逆元。
程序中给出了两种求模乘逆元的方法,一种是穷举法,此法对于模值较大时是一种低效率的方法;另外一种是利用扩展欧几里得算法求逆元。
求得模乘逆元之后,就可以算出同余方程组的解。
“物不知数”问题的一个解为23(最小值解),其解系为x=23+105k(k>=0)(105=3*5*7)。
#include <stdio.h>
// 递推法实现扩展欧几里德算法
long exgcd(long a, long b, long *x, long *y)
{
long x0=1, y0=0, x1=0, y1=1;
long r, q;
*x=0;
*y=1;
r = a % b;
q = (a - r) / b;
while(r)
{
*x = x0 - q * x1;
*y = y0 - q * y1;
x0 = x1;
y0 = y1;
x1 = *x;
y1 = *y;
a = b;
b = r;
r = a % b;
q = (a - r) / b;
}
return b;
}
// 扩展欧几里德算法求逆元
long minv(long a, long p)
{
long x, y;
exgcd(a, p, &x, &y);
return x<0 ? x+p : x;
}
// 试探法求逆元
long minv2(long a, long p)
{
long y=1, t;
int i;
if(a < 0) {
a = a % p;
a += p;
}
for(i=1; i<p; i++) {
t = a * i;
if(t % p == 1) {
y = i;
return y;
}
}
return y;
}
int main(void)
{
printf("a=%d m=%d x=%ld x=%ld\n", 65, 83, minv(65, 83), minv2(65, 83));
printf("a=%d m=%d x=%ld x=%ld\n", 11663, 103, minv(11663, 103), minv2(11663, 103));
printf("a=%d m=%d x=%ld x=%ld\n", 11227, 107, minv(11227, 107), minv2(11227, 107));
printf("a=%d m=%d x=%ld x=%ld\n", 11021, 103, minv(11021, 109), minv2(11021, 109));
// 孙子算经
long a[] = {2, 3, 2};
long m[] = {3, 5, 7};
int i, size= sizeof(a)/sizeof(long);
long bm=1, subm[size], x=0;
for(i=0; i<size; i++)
bm *= m[i];
for(i=0; i<size; i++)
subm[i] = bm / m[i];
for(i=0; i<size; i++) {
x += subm[i] * minv(subm[i], m[i]) * a[i];
x %= bm;
}
printf("x=%ld\n", x);
return 0;
}
程序运行结果:
a=65 m=83 x=23 x=23
a=11663 m=103 x=73 x=73
a=11227 m=107 x=40 x=40
a=11021 m=103 x=100 x=100
x=23