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 <d...@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 _______________________________________________ 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