[libmath-prime-util-perl] 07/55: Add last sequences to numseqs examples
Partha P. Mukherjee
ppm-guest at moszumanska.debian.org
Thu May 21 18:53:39 UTC 2015
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 c6dd262929c7f6b0d6b85479b433832a905a2d4a
Author: Dana Jacobsen <dana at acm.org>
Date: Mon Apr 28 14:57:18 2014 -0700
Add last sequences to numseqs examples
---
TODO | 2 +
examples/numseqs.pl | 205 ++++++++++++++++++++++++++++++++++++++++++++++------
2 files changed, 186 insertions(+), 21 deletions(-)
diff --git a/TODO b/TODO
index cc53f13..180aa39 100644
--- a/TODO
+++ b/TODO
@@ -74,3 +74,5 @@
(10**13,10**5) takes 2.5x longer, albeit with 6x less memory.
- lucas_sequence with n = 0 or 1
+
+- valuation
diff --git a/examples/numseqs.pl b/examples/numseqs.pl
old mode 100644
new mode 100755
index cb8f1ea..a35d731
--- a/examples/numseqs.pl
+++ b/examples/numseqs.pl
@@ -2,15 +2,21 @@
use warnings;
use strict;
use Math::Prime::Util qw/:all/;
-use List::Util qw/sum/;
+use List::Util qw/sum max/;
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.
-#
+
+# In general, these will work just fine for values up to 2^64, and typically
+# quite well beyond that. This is in contrast to most Math::NumSeq sequences
+# which limit themselves to 2^32 because Math::Factor::XS and Math::Prime::XS
+# use naive implementations which do not scale well.
+
+# The argument method is crippled here, just used for quick examples.
+
# 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.
@@ -18,14 +24,38 @@ use Math::BigInt try=>"GMP";
# The lucas_sequence function covers about 45 different OEIS sequences,
# including Fibonacci, Lucas, Pell, Jacobsthal, Jacobsthal-Lucas, etc.
+# These use the simple method of joining the results. For very large counts
+# this consumes a lot of memory, but is purely for the printing.
+
my $type = shift || 'AllPrimeFactors';
my $count = shift || 100;
my $arg = shift; $arg = '' unless defined $arg;
my @n;
if ($type eq 'Abundant') {
-# TODO
-
+ my $i = 1;
+ if ($arg eq 'deficient') {
+ while (@n < $count) {
+ $i++ while divisor_sum($i)-$i >= $i;
+ push @n, $i++;
+ }
+ } elsif ($arg eq 'primitive') {
+ while (@n < $count) {
+ $i++ while divisor_sum($i)-$i <= $i || abundant_divisors($i);
+ push @n, $i++;
+ }
+ } elsif ($arg eq 'non-primitive') {
+ while (@n < $count) {
+ $i++ while divisor_sum($i)-$i <= $i || !abundant_divisors($i);
+ push @n, $i++;
+ }
+ } else {
+ while (@n < $count) {
+ $i++ while divisor_sum($i)-$i <= $i;
+ push @n, $i++;
+ }
+ }
+ print join " ", @n;
} elsif ($type eq 'All') {
print join " ", 1..$count;
} elsif ($type eq 'AllPrimeFactors') {
@@ -56,6 +86,15 @@ if ($type eq 'Abundant') {
} elsif ($type eq 'DedekindPsiCumulative') {
my $c = 0;
print join " ", map { $c += psi($_) } 1..$count;
+} elsif ($type eq 'DedekindPsiSteps') {
+ print join " ", map { dedekind_psi_steps($_) } 1..$count;
+} elsif ($type eq 'DeletablePrimes') {
+ my $i = 0;
+ while (@n < $count) {
+ $i++ while !is_deletable_prime($i);
+ push @n, $i++;
+ }
+ print join " ", @n;
} elsif ($type eq 'DivisorCount') {
print join " ", map { scalar divisors($_) } 1..$count;
} elsif ($type eq 'DuffinianNumbers') {
@@ -65,14 +104,6 @@ if ($type eq 'Abundant') {
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) {
@@ -81,11 +112,28 @@ if ($type eq 'Abundant') {
$i = next_prime($i);
}
print join " ", @n;
+} elsif ($type eq 'ErdosSelfridgeClass') {
+ if ($arg eq 'primes') {
+ # Note we wouldn't have problems doing ith, as we have a fast nth_prime.
+ print "1" if $count >= 1;
+ forprimes {
+ print " ", erdos_selfridge_class($_);
+ } 3, nth_prime($count);
+ } else {
+ $arg = 1 unless $arg =~ /^-?\d+$/;
+ print join " ", map { erdos_selfridge_class($_,$arg) } 1..$count;
+ }
} 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 'GoldbachCount') {
+ if ($arg eq 'even') {
+ print join " ", map { goldbach_count($_<<1) } 1..$count;
+ } else {
+ print join " ", map { goldbach_count($_) } 1..$count;
+ }
} elsif ($type eq 'LemoineCount') {
print join " ", map { lemoine_count($_) } 1..$count;
} elsif ($type eq 'LiouvilleFunction') {
@@ -109,8 +157,42 @@ if ($type eq 'Abundant') {
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 '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 'PowerFlip') {
print join " ", map { powerflip($_) } 1..$count;
+} elsif ($type eq 'Powerful') {
+ my($which,$power) = ($arg =~ /^(all|some)?(\d+)?$/);
+ $which = 'some' unless defined $which;
+ $power = 2 unless defined $power;
+ my $i = 1;
+ if ($which eq 'some' && $power == 2) {
+ while (@n < $count) {
+ $i++ while moebius($i);
+ push @n, $i++;
+ }
+ } else {
+ my(@pe,$nmore);
+ $i = 0;
+ while (@n < $count) {
+ do {
+ @pe = factor_exp(++$i);
+ $nmore = scalar grep { $_->[1] >= $power } @pe;
+ } while ($which eq 'some' && $nmore == 0)
+ || ($which eq 'all' && $nmore != scalar @pe);
+ push @n, $i;
+ }
+ }
+ print join " ", @n;
+} elsif ($type eq 'PowerPart') {
+ $arg = 2 unless $arg =~ /^\d+$/;
+ print join " ", map { power_part($_,$arg) } 1..$count;
} elsif ($type eq 'Primes') {
print join " ", @{primes($count)};
} elsif ($type eq 'PrimeFactorCount') {
@@ -122,6 +204,15 @@ if ($type eq 'Abundant') {
} elsif ($type eq 'PrimeIndexPrimes') {
$arg = 2 unless $arg =~ /^\d+$/;
print join " ", map { primeindexprime($_,$arg) } 1..$count;
+} elsif ($type eq 'PrimeIndexOrder') {
+ if ($arg eq 'primes') {
+ print "1" if $count >= 1;
+ forprimes {
+ print " ", prime_index_order($_);
+ } 3, nth_prime($count);
+ } else {
+ print join " ", map { prime_index_order($_) } 1..$count;
+ }
} elsif ($type eq 'Primorials') {
print join " ", map { pn_primorial($_) } 0..$count-1;
} elsif ($type eq 'SophieGermainPrimes') {
@@ -165,7 +256,6 @@ if ($type eq 'Abundant') {
# 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
@@ -178,8 +268,6 @@ if ($type eq 'Abundant') {
# CollatzSteps
# ConcatNumbers
# CullenNumbers
-# DedekindPsiSteps TODO
-# DeleteablePrimes TODO
# DigitCount
# DigitCountHigh
# DigitCountLow
@@ -189,7 +277,6 @@ if ($type eq 'Abundant') {
# DigitProductSteps
# DigitSum
# DigitSumModulo
-# ErdosSelfridgeClass TODO
# Even
# Expression
# Factorials
@@ -201,7 +288,6 @@ if ($type eq 'Abundant') {
# FractionDigits
# GolayRudinShapiro
# GolayRudinShapiroCumulative
-# GoldbachCount TODO
# GolombSequence
# HafermanCarpet
# HappyNumbers
@@ -226,9 +312,6 @@ if ($type eq 'Abundant') {
# PisanoPeriod
# PisanoPeriodSteps
# Polygonal
-# PowerPart TODO
-# Powerful TODO
-# PrimeIndexOrder TODO
# Pronic
# ProthNumbers
# PythagoranHypots
@@ -265,6 +348,19 @@ exit(0);
# DedekindPsi
sub psi { jordan_totient(2,$_[0])/jordan_totient(1,$_[0]) }
+sub dedekind_psi_steps {
+ my $n = shift;
+ my $class = 0;
+ while (1) {
+ return $class if $n < 5;
+ my @pe = factor_exp($n);
+ return $class if scalar @pe == 1 && ($pe[0]->[0] == 2 || $pe[0]->[0] == 3);
+ return $class if scalar @pe == 2 && $pe[0]->[0] == 2 && $pe[1]->[0] == 3;
+ $class++;
+ $n = jordan_totient(2,$n)/jordan_totient(1,$n); # psi($n)
+ }
+}
+
sub is_duffinian {
my $n = shift;
return 0 if $n < 4 || is_prime($n);
@@ -323,6 +419,11 @@ sub primeindexprime {
$n;
}
+sub prime_index_order {
+ my $n = shift;
+ return is_prime($n) ? 1+prime_index_order(prime_count($n)) : 0;
+}
+
# TotientSteps
sub totient_steps {
my($n, $count) = (shift,0);
@@ -361,3 +462,65 @@ sub sg_upper_bound {
return $estimate;
}
+
+sub erdos_selfridge_class {
+ my($n,$add) = @_;
+ return 0 unless is_prime($n);
+ $n += (defined $add) ? $add : 1;
+ my $class = 1;
+ foreach my $pe (factor_exp($n)) {
+ next if $pe->[0] == 2 || $pe->[0] == 3;
+ my $nc = 1+erdos_selfridge_class($pe->[0],$add);
+ $class = $nc if $class < $nc;
+ }
+ $class;
+}
+
+sub abundant_divisors {
+ my($n,$is_abundant) = (shift, 0);
+ fordivisors {
+ $is_abundant = 1 if $_ > 1 && $_ < $n && divisor_sum($_)-$_ > $_;
+ } $n;
+ $is_abundant;
+}
+
+sub is_deletable_prime {
+ my $n = shift;
+ # Not deletable prime if n isn't itself prime
+ return 0 unless is_prime($n);
+ my $len = length($n);
+ # Length 1, return 1 because n is a prime
+ return 1 if $len == 1;
+ # Leading zeros aren't allowed, so check pos 1 specially.
+ return 1 if substr($n,1,1) != "0" && is_deletable_prime(substr($n,1));
+ # Now check deleting each other position.
+ foreach my $pos (1 .. $len-1) {
+ return 1 if is_deletable_prime(substr($n,0,$pos) . substr($n,$pos+1));
+ }
+ 0;
+}
+
+sub power_part {
+ my($n, $power) = @_;
+ return 1 if $power == 2 && moebius($n);
+ foreach my $d (reverse divisors($n)) {
+ if (is_power($d,$power)) {
+ return ($power == 2) ? int(sqrt($d))
+ : int($d ** (1.0/$power) + 0.000000001);
+# : int(Math::BigInt->new("$d")->broot($power)->bstr);
+ }
+ }
+ 1;
+}
+
+# This is simple and low memory, but not as fast as can be done with a prime
+# list. See Data::BitStream::Code::Additive for example.
+sub goldbach_count {
+ my $n = shift;
+ return is_prime($n-2) ? 1 : 0 if $n & 1;
+ my $count = 0;
+ forprimes {
+ $count++ if is_prime($n-$_);
+ } int($n/2);
+ $count;
+}
--
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