模乘逆元与孙子定理

孙子定理也称为中国剩余定理。

《孙子算经》卷下第二十六题（“物不知数”问题）：有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

孙子定理讲的是求解一元线性同余方程组的方法。

有人指出，孙子定理有以下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