12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788 |
- #!/usr/bin/perl
- # Daniel "Trizen" Șuteu
- # Date: 10 March 2021
- # https://github.com/trizen
- # Let's consider the following function:
- # a(n,v) = Sum_{k=1..n} (v mod k)
- # The goal is to compute a(n,v) in sublinear time with respect to v.
- # Formula:
- # a(n,v) = n*v - A024916(v) + Sum_{k=n+1..v} k*floor(v/k).
- # Formula derived from:
- # a(n,v) = Sum_{k=1..n} (v - k*floor(v/k))
- # = n*v - Sum_{k=1..n} k*floor(v/k)
- # = n*v - Sum_{k=1..v} k*floor(v/k) + Sum_{k=n+1..v} k*floor(v/k)
- # Related problem:
- # Is there a sublinear formula for computing: Sum_{1<=k<=n, gcd(k,n)=1} k*floor(n/k) ?
- # See also:
- # https://oeis.org/A099726 -- Sum of remainders of the n-th prime mod k, for k = 1,2,3,...,n.
- # https://oeis.org/A340976 -- Sum_{1 < k < n} sigma(n) mod k, where sigma = A000203.
- # https://oeis.org/A340180 -- a(n) = Sum_{x in C(n)} (sigma(n) mod x), where C(n) is the set of numbers < n coprime to n, and sigma = A000203.
- use 5.020;
- use strict;
- use warnings;
- use ntheory qw(:all);
- use experimental qw(signatures);
- sub triangular ($n) { # Sum_{k=1..n} k = n-th triangular number
- divint(mulint($n, addint($n, 1)), 2);
- }
- sub sum_of_sigma ($n) { # A024916(n) = Sum_{k=1..n} sigma(k) = Sum_{k=1..n} k*floor(n/k)
- my $T = 0;
- my $s = sqrtint($n);
- foreach my $k (1 .. $s) {
- my $t = divint($n, $k);
- $T = vecsum($T, triangular($t), mulint($k, $t));
- }
- subint($T, mulint(triangular($s), $s));
- }
- sub g ($a, $b) { # g(a,b) = Sum_{k=a..b} k*floor(b/k)
- my $T = 0;
- while ($a <= $b) {
- my $t = divint($b, $a);
- my $u = divint($b, $t);
- $T = addint($T, mulint($t, subint(triangular($u), triangular(subint($a, 1)))));
- $a = addint($u, 1);
- }
- return $T;
- }
- sub sum_remainders ($n, $v) { # sub-linear formula
- addint(subint(mulint($n, $v), sum_of_sigma($v)), g(addint($n, 1), $v));
- }
- say sprintf "[%s]", join(', ', map { sum_remainders($_, nth_prime($_)) } 1 .. 20); #=> A099726
- say sprintf "[%s]", join(', ', map { sum_remainders($_ - 1, divisor_sum($_)) } 1 .. 20); #=> A340976
- foreach my $k (1 .. 8) {
- say("A099726(10^$k) = ", sum_remainders(powint(10, $k), nth_prime(powint(10, $k))));
- }
- __END__
- A099726(10^1) = 30
- A099726(10^2) = 2443
- A099726(10^3) = 248372
- A099726(10^4) = 25372801
- A099726(10^5) = 2437160078
- A099726(10^6) = 252670261459
- A099726(10^7) = 24690625139657
- A099726(10^8) = 2516604108737704
|