孙子定理也称为中国剩余定理。
《孙子算经》卷下第二十六题(“物不知数”问题):有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二。问物几何?
孙子定理讲的是求解一元线性同余方程组的方法。
有人指出,孙子定理有以下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 ", 65, 83, minv(65, 83), minv2(65, 83)); printf("a=%d m=%d x=%ld x=%ld ", 11663, 103, minv(11663, 103), minv2(11663, 103)); printf("a=%d m=%d x=%ld x=%ld ", 11227, 107, minv(11227, 107), minv2(11227, 107)); printf("a=%d m=%d x=%ld x=%ld ", 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 ", 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