DLP

大数运算

大数运算实现先前的实验报告已经有过，不再赘述。

但是在 Pollard-Rho 中有大数的模 3 运算，在先前的大数库中没有实现。这里再在 Field::Int 类中增加模 3 运算的方法：

int mod3() const {
  __uint128_t res = 0;
  for(int i = F.N - 1;i >= 0; i--){
    res = (res * (static_cast<__uint128_t>(1) << 64) + a[i]) % 3;
  }
  return static_cast<int> (res);
}

大数字符串输入输出处理

本次实验不像先前的实验都是二进制输入输出，而是十进制字符串输入输出。

对于输出，由于输出的离散对数是有限域 \(n\) 以内的，而 \(n\) 在本次实验限定为 64 位以内，所以可以直接使用 C++ 本身自带的 64 位数输出。

对于输入，包含了大数，我在自己的大数库已经有实现，为了更高效的性能，我使用查表的方式实现。我的大数库最高支持 2048 bits 的大数，使用 uint64_t[32] 存储，对应到十进制是 618 位，所以建立的 LookUpTable 可以以 618 个十进制位、1 到 9 这 9 个十进制数作为索引，查询结果就是对应的大数的二进制表示。建表如下：

auto gen_TABLE(){
  array<array<array<uint64_t, 32>, 9>, 617> table = {0};
  uint64_t a[32] = {0};
  uint64_t a1[32] = {0};
  a[0] = 1;
  for(int i = 0; i < 617; ++i){
    for(int j = 0; j < 32; ++j){
      a1[j] = a[j];
    }
    int t = i ** 616 ? 3 : 9;
    for(int j = 0; j < t; ++j){
      uint64_t carry = 0;
      for (int k = 0; k < 32; ++k) {
        table[i][j][k] = a[k];
        uint64_t sum = a[k] + a1[k] + carry;
        if (sum < a[k] || sum < a1[k]) {
          carry = 1;
        } else {
          carry = 0;
        }
        a[k] = sum;
      }
    }
  }
  return table;
}

这里我之所以使用 array<array<array<uint64_t, 32>, 9>, 617> 存表，是因为我原本希望在编译时 constexpr 建表，提升运行效率，但是在 OJ 系统上会导致超出编译内存。

/*
not in constexpr, bcz adding it will exceed memory limit on OJ system.
*/
const auto TABLE = gen_TABLE();

输入就通过查表实现，根据输入的每一个十进制位查刚刚建立的表，查询结果加到最终结果之中：

// return the sum of used limbs
inline int string2limbs(const string& s, uint64_t * a){
  int l = s.length(), k = 1;
  bool flag = l > 19;
  a[0] = stoull(
    s.substr(
      flag ? l - 19 : 0,
      flag ? 19 : l
    )
  );
  for (int i = 19; i < l; ++i) {
    int j = s[l - i - 1] - '0';
    if(j ** 0) continue;
    uint64_t sum, carry = 0;
    k = 0;
    j--;
    while(TABLE[i][j][k] ** 0)k++;
    while(k < 32 && (TABLE[i][j][k] != 0 || carry != 0)){
      sum = a[k] + TABLE[i][j][k] + carry;
      if (sum < a[k] + carry || sum < TABLE[i][j][k] + carry) {
        carry = 1;
      } else {
        carry = 0;
      }
      a[k] = sum;
      k++;
    }
  }
  return k;
}

我的大数库采用面向对象的设计，所以输入也封装到了友元 >> 之中：

// 输入有限域（p）
inline friend std::istream& operator>>(std::istream& is, Field& f) {
  string s;
  is >> s;
  f.N = string2limbs(s, p);
  initField(f.p, f.p_2, f.R2modp, f.g, f.N);  // 有限域预计算处理
  return is;
}

// 输入整数
inline friend std::istream& operator>>(std::istream& is, Int& a) {
  string s;
  is >> s;
  string2limbs(s, a.a);
  return is;
}

Pollard-Rho

Pollard-Rho 的算法我也采用了面向对象的设计，将一个 DLP 问题的条件作为成员变量封装在里面：

• \(p\) 定义在 Field f，同时也是整个 Pollard-Rho 算法大数运算时所基于的有限域。
• uint64_t n 即为 \(n\)，因为输入限制在 64 位以内，所以不必使用大数。
• Int alpha = Int(f) 即为 \(\alpha\)，属于有限域 \(p\)。
• Int beta = Int(f) 即为 \(\beta\)，属于有限域 \(p\)。

