RogerBW's Blog

I’ve been doing the Weekly Challenges. The latest involved prime-counting and fraction trees. (Note that this is open until 9 January 2022.)

Task 1: 10001st Prime Number

Write a script to generate the 10001st prime number.

The fast way of generating small primes, the sieve of Eratosthenes, requires an upper bound, so clearly what I need here is some sort of prime-counting function – i.e. something to give me a number that is guaranteed to be no lower than the 10001st prime. I thought about implementing the Meissel-Lehmer algorithm but it seemed too much like work, so instead I tweaked my existing sieve of Eratosthenes generator for better performance, and used a relation from that Wikipedia page:

P(sub n) < n log (n log n) where n ≥ 6

This isn't entirely true; it's proved up to n ≤ e^95, which is plenty for this purpose. The result originates in Rosser, Barkley (1941). "Explicit bounds for some functions of prime numbers". American Journal of Mathematics. 63 (1): 211–232. doi:10.2307/2371291. JSTOR 2371291. Which of course I can't legally read, because it's vitally important that everyone who might want to look at an eighty-year-old mathematical paper should first pay their two obols to Johns Hopkins. Well, at least it's not Elsevier.

Anyway, I'll do that, then use the result as an upper bound for my existing prime generator. Raku:

sub nthprime($n) {
    my $m=15;
    if ($n >= 6) {
        $m=floor(1+$n*log($n*log($n)));
    }
    my @primes=genprimes($m);
    return @primes[$n-1];
}

The prime generator is slightly tweaked from the version of two weeks ago to use the 6n±1 progression of prime candidates, which means a short FIFO queue filled as necessary. (I suspect in Rust this could be done with a custom iterator, but I just used VecDeque.)

sub genprimes($mx) {
    my @primes;
    {
        my $primesh=(2,3).SetHash;
        loop (my $i=6;$i <= $mx+1; $i += 6) {
            for ($i-1,$i+1) -> $j {
                if ($j <= $mx) {
                    $primesh{$j}=True;
                }
            }
        }
        my $p=2;
        my @q=[2,3,5,7];
        my $mr=floor(sqrt($mx));
        while ($p <= $mr) {
            if ($primesh{$p}:exists) {
                my $i=$p*$p;
                while ($i <= $mx) {
                    $primesh{$i}:delete;
                    $i += $p;
                }
            }
            if (@q.elems < 2) {
                @q.push(@q[*-1]+4);
                @q.push(@q[*-1]+2);
            }
            $p=@q.shift;
        }
        @primes=$primesh.keys.sort;
    }
    return @primes;
}

I have code for PostScript that ought to work, but ran off the end of GhostScript's ability to manage memory.

Task 2: Curious Fraction Tree

Consider the following Curious Fraction Tree:

(image elided, see below)

You are given a fraction, member of the tree created similar to the above sample.

Write a script to find out the parent and grandparent of the given member.

This is a Calkin-Wilf tree rather than the Stern-Brocot that puts all its values in order. Every language I'm using except for PostScript has something like object orientation, so in each of them I defined a Fraction class, with methods set_from_string to set up the fraction from a string value, stringify to convert it to a string value, and get_parent to return a new Fraction object that's a parent of the current one. If I were doing more with this, it would be easy to extend it: get_left_child, get_right_child, reduce, arithmetic operations, and so on. Raku does actually offer a Rat class, but it seemed easier to write my own than to extend theirs.

In Lua:

Fraction = {n = 1, b = 1}

function Fraction:new (o)
   o = o or {}
   self.__index = self
   setmetatable(o,self)
   return o
end

function Fraction:get_parent ()
   p=Fraction:new()
   p.n = self.n
   p.d = self.d
   if p.n < p.d then
      p.d = p.d - p.n
   else
      p.n = p.n - p.d
   end
   return p
end

function Fraction:stringify ()
   return string.format("%d/%d",self.n,self.d)
end

function Fraction:set_from_string (s)
   for n, d in string.gmatch(s,"(%d+)/(%d+)") do
      self.n=tonumber(n)
      self.d=tonumber(d)
   end
   return self
end

function p_gp(s)
   f=Fraction:new():set_from_string(s)
   out={}
   for i = 1,2 do
      f=f:get_parent()
      table.insert(out,f:stringify())
   end
   return out
end

In PostScript I move the set_from_string logic down inside p_gp and just work with a pair of numbers on the stack. (I have library functions to do things like pushing to an undefined-length array or string (apush and spush), or to convert an integer to a string without predefining the string's length (i2s).

/stringify {
    exch i2s 0 string exch spush
    (/) spush
    exch i2s spush
} bind def

/get_parent {
    2 copy
    gt {
        exch 1 index sub exch
    } {
        1 index sub
    } ifelse
} bind def

/p_gp {
    /out 0 array def
    (/) search {
        cvi
        exch pop
        exch cvi
        1 1 2 {
            pop get_parent 2 copy stringify out exch apush /out exch def
        } for
        out
    } {
        pop out
    } ifelse
} bind def

While it's not relevant here, the Fusc sequence that we saw in challenge #104 is the sequence of numerators (or denominators) in a breadth-first traversal of the tree.

And I'm adding JavaScript to the array of languages I'm doing this in. (Running under Node.js for command-line testing, but I'm trying to avoid anything that wouldn't run as easily in a browser.) Where will this madness end? They each have a reason for me to learn to use them, though, and solving these problems in multiple languages has definitely helped get me loose from the regexes-and-hashes approach that doing everything in Perl encourages.

Full code on github.