『エラトステネスの篩』各言語記述例まとめ

タイトル通りです．筆者過去まとめ記事と趣旨は同じなのですが，一説に『ネット上のエラトステネスの篩のプログラムの大半がエラトステネスの篩じゃない』というのがあるらしく，その割合を少しでも減らす，という目論見もあったりします．なにそれ．

初版で用意した各言語記述例については，とりあえずWikipediaの説明に沿っているつもりですが，上記の記事の趣旨から，それちげーよというツッコミ大歓迎．また，最初にPython版を作成してそれを基に各言語に書き換えたという経緯から，各言語の特徴を活かした，よりシンプルでわかりやすい記述例も大募集．ただし，次の条件を満たすプログラムとします．

• アルゴリズムが『エラトステネスの篩』であること．
• 1から一億までの範囲の素数を一旦求め，その中で最大の素数（99999989）を出力すること（別館除く）．

差し替える際には，当方の次の環境で実行を確認することにしています．なお，非力なハードウェアなのはわざとです(^^;)．この環境で実行できる限り，他のプログラミング言語の記述例も募集中．

• ASUS Zenfone 5: Qualcomm Snapdragon 636 1.8GHz, 6GBメモリ
• Android 9 + Termux + AnLinux(Debian GNU/Linux 10) + 各言語処理系

$ uname -a
Linux localhost 4.4.141-perf+ #1 SMP PREEMPT Fri Jul 10 17:48:58 CST 2020 aarch64 GNU/Linux

※追記事項

【2022-01-12】【2022-02-18】
現在，次のPython記述相当のアルゴリズムへの統一および各言語版記事への個別化を進めています．個別記事化が終了次第，本記事をアルゴリズム解説を中心としたまとめ記事として改修する予定です．なお，既に個別記事としている旧『別館』の各言語版は修正済です．

x = int(input())     # 上限の整数値xを入力
a = [True] * (x + 1) # 素数判定用の配列a[1]～a[x]相当を全て真の値で初期化
a[1] = False         # a[1]に，iが素数ではないことを示す偽の値を格納
r = []               # 素数を格納する配列相当を初期化

i = 1
while i <= x:             # 整数i = 1, 2, ..., xそれぞれについて，
    if a[i]:              # もしiが真＝素数であり，
        if i <= x ** 0.5: # そのiがxの平方根の値以下であれば，
            j = i * i
            while j <= x: # i*iからxまでのiの倍数は素数ではないとする
                a[j] = False
                j = j + i
        r += [i]          # 素数iをrに格納
    i = i + 1

print(len(r)) # 素数の個数を表示
print(r[-1])  # 最も大きい素数を表示

# 配列の添字が1から始まる場合は，aの領域をx個確保し，1から始める．
# 偽による配列の初期化の方が容易であれば，配列aの真偽値を逆に扱う．
# 配列の確保が難しい場合は，rに最大の素数を順次格納して最後に表示する．
# xの平方根を求めることが難しい場合は，i <= x ** 0.5の判定を省略する．

また，次の試し割り法に基づく素因数分解を用いた素数判定とは異なるアルゴリズムとしています．

x = int(input()) # 上限の整数値xを入力
r = []           # 素数を格納する配列相当を初期化

i = 1
while i <= x:          # 整数i = 1, 2, ..., xそれぞれについて，
    n, p, f = i, 2, []
    while n > 1:       # 試し割り法でiの素因数分解を行い，
        if n % p == 0:
            f += [p]
            n = n // p
        else:
            p = p + 1
    if len(f) == 1:    # もしiの素因数が見つからなければ，
        r += [i]       # iを素数としてrに格納
    i = i + 1

print(len(r)) # 素数の個数を表示
print(r[-1])  # 最も大きい素数を表示

【2022-01-02】従来『別館』としていた次の各言語版は個別記事としました．

• シェルスクリプト版
• bash版
• awk版
• bc版
• Scratch版

各言語記述例

Python（個別記事移行）

Ruby

sieve.rb

x = gets.to_i
a = Array.new(x+1, true)
a[1] = false

