project-euler/Perl/708 Twos are all you need.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 22 May 2020
# https://github.com/trizen

# Twos are all you need
# https://projecteuler.net/problem=708

# Solution by counting the number of k-almost primes <= n.

# Let:
#   pi_k(n) = the number of k-almost primes <= n.

# Then:
#   S(n) = Sum_{k=1..n} 2^bigomega(k)
#        = Sum_{k=1..floor(log_2(n))} 2^k * pi_k(n)

# Runtime: 2 min, 23 sec

use 5.020;

use integer;
use ntheory qw(:all);
use experimental qw(signatures);

my @pi_table;
my $table_len = 1e5;

{
    my $count = 0;
    foreach my $k(0..$table_len) {
        if (is_prime($k)) {
            ++$count;
        }
        $pi_table[$k] = $count;
    }
}

sub prime_count($n) {
    ($n <= $table_len) ? $pi_table[$n] : ntheory::prime_count($n);
}

sub k_prime_count ($n, $k, $i = 14) {

    if ($k == 1) {
        return [0, 4, 25, 168, 1229, 9592, 78498, 664579, 5761455, 50847534, 455052511, 4118054813, 37607912018, 346065536839, 3204941750802, 29844570422669, 279238341033925, 2623557157654233, 24739954287740860, 234057667276344607, 2220819602560918840, 21127269486018731928, 201467286689315906290, 1925320391606803968923, 18435599767349200867866, 176846309399143769411680, 1699246750872437141327603, 16352460426841680446427399]->[$i];
    }

    if ($k == 2) {
        return [0, 4, 34, 299, 2625, 23378, 210035, 1904324, 17427258, 160788536, 1493776443, 13959990342, 131126017178, 1237088048653, 11715902308080, 111329817298881, 1061057292827269, 10139482913717352, 97123037685177087, 932300026230174178, 8966605849641219022, 86389956293761485464]->[$i];
    }

    if ($k == 3) {
        return [0, 1, 22, 247, 2569, 25556, 250853, 2444359, 23727305, 229924367, 2227121996, 21578747909, 209214982913, 2030133769624, 19717814526785, 191693417109381, 1865380637252270, 18168907486812690, 177123437184971927, 1728190923820610000]->[$i];
    }

    if ($k == 4) {
        return [    0, 0, 12, 149, 1712, 18744, 198062, 2050696, 20959322, 212385942, 2139236881, 21454599814, 214499908019, 2139634739326, 21306682904040, 211905511283590, 2105504493045818, 20905484578206982]->[$i];
    }

    if ($k == 5) {
        return [0, 0, 4, 76, 963, 11185, 124465, 1349779, 14371023, 150982388, 1570678136, 16218372618, 166497674684, 1701439985694, 17323079621014]->[$i];
    }

    if ($k == 6) {
        return [0, 0, 2, 37, 485, 5933, 68963, 774078, 8493366, 91683887, 977694273, 10327249593, 108264085934, 1128049914377, 11694704489580]->[$i];
    }

    if ($k == 7) {
        return [    0, 0, 0, 14, 231, 2973, 35585, 409849, 4600247, 50678212, 550454756, 5913771637, 62981797962, 665997804082, 7001087934965]->[$i];
    }

    if ($k == 8) {
        return [    0, 0, 0, 7, 105, 1418, 17572, 207207, 2367507, 26483012, 291646797, 3173159326, 34192782745, 365561221293, 3882841742380]->[$i];
    }

    if ($k == 9) {
        return [    0, 0, 0, 2, 47, 671, 8491, 101787, 1180751, 13377156, 148930536, 1636170477, 17787688377, 191742524399, 2052389350029]->[$i];
    }

    if ($k == 10) {
        return [0, 0, 0, 0, 22, 306, 4016, 49163, 578154, 6618221, 74342563, 823164388, 9011965866, 97765974368, 1052666075366]->[$i];
    }

    if ($k == 11) {
        return [0, 0, 0, 0, 7, 138, 1878, 23448, 279286, 3230577, 36585097, 407818620, 4490844534, 48972151631, 529781669333]->[$i];
    }

    if ($k == 12) {
        return [0, 0, 0, 0, 3, 63, 865, 11068, 133862, 1563465, 17836903, 200051717, 2214357712, 24255601105, 263439785143]->[$i];
    }

    my $count = 0;

    sub ($m, $p, $r) {

        my $s = rootint(divint($n, $m), $r);

        if ($r == 2) {
            my $j = prime_count($p) - 2;

            forprimes {
                $count += (prime_count(divint($n, $m * $_)) - ++$j);
            } $p, $s;

            return;
        }

        for (my $q = $p ; $q <= $s ; $q = next_prime($q)) {
            __SUB__->($m * $q, $q, $r - 1);
        }
    }->(1, 2, $k);

    return $count;
}

sub S($i) {
    my $t = 1;
    my $n = powint(10, $i);

    foreach my $k(1..logint($n, 2)) {
        say "Computing: $k";
        $t += powint(2, $k) * k_prime_count($n, $k, $i);
    }

    return $t;
}

say "S(10^8)  = ", S(8);
say "S(10^14) = ", S(14);