perl-scripts/Math/factorial_expansion_of_reciprocals.pl

#!/usr/bin/perl

# Daniel "Trizen" Șuteu
# Date: 20 December 2018
# https://github.com/trizen

# Factorial expansion of reciprocals of natural numbers.

# For n>1, the length of the factorial expansion of 1/n is the n-th Kempner number `k`, such that:
#   1/n = Sum_{j=0..k} f(j) / j!

# See also:
#   https://oeis.org/A002034
#   https://oeis.org/A122416
#   https://en.wikipedia.org/wiki/Kempner_function

use 5.020;
use strict;
use warnings;

use ntheory qw(factor_exp);
use experimental qw(signatures);
use Math::AnyNum qw(sum factorial bsearch_le);

sub f ($n, $x) {
    return $x->floor if ($n == 0);
    ($x * factorial($n))->floor - (($x * factorial($n - 1))->floor * $n);
}

sub kempner_function ($n) {

    return 0 if ($n == 1);

    my @f = factor_exp($n);
    my $x = Math::AnyNum->new_q(1, $n);

    my $min = 2;
    my $max = $f[-1][0] * $f[-1][1];

    my %seen;

    for (; ;) {

        my $k = bsearch_le(
            $min, $max,
            sub {
                0 <=> f($_, $x);
            }
        );

        if ($seen{$k}++) {
            ++$min;
            ++$max;
        }

        if (   f($k + 0, $x) != 0
            or f($k + 1, $x) != 0
            or f($k + 2, $x) != 0
            or f($k + 3, $x) != 0
        ) {
            $min = $k + 1;
            next;
        }

        if (f($k - 1, $x) == 0) {
            $max = $k;
            next;
        }

        return $k - 1;
    }
}

foreach my $n (1 .. 50) {

    my $x = Math::AnyNum->new_q(1, $n);
    my $k = kempner_function($n);

    my @a = map { f($_, $x) } 0 .. $k;
    my $r = sum(map { $a[$_] / factorial($_) } 0 .. $k);

    say "F($r) = [", join(', ', @a), "]";
    die "error: $r != $x" if ($r != $x);
}

say "\n[2..100] Kempner numbers: ", join(', ', map { kempner_function($_) } 2 .. 100);

