``````#!/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
# JugglerSteps
# Kolakoski
# LuckyNumbers
# MaxDigitCount
# Modulo
# Multiples
# NumAronson
# OEIS
# OEIS::Catalogue
# OEIS::Catalogue::Plugin
# Odd
# Palindromes
# Perrin
# PisanoPeriodSteps
# Polygonal
# Pronic
# ReReplace
# ReRound
# RepdigitAny
# Repdigits
# 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 {
return 0 unless is_prime(\$n);
my \$class = 1;
foreach my \$pe (factor_exp(\$n)) {
next if \$pe->[0] == 2 || \$pe->[0] == 3;
\$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;
}
``````