perl-scripts/Math/binary_prime_sieve_mpz.pl at master · trizen/perl-scripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# License: GPLv3
# Date: 18 May 2017
# https://github.com/trizen

# A binary sieve for prime numbers.

# Useful only when memory is very restricted.

use 5.010;
use strict;
use warnings;

use Math::GMPz;

sub binary_prime_sieve {
    my ($n) = @_;

    my $t = Math::GMPz::Rmpz_init_set_ui(1);
    my $c = Math::GMPz::Rmpz_init_set_ui(1);

    Math::GMPz::Rmpz_setbit($c, $n);

    foreach my $i (2 .. sqrt($n)) {
        Math::GMPz::Rmpz_mul_2exp($t, $t, $n - $i**2);

        for (my $j = $i**2 ; $j <= $n ; $j += $i) {
            Math::GMPz::Rmpz_ior($c, $c, $t);
            Math::GMPz::Rmpz_div_2exp($t, $t, $i);
        }

        Math::GMPz::Rmpz_set_ui($t, 1);
    }

    my $bin = Math::GMPz::Rmpz_get_str($c, 2);

    my @primes;
    foreach my $p (2 .. $n) {
        substr($bin, $p, 1) || push(@primes, $p);
    }
    return @primes;
}

my $n = shift(@ARGV) // 100;
my @primes = binary_prime_sieve($n);
say join(' ', @primes);
say "PI($n) = ", scalar(@primes);