123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132 |
- #!/usr/bin/perl
- # Author: Daniel "Trizen" Șuteu
- # Date: 22 November 2018
- # https://github.com/trizen
- # A new algorithm for computing the partial-sums of the Dedekind psi function `ψ_m(k)`, for `1 <= k <= n`:
- #
- # Sum_{k=1..n} ψ_m(k)
- #
- # for any fixed integer m >= 1.
- # Based on the formula:
- # Sum_{k=1..n} ψ_m(k) = Sum_{k=1..n} moebius(k)^2 * F(m, floor(n/k))
- #
- # where F(n,x) is Faulhaber's formula for `Sum_{k=1..x} k^n`, defined in terms of Bernoulli polynomials as:
- # F(n, x) = (Bernoulli(n+1, x+1) - Bernoulli(n+1, 1)) / (n+1)
- # Example for a(n) = Sum_{k=1..n} ψ_2(k):
- # a(10^1) = 462
- # a(10^2) = 400576
- # a(10^3) = 394504950
- # a(10^4) = 393921912410
- # a(10^5) = 393861539651230
- # a(10^6) = 393855661025817568
- # a(10^7) = 393855049001462029696
- # a(10^8) = 393854989687473892017612
- # a(10^9) = 393854983651633712634417940
- # a(10^10) = 393854983070527507612754907046
- # For m=1..3, we have the following asymptotic formulas:
- # Sum_{k=1..n} ψ_1(k) ~ n^2 * zeta(2) / (2*zeta(4))
- # Sum_{k=1..n} ψ_2(k) ~ n^3 * zeta(3) / (3*zeta(6))
- # Sum_{k=1..n} ψ_3(k) ~ n^4 * zeta(4) / (4*zeta(8))
- # In general, for m>=1, we have:
- # Sum_{k=1..n} ψ_m(k) ~ n^(m+1) * zeta(m+1) / ((m+1) * zeta(2*(m+1)))
- # See also:
- # https://oeis.org/A173290
- # https://en.wikipedia.org/wiki/M%C3%B6bius_function
- # https://en.wikipedia.org/wiki/Dedekind_psi_function
- # https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html
- use 5.020;
- use strict;
- use warnings;
- use experimental qw(signatures);
- use Math::AnyNum qw(ipow faulhaber_sum);
- use ntheory qw(jordan_totient moebius vecsum sqrtint forsquarefree is_square_free);
- sub squarefree_count {
- my ($n) = @_;
- my $k = 0;
- my $count = 0;
- foreach my $m (moebius(1, sqrtint($n))) {
- ++$k; $count += $m * int($n / $k / $k);
- }
- return $count;
- }
- sub dedekind_psi_partial_sum ($n, $m) { # O(sqrt(n)) complexity
- my $total = 0;
- my $s = sqrtint($n);
- my $u = int($n / ($s + 1));
- my $prev = squarefree_count($n);
- for my $k (1 .. $s) {
- my $curr = squarefree_count(int($n / ($k + 1)));
- $total += ($prev - $curr) * faulhaber_sum($k, $m);
- $prev = $curr;
- }
- forsquarefree {
- $total += faulhaber_sum(int($n / $_), $m);
- } $u;
- return $total;
- }
- sub dedekind_psi_partial_sum_2 ($n, $m) { # O(sqrt(n)) complexity
- my $total = 0;
- my $s = sqrtint($n);
- for my $k (1 .. $s) {
- $total += ipow($k, $m) * squarefree_count(int($n/$k));
- $total += faulhaber_sum(int($n/$k), $m) if is_square_free($k);
- }
- $total -= squarefree_count($s) * faulhaber_sum($s, $m);
- return $total;
- }
- sub dedekind_psi_partial_sum_test ($n, $m) { # just for testing
- vecsum(map { jordan_totient(2*$m, $_) / jordan_totient($m, $_) } 1 .. $n);
- }
- for my $m (1 .. 10) {
- my $n = int rand 1000;
- my $t1 = dedekind_psi_partial_sum($n, $m);
- my $t2 = dedekind_psi_partial_sum_2($n, $m);
- my $t3 = dedekind_psi_partial_sum_test($n, $m);
- die "error: $t1 != $t2" if ($t1 != $t2);
- die "error: $t1 != $t3" if ($t1 != $t3);
- say "Sum_{k=1..$n} psi_$m(k) = $t1";
- }
- __END__
- Sum_{k=1..626} psi_1(k) = 298020
- Sum_{k=1..203} psi_2(k) = 3314412
- Sum_{k=1..527} psi_3(k) = 20858324486
- Sum_{k=1..912} psi_4(k) = 131086192304600
- Sum_{k=1..221} psi_5(k) = 20014030184914
- Sum_{k=1..980} psi_6(k) = 125495875567427222916
- Sum_{k=1..892} psi_7(k) = 50529225624273249380976
- Sum_{k=1..831} psi_8(k) = 21153451972416324344508126
- Sum_{k=1..384} psi_9(k) = 7069511971715257063270976
- Sum_{k=1..434} psi_10(k) = 9477667039001209551910807864
|