#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util qw/:all/;
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/iths.
#
# For comparison, we can use something like:
# perl -MMath::NumSeq::Emirps -E 'my $seq = Math::NumSeq::Emirps->new; say 0+($seq->next)[1] for 1..1000'
# perl -MMath::NumSeq::Factorials -E 'my $seq = Math::NumSeq::Factorials->new; say join(" ",map { ($seq->next)[1] } 1..1000)' | md5sum
# In general, these will work just fine for values up to 2^64, and typically
# quite well beyond that. This is in contrast to many Math::NumSeq sequences
# which limit themselves to 2^32 because Math::Factor::XS and Math::Prime::XS
# do not scale well. Some other sequences such as Factorials and LucasNumbers
# are implemented well in Math::NumSeq.
# The argument method is really simple -- this is just to show code.
# 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.
# 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') {
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') {
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 'Catalan') {
# Done via ith. Much faster than MNS ith, but much slower than iterator
@n = map { binomial( $_<<1, $_) / ($_+1) } 0 .. $count-1;
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 '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') {
my $i = 0;
while (@n < $count) {
$i++ while !is_duffinian($i);
push @n, $i++;
}
print join " ", @n;
} elsif ($type eq 'Emirps') {
# About 15x faster until 200k or so, then exponentially faster.
my($i, $inc) = (13, 100+10*$count);
while (@n < $count) {
forprimes {
push @n, $_ if is_prime(reverse $_) && $_ ne reverse($_)
} $i, $i+$inc-1;
($i, $inc) = ($i+$inc, int($inc * 1.03) + 1000);
}
splice @n, $count;
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 'Factorials') {
print join " ", map { factorial($_) } 0..$count-1;
} elsif ($type eq 'Fibonacci') {
print join " ", map { lucasu(1, -1, $_) } 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') {
print join " ", map { liouville($_) } 1..$count;
} elsif ($type eq 'LucasNumbers') {
# Note the different starting point
print join " ", map { lucasv(1, -1, $_) } 1..$count;
} elsif ($type eq 'MephistoWaltz') {
print join " ", map { mephisto_waltz($_) } 0..$count-1;
} 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') {
print join " ", map { lucasu(2, -1, $_) } 0..$count-1;
} elsif ($type eq 'PisanoPeriod') {
print join " ", map { pisano($_) } 1..$count;
} 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') {
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 '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 'ProthNumbers') {
# The pred is faster and far simpler than MNS's pred, but slow as a sequence.
my $i = 0;
while (@n < $count) {
$i++ while !is_proth($i);
push @n, $i++;
}
print join " ", @n;
} elsif ($type eq 'PythagoreanHypots') {
my $i = 2;
if ($arg eq 'primitive') {
while (@n < $count) {
$i++ while scalar grep { 0 != ($_-1) % 4 } factor($i);
push @n, $i++;
}
} else {
while (@n < $count) {
$i++ while !scalar grep { 0 == ($_-1) % 4 } factor($i);
push @n, $i++;
}
}
print join " ", @n;
} 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 'SternDiatomic') {
# Slow direct way for ith value:
# vecsum( map { binomial($i-$_-1,$_) % 2 } 0..(($i-1)>>1) );
# Bitwise method described in MNS documentation:
print join " ", map { stern_diatomic($_) } 0..$count-1;
} elsif ($type eq 'Totient') {
print join " ", euler_phi(1,$count);
} elsif ($type eq 'TotientCumulative') {
# ith: vecsum(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.
# AlgebraicContinued
# AllDigits
# AsciiSelf
# BalancedBinary
# Base::IterateIth
# Base::IteratePred
# BaumSweet
# Beastly
# CollatzSteps
# ConcatNumbers
# CullenNumbers
# DigitCount
# DigitCountHigh
# DigitCountLow
# DigitLength
# DigitLengthCumulative
# DigitProduct
# DigitProductSteps
# DigitSum
# DigitSumModulo
# Even
# Expression
# Fibbinary
# FibbinaryBitCount
# FibonacciRepresentations
# FibonacciWord
# File
# FractionDigits
# GolayRudinShapiro
# GolayRudinShapiroCumulative
# GolombSequence
# HafermanCarpet
# HappyNumbers
# HappySteps
# HarshadNumbers
# HofstadterFigure
# JugglerSteps
# KlarnerRado
# Kolakoski
# LuckyNumbers
# MaxDigitCount
# Modulo
# Multiples
# NumAronson
# OEIS
# OEIS::Catalogue
# OEIS::Catalogue::Plugin
# Odd
# Palindromes
# Perrin
# PisanoPeriodSteps
# Polygonal
# Pronic
# RadixConversion
# RadixWithoutDigit
# ReReplace
# ReRound
# RepdigitAny
# RepdigitRadix
# Repdigits
# ReverseAdd
# ReverseAddSteps
# Runs
# SelfLengthCumulative
# SpiroFibonacci
# SqrtContinued
# SqrtContinuedPeriod
# SqrtDigits
# SqrtEngel
# StarNumbers
# Tetrahedral
# Triangular -stirling($_+1,$_) is a complicated solution
# UlamSequence
# UndulatingNumbers
# WoodallNumbers
# Xenodromes
die "sequence '$type' is not implemented here\n";
}
print "\n";
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);
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;
}
sub is_proth {
my $v = $_[0] - 1;
my $n2 = 1 << valuation($v,2);
$v/$n2 < $n2 && $v > 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;
}
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);
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;
}
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,\my $root)) {
return $root;
}
}
1;
}
# This isn't faster, but it was interesting.
sub mephisto_waltz {
my($n,$i) = (shift, 0);
while ($n > 1) {
$n /= 3**valuation($n,3);
$i++ if 2 == $n % 3;
$n = int($n/3);
}
$i % 2;
}
# 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;
}
sub pisano {
my $i = shift;
my @pe = factor_exp($i);
my @periods = (1);
foreach my $pe (@pe) {
my $period = $pe->[0] ** ($pe->[1] - 1);
my $modulus = $pe->[0];
{
my($f0,$f1,$per) = (0,1,1);
for ($per = 0; $f0 != 0 || $f1 != 1 || !$per; $per++) {
($f0,$f1) = ($f1, ($f0+$f1) % $modulus);
}
$period *= $per;
}
push @periods, $period;
}
lcm(@periods);
}
sub stern_diatomic {
my ($p,$q,$i) = (0,1,shift);
while ($i) {
if ($i & 1) { $p += $q; } else { $q += $p; }
$i >>= 1;
}
$p;
}