class PollardRho {
public:
  PollardRho(){};

private:
  void func(Int const & x, __uint128_t const & a, __uint128_t const & b, Int & x_, __uint128_t & a_, __uint128_t & b_);

  void func(Int & x, __uint128_t & a, __uint128_t & b);

  __uint128_t gcd(__uint128_t a, __uint128_t b);  // 返回 a,b 的 gcd

  __uint128_t mod_inverse(__uint128_t a, __uint128_t p);  // 返回 a^(-1) mod p

  // 值为 1 的大数
  string _One = string("1");
  const Int One = Int(f, _One);

public:
  Field f;
  Int alpha = Int(f);
  Int beta = Int(f);
  uint64_t n;

  uint64_t solve();
}

具体的算法流程如下：

\[ \begin{aligned} &(x1,a1,b1)\gets f(1,0,0)\\ &(x2,a2,b2)\gets f(x1,a1,b1)\\ &\bold{\text{while}}\ x1\neq x2\\ &\ \ \ \ (x1,a1,b1)\gets f(x1,a1,b1)\\ &\ \ \ \ (x2,a2,b2)\gets f(x2,a2,b2)\\ &\ \ \ \ (x2,a2,b2)\gets f(x2,a2,b2)\\ &\text{解线性同余方程:}\\ &\ \ \ \ (b2-b1)x=a1-a2\ \ \ \ (\bmod n) \end{aligned} \]

其中 \(f\) 的计算如下：

\[ f(x,a,b)= \begin{cases} (\beta x ,& a,& (b+1)\bmod n &), && x\bmod3=1\\ (x^2 ,& 2a ,& 2b\bmod n &), && x\bmod3=0\\ (\alpha x ,& (a+1)\bmod n ,& b &), && x\bmod3=2 \end{cases} \]

这里我对输入和输出是否为为同一变量分别作了实现避免 C++ 未定义行为（我真的怕极了，不想在遇到 undefined behavior 了）：

void func(Int const & x, __uint128_t const & a, __uint128_t const & b, Int & x_, __uint128_t & a_, __uint128_t & b_){
    switch(x.mod3()){
        case 1: {
            x_ = beta * x;
            a_ = a;
            b_ = (b + 1) % n;
        } break;
        case 0: {
            x_ = x * x;
            a_ = (a * 2) % n;
            b_ = (b * 2) % n;
        } break;
        default: {  // 2
            x_ = alpha * x;
            a_ = (a + 1) % n;
            b_ = b;
        }
    }
}

void func(Int & x, __uint128_t & a, __uint128_t & b){
    switch(x.mod3()){
        case 1: {
            x = beta * x;
            a = a;
            b = (b + 1) % n;
        } break;
        case 0: {
            x = x * x;
            a = (a * 2) % n;
            b = (b * 2) % n;
        } break;
        default: {  // 2
            x = alpha * x;
            a = (a + 1) % n;
            b = b;
        }
    }
}

然后求最终 \(b1,b2,a1,a2\) 的过程如下：

uint64_t solve(){
    Int x1(f), x2(f);
    __uint128_t a1, b1, a2, b2;
    func(One, 0, 0, x1, a1, b1);
    func(x1, a1, b1, x2, a2, b2);
    while(x1 != x2){
        func(x1, a1, b1);
        func(x2, a2, b2);
        func(x2, a2, b2);
    }
    // 解线性同余方程...
}

最后就是求解线性同余方程了，主要参考 OIWiki 中用逆元求解的思路，这里我的方程如下：

\[ (b2-b1)x=a1-a2\ \ \ \ (\bmod n) \]

当然这里我先要保证 \(b2-b1\) 和 \(a1-a2\) 都是模 \(n\) 下的正数，上述 Pollard-Rho 得到的结果并不保证 \(b2>b1\) 和 \(a1>a2\)，所以稍微处理一下：

while(b2 < b1) b2 += static_cast<__uint128_t>(n);
while(a1 < a2) a1 += static_cast<__uint128_t>(n);

接下来先求 \(g=gcd(b1-b2,n)\)，通过如下函数实现：

__uint128_t gcd(__uint128_t a, __uint128_t b){
    while (b != 0) {
        __uint128_t temp = b;
        b = a % b;
        a = temp;
    }
    return a;
}

当且仅当 \((b2-b1)\bmod g\neq0\) 时是无解的，反之则可以在 \((b2-b1)x=a1-a2\) 两侧同时除以 \(g\)，得到新方程：

