[libmath-prime-util-perl] 23/54: Add two more PE examples

[Partha P. Mukherjee]

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 684133940c30235ef6d09def31dbab0253725d56
Author: Dana Jacobsen <dana at acm.org>
Date:   Tue Feb 11 17:21:39 2014 -0800

    Add two more PE examples
---
 MANIFEST                      |   2 +
 examples/project_euler_214.pl |  21 +++++++++
 examples/project_euler_342.pl | 100 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 123 insertions(+)

diff --git a/MANIFEST b/MANIFEST
index 2a680d3..5b0e3b1 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -77,6 +77,8 @@ examples/project_euler_069.pl
 examples/project_euler_070.pl
 examples/project_euler_095.pl
 examples/project_euler_211.pl
+examples/project_euler_214.pl
+examples/project_euler_342.pl
 examples/project_euler_357.pl
 bin/primes.pl
 bin/factor.pl
diff --git a/examples/project_euler_214.pl b/examples/project_euler_214.pl
new file mode 100644
index 0000000..720ca6e
--- /dev/null
+++ b/examples/project_euler_214.pl
@@ -0,0 +1,21 @@
+#!/usr/bin/env perl
+use warnings;
+use strict;
+use Math::Prime::Util qw/forprimes euler_phi/;
+
+my $limit = shift || 40000000;
+my $cl = shift || 25;
+
+my @c;
+sub totchainlen {
+  my $n = shift;
+  return $n if $n <= 2;
+  $c[$n] //= 1 + totchainlen(euler_phi($n));
+  return $c[$n];
+}
+
+my $sum = 0;
+forprimes {
+  $sum += $_ if totchainlen($_) == $cl;
+} $limit;
+print "$sum\n";
diff --git a/examples/project_euler_342.pl b/examples/project_euler_342.pl
new file mode 100644
index 0000000..42051f9
--- /dev/null
+++ b/examples/project_euler_342.pl
@@ -0,0 +1,100 @@
+#!/usr/bin/env perl
+use warnings;
+use strict;
+use Math::Prime::Util qw/:all/;
+use Math::GMPz;
+
+# Sum of all n where is_power(euler_phi(n^2),3) = 1
+
+# Simple but very slow way.  The brute force method later in this file is
+# basically the same thing, but using the more efficient ranged moebius and
+# totient calls over intervals.
+#
+#  Pari:
+#    s=0; for(n=2,limit,if(ispower(n*eulerphi(n),3),s=s+n)); print(s)
+#  Perl/MPU:
+#    my $s=0;
+#    for my $n (2..$limit) { $s += $n if is_power($n*euler_phi($n),3); }
+#    say $s;
+#
+# TIMING:
+#              10^7    2*10^7   10^8     10^10
+# Clever       0.07s    0.09s   0.24s    5s
+# Brute        5.8s    11.0s   53.9s     5 hours
+# Simple MPU  13.9s    29.7s  192s       1-2 days?
+# Simple Pari 13.6s    33.4s  277s       5 days?
+#   
+
+my $limit = shift || 10**10-1;
+my $method = lc(shift || 'clever');
+die "Method must be 'clever' or 'brute'\n" unless $method =~ /^(clever|brute)$/;
+my $sum = 0;
+
+
+if ($method eq 'clever') {
+  # About 5 seconds for 10^10-1
+  my $sqlimit = Math::GMPz->new($limit) * $limit;
+  my $cblimit = int( ($limit*$limit) ** 0.3334 + 0.01 );
+  foreach my $k (2 .. $cblimit) {
+    next if $k & 1;
+    my($p, $e) = @{ (factor_exp($k))[-1] };
+    $e *= 3;
+    next unless $e & 1;
+    my $m = int($k / ($p ** int($e/3)));
+    $m **= 3;
+    next if $m % ($p-1);
+    $m = int($m / ($p-1));
+    my $n = $p ** (($e+1) >> 1);
+    next if $n >= $limit;
+    while ($m > 1) {
+      my ($p,$e) = @{ (factor_exp($m))[-1] };
+      last unless $e & 1;
+      last if $m % ($p-1);
+      $n *= $p ** (($e+1) >> 1);
+      last if $n >= $limit;
+      $m = int($m / ( ($p-1) * ($p**$e) ) );
+    }
+    if ($m == 1) {
+      #print "$n\n";
+      $sum += $n;
+    }
+  }
+} else {
+  # About 5 hours for 10^10-1
+  my $interval = 10_000_000;   # Window size for moebius/totient
+  #prime_precalc(10**9);        # Slightly faster ranged phi
+  my($beg,$end) = (0,0);
+  while ($beg < $limit) {
+    $end = $beg + $interval - 1;
+    $end = $limit if $end > $limit;
+    my $start = ($beg<2)?2:$beg;
+
+    my $glim = int(~0 / $end);
+    my @m = moebius($beg, $end);
+    my @t = euler_phi($beg, $end);
+
+    if ($end <= $glim) {   # Totient($n) * $n will always be < ~0
+      foreach my $n ($start .. $end) {
+        next unless $m[$n-$beg] == 0;
+        my $totn2 = $n * $t[$n-$beg];
+        if (is_power($totn2,3)) {
+          # print "$n\n";
+          $sum += $n
+        }
+      }
+    } else {
+      foreach my $n ($start .. $end) {
+        next unless $m[$n-$beg] == 0;
+        my $tot = $t[$n-$beg];
+        if ($tot <= $glim) {
+          print "$n\n" if is_power($n * $tot, 3);
+        } else {
+          $tot = Math::GMPz->new($n) * $tot;
+          print "$n\n" if Math::GMPz::Rmpz_perfect_power_p($tot) && is_power($tot,3);
+        }
+      }
+    }
+    $beg = $end+1;
+  }
+}
+print "$sum\n";

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