new file: Math/count_of_perfect_powers.pl

new file:   Math/sum_of_perfect_powers.pl
new file:   Math/sum_of_prime_powers.pl

Author: trizen

Math/brilliant_numbers_count.pl Math/count_of_brilliant_numbers.pl

@@ -0,0 +1,39 @@
+#!/usr/bin/perl
+
+# Efficient formula for counting the numbers of perfect powers <= n.
+
+# Formula:
+#   a(n) = n - Sum_{1..floor(log_2(n))} mu(k) * (floor(n^(1/k)) - 1)
+#        = 1 - Sum_{2..floor(log_2(n))} mu(k) * (floor(n^(1/k)) - 1)
+
+# See also:
+#   https://oeis.org/A069623
+
+use 5.036;
+use ntheory qw(logint rootint moebius vecsum);
+
+sub perfect_power_count ($n) {
+    1 - vecsum(map { moebius($_) * (rootint($n, $_) - 1) } 2 .. logint($n, 2));
+}
+
+foreach my $n (0 .. 15) {
+    printf("a(10^%d) = %s\n", $n, perfect_power_count(10**$n));
+}
+
+__END__
+a(10^0) = 1
+a(10^1) = 4
+a(10^2) = 13
+a(10^3) = 41
+a(10^4) = 125
+a(10^5) = 367
+a(10^6) = 1111
+a(10^7) = 3395
+a(10^8) = 10491
+a(10^9) = 32670
+a(10^10) = 102231
+a(10^11) = 320990
+a(10^12) = 1010196
+a(10^13) = 3184138
+a(10^14) = 10046921
+a(10^15) = 31723592

Math/count_of_perfect_powers.pl

@@ -16,6 +16,8 @@
 use strict;
 use warnings;
 
+no warnings 'recursion';
+
 use Memoize qw(memoize);
 use Math::AnyNum qw(:overload floor ceil);
 

Math/prime_power_count.pl Math/count_of_prime_power.pl

@@ -0,0 +1,43 @@
+#!/usr/bin/perl
+
+# Efficient formula for computing the sum of perfect powers <= n.
+
+# Formula:
+#   a(n) = faulhaber(n,1) - Sum_{1..floor(log_2(n))} mu(k) * (faulhaber(floor(n^(1/k)), k) - 1)
+#        = 1 - Sum_{2..floor(log_2(n))} mu(k) * (faulhaber(floor(n^(1/k)), k) - 1)
+#
+# where:
+#   faulhaber(n,k) = Sum_{j=1..n} j^k.
+
+# See also:
+#   https://oeis.org/A069623
+
+use 5.036;
+use ntheory      qw(moebius);
+use Math::AnyNum qw(faulhaber_sum sum ipow iroot ilog2);
+
+sub perfect_power_sum ($n) {
+    1 - sum(map { moebius($_) * (faulhaber_sum(iroot($n, $_), $_) - 1) } 2 .. ilog2($n));
+}
+
+foreach my $n (0 .. 15) {
+    printf("a(10^%d) = %s\n", $n, perfect_power_sum(ipow(10, $n)));
+}
+
+__END__
+a(10^0) = 1
+a(10^1) = 22
+a(10^2) = 452
+a(10^3) = 13050
+a(10^4) = 410552
+a(10^5) = 11888199
+a(10^6) = 361590619
+a(10^7) = 11120063109
+a(10^8) = 345454923761
+a(10^9) = 10800726331772
+a(10^10) = 338846269199225
+a(10^11) = 10659098451968490
+a(10^12) = 335867724220740686
+a(10^13) = 10595345580446344714
+a(10^14) = 334502268562161605300
+a(10^15) = 10566065095217905939231

Math/squarefree_count.pl Math/count_of_squarefree_numbers.pl

