inverse_totient

[hvds]

Is there a link to a description of Max Alekseyev's algorithm somewhere? I can attempt to reverse-engineer it from the code, but a description would be better.

Browsing the code I see that the PP code does an unnecessary (and costly) grep to free memory, I suggest:

-    delete @r{ grep { $_ != $n } keys %r };  # Delete all intermediate results
-    return Mvecsort($r{$n});
+    my $rn = $r{$n};  %r = (); # Delete all intermediate results
+    return Mvecsort($rn);

The C code has this early assert:

  MPUassert(n <= UV_MAX/7.5, "inverse_totient_list n too large");

.. which seems a shame - the Perl code has no such restriction. I guess the issue is that some results in the list can exceed UV_MAX in this case - the assert message is not very clear about that. (53# > 2^64, \prod_{p in P, p < 53}{p/(p-1)} ~= 7.21)

[trizen]

For Max Alexeyev's algorithm, see:

• https://github.com/maxale/gpscripts/blob/main/invphi.gp
• https://cs.uwaterloo.ca/journals/JIS/VOL19/Alekseyev/alek5.html

[hvds]

For Max Alexeyev's algorithm, see:

* https://github.com/maxale/gpscripts/blob/main/invphi.gp

* https://cs.uwaterloo.ca/journals/JIS/VOL19/Alekseyev/alek5.html

Thanks! I can't read the gp code, but hopefully the combination of MPU code and the worked example of \phi^-1 in the paper should be enough for me to get a feel for how and why it works.

[trizen]

Here is a simplified version of the algorithm:

use utf8;
use 5.036;
use ntheory qw(is_prime divisors valuation);

binmode(STDOUT, ':utf8');

sub inverse_euler_phi($n) {

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

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

        is_prime($d + 1) || next;

        my %temp;
        foreach my $k (1 .. (valuation($n, $d + 1) + 1)) {

            my $u = $d * ($d + 1)**($k - 1);
            my $v = ($d + 1)**$k;

            foreach my $f (divisors($n / $u)) {
                if (exists $r{$f}) {
                    push @{$temp{$f * $u}}, map { $v * $_ } @{$r{$f}};
                }
            }
        }

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

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

foreach my $n(1..70) {
    if (my @inv = inverse_euler_phi($n)) {
        say "φ−¹($n) = [", join(', ', @inv), "]";
    }
}

Alternatively:

use utf8;
use 5.036;
use ntheory qw(:all);
binmode(STDOUT, ':utf8');

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}) {
                    my @list = @{$r{$d}};
                    push @{$t{mulint($x, $d)}}, map { mulint($_, $y) } @list;
                }
            }
        }
        while (my ($k, $v) = each %t) {
            push @{$r{$k}}, @$v;
        }
    }

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

sub cook_phi ($N) {
    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 { [mulint($d, powint($p, $_ - 1)), powint($p, $_)] } 1 .. $v + 1;
    }

    [values %L];
}

sub inverse_phi ($N) {
    dynamicPreimage($N, cook_phi($N));
}

foreach my $n (1 .. 70) {
    if (my @inv = inverse_phi($n)) {
        say "φ−¹($n) = [", join(', ', @inv), "]";
    }
}

[hvds]

Here is a simplified version of the algorithm

Thank you, I was able to work through it.

I observe that some of the work done, and data kept, is unnecessary. Using variable names from the PP code, if it is desirable to minimize peak memory usage we can discard $r{$n/$d} after processing $d for any $d > 1 instead of discarding them all at the end. Avoiding insertion of new useless values is slightly trickier - if I have it right, we can skip pushing $d2 * $dp when $d2 * $dp < $n && $d2 * $dp * $d > $n.

It isn't clear to me whether that would also give a speed improvement - probably yes for the C code, more debatable for PP.

[danaj]