1.upto(x**0.5) { |i|
  if a[i]
    (i*i).step(x, i) { |j|
      a[j] = false
    }
  end
}

r = []
1.upto(x) { |i|
  if a[i] then r.push(i) end
}
p r[-1]

Ruby2.5.5

$ time ruby sieve.rb <<< 100000000
99999989

real    0m55.606s
user    0m53.630s
sys     0m0.870s

JavaScript

sieve.js

let x;
if (process.argv.length < 3)
  x = 10000000;
else
  x = Number(process.argv[2]);
let a = new Array(x+1).fill(true);
a[1] = false;

for (i = 1; i*i <= x; i++)
  if (a[i])
    for (j = i*i; j <= x; j += i)
      a[j] = false;

let r = [];
for (i = 1; i <= x; i++) if (a[i]) r.push(i);
console.log(r.slice(-1)[0]);

Node.js10.24.0

$ time node sieve.js 100000000
99999989

real    0m35.230s
user    0m33.300s
sys     0m1.310s

C言語

sieve.c

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

int main(int argc, char *argv[])
{
  int x = 100000000;
  if (argc > 1) x = atoi(argv[1]);
  if (x < 1) return 0;
  bool *a = calloc(sizeof(bool), x + 1);
  a[1] = true;

  for (int i = 1; i * i <= x; i++)
    if (!a[i])
      for (int j = i * i; j <= x; j += i)
        a[j] = true;

  unsigned int c = 0;
  for (int i = 1; i <= x; i++) if (!a[i]) c++;
  unsigned int *r = calloc(sizeof(unsigned int), c + 1);

  for (int i = 1; i <= x; i++)
    if (!a[i]) { r[c] = i; c--; }
  printf("%u\n", r[1]);

  return 0;
}

gcc8.3.0

$ gcc -O2 -o sieve sieve.c
$ time ./sieve 100000000
99999989

real    0m4.061s
user    0m3.920s
sys     0m0.090s

Scheme

sieve.scm

