project-euler / project-euler / 152 / euler-152-take2.pl

#!/usr/bin/perl

use strict;
use warnings;

use v5.16;

use Math::BigRat only => 'GMP';
use List::MoreUtils qw(any all uniq);
use Math::BigInt (only => 'GMP');

sub slow_sq_frac
{
    my ($n) = @_;

    return Math::BigRat->new('1/' . ($n*$n));
}

sub is_prime
{
    my ($n) = @_;

    if ($n <= 1)
    {
        return 0;
    }

    my $top = int(sqrt($n));

    for my $i (2 .. $top)
    {
        if ($n % $i == 0)
        {
            return 0;
        }
    }

    return 1;
}

sub factorize
{
    my ($n) = @_;
    my @ret;

    my $factor = 2;
    while ($n > 1)
    {
        if ($n % $factor == 0)
        {
            push @ret, $factor;
            $n /= $factor;
        }
        else
        {
            $factor++;
        }
    }
    return \@ret;
}

my $target = Math::BigRat->new('1/2');
my $limit = 80;

my $found_count = 0;

# We exclude 2 because the target is divided by it.
my @primes = (grep { is_prime($_) } (3 .. $limit));

my @sq_fracs = (map { $_ ? slow_sq_frac($_) : $_ } (0 .. $limit));

my %primes_lookup = (map { $_ => 1 } @primes);

my @remaining_sums;

$remaining_sums[$limit+1] = 0;

for my $n (reverse(2 .. $limit))
{
    $remaining_sums[$n] = $remaining_sums[$n+1] + $sq_fracs[$n];
}

my @end_at;
$end_at[2] = $limit+1;
for my $p (@primes)
{
    $end_at[$p] = int($limit/$p) * $p + 1;
}

my @valid_primes_multiples_sets;

