The Perl Toolchain Summit needs more sponsors. If your company depends on Perl, please support this very important event.
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 );
    my($l0, $l1);
    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;
}