perl-scripts/Math/partial_sums_of_exponential_prime_omega_functions.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 16 March 2021
# https://github.com/trizen

# Compute partial sums of the following three functions in sublinear time:
#   S1(n) = Sum_{k=1..n} v^bigomega(k)
#   S2(n) = Sum_{k=1..n} v^omega(k)
#   S3(n) = Sum_{k=1..n} v^omega(k) * mu(k)^2

use 5.020;
use warnings;

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

sub squarefree_almost_prime_count ($k, $n) {

    if ($k == 0) {
        return (($n <= 0) ? 0 : 1);
    }

    if ($k == 1) {
        return prime_count($n);
    }

    my $count = 0;

    sub ($m, $p, $k, $j = 0) {

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

        if ($k == 2) {

            foreach my $q (@{primes($p, $s)}) {

                ++$j;

                if (modint($m, $q) != 0) {
                    $count += prime_count(divint($n, mulint($m, $q))) - $j;
                }
            }

            return;
        }

        foreach my $p (@{primes($p, $s)}) {
            if (modint($m, $p) != 0) {
                __SUB__->(mulint($m, $p), $p, $k - 1, $j);
            }
            ++$j;
        }
    }->(1, 2, $k);

    return $count;
}

sub S1 ($n, $v = 2) {    # Sum_{k=1..n} v^bigomega(k)
    vecsum(map { mulint(powint($v, $_), almost_prime_count($_, $n)) } 0 .. logint($n, 2));
}

sub S2 ($n, $v = 2) {    # Sum_{k=1..n} v^omega(k)
    vecsum(map { mulint(powint($v, $_), omega_prime_count($_, $n)) } 0 .. logint($n, 2));
}

sub S3 ($n, $v = 2) {    # Sum_{k=1..n} v^omega(k) * mu(k)^2
    vecsum(map { mulint(powint($v, $_), squarefree_almost_prime_count($_, $n)) } 0 .. logint($n, 2));
}

say join ', ', map { S1($_) } 1 .. 20;  #=> A069205: [1, 3, 5, 9, 11, 15, 17, 25, 29, 33, 35, 43, 45, 49, 53, 69, 71, 79, 81]
say join ', ', map { S2($_) } 1 .. 20;  #=> A064608: [1, 3, 5, 7, 9, 13, 15, 17, 19, 23, 25, 29, 31, 35, 39, 41, 43, 47, 49]
say join ', ', map { S3($_) } 1 .. 20;  #=> A069201: [1, 3, 5, 5, 7, 11, 13, 13, 13, 17, 19, 19, 21, 25, 29, 29, 31, 31, 33]

say '';

say join ', ', map { S1($_, -1) } 1 .. 20;  #=> A002819: [1, 0, -1, 0, -1, 0, -1, -2, -1, 0, -1, -2, -3, -2, -1, 0, -1, -2, -3, -4]
say join ', ', map { S2($_, -1) } 1 .. 20;  #=> A174863: [1, 0, -1, -2, -3, -2, -3, -4, -5, -4, -5, -4, -5, -4, -3, -4, -5, -4, -5, -4]
say join ', ', map { S3($_, -1) } 1 .. 20;  #=> A002321: [1, 0, -1, -1, -2, -1, -2, -2, -2, -1, -2, -2, -3, -2, -1, -1, -2, -2, -3, -3]

__END__

# A069205(n) = Sum_{k=1..n} 2^bigomega(k)

A069205(10^1)  = 33
A069205(10^2)  = 811
A069205(10^3)  = 15301
A069205(10^4)  = 260615
A069205(10^5)  = 3942969
A069205(10^6)  = 55282297
A069205(10^7)  = 746263855
A069205(10^8)  = 9613563919
A069205(10^9)  = 120954854741
A069205(10^10) = 1491898574939
A069205(10^11) = 17944730372827
A069205(10^12) = 212986333467973
A069205(10^13) = 2498962573520227
A069205(10^14) = 28874142998632109

# A002819(n) = Sum_{k=1..n} (-1)^bigomega(k)
# See also: A090410

A002819(10^1) = 0
A002819(10^2) = -2
A002819(10^3) = -14
A002819(10^4) = -94
A002819(10^5) = -288
A002819(10^6) = -530
A002819(10^7) = -842
A002819(10^8) = -3884
A002819(10^9) = -25216
A002819(10^10) = -116026
A002819(10^11) = -342224
A002819(10^12) = -522626
A002819(10^13) = -966578
A002819(10^14) = -7424752

# A064608(n) = Sum_{k=1..n} 2^omega(k)
# See also: A180361

A064608(10^1)  = 23
A064608(10^2)  = 359
A064608(10^3)  = 4987
A064608(10^4)  = 63869
A064608(10^5)  = 778581
A064608(10^6)  = 9185685
A064608(10^7)  = 105854997
A064608(10^8)  = 1198530315
A064608(10^9)  = 13385107495
A064608(10^10) = 147849112851
A064608(10^11) = 1618471517571
A064608(10^12) = 17584519050293

# A174863(n) = Sum_{k=1..n} (-1)^omega(k)

A174863(10^1)  = -4
A174863(10^2)  = 14
A174863(10^3)  = 64
A174863(10^4)  = -16
A174863(10^5)  = -720
A174863(10^6)  = -1908
A174863(10^7)  = -1650
A174863(10^8)  = 10734
A174863(10^9)  = 53740
A174863(10^10) = 108654
A174863(10^11) = 195702
A174863(10^12) = 27158

# A069201(n) = Sum_{k=1..n} mu(k)^2 * 2^omega(k)

A069201(10^1)  = 17
A069201(10^2)  = 211
A069201(10^3)  = 2825
A069201(10^4)  = 34891
A069201(10^5)  = 414813
A069201(10^6)  = 4808081
A069201(10^7)  = 54684335
A069201(10^8)  = 612868643
A069201(10^9)  = 6788951097
A069201(10^10) = 74492096539
A069201(10^11) = 810947010335
A069201(10^12) = 8769730440341

# A002321(n) = Sum_{k=1..n} (-1)^omega(k) * mu(k)^2 = Sum_{k=1..n} mu(k)
# See also: A084237

A002321(10^1) = -1
A002321(10^2) = 1
A002321(10^3) = 2
A002321(10^4) = -23
A002321(10^5) = -48
A002321(10^6) = 212
A002321(10^7) = 1037
A002321(10^8) = 1928
A002321(10^9) = -222
A002321(10^10) = -33722
A002321(10^11) = -87856
A002321(10^12) = 62366

# A013928(n) = Sum_{k=1..n} mu(k)^2
# See also: A071172

A013928(10^1)  = 7
A013928(10^2)  = 61
A013928(10^3)  = 608
A013928(10^4)  = 6083
A013928(10^5)  = 60794
A013928(10^6)  = 607926
A013928(10^7)  = 6079291
A013928(10^8)  = 60792694
A013928(10^9)  = 607927124
A013928(10^10) = 6079270942
A013928(10^11) = 60792710280
A013928(10^12) = 607927102274