123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 05 April 2022
- # https://github.com/trizen
- # A sublinear algorithm for computing the partial sums of the Euler totient function times k.
- # The partial sums of the Euler totient function is defined as:
- #
- # a(n,m) = Sum_{k=1..n} k * phi(k)
- #
- # where phi(k) is the Euler totient function.
- # Example:
- # a(10^1) = 217
- # a(10^2) = 203085
- # a(10^3) = 202870719
- # a(10^4) = 202653667159
- # a(10^5) = 202643891472849
- # a(10^6) = 202642368741515819
- # a(10^7) = 202642380629476099463
- # a(10^8) = 202642367994273571457613
- # a(10^9) = 202642367530671221417109931
- # a(10^10) = 202642367286524384080814204093
- # General asymptotic formula:
- #
- # Sum_{k=1..n} k^m * phi(k) ~ F_(m+1)(n) / zeta(2).
- #
- # where F_m(n) are the Faulhaber polynomials.
- # OEIS sequences:
- # https://oeis.org/A011755 -- Sum_{k=1..n} k*phi(k).
- # https://oeis.org/A002088 -- Sum of totient function: a(n) = Sum_{k=1..n} phi(k).
- # https://oeis.org/A064018 -- Sum of the Euler totients phi for 10^n.
- # https://oeis.org/A272718 -- Partial sums of gcd-sum sequence A018804.
- # See also:
- # https://en.wikipedia.org/wiki/Faulhaber's_formula
- # https://en.wikipedia.org/wiki/Dirichlet_hyperbola_method
- # https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use ntheory qw(:all);
- sub triangular ($n) {
- divint(mulint($n, $n + 1), 2);
- }
- sub square_pyramidal ($n) {
- divint(vecprod($n, $n + 1, mulint(2, $n) + 1), 6);
- }
- sub partial_sums_of_euler_totient ($n) {
- my $s = sqrtint($n);
- my @euler_sum_lookup = (0);
- my $lookup_size = int(2 * rootint($n, 3)**2);
- my @euler_phi = euler_phi(0, $lookup_size);
- foreach my $i (1 .. $lookup_size) {
- $euler_sum_lookup[$i] = addint($euler_sum_lookup[$i - 1], mulint($i, $euler_phi[$i]));
- }
- my %seen;
- sub ($n) {
- if ($n <= $lookup_size) {
- return $euler_sum_lookup[$n];
- }
- if (exists $seen{$n}) {
- return $seen{$n};
- }
- my $s = sqrtint($n);
- my $T = square_pyramidal($n);
- foreach my $k (2 .. divint($n, $s + 1)) {
- $T = subint($T, mulint($k, __SUB__->(divint($n, $k))));
- }
- my $prev = triangular($n);
- foreach my $k (1 .. $s) {
- my $curr = triangular(divint($n, $k + 1));
- $T = subint($T, mulint(subint($prev, $curr), __SUB__->($k)));
- $prev = $curr;
- }
- $seen{$n} = $T;
- }->($n);
- }
- foreach my $n (1 .. 8) { # takes ~5 seconds
- say "a(10^$n) = ", partial_sums_of_euler_totient(powint(10, $n));
- }
|