partial_sums_of_n_over_k-almost_prime_divisors.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 20 July 2020
# https://github.com/trizen

# Sublinear algorithm for computing the following partial sum:
#   S(n) = Sum_{k=1..n} Sum_{d|k, d is r-almost prime} (k/d)^m

# Equivalently:
#   S(n) = Sum_{t is r-almost prime <= n} F_m(floor(n/t))
# where F_m(x) are the Faulhaber polynomials.

use 5.020;
use strict;
use warnings;

use ntheory qw(:all);
use experimental qw(signatures);
use Math::AnyNum qw(faulhaber_sum ipow);

sub f($n, $r = 1, $m = 0) {

    my $total = 0;
    my $s = sqrtint($n);

    for my $k (1 .. $s) {
        $total += ipow($k, $m) * almost_prime_count($r, int($n/$k));
        $total += faulhaber_sum(int($n/$k), $m) if is_almost_prime($r, $k);
    }

    $total -= faulhaber_sum($s, $m) * almost_prime_count($r, $s);
    $total;
}

my $n = 100;

say f($n, 1, 0);      #=> Sum_{p     <= n} floor(n/p)       = Sum_{k=1..n} omega(k)
say f($n, 2, 0);      #=> Sum_{p*q   <= n} floor(n/(p*q))   = Sum_{k=1..n} (number of semiprime divisors of k)
say f($n, 3, 0);      #=> Sum_{p*q*r <= n} floor(n/(p*q*r)) = Sum_{k=1..n} (number of 3-almost prime divisors of k)

say '';

say f($n, 1, 1);      #=> Sum_{k=1..n} Sum_{d|k, d is prime} k/d
say f($n, 2, 1);      #=> Sum_{k=1..n} Sum_{d|k, d is semiprime} k/d
say f($n, 3, 1);      #=> Sum_{k=1..n} Sum_{d|k, d is 3-almost prime} k/d

say '';

say f($n, 1, 2);      #=> Sum_{k=1..n} Sum_{d|k, d is prime} (k/d)^2
say f($n, 2, 2);      #=> Sum_{k=1..n} Sum_{d|k, d is semiprime} (k/d)^2
say f($n, 3, 2);      #=> Sum_{k=1..n} Sum_{d|k, d is 3-almost prime} (k/d)^2

say "=> Sum_{k=1..10^n} (number of r-almost prime divisors of k)";

foreach my $r(1..10) {
    say "r = $r: {", join(', ', map{ f(powint(10, $_), $r, 0) } 1..10), "}";
}

say "\n=> Sum_{k=1..10^n} Sum_{d|k, d is r-almost prime} k/d";

foreach my $r(1..10) {
    say "r = $r: {", join(', ', map{ f(powint(10, $_), $r, 1) } 1..10), "}";
}

__END__
=> Sum_{k=1..10^n} (number of r-almost prime divisors of k)
r = 1: {11, 171, 2126, 24300, 266400, 2853708, 30130317, 315037281, 3271067968, 33787242719}
r = 2: {5, 122, 1913, 25368, 309107, 3587501, 40365331, 444658798, 4824183366, 51743978073}
r = 3: {1, 58, 1133, 17179, 230719, 2887977, 34547708, 400531419, 4538949470, 50558632114}
r = 4: {0, 22, 540, 9233, 134679, 1797417, 22659565, 274626874, 3233939674, 37258074465}
r = 5: {0, 7, 227, 4370, 68530, 965003, 12701142, 159627891, 1939960994, 22982979719}
r = 6: {0, 2, 87, 1916, 32224, 475757, 6492864, 84065469, 1048002136, 12697321609}
r = 7: {0, 0, 31, 798, 14434, 222925, 3142601, 41737061, 531430463, 6557159407}
r = 8: {0, 0, 10, 320, 6254, 101133, 1470682, 19990495, 259291249, 3249251063}
r = 9: {0, 0, 2, 123, 2636, 44843, 673192, 9358736, 123499047, 1569291893}
r = 10: {0, 0, 0, 43, 1082, 19518, 303259, 4314150, 57902495, 745552461}

=> Sum_{k=1..10^n} Sum_{d|k, d is r-almost prime} k/d
r = 1: {25, 2298, 226342, 22616110, 2261266482, 226124236118, 22612374197143, 2261237139656553, 226123710243814636, 22612371006991736766}
r = 2: {6, 708, 70451, 7039258, 703809052, 70380387011, 7038023049102, 703802183270761, 70380217285372212, 7038021718888470558}
r = 3: {1, 185, 19261, 1926267, 192581190, 19258134188, 1925810130677, 192580966614994, 19258096515198495, 1925809649512680144}
r = 4: {0, 45, 4923, 500170, 50040884, 5004660706, 500471363203, 50047175747701, 5004718038062777, 500471809568738447}
r = 5: {0, 11, 1223, 126815, 12721482, 1272501930, 127253328013, 12725377777502, 1272538042723713, 127253807917463043}
r = 6: {0, 2, 294, 31833, 3202085, 320487410, 32051378868, 3205166314991, 320516898071185, 32051692261591786}
r = 7: {0, 0, 71, 7961, 802623, 80380033, 8039296889, 803941592045, 80394302031247, 8039431476576389}
r = 8: {0, 0, 14, 1987, 200573, 20122035, 2012708079, 201279547587, 20128037882005, 2012804711838236}
r = 9: {0, 0, 2, 478, 50020, 5033105, 503486440, 50352373220, 5035281352929, 503528648179002}
r = 10: {0, 0, 0, 106, 12431, 1257575, 125898801, 12591617913, 1259181979675, 125918535892823}