A090252

#!/opt/maths/bin/perl
use strict;
use warnings;
use Math::Prime::Util qw{ next_prime factor_exp };
use Math::GMP;

$| = 1;
my $style = shift(@ARGV) // 'default';

# Take care when multiplying or squaring that we transition to Math::GMP
# objects for any result that exceeds Perl's native integer precision.
my $MAX_INT = ~0;
my $MAX_SQ = int(sqrt($MAX_INT));

=head1 A090252

A090252 is the lexically first sequence of positive integers such that
each element a(n) is coprime to the next n elements.

The vast majority of the elements are the primes, in order; we therefore
record only the composite values.

Let m(n) be max(a(i): 1 <= i <= n). The list of disallowed primes is then
expected to consist of a) a handful of small primes (rarely if ever more
than 3 of them), b) a contiguous range of primes greater than sqrt(m(n)),
which have never been used in a composite value.

=cut

my %comp_index; # n => a(n) where a(n) composite
my %comp_value; # a(n) => 1 where a(n) composite
my $n = 1;
# The next prime that will be used
my $next_used = 2;
# The next prime that will become free
my $next_free = 2;

# The available small primes
my @small;
# Start and end of p > sqrt(m(n)) range, or undef if range is empty;
# these are all primes that have been used individually, but never as
# part of any composite.
my($min, $max);

sub NOVAL () { 0 }
sub TOOBIG () { 1 }
sub REUSED () { 2 }
sub GOOD () { 3 }

# Flags marking what we're about to use: S for one or more small primes,
# R for the first element of the range; U for the next unused prime;
sub S () { 1 } sub R () { 2 } sub U () { 4 }

my $best;

my $out = {
    default => sub {
        my($v, $freed) = @_;
        printf "%s\t%s\t%s\t%s\n", $n, $freed, $v, join(' ',
            grep length, join(',', @small), ($min ? "R$min..$max" : "")
        );
        return;
    },
    a090252 => sub {
        print "$n $_[0]\n";
        return;
    },
    a354164 => sub {
        my $r = $best->{free};
        return if $r & U;
        my $v = $_[0];
        my $c = (!($r & S) || (!($r & R) && @{$best->{small}} == 1))
            ? 1 : 2;
        print "$n $v $c\n";
        return;
    },
    lim => sub {
        my $v = $_[0];
        exit(0) if $v > 1e9;
        print "$n $_[0]\n" unless $best->{free} & U;
        return;
    },
    S => sub {
        use feature 'state';
        state $s = shift(@ARGV) // die "Usage: $0 S <n>";
        my $v = $_[0];
        return if $v % $s;
        print "$n $v\n";
        return;
    }
}->{$style} // die "Unknown style '$style'";


while (1) {
    ++$n;
    my $freed = ''; # for diagnostics
    if ($n & 1) {
        # factors of a(n/2) newly become available
        if (my $free = delete $comp_index{$n >> 1}) {
            # it's a composite
            my @f = map $_->[0], factor_exp($free);
            @small = sort { $a <=> $b } (@small, @f);
            $freed = join ',', @f;
        } elsif ($n > 3) {
            $min //= $next_free;
            $max = $next_free;
            $next_free = next_prime($next_free);
        }
    }

    # Given s < s' ... < sqrt(u) < m < u, we want the least unused of
    # s^k ...; ss' ...; sm ...; m^2; u and fail on s^2s'
    $best = { value => $next_used, free => U };

    try(square($min), R), goto done if @small == 0;
    trypow($small[$_], S, [ $_ ]) for 0 .. $#small;
    try(prod($small[0], $min), S | R, [ 0 ]);
    goto done if @small == 1;

    I: for my $i (0 .. $#small - 1) {
        for my $j ($i + 1 .. $#small) {
            my $t = try(prod($small[$i], $small[$j]), S, [ $i, $j ]);
            if ($t == TOOBIG || $t == GOOD) {
                next I if $j > $i + 1;
                last I;
            }
        }
    }

    if (try(prod($small[0], prod($small[0], $small[1])), S, [ 0, 1 ]) == GOOD) {
        fail("s^2s' $small[0], $small[1]");
    }

  done:
    # now act on best
    my $r = $best->{free};
    if ($r & R) {
        if ($min == $max) {
            ($min, $max) = (undef, undef);
        } else {
            $min = next_prime($min);
        }
    }
    if ($r & S) {
        splice(@small, $_, 1) for reverse @{ $best->{small} };
    }

    my $v = $best->{value};
    if ($r & U) {
        $next_used = next_prime($next_used) if $r & U;
    } else {
        fail("attempt to reuse comp $v") if $comp_value{$v};
        $comp_value{$v} = 1;
        $comp_index{$n} = $v;
    }
    $out->($v, $freed);
}

sub try {
    my($value, $free, $small) = @_;
    return NOVAL unless $value;
    return TOOBIG unless $value < $best->{value};
    return REUSED if $comp_value{$value};
    %$best = (value => $value, free => $free, small => $small);
    return GOOD;
}

sub fail {
    my($s) = @_;
    die sprintf "At %s, %s with small=(%s), range=%s\n",
            $n, $s, (join ',', @small), $min ? "$min..$max" : "";
}

sub square {
    my($n) = @_;
    return undef unless defined $n;
    return $n * $n if ref $n;
    return $n * $n if $n <= $MAX_SQ;
    return Math::GMP->new($n) * $n;
}

sub prod {
    my($x, $y) = @_;
    return undef unless defined($x) && defined($y);
    return $x * $y if ref($x) || ref($y);
    return $x * $y if $MAX_INT / $x > $y;
    return Math::GMP->new($x) * $y;
}

sub trypow {
    my($p, $free, $small) = @_;
    my $ilim = $MAX_INT / $p;
    my $n = $p;
    while (1) {
        $n = (($n < $ilim || $n > $MAX_INT) ? $n : Math::GMP->new($n)) * $p;
        my $t = try($n, $free, $small);
        return $t unless $t == REUSED;
    }
}