bisect-g

#!/opt/maths/bin/perl
use strict;
use warnings;
use lib 'lib';
use List::Util qw{ sum };
use Type;
use Seq::Db;

my $typename = 't';
my($opt_g) = (0);
my @opt_m;
my $opt_M;

while (@ARGV && $ARGV[0] =~ /^-/) {
    my $arg = shift @ARGV;
    last if $arg eq '--';
    $typename = $1, next if $arg =~ /^-y(.*)/;
    $opt_g = 1, next if $arg eq '-g';
    push(@opt_m, $arg), next if $arg =~ /^-m/;
    $opt_M = $1, next if $arg =~ /^-M(\d+)$/;
    die "Unknown option '$arg'";
}
my $type = Type->new($typename);
option_g(@ARGV) if $opt_g;  # does not return
@ARGV == 5 or die "500 Usage: $0 n ming maxg d optc";
my($n, $ming, $maxg, $d, $optc) = @ARGV;

=head1 bisect-g: find upper bound for g(n)

Given inputs I<n>, I<ming>, I<maxg>, I<d>, I<optc>, attempts to refine
the upper bound I<maxg> for C<g(n)> by bisection.

I<ming> should be the best known lower bound for C<g(n)>; I<maxg> should
be the current best known upper bound for C<g(n)>; I<d> should be the
greatest difference checked to; and I<optc> is the requested number of
values to check.

The program will invoke C<gtauseq> to perform a binary chop, using options
C<< -c I<optc> -cp 20 -D >>, and report the new best upper bound.

=cut

my($good, $bad) = ($ming, $maxg + 1);
my @M = $opt_M ? (map "-m${opt_M}=$_", 0 .. $opt_M - 1) : (undef);
while ($good + 1 < $bad) {
    my $mid = ($good + $bad) >> 1;
    my @nextM = grep {
        local $opt_m[@opt_m] = $_ if defined $_;
        !bad($type, $mid);
    } @M;
    if (@nextM) {
        $good = $mid;
        @M = @nextM;
    } else {
        $bad = $mid;
    }
}
my $time = sum(times());
printf "200 g(%s) <= %s (%.2fs)\n",
        $n, $good, $time;
exit 0;

sub bad {
    my($type, $k) = @_;
    my $log = sprintf "%s/bg%s.%s-%s", $type->logpath, $n, $k, $$;
    my @args = (
        '-n', $d, '-x', $d, '-c', $optc, '-cp', 20, @opt_m, '-D',
        # force to avoid premature exit when dependent
        '-f',
        $n, $k,
    );
    my $lines = $type->invoked('gtauseq', "bgt($n, $k)", \@args, $log);
    # optionally parse '302' and/or '303' lines to refine -cp value
    if ($lines->{402}) {
        # 402 Error: all values (mod 2) disallowed (4.680s)
        return 1;
    }
    if ($lines->{403}) {
        # 403 Error: f(243) > 4 known impossible by exception (4.680s)
        return 1;
    }
    if ($lines->{404}) {
        # 404 Error: n + ${k}d must be divisible by n
        return 1;
    }
    if ($lines->{405}) {
        # 405 Error: Fixed power v_$k = ${x}y^$z is non-residue $v (mod $m)
        return 1;
    }
    if ($lines->{502}) {
        # Error: fixed 18 not available in tau 48
        return 1;
    }
    if ($lines->{309}) {
        # Prep finished, frequency 36.91 (184.980s)
        return 0;
    }
    # Not an expected result.
    die "501 Error parsing logs from '$log'\n";
}

# -g is a request to use database detail to decide parameters, and
# update database with results. We create a Run object to mediate
# the process, then reinvoke ourselves through it.
sub option_g {
    die "-m not supported with -g, only complete runs permitted"
            if @opt_m;
    @_ == 2 or die "Usage: $0 -g [options, ...] <n> <optc>\n";
    my($n, $optc) = @_;
    $type->bind_owner('harness');
    $type->bind($n);
    my $db = Seq::Db->new($type, 0);
    Seq::TauG->genTo($db, $n);
    my $g = $db->resultset('TauG')->find({ n => $n })
            // die "cannot find TauG entry for n=$n";
    my $f = $g->fnext($db)
            // die "cannot find next TauF entry for n=$n";
    my @opts = (
        ($opt_M ? "-M$opt_M" : ()),
        (map "-m$_", @{ $f->optm }),
    );
    my $r = Seq::Run::BisectG->new($type, $g, undef, $optc, \@opts);
    my @t0 = times;
    my $pid = $r->run;
    waitpid($pid, 0);
    my @t1 = times;
    my $time = $t1[2] + $t1[3] - $t0[2] - $t0[3];
    my $log = $r->logpath;
    system("cat $log");
    $r->finalize($db, $time);
    exit(0);
}