@@ -0,0 +1,123 @@
+#!/usr/bin/perl
+
+# Three sublinear algorithms for computing the sum of prime powers <= n,
+# based on the sublinear algorithm for computing the sum of primes <= n.
+
+# See also:
+#   https://oeis.org/A074793
+
+use 5.036;
+use Math::GMPz;
+use ntheory                qw(:all);
+use Math::Prime::Util::GMP qw(faulhaber_sum);
+
+sub sum_of_primes ($n, $k = 1) {    # Sum_{p prime <= n} p^k
+
+    return sum_primes($n) if ($k == 1);    # optimization
+
+    $n > ~0 and return undef;
+    $n <= 1 and return 0;
+
+    my $r = sqrtint($n);
+    my @V = map { divint($n, $_) } 1 .. $r;
+    push @V, CORE::reverse(1 .. $V[-1] - 1);
+
+    my $t = Math::GMPz::Rmpz_init_set_ui(0);
+    my $u = Math::GMPz::Rmpz_init();
+
+    my %S;
+    @S{@V} = map { Math::GMPz::Rmpz_init_set_str(faulhaber_sum($_, $k), 10) } @V;
+
+    foreach my $p (2 .. $r) {
+        if ($S{$p} > $S{$p - 1}) {
+            my $cp = $S{$p - 1};
+            my $p2 = $p * $p;
+            Math::GMPz::Rmpz_ui_pow_ui($t, $p, $k);
+            foreach my $v (@V) {
+                last if ($v < $p2);
+                Math::GMPz::Rmpz_sub($u, $S{divint($v, $p)}, $cp);
+                Math::GMPz::Rmpz_submul($S{$v}, $u, $t);
+            }
+        }
+    }
+
+    $S{$n} - 1;
+}
+
+sub sum_of_prime_powers ($n) {
+
+    # a(n) = Sum_{p prime <= n} p
+    # b(n) = Sum_{p prime <= n^(1/2)} p^2
+    # c(n) = Sum_{p prime <= n^(1/3)} f(p)
+
+    # sum_of_prime_powers(n) = a(n) + b(n) + c(n)
+
+    my $ps1 = sum_of_primes($n);
+    my $ps2 = sum_of_primes(sqrtint($n), 2);
+
+    # f(p) = (Sum_{k=1..floor(log_p(n))} p^k) - p^2 - p
+    #      = (p^(1+floor(log_p(n))) - 1)/(p-1) - p^2 - p - 1
+
+    my $ps3 = 0;
+    foreach my $p (@{primes(rootint($n, 3))}) {
+        $ps3 += divint(powint($p, logint($n, $p) + 1) - 1, $p - 1) - $p * $p - $p - 1;
+    }
+
+    return vecsum($ps1, $ps2, $ps3);
+}
+
+sub sum_of_prime_powers_2 ($n) {
+
+    # a(n) = Sum_{p prime <= n} p
+    # b(n) = Sum_{p prime <= n^(1/2)} f(p)
+
+    # sum_of_prime_powers(n) = a(n) + b(n)
+
+    my $ps1 = sum_of_primes($n);
+
+    # f(p) = (Sum_{k=1..floor(log_p(n))} p^k) - p
+    #      = (p^(1+floor(log_p(n))) - 1)/(p-1) - p - 1
+
+    my $ps2 = 0;
+    forprimes {
+        $ps2 += divint(powint($_, logint($n, $_) + 1) - 1, $_ - 1) - $_ - 1;
+    } sqrtint($n);
+
+    return vecsum($ps1, $ps2);
+}
+
+sub sum_of_prime_powers_3 ($n) {
+
+    # a(n) = Sum_{k=1..floor(log_2(n))} Sum_{p prime <= n^(1/k)} p^k.
+    vecsum(map { sum_of_primes(rootint($n, $_), $_) } 1 .. logint($n, 2));
+}
+
+foreach my $n (0 .. 10) {
+    say "a(10^$n) = ", sum_of_prime_powers(powint(10, $n));
+}
+
+foreach my $k (1 .. 100) {
+    my $n = int(rand(1e3)) + 1;
+
+    my $x = sum_of_prime_powers($n);
+    my $y = sum_of_prime_powers_2($n);
+    my $z = sum_of_prime_powers_3($n);
+
+    $x == $y or die "error";
+    $x == $z or die "error";
+}
+
+__END__
+a(10^0) = 0
+a(10^1) = 38
+a(10^2) = 1375
+a(10^3) = 82674
+a(10^4) = 5850315
+a(10^5) = 457028152
+a(10^6) = 37610438089
+a(10^7) = 3204814813355
+a(10^8) = 279250347324393
+a(10^9) = 24740607755154524
+a(10^10) = 2220853189506845580
+a(10^11) = 201467948093608962539
+a(10^12) = 18435613572072500152927

Math/partitions_count.pl

@@ -502,7 +502,6 @@ A nice collection of day-to-day Perl scripts.
     * [BPSW primality test](./Math/BPSW_primality_test.pl)
     * [BPSW primality test mpz](./Math/BPSW_primality_test_mpz.pl)
     * [Brazilian primes constant](./Math/brazilian_primes_constant.pl)
-    * [Brilliant numbers count](./Math/brilliant_numbers_count.pl)
     * [Brown numbers](./Math/brown_numbers.pl)
     * [Carmichael factorization method](./Math/carmichael_factorization_method.pl)
     * [Carmichael factorization method generalized](./Math/carmichael_factorization_method_generalized.pl)