I was using the Pari code as a guideline. [Completely different Pari code than the one linked: I used this from 2012: https://drememi.ludost.net/invphi2.gp]

I'm looking into being more exact with the limit in the XS code. If I remember right, memory was an issue in the PP code (at least for me with a 16GB machine so anything over 8GB will start swapping). It sounds like a good idea to look into being more aggressive with it, especially your ideas for proactively avoiding adding things.

[danaj]

The C code has this early assert:

  MPUassert(n <= UV_MAX/7.5, "inverse_totient_list n too large");

.. which seems a shame - the Perl code has no such restriction. I guess the issue is that some results in the list can exceed UV_MAX in this case - the assert message is not very clear about that. (53# > 2^64, \prod_{p in P, p < 53}{p/(p-1)} ~= 7.21)

The assert in the inverse_totient code ties into the XS function that shouldn't call the function if n is too large.

It is indeed an overflow concern.

(1) It looks like the limit it lower than it needs to be. In 32-bit it's much too low.

(2) We could test for overflow and return an overflow indicator only in the cases that actually overflow. A small amount more overhead. I'll do some benchmarks.

(3) The PP code is heavy handed in turning anything over 2^49 into a bigint, which will slow everything down. A fairly recent change is the ability to use Math::GMP, Math::GMPz, etc. instead of the default Math::BigInt -- that should make it run faster. I will also try using mulint instead of native multiply, which could turn out to be better. Maybe.

[danaj]

See https://arxiv.org/pdf/1002.1998 and https://arxiv.org/pdf/1207.4446 for more theoretical inequalities.

mpu '@a=inverse_totient(716636160); say $a[-1];'
4306662360
# This value is larger than 2^32

Given this input (2^16 * 3^7 * 5), the largest result overflows on a 32-bit machine. That's about UV_MAX/6, and seems to be the first value that triggers an overflow. Given this fairly small size for 32-bit, I can exhaustively check for the first value. If the guards are removed then as expected, a 32-bit build gives the wrong answer showing it overflows.

We can raise that limit slightly in most cases, as 1145325184 is the first 32-bit value that overflows when we have 16 or more divisors. There's probably more that could be done here.

For 64-bit, the first value I see is 2888172130711633920 which is ~ UV_MAX/6.4, but I'm sure there are smaller values (therefore a larger multiplier).

Inside setlist_addlist (keyval.h), we multiply the list by the given multiplier and that's where it overflows. That's the map { $_ * $pp } from the PP inner loop.

[danaj]

For inputs with many divisors, as hvds says, it uses lots of memory. Reducing that would be very helpful. Changing the code from using bigints throughout if n > 2^49 to instead using the addint/mulint/divint functions was a big improvement for me. Those smart operators keep most values as native integers, only using bigints when needed. That makes the operations faster as well.

The result is similar to what trizen's second code obtains, since it uses those operators also.

I did another small optimization, quite effective on some examples such as 2888172130711633920, where we do the divisors in reverse order, then on the final loop (divisor 1) only add values to @T if it will be part of the final result.

Unfortunately, inputs such as 2459565884805120000 or 2930663548728115200 will still blow it up to 50+GB.

[danaj]

@hvds I just pushed a change to inverse_totient. The XS code will return an overflow flag. That isn't a huge change but it's better than what was there before.

The big change is the PP code. It's not as simple as before, but it uses a lot less memory, though I wish we could trim it even more. As an example, for the input 2888172130711633920 we get a list with 23688815 values. Trizen's second program (with dynamicPreimage) uses 18GB and takes 5 minutes (in theory it should be closer to 2 minutes if I had more RAM). The updated code uses 3GB and takes 9 seconds. For the input 2772864899481600000, the new code uses 7GB and takes 22 seconds, while we used to use 32+GB. A few less obnoxious inputs are 8-9 times faster.

I have not changed the C code (yet). It's worth looking into.

[hvds]

Cool; I was thinking more on this earlier today, and I think that while the absolute numbers are interesting, it would better describe what is going on if we relate the peak memory use to the size of the results list (which is an unavoidable cost), perhaps by describing both in terms of how many values are involved.

For example, I think you're saying that inverse_totient_list(2888172130711633920) returns x=23688815 values, and has peak memory use of 3GB. How many values make up that 3GB? If it is of the approximate order of x (eg noticeably less than $x \ln x$) then that would suggest there is little more to gain without making it much slower.

(Note that I have not yet written any code actually using these functions - I was browsing the code during the process of thinking about a problem, and just wanted to capture my observations rather than forget them.)

(Note 2: omg, github has added tex support while I wasn't looking - I just put dollars around the x ln x to make it stand out, but it actually rendered it. That's going to make life so much easier.:)

[danaj]

I'll only use one example N value, the others are basically the same ratios.

N = 4311014400
has 135519 inverse totient values returned. All native integers on 64-bit.

Using trizen second program, should be identical results to what we had before:

1,841,396 intermediate values in %r (all of these to be deleted)
135,519 values in $r{$N}

For the PP code after my pushed changes:

118,007 intermediate values in %r (all of these to be deleted)
135,519 values in $r{$N}

I found a reasonably efficient way to remove the intermediate values after multiplying through at each level. When we finish the loop only the final list remains. It makes a meaningful difference in memory use on my system for a larger example (2772864768682229760 with 98,628,844 values returned).

[danaj]

(Note that I have not yet written any code actually using these functions - I was browsing the code during the process of thinking about a problem, and just wanted to capture my observations rather than forget them.)

I greatly appreciate it.

[hvds]

For the PP code after my pushed changes:

118,007 intermediate values in %r (all of these to be deleted) 135,519 values in $r{$N}

That seems like a perfectly reasonable state now.

[hvds]

In case it was of interest, here is a first pass at the totient-related problem I was thinking about. I plan to make a C-based version over the holiday period to see if I can find $n: a(n)=5$.

[danaj]

Running the program as-is (other than replacing H::Diag with printing), with the current code, I get to t=100M while under 20MB memory used. Unfortunately nothing interesting found after an hour getting to t=320M with about 25MB memory used. Only the same ones you've seen:

4 [2601630 2613510]
4 [32324236 32365816]
4 [61376256 61454592]

I'm not sure why you're getting OOM at only 32M.