\[ b'x'=a'\ \ \ \ (\bmod n')\ \ \ \ \begin{cases} b'=(b2-b1)\div g\\ a'=(a1-a2)\div g\\ n'=n\div g\\ \end{cases} \]

就可以通过求逆得到一个基础解 \(x'=a'\cdot b'^{-1}\bmod n\)。求逆代码实现如下：

__uint128_t mod_inverse(__uint128_t a, __uint128_t p) {
    __uint128_t u = a, v = p, x1 = 1, x2 = 0, q, r, x;
    while(u != 1){
        q = v / u;
        r = v - q * u;
        while(x2 < q * x1) x2 += p;
        x = x2 - q * x1;
        v = u;
        u = r;
        x2 = x1;
        x1 = x;
    }
    return x1 % p;
}

但这个解不一定是 DLP 的解，因为线性同余方程有 \(g=gcd(b1-b2,n)\) 个解，只有一个才是我们想要的，那就只好尝试 \(g\) 次试出最终解了，这 \(g\) 个解通过 \(x'\) 构造如下：

\[ x_i=(x'+i\cdot n')\bmod n\ ,\ \ \ 0\leq i\lt g \]

最终求解原方程代码如下：

if ((b2 - b1) % g != 0){
    // 无解
} else {
    __uint128_t n_ = n / g;
    __uint128_t a_ = (a1 - a2) / g;
    __uint128_t b_ = (b2 - b1) / g;
    __uint128_t b_inv = mod_inverse(b_, n_);
    __uint128_t x_ = (a_ * b_inv) % n_;
    Int t(f);
    for(uint64_t i = 0; i < g; ++i){
        uint64_t tmp = static_cast<uint64_t>(x_);
        t.a[0] = tmp;
        if((alpha ^ t) ** beta){
            return tmp;
        }
        x_ = (x_ + n_) % n;
    }
}
return 0xFFFF'FFFF'FFFF'FFFF;  // 无解，我返回一个特殊的数

我学号的 DLP

相关代码进行计时，输入输出：

#ifndef ONLINE_JUDGE
#include <chrono>
#include <iomanip>
FILE *output, *input;
chrono::time_point<chrono::high_resolution_clock> start_time, end_time;
#endif

int main(){
  cin.tie(0);
  ios::sync_with_stdio(false);
  cout.tie(0);

  #ifndef ONLINE_JUDGE
  output = freopen("output", "w", stdout);
  input = freopen("input", "r", stdin);
  start_time = chrono::high_resolution_clock::now();
  #endif

  PollardRho solution;
  string s;
  cin >> solution.f >> solution.n >> solution.alpha >> solution.beta;
  cout << solution.solve();
  cout.flush();

  #ifndef ONLINE_JUDGE
  output = freopen("time", "w", stdout);
  end_time = chrono::high_resolution_clock::now();
  auto ns = chrono::duration_cast<chrono::nanoseconds>(end_time - start_time);
  auto us = chrono::duration_cast<chrono::microseconds>(end_time - start_time);
  auto ms = chrono::duration_cast<chrono::milliseconds>(end_time - start_time);
  cout << ns.count() << " ns" << endl
  fclose(input);
  fclose(output);
  #endif
}

我的学号 \(22342043\) 所确定的 DLP 问题：

\[ \begin{aligned} p&=3768901521908407201157691198029711972876087647970824596533\\ (&p−1=2×2×23×8783×2419781956425763×192888768642311611×9993115456385501509)\\ n&=9993115456385501509\\ \alpha&=3107382411142271813235322646657672922264748410711464860476^{2×2×23×8783×2419781956425763×192888768642311611×22342043}\\ &=1672705575588793411750559287020484713881877693514202546508\\ \beta&=2120553873612439845419858696451540936395844505496867133711 \end{aligned} \]

输入文件 input 内容如下：

3768901521908407201157691198029711972876087647970824596533
9993115456385501509
1672705575588793411750559287020484713881877693514202546508
2120553873612439845419858696451540936395844505496867133711

我在 Linux 虚拟机环境下计算，配置：

• CPU：Intel(R) Xeon(R) CPU E7-4830 v3 @ 2.10GHz 总共 8 核
• 内存： 16.00 GiB
• OS：Ubuntu 22.04

(fish) 开始计算：

nohup bash ./my_program > output.log 2>&1 & set pid $last_pid; echo $pid > pid

计算结果：

\[ 4823944491583471581 \]

计算时间：