sub recurse
{
    my ($to_check, $so_far, $sum, $denom_factors_so_far) = @_;

    # print "Checking: ToCheck=@$to_check ; $sum+[@$so_far]\n";

    if ($sum == $target)
    {
        $found_count++;
        print "Found {", join(',', @$so_far), "}\n";
        print "Found so far: $found_count\n";
        return;
    }

    if (! @$to_check)
    {
        return;
    }

    my $start_from = $to_check->[0];

    if ($sum + $remaining_sums[$start_from] < $target)
    {
        # print "Remaining sum prune\n";
        return;
    }

    if (any { $_ > 2 && $end_at[$_] <= $start_from } @{factorize($sum->denominator())})
    {
        return;
    }

    FIRST_LOOP:
    foreach my $first_idx (keys(@$to_check))
    {
        my $first = $to_check->[$first_idx];
        my $new_sum = ($sum + $sq_fracs[$first])->bnorm();
        if ($new_sum > $target)
        {
            next FIRST_LOOP;
        }

        my @new_factors = grep { $_ > 2 } uniq(@{ factorize($new_sum->denominator()) });
        my $new_to_check = [@$to_check[$first_idx+1..$#$to_check]];

        if (! @new_factors)
        {
            recurse(
                $new_to_check,
                [sort { $a <=> $b } (@$so_far, $first)],
                $new_sum,
                $denom_factors_so_far,
            );
            next FIRST_LOOP;
        }

        my @new_factors_contains = (map {
            my $new_factor = $_;
            [
                grep
                { $new_to_check->[$_] % $new_factor == 0 }
                keys @$new_to_check
            ]
            }
            @new_factors
        );

        if (all { scalar(@$_) } @new_factors_contains)
        {
            my %new_factors_contains_lookup =
            (map { map { $_ => 1 } @$_ } @new_factors_contains);

            my @factors_not_contains = (grep
                { !exists($new_factors_contains_lookup{$_}) }
                keys @$new_to_check
            );

            my %encountered_factors;

            # my $sum_threshold = $target - $remaining_sums[$new_to_check->[0]];
            my $sum_threshold = $target;

            my $iter_factors_recurse = sub {
                my ($idx, $factor_idx, $factors_href, $new_new_sum, $new_denom_factors_so_far) = @_;

                if ($new_new_sum > $sum_threshold)
                {
                    return;
                }
                if ($idx == @new_factors)
                {
                    my @factors = sort {$a <=> $b } keys(%$factors_href);
                    if (! $encountered_factors{join(',',@factors)}++)
                    {
                        recurse([@$new_to_check[@factors_not_contains]],
                            [sort { $a <=> $b } (@$so_far, $first, @factors)],
                            $new_new_sum->bnorm(),
                            $new_denom_factors_so_far
                        );
                    }
                    return;
                }


                if ($factor_idx == @{$new_factors_contains[$idx]})
                {
                    my $denom = $new_new_sum->bnorm->denominator();
                    my $next_factors_so_far = {
                        %$new_denom_factors_so_far, $new_factors[$idx] => 1
                    };

                    # if (any { exists($next_factors_so_far->{$_}) } @{factorize($denom)})
                    if (any { $denom % $_ == 0 } keys(%$next_factors_so_far))
                    {
                        return;
                    }
                    else
                    {
                        return __SUB__->(
                            $idx+1, 0, $factors_href, $new_new_sum, $next_factors_so_far,
                        );
                    }
                }

                my $factor = $new_to_check->[$new_factors_contains[$idx][$factor_idx]];

                if (!exists($factors_href->{$factor}))
                {
                    __SUB__->($idx, $factor_idx+1, {%$factors_href, $factor => 1}, $new_new_sum + $sq_fracs[$factor], $new_denom_factors_so_far);
                }
                __SUB__->($idx, $factor_idx+1, $factors_href, $new_new_sum, $new_denom_factors_so_far);

                return;
            };

            $iter_factors_recurse->(0, 0, +{}, $new_sum, $denom_factors_so_far);
        }

        # recurse($first+1, [@$so_far, $first], $new_sum);
    }

    return;
}

# Filter out the large primes and the primes which only have 2*p in the limit
# and their product.
my @init_to_check = (grep { (!(exists($primes_lookup{$_}) || ((($_&0x1) == 0) && exists($primes_lookup{$_>>1})))) || $_ < $limit/3 } (2 .. $limit));

foreach my $p (@init_to_check)
{
    if ($p > 2 && (exists($primes_lookup{$p})))
    {
        print "Processing $p\n";
        my @products = (map { $p * $_ } (1 .. int($limit/$p)));

        my $lcm = Math::BigInt::blcm(map { $_ * $_ } @products);
        print "LCM=$lcm\n";
        my $mod = $lcm / $p;

        my @numers = map { $lcm / $_ } @products;
        my @combos;

        my @rev_sums;
        {
            my $s = 0;
            foreach my $n (reverse @numers)
            {
                $s += $n;
                unshift(@rev_sums, $s);
            }
        }

        my $recurse = sub {
            my ($i, $sum, $combo_so_far) = @_;
            # print "Combo=@$combo_so_far\n";
            if ($i == @products)
            {
                return;
            }
            my $m = $sum % $mod;
            if ($sum and $m == 0)
            {
                print "Found [@$combo_so_far]\n for $p\n";
                # We can combine distinct combos later.
                push @combos, $combo_so_far;
                return;
            }
            if ($rev_sums[$i+1] < $mod-$m)
            {
                # print "Flavoo = {$i}\n";
                return;
            }
            __SUB__->($i+1, $sum, $combo_so_far);
            __SUB__->($i+1, $sum+$numers[$i], [@$combo_so_far, $i]);

            return;
        };
        $recurse->(0, 0, []);

        $valid_primes_multiples_sets[$p] = \@combos;
    }
}

recurse([@init_to_check], [], Math::BigRat->new('0/1'), +{});



=begin removed

        foreach my $mask (1 .. ((1 << @products)-1))
        {
            # print "M=$mask\n";
            my $sum = 0;

            my @combo;

            foreach my $i (keys(@products))
            {
                if (($mask>>$i)&0x1)
                {
                    my $prod = $products[$i];
                    push @combo, $prod;
                    $sum += $sq_fracs[$prod];
                }
            }

            if ($sum->bnorm()->denominator() % $p)
            {
                push @combos, {map { $_ => 1 } @combo};
            }
        }

=end removed

=cut
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.