This is an automated email from the git hooks/post-receive script.

ppm-guest pushed a commit to annotated tag v0.38
in repository libmath-prime-util-perl.

commit dd68fe4c47db0cb5330af8138e076f2a88b55a42
Author: Dana Jacobsen <d...@acm.org>
Date:   Mon Feb 3 13:19:24 2014 -0800

    Add inverse totient example
---
 MANIFEST                      |  1 +
 examples/README               | 12 ++++++
 examples/inverse_totient.pl   | 92 +++++++++++++++++++++++++++++++++++++++++++
 examples/project_euler_069.pl |  4 +-
 examples/project_euler_070.pl |  6 +--
 5 files changed, 110 insertions(+), 5 deletions(-)

diff --git a/MANIFEST b/MANIFEST
index 20849e5..2a680d3 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -62,6 +62,7 @@ examples/sophie_germain.pl
 examples/twin_primes.pl
 examples/abundant.pl
 examples/find_mr_bases.pl
+examples/inverse_totient.pl
 examples/parallel_fibprime.pl
 examples/porter.pl
 examples/verify-gmp-ecpp-cert.pl
diff --git a/examples/README b/examples/README
index 49d6706..d58f90b 100644
--- a/examples/README
+++ b/examples/README
@@ -43,6 +43,18 @@ porter.pl
     a(n) = m s.t. sigma(m) + sigma(m+1) + ... + sigma(m+n-1) is prime.
   Includes comparison to Pari/GP.
 
+inverse_totient.pl
+
+  Computes the image of phi(n) for a given m.  That is, given a number m,
+  the function computes all n where euler_phi(n) = m.  It returns just the
+  count in scalar context (which can be faster and lower memory for inputs
+  such as factorials that have huge images).
+
+project_euler_*.pl
+
+  Example solutions for some Project Euler problems.  If you participate
+  in PE, you really should solve the problems yourself first.  These
+  provide good examples how how to use some of the module functionality.
 
 
 verify-cert.pl
diff --git a/examples/inverse_totient.pl b/examples/inverse_totient.pl
new file mode 100644
index 0000000..f984707
--- /dev/null
+++ b/examples/inverse_totient.pl
@@ -0,0 +1,92 @@
+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;
+  }
+
+  # Based on invphi.gp v1.3 by Max Alekseyev
+  my @L;
+  foreach my $n (divisors($N)) {
+    $n = $do_bigint->new("$n") if $do_bigint;
+    my $p = $n+1;
+    next unless is_prime($p);
+    if ( ($N % $p) != 0 ) {
+      push @L, [ [$n, $p] ];
+    } else {
+      my ($v, $t) = (2, int($N/$p));
+      while (!($t % $p)) { $v++; $t = int($t/$p); }
+      push @L, [ [$n,$p], map { [$n*$p**($_-1), $p**$_] } 2..$v ];
+    }
+  }
+
+  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;
+        foreach my $n (divisors($N / $l0)) {
+          $t{$n*$l0} += $r{$n} if defined $r{$n};
+        }
+      }
+      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;
+}
diff --git a/examples/project_euler_069.pl b/examples/project_euler_069.pl
index 855897d..b126ed3 100644
--- a/examples/project_euler_069.pl
+++ b/examples/project_euler_069.pl
@@ -9,9 +9,9 @@ $n++ while pn_primorial($n+1) < 1000000;
 print pn_primorial($n), "\n";
 
 # Brute force
-my ($maxn, $maxratio) = (0, 0);
+my ($maxn, $maxratio, $ratio) = (0, 0);
 foreach my $n (1 .. 1000000) {
-  my $ratio = $n / euler_phi($n);
+  $ratio = $n / euler_phi($n);
   ($maxn, $maxratio) = ($n, $ratio) if $ratio > $maxratio;
 }
 print "$maxn  $maxratio\n";
diff --git a/examples/project_euler_070.pl b/examples/project_euler_070.pl
index f4a2636..d5a04c4 100644
--- a/examples/project_euler_070.pl
+++ b/examples/project_euler_070.pl
@@ -9,10 +9,10 @@ sub is_perm {
          join("",sort split(//,$a)) eq join("",sort split(//,$b));
 }
 
-my ($maxn, $minratio) = (0, 1000000);
+my ($maxn, $minratio, $totient, $ratio) = (0, 1000000);
 foreach my $n (2 .. 10_000_000) {
-  my $totient = euler_phi($n);
-  my $ratio = $n / $totient;
+  $totient = euler_phi($n);
+  $ratio = $n / $totient;
   ($maxn, $minratio) = ($n, $ratio) if $ratio < $minratio && is_perm($totient, 
$n);
 }
 print "$maxn  $minratio\n";

-- 
Alioth's /usr/local/bin/git-commit-notice on 
/srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git

_______________________________________________
Pkg-perl-cvs-commits mailing list
Pkg-perl-cvs-commits@lists.alioth.debian.org
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits

Reply via email to