view release on metacpan or  search on metacpan

examples/partial_sums_of_sigma_function.pl  view on Meta::CPAN

#!/usr/bin/perl

# Author: Daniel "Trizen" Șuteu
# Date: 10 Novermber 2018
# https://github.com/trizen

# A generalized formula with O(sqrt(n)) complexity for computing the partial-sum of `k^m * sigma_j(k)`, for `1 <= k <= n`:
#
#   Sum_{k=1..n} k^m * sigma_j(k)
#
# for any fixed integers m >= 0 and j >= 0.

# Formula:
#   Sum_{k=1..n} k^m * sigma_j(k) =   Sum_{k=1..floor(sqrt(n))} F(m, k) * (F(m+j, floor(n/k)) - F(m+j, floor(n/(k+1))))
#                                   + Sum_{k=1..floor(n/(floor(sqrt(n))+1))} k^(m+j) * 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)) / (n+1)
#
# where Bernoulli(n,x) are the Bernoulli polynomials and Bernoulli(n) is the n-th Bernoulli number.

# See also:
#   https://en.wikipedia.org/wiki/Divisor_function
#   https://en.wikipedia.org/wiki/Faulhaber%27s_formula
#   https://en.wikipedia.org/wiki/Bernoulli_polynomials
#   https://trizenx.blogspot.com/2018/11/partial-sums-of-arithmetical-functions.html

use 5.020;
use strict;
use warnings;

use lib qw(../lib);
use experimental qw(signatures);
use Math::AnyNum qw(isqrt idiv ipow bernoulli faulhaber_sum dirichlet_sum);

sub faulhaber_partial_sum_of_sigma($n, $m, $j) {      # using Faulhaber's formula

    my $total = 0;

    my $s = isqrt($n);
    my $u = int($n / ($s + 1));

    for my $k (1 .. $s) {
        $total += faulhaber_sum($k, $m) * (
                  faulhaber_sum(int($n/$k),     $m+$j)
                - faulhaber_sum(int($n/($k+1)), $m+$j)
            );
    }

    for my $k (1 .. $u) {
        $total += ipow($k, $m+$j) * faulhaber_sum(int($n/$k), $m);
    }

    return $total;
}

sub bernoulli_partial_sum_of_sigma($n, $m, $j) {      # using Bernoulli polynomials

    my $total = 0;

    my $s = isqrt($n);
    my $u = int($n / ($s + 1));

    for my $k (1 .. $s) {
        $total += (
                      bernoulli($m+1, $k+1)
                    - bernoulli($m+1)
                  )/($m+1)
                * (
                      bernoulli($m+$j+1, 1+int($n/$k))
                    - bernoulli($m+$j+1, 1+int($n/($k+1)))
                  )/($m+$j+1);
    }

    for my $k (1 .. $u) {
        $total += ipow($k, $m+$j) * (bernoulli($m+1, 1+int($n/$k)) - bernoulli($m+1)) / ($m+1);
    }

    return $total;
}

sub dirichlet_partial_sum_of_sigma ($n, $m, $j) {    # using Dirichlet's hyperbola method

    my $total = 0;
    my $s = isqrt($n);

    for my $k (1 .. $s) {
        $total += ipow($k, $m)    * faulhaber_sum(int($n/$k), $m+$j);
        $total += ipow($k, $m+$j) * faulhaber_sum(int($n/$k), $m);
    }

    $total -= faulhaber_sum($s, $m) * faulhaber_sum($s, $j+$m);

    return $total;
}