``````use warnings;
use strict;
use Math::Prime::Util qw/:all/;
use Getopt::Long;

my %opts;
GetOptions(\%opts,
'count',
'help',
) || die_usage();
die_usage() if exists \$opts{'help'};
my \$n = shift;
die_usage() unless defined \$n && length(\$n) > 0 && \$n !~ tr/0123456789//c;
if (exists \$opts{'count'}) {
print scalar inverse_euler_phi(\$n), "\n";
} else {
print join("\n", inverse_euler_phi(\$n)), "\n";
}

sub die_usage {
die "Usage: \$0 [-count] <m>\n\nPrint all n such that euler_phi(n) = m.\nIf -count is used, just prints the number of such n.\n";
}

sub inverse_euler_phi {
my \$N = shift;

my \$do_bigint = (\$N > 2**49);
if (\$do_bigint) {
# Math::GMPz and Math::GMP are fast.  Math::BigInt::GMP is 10x slower.
eval { use Math::GMPz; \$do_bigint = "Math::GMPz"; 1; } ||
eval { use Math::GMP; \$do_bigint = "Math::GMP"; 1; } ||
eval { use Math::BigInt try=>"GMP,Pari"; \$do_bigint = "Math::BigInt"; 1; };
\$N = \$do_bigint->new("\$N");
}

return wantarray ? (1,2) : 2 if \$N == 1;
return wantarray ? () : 0 if \$N < 1 || (\$N & 1);
if (is_prime(\$N >> 1)) {   # Coleman Remark 3.3 (Thm 3.1) and Prop 6.2
return wantarray ? () : 0             if !is_prime(\$N+1);
return wantarray ? (\$N+1, 2*\$N+2) : 2 if \$N >= 10;
}

#if (!wantarray) { return a014197(\$N) }

# Based on invphi.gp v1.3 by Max Alekseyev
my @L;
fordivisors { \$n=\$_;
\$n = \$do_bigint->new("\$n") if \$do_bigint;
my \$p = \$n+1;
if (is_prime(\$p)) {
if ( (\$N % \$p) != 0 ) {
push @L, [ [\$n, \$p] ];
} else {
my \$v = valuation(\$N, \$p);
my \$t = \$N / \$p**\$v;
push @L, [ [\$n,\$p], map { [\$n*\$p**(\$_-1), \$p**\$_] } 2..\$v+1 ];
}
}
} \$N;

if (!wantarray) {   # Just count.  Much less memory.
my %r = ( 1 => 1 );
foreach my \$Li (@L) {
my %t;
foreach my \$Lij (@\$Li) {
my(\$l0, \$l1) = @\$Lij;
fordivisors {
\$t{\$_*\$l0} += \$r{\$_} if defined \$r{\$_};
} \$N / \$l0;
}
while (my(\$i,\$vec) = each(%t)) { \$r{\$i} += \$t{\$i}; }
}
return (defined \$r{\$N}) ? \$r{\$N} : 0;
}

my %r = ( 1 => [1] );
my(\$l0, \$l1);
foreach my \$Li (@L) {
my %t;
foreach my \$Lij (@\$Li) {
my(\$l0, \$l1) = @\$Lij;
foreach my \$n (divisors(\$N / \$l0)) {
push @{ \$t{\$n*\$l0} }, map { \$_ * \$l1 } @{ \$r{\$n} }
if defined \$r{\$n};
}
}
while (my(\$i,\$vec) = each(%t)) { push @{\$r{\$i}}, @\$vec; }
}
return () unless defined \$r{\$N};
delete @r{ grep { \$_ != \$N } keys %r };  # Delete all intermediate results
my @result = sort { \$a <=> \$b } @{\$r{\$N}};
return @result;
}

# Simple recursive count translated from Pari.
sub a014197 {
my(\$n,\$m) = @_;
\$m=1 unless defined \$m;
return 1+(\$m<2) if \$n == 1;
# TODO: divisor_sum with sub ought to be faster
#divisor_sum( \$n, sub { my \$d=shift;
#  return 0 if \$d < \$m || !is_prime(\$d+1);
#  my(\$p, \$q) = (\$d+1, \$n/\$d);
#  vecsum( map { a014197(\$q/(\$p**\$_), \$p) } 0 .. valuation(\$q,\$p) );
#} );
my(\$sum,\$p,\$q) = (0);
fordivisors {
if (\$_ >= \$m && is_prime(\$_+1)) {
(\$p,\$q)=(\$_+1,\$n/\$_);
\$sum += vecsum( map { a014197(\$q/(\$p**\$_), \$p) } 0 .. valuation(\$q,\$p) );
}
} \$n;
\$sum;
}
``````