素数を数えて落ち着くんだ･･･

# はじめに
2024年10月21日に「最大の素数発見」の記事が流れました。
メルセンヌ素数を候補として、スパコンで探したらしいです。発見したメルセンヌ素数は、
```math
2^{136279841} - 1
```
約4100万桁

　素数は、近年の通信の暗号化などでよく引き合いに出され、暗号解読には計算量が指数関数的に増える特徴を踏まえています。

円周率や素数は計算例題の定番なので、調べて、一般的なパソコンで動かしてみます。

大きな値は扱いにくいので、10桁程度の計算をしてみて、巨大素数の計算がいかに大変であるかを感じてみます。

# 前知識
円周率にも言えることだが、複数の計算方法があり、計算後に本当に正しいかどうかの検証が必要となる。
素数計算には、判定式というものがあり、確実に判定できる方法と、確率的に判定する方法がある。


# 環境
Mac OSX 13.4.1
Apple M1 Pro
Processing 4.2

# 総当たり
```java:all
void setup() {
  for (int n=2; n<100000; n++) {
    if (prime(n)) {
      //println(n);
    }
  }
  println(millis());
}

boolean prime(int n) {
  int i = 2;
  while (i < n) {
    if (n%i == 0) {
      return false;
    }
    i++;
  }
  return true;
}
```

# 奇数のみ
```java:odd
void setup() {
  for (int n=3; n<1000000; n++) {
    if (prime(n)) {
      //println(n);
    }
  }
  println(millis());
}

boolean prime(int n) {
  if(n%2==0){
    return false;
  }
    
  int i = 3;
  while (i < n) {
    if (n%i == 0) {
      return false;
    }
    i+=2;
  }
  return true;
}
```

# sqrt(n)
```java:sqrt
void setup() {
  for (int n=3; n<100000000; n++) {
    if (prime(n)) {
      //println(n);
    }
  }
  println(millis());
}

boolean prime(int n) {
  if(n%2==0)
    return false;
  
  int i = 3;
  while (i <= sqrt(n)) {
    if (n%i == 0)
      return false;
    i+=2;
  }
  return true;
}
```
# wheel factorization
```java:wheel
void setup() {
  for (int n=7; n<10000000; n++) {
    if (prime(n)) {
      //println(n);
    }
  }
  println(millis());
}

boolean prime(int n) {
  if(n%2==0 || n%3==0 || n%5==0)
    return false;
  
  int i = 7;
  int c = 0;
  int [] wheel = {4,2,4,2,4,6,2,6}; 
  while (i < sqrt(n)+1) {
    if (n%i == 0)
      return false;
    i+=wheel[c%8];
    c+=1;
  }
  return true;
}
```

# ミラーラビン素数判定法
```java
int num=1;

void setup() {
  //println(2);
  for (int n=2; n<10000000; n++) {
    if (prime(n)) {
      //println(n);
      num++;
    }
  }
  println(millis()+":"+num);
}

boolean prime(int n) {
  if (n%2==0)
    return false;

  int d = n - 1;
  int s=0;
  while (d%2==0) {
    d >>= 1;
    s++;
  }

  for (int i=0; i<20; i++) {
    boolean isPP = false;
    int a = int(random(2, n-1));
    int x = modPow(a, d, n);
    if (x==1 || x==n-1) {
      isPP = true;
    }

    for (int j=0; j<s; j++) {
      x=modPow(x, 2, n);
      if (x==n-1) {
        isPP = true;
      }
    }
    if (!isPP)
      return false;
  }
  return true;
}

int modPow(long x, int k, int m) {
  if (k == 0) {
    return 1;
  }
  if (k % 2 == 0) {
    return modPow(x * x % m, k / 2, m);
  } else {
    return (int) (modPow(x, k - 1, m) * x % m);
  }
}
```
## メモ
pythonは`pow(x, 2, n)`と書いて、(x^2)%nが計算でき、かつ、桁数を気にする必要がない。
しかし、processing(java)は桁数を気にする必要があるので、`pow(x,2)%n`と書いたとしても、floatの32bit精度(7桁程度)になってしまう問題がある。`int(pow(x, 2))%n`と書いても、だめ。
本記事ではmodPowは以下のサイトを参考にした。
このアルゴリズムは「繰返し二乗法」というらしい。

