RabbitFarm

2022-01-09

Sieve of Atkin / Curious Fraction Tree

The examples used here are from The Weekly Challenge problem statement and demonstrate the working solution.

Part 1

Write a script to generate the 10001st prime number.

Solution


use strict;
use warnings;

use boolean; 
use Getopt::Long;
use LWP::UserAgent;

use constant N => 10_001;   
use constant PRIME_URL => "http://primes.utm.edu/lists/small/100000.txt";

sub get_primes{
    my @primes;
    my $ua = new LWP::UserAgent(
        ssl_opts => {verify_hostname => 0}
    );
    my $response = $ua->get(PRIME_URL);
    my @lines = split(/\n/,$response->decoded_content);
    foreach my $line (@lines){
        my @p = split(/\s+/, $line);
        unless(@p < 10){
            push @primes, @p[1..(@p - 1)];
        }
    }
    return @primes;
}

sub sieve_atkin{
    my($n) = @_;
    my @primes = (2, 3, 5);
    my $upper_bound = int($n * log($n) + $n * log(log($n)));
    my @atkin = (false) x $upper_bound;    
    my @sieve = (1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59);
    for my $x (1 .. sqrt($upper_bound)){
        for(my $y = 1; $y <= sqrt($upper_bound); $y+=2){
            my $m = (4 * $x ** 2) + ($y ** 2);
            my @remainders;  
            @remainders = grep {$m % 60 == $_} (1, 13, 17, 29, 37, 41, 49, 53) if $m <= $upper_bound; 
            $atkin[$m] = !$atkin[$m] if @remainders; 
        }          
    } 
    for(my $x = 1; $x <= sqrt($upper_bound); $x += 2){
        for(my $y = 2; $y <= sqrt($upper_bound); $y += 2){
            my $m = (3 * $x ** 2) + ($y ** 2);
            my @remainders;  
            @remainders = grep {$m % 60 == $_} (7, 19, 31, 43) if $m <= $upper_bound; 
            $atkin[$m] = !$atkin[$m] if @remainders; 
        }          
    }   
    for(my $x = 2; $x <= sqrt($upper_bound); $x++){
        for(my $y = $x - 1; $y >= 1; $y -= 2){
            my $m = (3 * $x ** 2) - ($y ** 2);
            my @remainders;  
            @remainders = grep {$m % 60 == $_} (11, 23, 47, 59) if $m <= $upper_bound; 
            $atkin[$m] = !$atkin[$m] if @remainders; 
        }          
    } 
    my @m;
    for my $w (0 .. ($upper_bound / 60)){
        for my $s (@sieve){
            push @m, 60 * $w + $s;  
        }
    }
    for my $m (@m){
        last if $upper_bound < ($m ** 2);
        my $mm = $m ** 2;
        if($atkin[$m]){
            for my $m2 (@m){
                my $c = $mm * $m2;
                last if $c > $upper_bound;
                $atkin[$c] = false;
            }
        }
    }
    map{ push @primes, $_ if $atkin[$_] } 0 .. @atkin - 1;
    return @primes; 
}

sub get_nth_prime{
    my($n, $generate) = @_; 
    my @primes;
    unless($generate){
        @primes = get_primes;
    }
    else{
        @primes = sieve_atkin($n);
    }
    return $primes[$n - 1]; 
}


MAIN:{
    my $n = N;
    my $generate = false;
    GetOptions("n=i" => \$n, generate => \$generate);
    print get_nth_prime($n, $generate) . "\n"; 
}

Sample Run


$ perl perl/ch-1.pl
104743
$ perl perl/ch-1.pl --generate
104743
$ perl perl/ch-1.pl --generate
104743
$ perl perl/ch-1.pl --generate --n 101
547
$ perl perl/ch-1.pl --generate --n 11
31
$ perl perl/ch-1.pl --n 10001
104743
$ perl perl/ch-1.pl --n 11
31

Notes

I've mentioned it before, but for anything that asks for or needs prime numbers I always ust grab them from one of several convenient online sources, rather than generate them myself.

This time around I figured it'd be sporting to generate them myself, but maybe in an interesting way. Here I implement a sieve method for determining prime numbers. This Sieve of Atkin_ has a claim to fame of being the most performant among prime number generating sieve techniques. The code is a bit convoluted looking, I will admit, but is a faithful Perl representation of the algorithm (follow the reference link for pseudocode). Also, rather than try and explain the algorithm myself anyone interested can find full in depth treatments elsewhere. A background in number theory helps for some of the details.

Since I have some existing code for getting the pre-computed primes I figured I would use that as a check and extra feature. Command line options allow for the default behavior (fetch pre-computed primes for an N of 10,001) to be overridden.

Part 2

Given a fraction return the parent and grandparent of the fraction from the Curious Fraction Tree.

Solution


use strict;
use warnings;

use Graph;
use constant ROOT => "1/1";
use constant SEPARATOR => "/";

sub initialize{
    my($member) = @_;
    my $graph = new Graph();
    $graph->add_vertex(ROOT);
    my @next = (ROOT);
    my @changes = ([0, 1], [1, 0]);
    my $level = 0;
    {
        my @temp_next;
        my @temp_changes;
        do{
            $level++;
            my $next = shift @next;
            my($top, $bottom) = split(/\//, $next);
            my $change_left = shift @changes;
            my $change_right = shift @changes;
            my $v_left = ($top + $change_left->[0]) . SEPARATOR . ($bottom + $change_left->[1]);
            my $v_right = ($top + $change_right->[0]) . SEPARATOR . ($bottom + $change_right->[1]);    
            $graph->add_edge($next, $v_left);
            $graph->add_edge($next, $v_right);
            push @temp_next, $v_left, $v_right;
            push @temp_changes, $change_left;
            push @temp_changes, [$level + 1, 0], [0, $level + 1];
            push @temp_changes, $change_right;
        }while(@next && !$graph->has_vertex($member));
        @next = @temp_next;
        @changes = @temp_changes; 
        redo if !$graph->has_vertex($member);
    }
    return $graph;
}

sub curious_fraction_tree{
    my($member) = @_;
    my $graph = initialize($member);
    my($parent) = $graph->predecessors($member);
    my($grandparent) = $graph->predecessors($parent);
    return ($parent, $grandparent);
}

MAIN:{
    my($member, $parent, $grandparent);
    $member = "3/5";
    ($parent, $grandparent) = curious_fraction_tree($member);
    print "member = '$member'\n";
    print "parent = '$parent' and grandparent = '$grandparent'\n";
    print "\n";
    $member = "4/3";
    ($parent, $grandparent) = curious_fraction_tree($member);
    print "member = '$member'\n";
    print "parent = '$parent' and grandparent = '$grandparent'\n";
}

Sample Run


$ perl perl/ch-2.pl
member = '3/5'
parent = '3/2' and grandparent = '1/2'

member = '4/3'
parent = '1/3' and grandparent = '1/2'

Notes

My thought process on this problem started somewhat backwards. After reading the problem statement I thought of the Graph module and remembered that it defines a function predecessors() which would be very useful for this. After convincing myself to use Graph; I then probably spent the majority of the time for this just getting my head around how to define new vertices at each level of the tree. Like all trees there is some recursiveness to the structure, but an iterative implementation still looks clean as well.

Once the graph is constructed the solution as required comes from calling predecessors() to get the parent and grandparent vertices.

References

Challenge 146

Sieve of Atkin

Prime Pages

Graph

posted at: 17:32 by: Adam Russell | path: /perl | permanent link to this entry