離散対数問題の準指数時間解法

# 離散対数問題
離散対数問題とは、DiffieHellman鍵交換方式で使われる理論で
$p$を素数、$g$を原始根とする時の、巡回群 $Z_p^\times$ の元$A$における以下を満たす指数$a$を求める問題。

```math
A=g^a \bmod p
```

# 計算量
$g$が原始根であるから巡回群 $Z_p^\times$ をすべて生成するのには$g$を$p-1$乗すれば良い。よって離散対数問題は、$O(p)$で解ける。
計算量は、$p$をbitで表現した時の長さに関する指数時間となるが、実際はもう少し緩く、準指数時間と呼ばれる計算量クラスに分類されている。 評価式は以下の資料を参照。
https://www.journal.ieice.org/conts/kaishi_wadainokiji/2008/200806.pdf
従って、$p$が1024bitなどで表現されると、高速には解けないところが暗号に利用できる根拠となる。

# 準指数時間での解法

素数を集め基底とする。ここでは4個を指定した。

```c++
vector<int>base;
base.push_back(3);
base.push_back(5);
base.push_back(7);
base.push_back(11);
```

$k$を適当に選び$g^k\bmod p$を上記の基底で素因数分解をする

```c++
for(int i=0;i<modp;i++){
   int tmp_g_kai=g_kai;
   vector<int>coefficient;
   for(int j=0;j<basesize;j++){
      int cnt=0;
      while(tmp_g_kai%base[j]==0){
         tmp_g_kai/=base[j];
         cnt++;
      }
      coefficient.push_back(cnt);
   }
   if(tmp_g_kai==1){
      //ok
      //処理内容は下で記載
   }
   //・・・
   g_kai*=g;
   g_kai%=modp;
}
```

$g^k\bmod p$に関する素因数分解の式の両辺の$log$をとる。指数に関しては、$Z_p^\times$の位数と同じ$\bmod p-1$となることに注意。

```math
g^k\bmod p=\prod_{q}q^{a_q}　ただし、qは素数を集めた基底\\
k=\sum_{q}a_q log_g(q) \bmod p-1
```
これを満たす$k$を、選んだ基底の個数分収集することで連立方程式を作り、これを解く。

```math
\left(
\begin{array}{ccc}
k_0\\
k_1\\
k_2\\
k_3
\end{array}
\right)
=
\left(
\begin{array}{ccc}
a_{q_{0_0}} & a_{q_{0_1}} & a_{q_{0_2}} & a_{q_{0_3}} \\
a_{q_{1_0}} & a_{q_{1_1}} & a_{q_{1_2}} & a_{q_{1_3}} \\
a_{q_{2_0}} & a_{q_{2_1}} & a_{q_{2_2}} & a_{q_{2_3}} \\
a_{q_{3_0}} & a_{q_{3_1}} & a_{q_{3_2}} & a_{q_{3_3}} \\
\end{array}
\right)

\left(
\begin{array}{ccc}
log (q_0)\\
log (q_1)\\
log (q_2)\\
log (q_3)
\end{array}
\right)
```

```c++
if(tmp_g_kai==1){
   for(int j=0;j<basesize;j++){
      a[basesize*j+n]=coefficient[j];
   }
   b[n]=i;
   n++;
   //・・・
}
```

連立方程式をとくには、IntelのMKLに入っているdgesv_()を用いるのが簡単。

```c++
#include <Accelerate/Accelerate.h>
//・・・
dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info);
```

最後に$g^r A\bmod p$を上記の基底で素因数分解して両辺の$log$をとると

```math
g^rA \bmod p=\prod_{q}q^{a_q}\\
r+a=\sum_{q}a_q log_g(q) \bmod p-1\\
a=\sum_{q}a_q log_g(q)-r \bmod p-1\\
```
$a$が求まる。

```c++
for(int i=0;i<modp;i++){
   int tmp_A=A;
   vector<int>coefficient;
   for(int j=0;j<basesize;j++){
      int cnt=0;
      while(tmp_A%base[j]==0){
         tmp_A/=base[j];
         cnt++;
      }
      coefficient.push_back(cnt);
   }
   if(tmp_A==1){
      //ok
      //処理内容は下で記載
   }
    A*=g;
    A%=modp;
}
```

dgesv_()では、連立方程式の解は配列bに入っているので

```c++
if(tmp_A==1){
   int ans=0;
   for(int j=0;j<basesize;j++){
      ans+=coefficient[j]*b[j];
      ans%=modp-1;
   }
   ans-=i;
   ans%=modp-1;
   cout << ans << endl;
   break;
}
```

全体のコードは以下

```c++
#include <iostream>
#include <Accelerate/Accelerate.h>
#include <vector>
using namespace std;

int main(int argc, const char * argv[]) {
    int modp=107;
    int g=2;
    vector<int>base;
    //base.push_back(2);
    base.push_back(3);
    base.push_back(5);
    base.push_back(7);
    base.push_back(11);
    
    int basesize=(int)base.size();
    //dgesv
    int n=0;
    int nrhs=1;
    int lda=basesize;
    int ldb=basesize;
    int info;
    int ipiv[4/*basesize*/];//Local arrays
    double a[16/*lda*basesize*/];
    double b[4/*ldb*nrhs*/];
    int g_kai=pow(2,10);
    g_kai%=modp;
    for(int i=0;i<modp;i++){
        int tmp_g_kai=g_kai;
        vector<int>coefficient;
        for(int j=0;j<basesize;j++){
            int cnt=0;
            while(tmp_g_kai%base[j]==0){
                tmp_g_kai/=base[j];
                cnt++;
            }
            coefficient.push_back(cnt);
        }
        
        if(tmp_g_kai==1){
            //rankの確認
            for(int j=0;j<basesize;j++){
                a[basesize*j+n]=coefficient[j];
            }
            b[n]=10+i;
            n++;
            if(n==basesize)break;
        }
        
        g_kai*=g;
        g_kai%=modp;
    }
    /* Solve the equations A*X = B */
    if(n==basesize){
        dgesv_(&n,&nrhs,a,&lda,ipiv,b,&ldb,&info);
    }
    int A=96;
    A%=modp-1;
    for(int i=0;i<modp;i++){
        int tmp_A=A;
        vector<int>coefficient;
        for(int j=0;j<basesize;j++){
            int cnt=0;
            while(tmp_A%base[j]==0){
                tmp_A/=base[j];
                cnt++;
            }
            coefficient.push_back(cnt);
        }
        
        if(tmp_A==1){
            int ans=0;
            for(int j=0;j<basesize;j++){
                ans+=coefficient[j]*b[j];
                ans%=modp-1;
            }
            ans-=i;
            ans%=modp-1;
            cout << ans << endl;
            break;
        }
        A*=g;
        A%=modp;
    }
    
    return 0;
}
```

dgesv_()は行列に対応する配列の使い方に癖があるので注意。
$p$を全部を回していないので、$p$をbitで表現した時の長さに関する指数時間より小さいはず。