https://qiita.com/mojaie/items/e306fa2b335948444455

# PHPのgmp_prob_prime()
PHPにはミラーラビンの関数が用意されていて「素数じゃない,多分素数,素数」を判定してくれる。
PHPの整数の最大値は64bit環境では、2^63-1=9223372036854775807(約19桁)
こりゃ便利だね。

# エラトステネスの篩
```java
int num = 1000000000;
num = (num - 3) / 2;
boolean[] pList = new boolean[num];

for (int i = 0; i < num; i++) {
  pList[i] = true; // 1:素数、0:非素数
}

int count = 1; // 素数2
for (int i = 0; i < num; i++) {
  if (pList[i] == true) {
    int pNum = i + i + 3;
    //println(pNum);
    count++;

    for (int j = i + pNum; j < num; j += pNum) {
      pList[j] = false;
    }
  }
}
println("pNumberCount -> " + count);

println(millis());
```

# 計算時間の比較
| n | all | odd | sqrt | wheel |ミラー|篩|素数個数|
|--:|--:|--:|--:|--:|--:|--:|--:|
| 100000  |   0.666 |  0.518 |       |     |||9592|
| 1000000 | 23.947 | 12.195  | 0.487 |     |||78498|
| 10000000|        |         | 1.250 | 1.057|3.090||664579|
|100000000|        |         |23.457 |17.505|27.625|0.662|5761455|
|1000000000|||||281.006|5.432|50847534|

単位：秒

# 考察
エラトステネスの篩はとても速いが、すべての素数を保持するためメモリを使う。nが10倍になるとメモリも10倍必要。
ミラーラビンはエラトステネスの篩ほどで速くはないが、nが10倍になると時間も10倍で済むので、桁数が増えれば増えるほど有利になる。

また、ミラーラビンはピンポイントで素数か否かを判定するので、分散処理が簡単になるはず。

エラトステネスの篩に必要なメモリは
```math
2^{32}=4294967296bit=536870912Byte=0.5GB
```
なので、javaの配列サイズ制限も、メモリの限界には近く、単純には劇的な改善は望めない。

# processing3と4を比較して、単純計算が速くなってる
allのコードで計算して、
processing3環境で、37.754sec
processing4環境で、23.947sec
1.57倍も速くなった。内部でjava8からjava17へ変わったことで、最適化が進み速くなっているのか？


# 参考
https://qiita.com/ppza53893/items/e0f464340d6f97760cd5

https://en.wikipedia.org/wiki/Mersenne_prime

JavaのBigIntegerにはmodPow()がある

https://2357.aimary.com/total_prime.html

https://qiita.com/hextomino/items/0de35974bb3d167e57bf

https://tools.m-bsys.com/calculators/prime_number.php

https://tools.m-bsys.com/calculators/primality_test_for_large_number.php



# python版 sqrt(n)
```py
def is_prime(n):
    if n < 2:
        return False

    for i in range(2, int(n**0.5) + 1):
        if n % i == 0:
            return False
    return True
    
    
c=0
for (i)in range(100000):
    a = is_prime(i) 
    if(a):
        #print(i)
        c+=1
print(c)#9592,78498
#println(millis())#4516
```
https://paiza.io/ で試すと、10万ならすぐ終わる。100万ならTimeOutで止まる。
processing.py で試すと10万ならすぐ終わる。100万なら4.5秒で終わる。
processing(java)で100万なら0.5秒なので、実行速度は期待できないのかな。


# スーパー素数(super prime)

https://en.wikipedia.org/wiki/Super-prime

素数の数列から、再度、素数番目の数を取り出した数列に、法則性はあるか？・・・ない。

| n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|p(n)|2|3|5|7|11|13|17|19|
|p(p(n))|3|5|11|17|31|41|59|67|

https://oeis.org/A006450