__END__
F(1) = [1]
F(1/2) = [0, 0, 1]
F(1/3) = [0, 0, 0, 2]
F(1/4) = [0, 0, 0, 1, 2]
F(1/5) = [0, 0, 0, 1, 0, 4]
F(1/6) = [0, 0, 0, 1]
F(1/7) = [0, 0, 0, 0, 3, 2, 0, 6]
F(1/8) = [0, 0, 0, 0, 3]
F(1/9) = [0, 0, 0, 0, 2, 3, 2]
F(1/10) = [0, 0, 0, 0, 2, 2]
F(1/11) = [0, 0, 0, 0, 2, 0, 5, 3, 1, 4, 0, 10]
F(1/12) = [0, 0, 0, 0, 2]
F(1/13) = [0, 0, 0, 0, 1, 4, 1, 2, 5, 4, 8, 5, 0, 12]
F(1/14) = [0, 0, 0, 0, 1, 3, 3, 3]
F(1/15) = [0, 0, 0, 0, 1, 3]
F(1/16) = [0, 0, 0, 0, 1, 2, 3]
F(1/17) = [0, 0, 0, 0, 1, 2, 0, 2, 3, 6, 8, 9, 0, 9, 2, 7, 0, 16]
F(1/18) = [0, 0, 0, 0, 1, 1, 4]
F(1/19) = [0, 0, 0, 0, 1, 1, 1, 6, 2, 0, 9, 5, 2, 6, 11, 11, 13, 8, 0, 18]
F(1/20) = [0, 0, 0, 0, 1, 1]
F(1/21) = [0, 0, 0, 0, 1, 0, 4, 2]
F(1/22) = [0, 0, 0, 0, 1, 0, 2, 5, 0, 6, 5, 5]
F(1/23) = [0, 0, 0, 0, 1, 0, 1, 2, 1, 0, 3, 10, 0, 6, 10, 14, 5, 9, 10, 18, 3, 10, 0, 22]
F(1/24) = [0, 0, 0, 0, 1]
F(1/25) = [0, 0, 0, 0, 0, 4, 4, 5, 4, 7, 2]
F(1/26) = [0, 0, 0, 0, 0, 4, 3, 4, 6, 6, 9, 2, 6, 6]
F(1/27) = [0, 0, 0, 0, 0, 4, 2, 4, 5, 3]
F(1/28) = [0, 0, 0, 0, 0, 4, 1, 5]
F(1/29) = [0, 0, 0, 0, 0, 4, 0, 5, 6, 3, 1, 0, 4, 7, 2, 6, 3, 5, 4, 18, 6, 18, 18, 4, 18, 5, 4, 13, 0, 28]
F(1/30) = [0, 0, 0, 0, 0, 4]
F(1/31) = [0, 0, 0, 0, 0, 3, 5, 1, 4, 5, 8, 0, 8, 6, 9, 14, 0, 8, 13, 17, 15, 10, 3, 12, 14, 17, 19, 7, 23, 14, 0, 30]
F(1/32) = [0, 0, 0, 0, 0, 3, 4, 3, 4]
F(1/33) = [0, 0, 0, 0, 0, 3, 3, 5, 5, 7, 3, 7]
F(1/34) = [0, 0, 0, 0, 0, 3, 3, 1, 1, 7, 9, 4, 6, 4, 8, 3, 8, 8]
F(1/35) = [0, 0, 0, 0, 0, 3, 2, 4]
F(1/36) = [0, 0, 0, 0, 0, 3, 2]
F(1/37) = [0, 0, 0, 0, 0, 3, 1, 3, 1, 6, 5, 7, 5, 2, 6, 6, 7, 13, 5, 15, 18, 7, 20, 18, 15, 14, 4, 24, 22, 20, 11, 10, 28, 17, 28, 17, 0, 36]
F(1/38) = [0, 0, 0, 0, 0, 3, 0, 6, 5, 0, 4, 8, 1, 3, 5, 13, 6, 12, 9, 9]
F(1/39) = [0, 0, 0, 0, 0, 3, 0, 3, 1, 7, 6, 1, 8, 4]
F(1/40) = [0, 0, 0, 0, 0, 3]
F(1/41) = [0, 0, 0, 0, 0, 2, 5, 3, 7, 3, 7, 3, 5, 11, 1, 5, 7, 13, 12, 5, 11, 4, 13, 9, 12, 21, 24, 19, 21, 4, 28, 16, 20, 9, 22, 13, 23, 26, 6, 19, 0, 40]
F(1/42) = [0, 0, 0, 0, 0, 2, 5, 1]
F(1/43) = [0, 0, 0, 0, 0, 2, 4, 5, 1, 6, 0, 7, 8, 1, 2, 13, 15, 4, 6, 5, 6, 0, 21, 11, 5, 14, 13, 24, 13, 19, 16, 23, 2, 7, 22, 32, 20, 3, 16, 30, 33, 20, 0, 42]
F(1/44) = [0, 0, 0, 0, 0, 2, 4, 2, 4, 3, 2, 8]
F(1/45) = [0, 0, 0, 0, 0, 2, 4]
F(1/46) = [0, 0, 0, 0, 0, 2, 3, 4, 4, 4, 6, 10, 6, 3, 5, 7, 2, 13, 5, 9, 1, 15, 11, 11]
F(1/47) = [0, 0, 0, 0, 0, 2, 3, 2, 1, 7, 8, 5, 7, 5, 3, 8, 9, 14, 14, 17, 15, 15, 14, 0, 23, 12, 6, 2, 8, 9, 26, 5, 8, 28, 2, 31, 9, 35, 16, 6, 25, 21, 33, 41, 7, 22, 0, 46]
F(1/48) = [0, 0, 0, 0, 0, 2, 3]
F(1/49) = [0, 0, 0, 0, 0, 2, 2, 4, 6, 7, 7, 1, 6, 11, 2]
F(1/50) = [0, 0, 0, 0, 0, 2, 2, 2, 6, 3, 6]

[2..100] Kempner numbers: 2, 3, 4, 5, 3, 7, 4, 6, 5, 11, 4, 13, 7, 5, 6, 17, 6, 19, 5, 7, 11, 23, 4, 10, 13, 9, 7, 29, 5, 31, 8, 11, 17, 7, 6, 37, 19, 13, 5, 41, 7, 43, 11, 6, 23, 47, 6, 14, 10, 17, 13, 53, 9, 11, 7, 19, 29, 59, 5, 61, 31, 7, 8, 13, 11, 67, 17, 23, 7, 71, 6, 73, 37, 10, 19, 11, 13, 79, 6, 9, 41, 83, 7, 17, 43, 29, 11, 89, 6, 13, 23, 31, 47, 19, 8, 97, 14, 11, 10