[libmath-prime-util-perl] 31/33: Add Porter sequence example

[Partha P. Mukherjee]

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

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

commit 575a4ceed7cb116d6bec10182e87d67fe2a8c7e2
Author: Dana Jacobsen <dana at acm.org>
Date:   Sun Jan 26 01:43:32 2014 -0800

    Add Porter sequence example
---
 MANIFEST           |   1 +
 examples/README    |  11 +++++
 examples/porter.pl | 116 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 128 insertions(+)

diff --git a/MANIFEST b/MANIFEST
index 253dcd5..bedbeb0 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -63,6 +63,7 @@ examples/twin_primes.pl
 examples/abundant.pl
 examples/find_mr_bases.pl
 examples/parallel_fibprime.pl
+examples/porter.pl
 examples/verify-gmp-ecpp-cert.pl
 examples/verify-sage-ecpp-cert.pl
 examples/verify-cert.pl
diff --git a/examples/README b/examples/README
index 1a40ac1..49d6706 100644
--- a/examples/README
+++ b/examples/README
@@ -7,12 +7,14 @@ abundant.pl
   perl abundant.pl 100  deficient
   perl abundant.pl 15   perfect
 
+
 sophie_germain.pl
 
   Prints the first N Sophie-Germain primes.  E.g.:
 
   perl sophia_germain.pl 100000
 
+
 twin_primes.pl
 
   Prints the first N twin-primes (first value of the pair).  E.g.:
@@ -34,6 +36,15 @@ parallel_fibprime.pl
   Find Fibonacci primes, in parallel.  You will want Math::Prime::Util::GMP
   installed, as these are many-thousand-digit numbers.
 
+
+porter.pl
+
+  Various ways of constructing a sequence suggested by Michael B. Porter:
+    a(n) = m s.t. sigma(m) + sigma(m+1) + ... + sigma(m+n-1) is prime.
+  Includes comparison to Pari/GP.
+
+
+
 verify-cert.pl
 
   Takes an MPU or Primo primality certificate and verifies it.  This is
diff --git a/examples/porter.pl b/examples/porter.pl
new file mode 100755
index 0000000..3401e74
--- /dev/null
+++ b/examples/porter.pl
@@ -0,0 +1,116 @@
+#!/usr/bin/env perl
+use warnings;
+use strict;
+use 5.14.0;
+use Math::Prime::Util qw/:all/;
+use List::Util qw/sum/;
+use Benchmark qw/:all/;
+
+my $lim = shift || 1000;
+
+# Michael B Porter proposed this OEIS sequence:
+#
+#   a(n) = m such that sigma(m) + sigma(m+1) + ... + sigma(m+n-1) is prime
+#
+# http://oeis.org/wiki/User:Michael_B._Porter
+#
+# Charles R Greathouse IV suggested this as an efficient computation:
+#   a(n)=my(t=sum(i=1,n,sigma(i)),k=1);while(!isprime(t),t-=sigma(k)-sigma(n+k);k++);k
+# which can be turned into a vector as:
+#   vector(1000,i,a(i))
+#
+# Pari does this for 10k elements in ~15 seconds.
+# Version opt2 does it in Perl in 3.0s.
+# For 20k it's 63s in Pari, 12s in Perl.
+# Of course Pari could be optimized as well.
+
+sub simple {
+  my $lim = shift;
+  my @list;
+  foreach my $n (1 .. $lim) {
+    my($m, $sum) = (1, 0);
+    while (!is_prime($sum)) {
+      $sum = 0;
+      $sum += divisor_sum($m+$_) for 0..$n-1;
+      $m++;
+    }
+    push @list, $m-1;
+  }
+  return @list;
+}
+# perl -MMath::Prime::Util=:all -E 'my @list; foreach my $n (1 .. 1000) { my ($m,$sum) = (1,0); while (!is_prime($sum)) { $sum = 0; $sum += divisor_sum($m+$_) for 0..$n-1; $m++; } push @list, $m-1; } say join ",", @list;'
+
+sub crg4 {
+  my $lim = shift;
+  my @list;
+  foreach my $n (1 .. $lim) {
+    my($k, $t) = (1,0);
+    $t += divisor_sum($_) for 1..$n;
+    while (!is_prime($t)) {
+      $t -= divisor_sum($k)-divisor_sum($n+$k);
+      $k++;
+    }
+    push @list,$k;
+  }
+  return @list;
+}
+# perl -MMath::Prime::Util=:all -E 'my @list; foreach my $n (1 .. 10000) { my($k,$t)=(1,0); $t += divisor_sum($_) for 1..$n; while (!is_prime($t)) { $t -= divisor_sum($k)-divisor_sum($n+$k); $k++; } push @list, $k; } say join ",", @list;'
+# 9.8s for 10k
+
+sub opt1 {
+  my $lim = shift;
+  my @list = map {
+    my($n,$t,$k) = ($_,0,1);
+    $t += divisor_sum($_) for 1..$n;
+    while (!is_prime($t)) {
+      $t -= divisor_sum($k) - divisor_sum($n+$k);
+      $k++;
+    }
+    $k;
+  } 1 .. $lim;
+  return @list;
+}
+# perl -MMath::Prime::Util=:all -E 'say join ",", map { my($n,$t,$k) = ($_,0,1); $t += divisor_sum($_) for 1..$n; while (!is_prime($t)) { $t -= divisor_sum($k) - divisor_sum($n+$k); $k++; } $k; } 1 .. 10000'
+# 9.5s for 10k
+
+sub opt2 {
+  my $lim = shift;
+  my @ds;
+  my @list = map {
+    my($n,$t,$k) = ($_,0,1);
+    $ds[$n] //= divisor_sum($n);
+    $t += $ds[$_] for 1..$n;
+    while (!is_prime($t)) {
+      $ds[$n+$k] //= divisor_sum($n+$k);
+      $t -= $ds[$k] - $ds[$n+$k];
+      $k++;
+    }
+    $k;
+  } 1 .. $lim;
+  return @list;
+}
+# perl -MMath::Prime::Util=:all -E '@ds = (1,1); say join ",", map { my($n,$t,$k) = ($_,0,1); $t += $ds[$_] for 1..$n; while (!is_prime($t)) { $ds[$n+$k] //= divisor_sum($n+$k); $t -= $ds[$k] - $ds[$n+$k]; $k++; } $k; } 1..10000'
+# 3.0s for 10k
+
+# Verify
+{
+  my $vlim = 100;
+  my @a1 = simple($vlim);
+  my @a2 = crg4($vlim);
+  my @a3 = opt1($vlim);
+  my @a4 = opt2($vlim);
+  foreach my $i (0 .. $vlim-1) {
+    die "Mismatch in crg4 at $i" unless $a1[$i] == $a2[$i];
+    die "Mismatch in opt1 at $i" unless $a1[$i] == $a3[$i];
+    die "Mismatch in opt2 at $i" unless $a1[$i] == $a4[$i];
+  }
+}
+
+cmpthese(-5, {
+  #'simple' => sub { simple($lim) },
+  'crg4'   => sub { crg4($lim) },
+  'opt1'   => sub { opt1($lim) },
+  'opt2'   => sub { opt2($lim) },
+});
+
+#say join ", ", opt1($lim);

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