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

はじめに

2024年10月21日に「最大の素数発見」の記事が流れました。
メルセンヌ素数を候補として、スパコンで探したらしいです。発見したメルセンヌ素数は、

2^{136279841} - 1

約4100万桁

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

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

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

前知識

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

環境

Mac OSX 13.4.1
Apple M1 Pro
Processing 4.2

総当たり

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;
}

奇数のみ

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)

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

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;
}

ミラーラビン素数判定法

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は以下のサイトを参考にした。
このアルゴリズムは「繰返し二乗法」というらしい。

PHPのgmp_prob_prime()

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

エラトステネスの篩

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倍で済むので、桁数が増えれば増えるほど有利になる。

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

エラトステネスの篩に必要なメモリは

2^{32}=4294967296bit=536870912Byte=0.5GB

なので、javaの配列サイズ制限も、メモリの限界には近く、単純には劇的な改善は望めない。

processing3と4を比較して、単純計算が速くなってる

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

参考

JavaのBigIntegerにはmodPow()がある

python版 sqrt(n)

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秒なので、実行速度は期待できないのかな。