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

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

commit 9099e0ee85e001c840ba333c362398c5c385d2bd
Author: Dana Jacobsen <d...@acm.org>
Date:   Fri Apr 25 20:22:14 2014 -0700

    Add numseqs example file, first version
---
 MANIFEST            |   1 +
 TODO                |   2 +
 examples/numseqs.pl | 363 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 366 insertions(+)

diff --git a/MANIFEST b/MANIFEST
index d7af73a..a5ed401 100644
--- a/MANIFEST
+++ b/MANIFEST
@@ -63,6 +63,7 @@ examples/twin_primes.pl
 examples/abundant.pl
 examples/find_mr_bases.pl
 examples/inverse_totient.pl
+examples/numseqs.pl
 examples/parallel_fibprime.pl
 examples/porter.pl
 examples/verify-gmp-ecpp-cert.pl
diff --git a/TODO b/TODO
index f8132c7..cc53f13 100644
--- a/TODO
+++ b/TODO
@@ -72,3 +72,5 @@
 - We don't use legendre_phi for other functions any more, but it'd be nice
   to speed it up using some ideas from the Ohana 2011 SAGE branch.  For example
   (10**13,10**5) takes 2.5x longer, albeit with 6x less memory.
+
+- lucas_sequence with n = 0 or 1
diff --git a/examples/numseqs.pl b/examples/numseqs.pl
new file mode 100644
index 0000000..cb8f1ea
--- /dev/null
+++ b/examples/numseqs.pl
@@ -0,0 +1,363 @@
+#!/usr/bin/env perl
+use warnings;
+use strict;
+use Math::Prime::Util qw/:all/;
+use List::Util qw/sum/;
+use Math::BigInt try=>"GMP";
+
+# This shows examples of many sequences from:
+#   https://metacpan.org/release/Math-NumSeq
+# Some of them are faster, some are much faster, a few are slower.
+# This usually shows up once past ~ 10k values, or for large preds.
+# These also do not have the limit of 2^32 of most Math::NumSeq 
implementations.
+#
+# Note that this completely lacks the framework of the module, and Math::NumSeq
+# often implements various options that aren't always here.  It's just
+# showing some examples of using MPU to solve these sort of problems.
+
+# The lucas_sequence function covers about 45 different OEIS sequences,
+# including Fibonacci, Lucas, Pell, Jacobsthal, Jacobsthal-Lucas, etc.
+
+my $type = shift || 'AllPrimeFactors';
+my $count = shift || 100;
+my $arg = shift;  $arg = '' unless defined $arg;
+my @n;
+
+if      ($type eq 'Abundant') {
+# TODO
+
+} elsif ($type eq 'All') {
+  print join " ", 1..$count;
+} elsif ($type eq 'AllPrimeFactors') {
+  my $i = 2;
+  if ($arg eq 'descending') {
+    push(@n, reverse factor($i++)) while scalar @n < $count;
+  } else {
+    push(@n, factor($i++)) while scalar @n < $count;
+  }
+  print join " ", @n[0..$count-1];
+} elsif ($type eq 'AlmostPrimes') {
+  $arg = 2 unless $arg =~ /^\d+$/;
+  my $i = 1;
+  while (@n < $count) {
+    # use factor_exp for distinct
+    $i++ while scalar factor($i) != $arg;
+    push @n, $i++;
+  }
+  print join " ", @n;
+} elsif ($type eq 'Cubes') {
+  # Done via pred to show use
+  my $i = 0;
+  while (@n < $count) {
+    $i++ while !is_power($i,3);
+    push @n, $i++;
+  }
+  print join " ", @n;
+} elsif ($type eq 'DedekindPsiCumulative') {
+  my $c = 0;
+  print join " ", map { $c += psi($_) } 1..$count;
+} elsif ($type eq 'DivisorCount') {
+  print join " ", map { scalar divisors($_) } 1..$count;
+} elsif ($type eq 'DuffinianNumbers') {
+  my $i = 0;
+  while (@n < $count) {
+    $i++ while !is_duffinian($i);
+    push @n, $i++;
+  }
+  print join " ", @n;
+} elsif ($type eq 'PolignacObstinate') {
+  my $i = 1;
+  while (@n < $count) {
+    $i += 2 while !is_polignac_obstinate($i);
+    push @n, $i;
+    $i += 2;
+  }
+  print join " ", @n;
+} elsif ($type eq 'Emirps') {
+  my $i = 13;
+  while (@n < $count) {
+    $i = next_prime($i) while !is_prime(reverse $i) || $i eq reverse($i);
+    push @n, $i;
+    $i = next_prime($i);
+  }
+  print join " ", @n;
+} elsif ($type eq 'Fibonacci') {
+  # This is not a good way to do it, but does show a use for the function.
+  my $lim = ~0;
+  $lim = Math::BigInt->new(2) ** $count  if $count > 70;
+  print join " ", map { (lucas_sequence($lim, 1, -1, $_))[0] } 0..$count-1;
+} elsif ($type eq 'LemoineCount') {
+  print join " ", map { lemoine_count($_) } 1..$count;
+} elsif ($type eq 'LiouvilleFunction') {
+  print join " ", map { liouville($_) } 1..$count;
+} elsif ($type eq 'LucasNumbers') {
+  # This is not a good way to do it, but does show a use for the function.
+  my $lim = ~0;
+  $lim = Math::BigInt->new(2) ** $count  if $count > 91;
+  print join " ", map { (lucas_sequence($lim, 1, -1, $_))[1] } 1..$count;
+} elsif ($type eq 'MobiusFunction') {
+  print join " ", moebius(1,$count);
+} elsif ($type eq 'MoranNumbers') {
+  my $i = 1;
+  while (@n < $count) {
+    $i++ while !is_moran($i);
+    push @n, $i++;
+  }
+  print join " ", @n;
+} elsif ($type eq 'Pell') {
+  # This is not a good way to do it, but does show a use for the function.
+  my $lim = ~0;
+  $lim = Math::BigInt->new(3) ** $count  if $count > 51;
+  print join " ", map { (lucas_sequence($lim, 2, -1, $_))[0] } 0..$count-1;
+} elsif ($type eq 'PowerFlip') {
+  print join " ", map { powerflip($_) } 1..$count;
+} elsif ($type eq 'Primes') {
+  print join " ", @{primes($count)};
+} elsif ($type eq 'PrimeFactorCount') {
+  if ($arg eq 'distinct') {
+    print join " ", map { scalar factor_exp($_) } 1..$count;
+  } else {
+    print join " ", map { scalar factor($_) } 1..$count;
+  }
+} elsif ($type eq 'PrimeIndexPrimes') {
+  $arg = 2 unless $arg =~ /^\d+$/;
+  print join " ", map { primeindexprime($_,$arg) } 1..$count;
+} elsif ($type eq 'Primorials') {
+  print join " ", map { pn_primorial($_) } 0..$count-1;
+} elsif ($type eq 'SophieGermainPrimes') {
+  my $estimate = sg_upper_bound($count);
+  my $numfound = 0;
+  forprimes {  push @n, $_ if is_prime(2*$_+1);  } $estimate;
+  print join " ", @n[0..$count-1];
+} elsif ($type eq 'Squares') {
+  # Done via pred to show use
+  my $i = 0;
+  while (@n < $count) {
+    $i++ while !is_power($i,2);
+    push @n, $i++;
+  }
+  print join " ", @n;
+} elsif ($type eq 'Totient') {
+  print join " ", euler_phi(1,$count);
+} elsif ($type eq 'TotientCumulative') {
+  # pred:   sum(euler_phi(0,$_[0]));
+  my $c = 0;
+  print join " ", map { $c += euler_phi($_) } 0..$count-1;
+} elsif ($type eq 'TotientPerfect') {
+  my $i = 1;
+  while (@n < $count) {
+    $i += 2 while $i != totient_steps_sum($i,0);
+    push @n, $i;
+    $i += 2;
+  }
+  print join " ", @n;
+} elsif ($type eq 'TotientSteps') {
+  print join " ", map { totient_steps($_) } 1..$count;
+} elsif ($type eq 'TotientStepsSum') {
+  print join " ", map { totient_steps_sum($_) } 1..$count;
+} elsif ($type eq 'TwinPrimes') {
+  my $l = 2;
+  my $upper = 400 + int(1.01 * nth_twin_prime_approx($count));
+  $l=2; forprimes { push @n, $l if $l+2==$_; $l=$_; } $upper;
+  print join " ", @n[0..$count-1];
+} else {
+
+# The following sequences, other than those marked TODO, do not exercise the
+# features of MPU, hence there is little point reproducing them here.
+
+# Abundant             TODO
+# AlgebraicContinued
+# AllDigits
+# AsciiSelf
+# BalancedBinary
+# Base::IteraeIth
+# Base::IteratePred
+# BaumSweet
+# Beastly
+# Catalan
+# CollatzSteps
+# ConcatNumbers
+# CullenNumbers
+# DedekindPsiSteps     TODO
+# DeleteablePrimes     TODO
+# DigitCount
+# DigitCountHigh
+# DigitCountLow
+# DigitLength
+# DigitLengthCumulative
+# DigitProduct
+# DigitProductSteps
+# DigitSum
+# DigitSumModulo
+# ErdosSelfridgeClass  TODO
+# Even
+# Expression
+# Factorials
+# Fibbinary
+# FibbinaryBitCount
+# FibonacciRepresentations
+# FibonacciWord
+# File
+# FractionDigits
+# GolayRudinShapiro
+# GolayRudinShapiroCumulative
+# GoldbachCount        TODO
+# GolombSequence
+# HafermanCarpet
+# HappyNumbers
+# HappySteps
+# HarshadNumbers
+# HofstadterFigure
+# JugglerSteps
+# KlarnerRado
+# Kolakoski
+# LuckyNumbers
+# MaxDigitCount
+# MephistoWaltz
+# Modulo
+# Multiples
+# NumAronson
+# OEIS
+# OEIS::Catalogue
+# OEIS::Catalogue::Plugin
+# Odd
+# Palindromes
+# Perrin
+# PisanoPeriod
+# PisanoPeriodSteps
+# Polygonal
+# PowerPart          TODO
+# Powerful           TODO
+# PrimeIndexOrder    TODO
+# Pronic
+# ProthNumbers
+# PythagoranHypots
+# RadixConversion
+# RadixWithoutDigit
+# ReReplace
+# ReRound
+# RepdigitAny
+# RepdigitRadix
+# Repdigits
+# ReverseAdd
+# ReverseAddSteps
+# Runs
+# SelfLengthCumulative
+# SpiroFibonacci
+# SqrtContinued
+# SqrtContinuedPeriod
+# SqrtDigits
+# SqrtEngel
+# StarNumbers
+# SternDiatomic
+# Tetrahedral
+# Triangular
+# UlamSequence
+# UndulatingNumbers
+# WoodallNumbers
+
+  die "sequence '$type' is not implemented here\n";
+}
+print "\n";
+exit(0);
+
+
+# DedekindPsi
+sub psi { jordan_totient(2,$_[0])/jordan_totient(1,$_[0]) }
+
+sub is_duffinian {
+  my $n = shift;
+  return 0 if $n < 4 || is_prime($n);
+  my $dsum = divisor_sum($n);
+  foreach my $d (divisors($n)) {
+    return 0 unless $d == 1 || $dsum % $d;
+  }
+  1;
+}
+
+sub is_moran {
+  my $n = shift;
+  my $digsum = sum(split('',$n));
+  return 0 if $n % $digsum;
+  return 0 unless is_prime($n/$digsum);
+  1;
+}
+
+sub is_polignac_obstinate {
+  my $n = shift;
+  return (0,1,0,0)[$n] if $n <= 3;
+  return 0 unless $n & 1;
+  my $k = 1;
+  while (($n >> $k) > 0) {
+    return 0 if is_prime($n - (1 << $k));
+    $k++;
+  }
+  1;
+}
+  
+
+# Lemoine Count (A046926)
+sub lemoine_count {
+  my($n, $count) = (shift, 0);
+  return is_prime(($n>>1)-1) ? 1 : 0 unless $n & 1;
+  forprimes { $count++ if is_prime($n-2*$_) } $n>>1;
+  $count;
+}
+
+sub powerflip {
+  my($n, $prod) = (shift, 1);
+  # The spiffy log solution for bigints taken from Math::NumSeq
+  my $log = 0;
+  foreach my $pe (factor_exp($n)) {
+    my ($p,$e) = @$pe;
+    $log += $p * log($e);
+    $e = Math::BigInt->new($e) if $log > 31;
+    $prod *= $e ** $p;
+  }
+  $prod;
+}
+
+sub primeindexprime {
+  my($n,$level) = @_;
+  $n = nth_prime($n) for 1..$level;
+  $n;
+}
+
+# TotientSteps
+sub totient_steps {
+  my($n, $count) = (shift,0);
+  while ($n > 1) {
+    $n = euler_phi($n);
+    $count++;
+  }
+  $count;
+}
+
+# TotientStepsSum
+sub totient_steps_sum {
+  my $n = shift;
+  my $sum = shift;  $sum = $n unless defined $sum;
+  while ($n > 1) {
+    $n = euler_phi($n);
+    $sum += $n;
+  }
+  $sum;
+}
+
+# Sophie-Germaine primes upper bound.  Messy.
+sub sg_upper_bound {
+  my $count = shift;
+  my $nth = nth_prime_upper($count);
+  # For lack of a better formula, do this step-wise estimate.
+  my $estimate = ($count <   5000) ? 150 + int( $nth * log($nth) * 1.2 )
+               : ($count <  19000) ? int( $nth * log($nth) * 1.135 )
+               : ($count <  45000) ? int( $nth * log($nth) * 1.10 )
+               : ($count < 100000) ? int( $nth * log($nth) * 1.08 )
+               : ($count < 165000) ? int( $nth * log($nth) * 1.06 )
+               : ($count < 360000) ? int( $nth * log($nth) * 1.05 )
+               : ($count < 750000) ? int( $nth * log($nth) * 1.04 )
+               : ($count <1700000) ? int( $nth * log($nth) * 1.03 )
+               :                     int( $nth * log($nth) * 1.02 );
+
+  return $estimate;
+}

-- 
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