new file: Math/sum_of_perfect_powers.sf

Author: trizen

Math/sum_of_perfect_powers.sf

@@ -0,0 +1,41 @@
+#!/usr/bin/ruby
+
+# 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)
+#
+# where:
+#   faulhaber(n,k) = Sum_{j=1..n} j^k.
+
+# See also:
+#   https://oeis.org/A069623
+
+func perfect_power_sum(n) {
+    n.faulhaber(1) - sum(1..n.ilog(2), {|k|
+        mu(k) * (n.iroot(k).faulhaber(k) - 1)
+    })
+}
+
+for n in (0..15) {
+    printf("a(10^%d) = %s\n", n, perfect_power_sum(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

README.md

@@ -667,6 +667,7 @@ A nice collection of day-to-day Sidef scripts.
     * [Sum of nth power digits](./Math/sum_of_nth_power_digits.sf)
     * [Sum of number of divisors of gcd x y](./Math/sum_of_number_of_divisors_of_gcd_x_y.sf)
     * [Sum of number of unitary divisors](./Math/sum_of_number_of_unitary_divisors.sf)
+    * [Sum of perfect powers](./Math/sum_of_perfect_powers.sf)
     * [Sum of polygonal numbers function recursive](./Math/sum_of_polygonal_numbers_function_recursive.sf)
     * [Sum of prime-power exponents of factorial](./Math/sum_of_prime-power_exponents_of_factorial.sf)
     * [Sum of prime-power exponents of product of binomials](./Math/sum_of_prime-power_exponents_of_product_of_binomials.sf)