大数有理数取模与扩展欧几里得算法实践

马迪姐

1. 有理数取模问题概述

在算法竞赛和密码学应用中，我们经常需要处理大数的模运算问题。有理数取模是一个经典问题，它要求我们将分数形式的有理数转换为模意义下的整数表示。具体来说，给定一个有理数c=a/b，我们需要找到整数x，使得bx≡a(mod p)，其中p是一个给定的质数（在本题中p=19260817）。

这个问题看似简单，但当a和b都是超大整数（比如长达10001位）时，传统的计算方法就会遇到瓶颈。我们需要一种既能处理大数运算，又能高效计算模逆元的解决方案。

2. 核心算法原理

2.1 模逆元的概念

模逆元是解决有理数取模问题的关键。在模p运算中，如果a和p互质，那么存在唯一的整数x（1≤x<p），使得a×x≡1(mod p)。这个x就是a在模p下的逆元，记作a⁻¹。

对于有理数a/b mod p，可以转化为a×b⁻¹ mod p，其中b⁻¹是b的模逆元。因此，问题的核心转化为如何高效计算b的模逆元。

2.2 扩展欧几里得算法

扩展欧几里得算法不仅能计算两个数的最大公约数，还能找到满足贝祖等式ax+by=gcd(a,b)的整数解x和y。当a和b互质时，这个等式就变成了ax+by=1，此时x就是a在模b下的逆元。

算法的时间复杂度是O(log min(a,b))，非常高效。这使得它成为计算模逆元的首选方法，尤其适合处理大数运算。

3. 实现细节解析

3.1 大数处理技巧

题目中a和b的长度可能达到10001位，直接处理这样的数字是不现实的。我们需要利用模运算的性质：

(a + b) mod p = [(a mod p) + (b mod p)] mod p
(a × b) mod p = [(a mod p) × (b mod p)] mod p

基于这个性质，我们可以边读入数字边取模：

cpp复制auto ToInt = [&](const char* s) {
    int ret = 0;
    for (int i = 0; i < strlen(s); i++) {
        ret = (ret * 10 + s[i] - '0') % m_iMod;
    }
    return ret;
};

这种方法将O(n)的大数运算转化为O(1)的整数运算，极大提高了效率。

3.2 扩展欧几里得实现

我们实现了一个Caxby1类来封装扩展欧几里得算法：

cpp复制class Caxby1 {
public:
    Caxby1(int a, int b) :G(m_iG) {
        m_iSignA = (a < 0) ? -1 : 1;
        m_iSignB = (b < 0) ? -1 : 1;
        m_a = abs(a);
        m_b = abs(b);
        m_iG = gcd(m_a, m_b);
    }
    
    tuple<int, int> Any() {
        if (0 == m_iG) { return { 0,0 }; }
        auto [x, y] = Any(m_a / m_iG, m_b / m_iG);
        return { m_iSignA * x,m_iSignB * y };
    }
    
    const int& G;
private:
    static pair<int, int> Any(long long a, long long b) {
        if (1 == a) { return { 1,0 }; }
        if (1 == b) { return { 0,1 }; }
        if (a > b) {
            auto c = a / b, d = a % b;
            if (0 == d) { d = b; c--; }
            auto [x1, y1] = Any(d, b);
            return { x1, y1 - c * x1 };
        }
        auto [x2, y2] = Any(b, a);
        return { y2,x2 };
    }

private:
    int m_a, m_b, m_iG, m_iSignA, m_iSignB;
};

这个实现有几个关键点：

处理了负数情况，保证内部计算使用正数
递归实现扩展欧几里得算法
返回解的符号与原始输入一致

3.3 解的存在性判断

在计算bx≡a(mod p)时，解存在的条件是gcd(b,p)能整除a。我们的实现中通过检查gcd(ab.G, ia) == ab.G来判断：

cpp复制if (gcd(ab.G, ia) != ab.G) { return -1; }

如果条件不满足，直接返回-1表示无解。

4. 完整解决方案

将上述组件组合起来，我们得到完整的解决方案：

cpp复制class Solution {
public:
    int Ans(const char* a, const char* b) {
        auto ToInt = [&](const char* s) {
            int ret = 0;
            for (int i = 0; i < strlen(s); i++) {
                ret = (ret * 10 + s[i] - '0') % m_iMod;
            }
            return ret;
        };
        
        int ia = ToInt(a);
        int ib = ToInt(b);        
        Caxby1 ab(ib, m_iMod);
        
        if (gcd(ab.G, ia) != ab.G) { return -1; }
        
        auto [x1,y1] = ab.Any();
        x1 = ((x1 % m_iMod) + m_iMod) % m_iMod;
        return (ia * (long long)x1) % m_iMod;
    }
    
    const int m_iMod = 19260817;
};

5. 关键注意事项

大数处理：必须边读入边取模，否则会溢出。即使使用long long类型，也无法直接存储10001位数。
解的存在性：不是所有有理数都能在模意义下表示。当gcd(b,p)不能整除a时，问题无解，应返回"Angry!"。
负数的处理：计算得到的逆元可能是负数，需要通过加模数再取模的方式转换为正数：
```
cpp复制x1 = ((x1 % m_iMod) + m_iMod) % m_iMod;
```
中间结果溢出：在计算ia * x1时，即使ia和x1都在int范围内，乘积也可能溢出，因此需要先转换为long long：
```
cpp复制return (ia * (long long)x1) % m_iMod;
```