(define (sieve x)
  (let ((a (make-vector (+ x 1) #t)) (s (sqrt x)))
    (vector-set! a 1 #f)
    (let loop1 ((i 1))
      (if (<= i s)
          (if (vector-ref a i)
              (let loop2 ((j (* i i)))
                (cond ((<= j x) (vector-set! a j #f)
                                (loop2 (+ j i)))
                      (else (loop1 (+ i 1)))))
              (loop1 (+ i 1)))))
    (let loop ((i 0) (r '()))
      (if (<= i x)
          (if (vector-ref a i)
              (loop (+ i 1) (cons i r))
              (loop (+ i 1) r))
          r))))

(display (car (sieve (read)))) (newline)

Gauche0.9.11-p1

$ time gosh sieve.scm <<< 100000000
99999989

real    0m32.206s
user    0m31.180s
sys     0m0.550s

Julia

sieve.jl

x = parse(Int, readline())
a = trues(x)
a[1] = false

for i in 2:isqrt(x)
    a[i] && (a[i*i:i:x] .= false)
end

println((1:x)[a][end])

Julia1.0.3

$ time julia sieve.jl <<< 100000000
99999989

real    0m3.354s
user    0m2.950s
sys     0m1.460s

R

sieve.r

f <- file("stdin", "r")
x <- as.numeric(readLines(f, 1))
a <- 1:x
a[1] <- 0

for (i in 2:floor(sqrt(x)))
    if (a[i])
        a[i:(x %/% i)*i] <- 0

r <- a[a > 0]
r[length(r)]

R3.5.2

$ time R -f sieve.r --slave <<< 100000000
[1] 99999989

real    0m18.204s
user    0m14.550s
sys     0m4.180s

Common Lisp

sieve.lisp

(defun sieve (x)
  (let ((a (make-array (+ x 1) :initial-element t))
        (s (sqrt x)))
    (setf (aref a 1) nil)
    (do ((i 1 (+ i 1)))
        ((> i s))
      (if (aref a i)
          (do ((j (* i i) (+ j i)))
              ((> j x))
            (setf (aref a j) nil))))
    (do ((i 0 (+ i 1))
         (r '() (if (aref a i) (cons i r) r)))
        ((> i x) r))))

(princ (car (sieve (read)))) (terpri)

SBCL1.4.16

$ time sbcl --script sieve.lisp <<< 100000000
99999989

real    0m21.725s
user    0m18.960s
sys     0m1.740s

Perl

sieve.pl

my $x = <STDIN>;
my @a;
my $i;
my $j;
$a[1] = 1;

for ($i = 1; $i * $i <= $x; $i++) {
  if (!$a[$i]) {
    for ($j = $i * $i; $j <= $x; $j += $i) {
      $a[$j] = 1;
    }
  }
}

my @r;
my $c = 0;
for ($i = 1; $i <= $x; $i++) {
  if (!$a[$i]) { $r[$c] = $i; $c++; }
}
print $r[$c-1], "\n";

Perl5.28.1

$ time perl sieve.pl <<< 100000000
99999989

real    2m10.029s
user    2m0.990s
sys     0m4.830s

PHP

sieve.php

<?php
declare(strict_types = 1);

$x = $argc > 1 ? intval($argv[1]) : 100000000;
if ($x < 1) exit;
$a = array_fill(1, $x, true);
$a[1] = false;

for ($i = 1; $i * $i <= $x; $i++) {
  if ($a[$i]) {
    for ($j = $i * $i; $j <= $x; $j += $i) {
      $a[$j] = false;
    }
  }
}

$r = [];
$c = 0;
for ($i = 1; $i <= $x; $i++) {
  if ($a[$i]) $r[$c++] = $i;
}
echo $r[$c - 1], PHP_EOL;

PHP7.3.29

$ time php sieve.php 100000000
99999989

real    0m56.325s
user    0m53.050s
sys     0m2.430s

Lua

※インタプリタ実行では一億までの処理が行えず強制終了してしまうため，バイトコードコンパイルで実行しています．

sieve.lua

x = tonumber(io.read())
a = {}
a[1] = true

for i = 1, math.sqrt(x) do
  if not a[i] then
    for j = i * i, x, i do
      a[j] = true
    end
  end
end

r = {}
for i = 1, x do
  if not a[i] then
    table.insert(r, i)
  end
end
print(r[#r])

Lua5.3.3

$ luac sieve.lua
$ time lua luac.out <<< 100000000
99999989

real    1m0.471s
user    0m54.910s
sys     0m4.410s

Java

sieve.java

import java.util.ArrayList;

class sieve
{
    public static void main(String[] args)
    {
        int x = args.length > 0 ? Integer.parseInt(args[0]) : 100000000;
        if (x < 1) return;
        boolean[] a = new boolean[x + 1];
        a[1] = true;

        for (int i = 1; i * i <= x; i++)
            if (!a[i])
                for (int j = i * i; j <= x; j += i)
                    a[j] = true;

        ArrayList<Integer> r = new ArrayList<>();
        for (int i = 1; i <= x; i++)
            if (!a[i]) r.add(i);
        System.out.println(r.get(r.size() - 1));
    }
}

OpenJDK11.0.12

$ javac sieve.java
$ time java sieve 100000000
99999989

real    0m5.017s
user    0m7.860s
sys     0m0.470s

Pascal

sieve.pp

program eratosthenes;

uses SysUtils;

var
  x: longint;
  a: array of boolean;
  i: longint;
  j: longint;
  c: longint;
  r: array of longint;

begin
  x := 100000000;
  if argc > 1 then
    x := StrToInt(argv[1]);
  if x < 1 then exit;
  SetLength(a, x + 1);
  a[1] := true;
  i := 1;
  while i * i <= x do begin
    if not a[i] then begin
      j := i * i;
      while j <= x do begin
        a[j] := true;
        j := j + i;
      end;
    end;
    i := i + 1;
  end;

  c := 0;
  for i := 1 to x do
    if not a[i] then c := c + 1;
  SetLength(r, c + 1);

  j := 0;
  for i := 1 to x do
    if not a[i] then begin
      r[j] := i;
      j := j + 1;
    end;
  writeln(r[c - 1]);
end.

fp-compiler3.0.4

$ fpc -O4 sieve
Free Pascal Compiler version 3.0.4+dfsg-22 [2019/01/24] for aarch64
Copyright (c) 1993-2017 by Florian Klaempfl and others
Target OS: Linux for AArch64
Compiling sieve.pp
Assembling eratosthenes
Linking sieve
45 lines compiled, 1.1 sec
$ time ./sieve 100000000
99999989

real    0m3.802s
user    0m3.680s
sys     0m0.070s

C++

sieve.cpp

#include <iostream>
#include <cstdlib>
#include <list>

int main(int argc, char *argv[])
{
    int x = 100000000;
    if (argc > 1) x = atoi(argv[1]);
    if (x < 1) return 0;
    bool *a = new bool[x + 1];
    a[1] = true;

    for (int i = 1; i * i <= x; i++)
        if (!a[i])
            for (int j = i * i; j <= x; j += i)
                a[j] = true;

    std::list<int> r;
    for (int i = 1; i <= x; i++)
        if (!a[i]) r.push_front(i);
    std::cout << r.front() << std::endl;
}

gcc8.3.0(GNUCompilerCollection)

$ g++ -O2 -o sieve sieve.cpp
$ time ./sieve 100000000
99999989

real    0m4.774s
user    0m4.180s
sys     0m0.330s

FORTRAN

sieve.f

! Sieve of Eratosthenes
implicit none
integer c, x, y, i, j
logical,allocatable :: a(:)
integer,allocatable :: p(:)

read(*,*) x
allocate(a(x))
a(1) = .true.
                                                       y = int(sqrt(real(x)))
do i = 1, y
        if (.not. a(i)) then
                do j = i * i, x, i
                        a(j) = .true.
                end do
        end if
end do

c = 0
do i = 1, x
        if (.not. a(i)) then
                c = c + 1
        end if
end do
print *, c
allocate(p(c))
c = 0
do i = 1, x
        if (.not. a(i)) then
                p(c) = i
                c = c + 1
        end if
end do
print *, p(c - 1)
end

GNUFortran8.3.0

$ gfortran -O3 -o sieve -ffree-form sieve.f
$ time ./sieve <<< 100000000
     5761455
    99999989

real    0m6.993s
user    0m6.540s
sys     0m0.320s

Go

sieve.go

package main
import "fmt"
func main() {
  var x int
  fmt.Scan(&x)
  a := make([]bool, x+1)
  a[1] = true

  for i := 1; i*i <= x; i++ {
    if !a[i] {
      for j := i*i; j <= x; j += i {
        a[j] = true
      }
    }
  }

  c := 0
  for i := 1; i <= x; i++ {
    if !a[i] { c++ }
  }
  r := make([]int, c+1)

  for i := 1; i <= x; i++ {
    if !a[i] {
      r[c] = i
      c--
    }
  }
  fmt.Println(r[1]);
}

go1.11.6

$ time go run sieve.go <<< 100000000
99999989

real    0m6.326s
user    0m4.890s
sys     0m1.010s

COBOL

sieve.cob

    *> Sieve of Eratosthenes
    program-id. sieve.

    data division.
    working-storage section.
    77 a occurs 0 to 100000000
        depending on a-size
        indexed by a-idx
        pic 9.
    77 a-size pic 9(9).
    77 p occurs 0 to 100000000
        depending on p-size
        pic 9(9).
    77 p-size pic 9(9).
    77 c pic 9(9).
    77 i pic 9(9).
    77 j pic 9(9).
    77 k pic 9(9).
    77 x pic 9(9).
    77 y pic 9(9).

    procedure division.
    accept x
    set a-size to x
    move 1 to a(1)

    compute y = function sqrt(x)
    perform varying i from 1 by 1 until i > y
        if (a(i) = 0) then
            compute k = i * i
            perform varying j from k by i until j > x
                move i to a(j)
            end-perform
        end-if
    end-perform

    move 0 to c
    perform varying i from 1 by 1 until i > x
        if (a(i) = 0) then
            add 1 to c
        end-if
    end-perform
    display c
    set p-size to c
    move 0 to c
    perform varying i from 1 by 1 until i > x
        if (a(i) = 0) then
            move i to p(c)
            add 1 to c
        end-if
    end-perform
    subtract 1 from c
    display p(c).

GnuCOBOL2.2.0

$ cobc -x -O3 -F sieve.cob
$ time ./sieve <<< 100000000
005761455
099999989

real    2m42.231s
user    2m41.230s
sys     0m0.150s

Rust

sieve.rs

use std::io;

fn main()
{
  let mut s = String::new();
  io::stdin().read_line(&mut s).ok();
  let x: usize = s.trim().parse().ok().unwrap();

  let mut a = vec![false];
  let mut i = 1;
  while i <= x { a.push(true); i += 1; }
  a[1] = false;

  let mut i = 1;
  while i*i <= x {
    if a[i] {
      let mut j = i*i;
      while j <= x {
        a[j] = false;
        j += i;
      }
    }
    i += 1;
  }

  let mut i = 1;
  let mut r = 0;
  let mut c = 0;
  while i <= x {
    if a[i] { r = i; c += 1; }
    i += 1;
  }

  println!("{}\n{}", c, r);
}

Rust1.41.1

$ rustc -O sieve.rs
$ time ./sieve <<< 100000000
5761455
99999989

real    0m4.234s
user    0m4.090s
sys     0m0.080s

D

sieve.d

import std.stdio, std.bitmanip, std.conv : to;

void main(string[] args)
{
    int x = args.length > 1 ? to!(int)(args[1]) : 100000000;
    BitArray a = BitArray([]);
    a.length = x + 1;

    for (int i = 2; i * i <= x; i++)
        if (!a[i])
            for (int j = i * i; j <= x; j += i)
                a[j] = true;

    int[] r;
    for (int i = 2; i <= x; i++)
        if (!a[i]) r ~= i;
    if (r.length > 0) writeln(r[$ - 1]);
}

LLVM_D_compiler1.12.0

$ ldc2 -O5 sieve.d
$ time ./sieve <<< 100000000
99999989

real    0m3.908s
user    0m3.830s
sys     0m0.020s

備考

更新履歴

• 2022-02-18：アルゴリズムのPython記述修正および試し割り法に基づくPython記述を追加
• 2022-02-18：移行に関する追記事項のアルゴリズムをPython記述＋コメントに変更
• 2022-01-17：Python版を個別記事に移行
• 2022-01-12：移行に関する追記事項を冒頭に追加
• 2022-01-06：Scratch版記事へのリンクを追加
• 2022-01-05：Rust版を追加
• 2022-01-05：COBOL版，D版を追加（コメントより）
• 2022-01-04：Go版を追加
• 2022-01-04：FORTRAN版を追加（コメントより）
• 2022-01-03：Pascal版の実行確認を修正（コメントより）
• 2022-01-02：PHP版，Java版，Pascal版，C++版を追加（コメントより）
• 2022-01-02：別館（POSIXシェル，bash，awk，bc）を個別記事に分離
• 2022-01-02：Lua版をバイトコードコンパイルの結果に置き換えて本館に移動
• 2022-01-02：別館を追加（POSIXシェル，bash，awk，bc，Lua）
• 2022-01-01：Perl版を追加
• 2022-01-01：Common Lisp版を追加
• 2021-12-31：各言語間で変数名を可能な限り統一
• 2021-12-31：C言語版の実行確認を修正（コメントより）
• 2021-12-31：R版を追加（コメントより）
• 2021-12-31：添字の扱い等について各言語版に修正
• 2021-12-31：実行確認を『1から一億までの範囲の素数』に変更
• 2021-12-31：Julia版を追加（コメントより）
• 2021-12-30：Python版を変更（コメントより）
• 2021-12-30：初版公開（Python，Ruby，JavaScript，C言語，Scheme）