[feature request] Image of sigma() function

[trizen]

Hi. I would like to suggest a new function, called inverse_sigma(n,k=1), which returns a list with all the solutions x that satisfy sigma_k(x) = n.

Similar in functionality to inverse_totient(n).

Algorithm in Perl (based on invphi.gp ver. 2.1 by Max Alekseyev):

use 5.020;
use strict;
use warnings;

use ntheory qw(:all);
use experimental qw(signatures);

sub dynamicPreimage ($N, $L) {

    my %r = (1 => [1]);

    foreach my $l (@$L) {
        my %t;

        foreach my $pair (@$l) {
            my ($x, $y) = @$pair;

            foreach my $d (divisors(divint($N, $x))) {
                if (exists $r{$d}) {
                    push @{$t{mulint($x, $d)}}, map { mulint($_, $y) } @{$r{$d}};
                }
            }
        }
        while (my ($k, $v) = each %t) {
            push @{$r{$k}}, @$v;
        }
    }

    return if not exists $r{$N};
    sort { $a <=> $b } @{$r{$N}};
}

sub cook_sigma ($N, $k) {
    my %L;

    foreach my $d (divisors($N)) {

        next if ($d == 1);

        foreach my $p (map { $_->[0] } factor_exp(subint($d, 1))) {

            my $q = addint(mulint($d, subint(powint($p, $k), 1)), 1);
            my $t = valuation($q, $p);

            next if ($t <= $k or ($t % $k) or $q != powint($p, $t));

            push @{$L{$p}}, [$d, powint($p, subint(divint($t, $k), 1))];
        }
    }

    [values %L];
}

sub inverse_sigma ($N, $k = 1) {
    ($N == 1) ? (1) : dynamicPreimage($N, cook_sigma($N, $k));
}

say "solutions for sigma_1(x) = 5040:  ", "[", join(", ", inverse_sigma(5040)), "]";
say "solutions for sigma_2(x) = 22100: ", "[", join(", ", inverse_sigma(22100, 2)), "]";

Extra: algorithm for computing the image of the usigma(n) function (OEIS: A034448):

use 5.020;
use strict;
use warnings;

use ntheory qw(:all);
use experimental qw(signatures);

sub inverse_usigma ($n) {

    my %r = (1 => [1]);

    foreach my $d (divisors($n)) {

        my $D = subint($d, 1);
        is_prime_power($D) || next;

        my %temp;

        foreach my $f (divisors(divint($n, $d))) {
            if (exists $r{$f}) {
                push @{$temp{mulint($f, $d)}}, map { mulint($D, $_) }
                  grep { gcd($D, $_) == 1 } @{$r{$f}};
            }
        }

        while (my ($key, $value) = each(%temp)) {
            push @{$r{$key}}, @$value;
        }
    }

    return if not exists $r{$n};
    return sort { $a <=> $b } @{$r{$n}};
}

say "Solutions for usigma(x) = 5040: ", "[", join(', ', inverse_usigma(5040)), "]";

[danaj]

I wish there was a good way to generalize this. I can see how it can be done in C, and it would be fast, but it looks like 90% of the code would be identical to the inverse totient code.

[trizen]

By reusing the dynamicPreimage function, we have:

sub cook_phi ($N) {   # Euler phi
    my %L;

    foreach my $d (divisors($N)) {

        my $p = addint($d, 1);
        is_prime($p) || next;

        my $v = valuation($N, $p);

        push @{$L{$p}}, map {
            # [(p-1)*p^(k-1), p^k]
            [mulint(subint($p, 1), powint($p, $_ - 1)), powint($p, $_)]
        } 1 .. $v + 1;
    }

    [values %L];
}

sub inverse_phi ($N) {    # Inverse of Euler phi
    dynamicPreimage($N, cook_phi($N));
}

Similarly for the inverse of the Dedekind psi function:

sub cook_psi ($N) {     # Dedekind psi
    my %L;

    foreach my $d (divisors($N)) {

        my $p = subint($d, 1);
        is_prime($p) || next;

        my $v = valuation($N, $p);

        push @{$L{$p}}, map {
            # [(p+1)*p^(k-1), p^k]
            [mulint(addint($p, 1), powint($p, $_ - 1)), powint($p, $_)]
        } 1 .. $v + 1;
    }

    [values %L];
}

sub inverse_psi ($N) {     # inverse of Dedekind psi
    dynamicPreimage($N, cook_psi($N));
}

[trizen]

Optionally, the dynamicPreimage can also be extended to compute the inverse of some unitary functions, like uphi and usigma:

use 5.020;
use strict;
use warnings;

use ntheory qw(:all);
use experimental qw(signatures);

sub dynamicPreimage ($N, $L, %opt) {

    my %r = (1 => [1]);

    foreach my $l (@$L) {
        my %t;

        foreach my $pair (@$l) {
            my ($x, $y) = @$pair;

            foreach my $d (divisors(divint($N, $x))) {
                if (exists $r{$d}) {

                    my @list = @{$r{$d}};

                    if ($opt{unitary}) {
                        @list = grep { gcd($y, $_) == 1 } @list;
                    }

                    push @{$t{mulint($x, $d)}}, map { mulint($_, $y) } @list;
                }
            }
        }
        while (my ($k, $v) = each %t) {
            push @{$r{$k}}, @$v;
        }
    }

    return if not exists $r{$N};
    sort { $a <=> $b } @{$r{$N}};
}

sub cook_usigma ($N) {
    my %L;

    foreach my $d (divisors($N)) {

        my $p = subint($d, 1);
        is_prime_power($p) || next;

        push @{$L{$p}}, [$d, $p];    # [p^k+1, p^k]
    }

    [values %L];
}

sub cook_uphi ($N) {
    my %L;

    foreach my $d (divisors($N)) {

        my $p = addint($d, 1);
        is_prime_power($p) || next;

        push @{$L{$p}}, [$d, $p];    # [p^k-1, p^k]
    }

    [values %L];
}

sub inverse_usigma ($N) {
    dynamicPreimage($N, cook_usigma($N), unitary => 1);
}

sub inverse_uphi ($N) {
    dynamicPreimage($N, cook_uphi($N), unitary => 1);
}

say "Solutions for usigma(x) = 5040: ", "[", join(', ', inverse_usigma(5040)), "]";
say "Solutions for uphi(x) = 5040:   ", "[", join(', ', inverse_uphi(5040)),   "]";