@@ -544,19 +543,23 @@ A nice collection of day-to-day Perl scripts.
     * [Continued fractions prime constant](./Math/continued_fractions_prime_constant.pl)
     * [Convergent series](./Math/convergent_series.pl)
     * [Cosmic calendar](./Math/cosmic_calendar.pl)
+    * [Count of brilliant numbers](./Math/count_of_brilliant_numbers.pl)
     * [Count of cube-full numbers](./Math/count_of_cube-full_numbers.pl)
     * [Count of integers with gpf of n equals p](./Math/count_of_integers_with_gpf_of_n_equals_p.pl)
     * [Count of integers with lpf of n equals p](./Math/count_of_integers_with_lpf_of_n_equals_p.pl)
     * [Count of k-almost primes](./Math/count_of_k-almost_primes.pl)
     * [Count of k-omega primes](./Math/count_of_k-omega_primes.pl)
     * [Count of k-powerfree numbers](./Math/count_of_k-powerfree_numbers.pl)
     * [Count of k-powerful numbers](./Math/count_of_k-powerful_numbers.pl)
+    * [Count of perfect powers](./Math/count_of_perfect_powers.pl)
+    * [Count of prime power](./Math/count_of_prime_power.pl)
     * [Count of rough numbers](./Math/count_of_rough_numbers.pl)
     * [Count of smooth numbers](./Math/count_of_smooth_numbers.pl)
     * [Count of smooth numbers mpz](./Math/count_of_smooth_numbers_mpz.pl)
     * [Count of smooth numbers mpz 2](./Math/count_of_smooth_numbers_mpz_2.pl)
     * [Count of smooth numbers with k factors](./Math/count_of_smooth_numbers_with_k_factors.pl)
     * [Count of squarefree k-almost primes](./Math/count_of_squarefree_k-almost_primes.pl)
+    * [Count of squarefree numbers](./Math/count_of_squarefree_numbers.pl)
     * [Count subtriangles](./Math/count_subtriangles.pl)
     * [Cube-full numbers](./Math/cube-full_numbers.pl)
     * [Cuboid](./Math/cuboid.pl)
@@ -868,7 +871,6 @@ A nice collection of day-to-day Perl scripts.
     * [Prime functions in terms of zeros of zeta](./Math/prime_functions_in_terms_of_zeros_of_zeta.pl)
     * [Prime numbers generator](./Math/prime_numbers_generator.pl)
     * [Prime omega function generalized](./Math/prime_omega_function_generalized.pl)
-    * [Prime power count](./Math/prime_power_count.pl)
     * [Prime quadratic polynomial analyzer](./Math/prime_quadratic_polynomial_analyzer.pl)
     * [Prime quadratic polynomials](./Math/prime_quadratic_polynomials.pl)
     * [Prime summation](./Math/prime_summation.pl)
@@ -956,7 +958,6 @@ A nice collection of day-to-day Perl scripts.
     * [Squarefree almost primes from factor list](./Math/squarefree_almost_primes_from_factor_list.pl)
     * [Squarefree almost primes in range](./Math/squarefree_almost_primes_in_range.pl)
     * [Squarefree almost primes in range from factor list](./Math/squarefree_almost_primes_in_range_from_factor_list.pl)
-    * [Squarefree count](./Math/squarefree_count.pl)
     * [Squarefree divisors](./Math/squarefree_divisors.pl)
     * [Squarefree fermat overpseudoprimes in range](./Math/squarefree_fermat_overpseudoprimes_in_range.pl)
     * [Squarefree fermat pseudoprimes in range](./Math/squarefree_fermat_pseudoprimes_in_range.pl)
@@ -975,8 +976,10 @@ A nice collection of day-to-day Perl scripts.
     * [Sum of digits subquadratic algorithm](./Math/sum_of_digits_subquadratic_algorithm.pl)
     * [Sum of digits subquadratic algorithm mpz](./Math/sum_of_digits_subquadratic_algorithm_mpz.pl)
     * [Sum of natural powers in constant base](./Math/sum_of_natural_powers_in_constant_base.pl)
+    * [Sum of perfect powers](./Math/sum_of_perfect_powers.pl)
     * [Sum of prime-power exponents of factorial](./Math/sum_of_prime-power_exponents_of_factorial.pl)
     * [Sum of prime-power exponents of product of binomials](./Math/sum_of_prime-power_exponents_of_product_of_binomials.pl)
+    * [Sum of prime powers](./Math/sum_of_prime_powers.pl)
     * [Sum of primes generalized](./Math/sum_of_primes_generalized.pl)
     * [Sum of sigma](./Math/sum_of_sigma.pl)
     * [Sum of sigma 2](./Math/sum_of_sigma_2.pl)