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

Shlomi Fish a5bd21c 




















































































































































































































































Shlomi Fish 5ea826c 
Shlomi Fish a5bd21c 
Shlomi Fish 5ea826c 
Shlomi Fish a5bd21c 
Shlomi Fish 5ea826c 




Shlomi Fish a5bd21c 
Shlomi Fish 5ea826c 


Shlomi Fish a5bd21c 

Shlomi Fish e45f0bd 









Shlomi Fish 5ea826c 

Shlomi Fish a5bd21c 


Shlomi Fish 5ea826c 
Shlomi Fish a5bd21c 


Shlomi Fish e45f0bd 

Shlomi Fish a5bd21c 





Shlomi Fish 5ea826c 
Shlomi Fish e45f0bd 



Shlomi Fish a5bd21c 










































#!/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 @sq = (map { $_ * $_ } @products);

        my $lcm = Math::BigInt::blcm(@sq);
        print "LCM=$lcm\n";
        my $bef_mod = $lcm;
        while ($bef_mod % $p == 0)
        {
            $bef_mod /= $p;
        }

        my $mod = $lcm / $bef_mod;

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

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

        my $half_lcmed_limit = ($lcm>>1);

        my $recurse = sub {
            my ($i, $sum, $combo_so_far) = @_;
            # print "Combo=@$combo_so_far\n";
            if ($i == @sq)
            {
                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 or $sum >= $half_lcmed_limit)
            {
                # 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.