[libmath-prime-util-perl] 17/54: Add inverse totient example
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:52:08 UTC 2015
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 <dana at 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
More information about the Pkg-perl-cvs-commits
mailing list