RogerBW's Blog

I’ve been doing the Weekly Challenges. The latest involved more number theory problems. (Note that this is open until 13 March 2022.)

Task 1: Fortunate Numbers

Write a script to produce first 8 Fortunate Numbers (unique and sorted).

An odd business, this. The fortunate number fN is the smallest number which, added to the product of the first N primes (pN#), produces a prime. The series isn't an ascending one, though there are some constraints that produce useful ending conditions. Specifically fN is always larger than the last prime in the sequence that makes up pn# - which gives me a place to start looking.

I borrow my prime-generating code from challenge #146, and primality-testing code from #154.

sub fortunate {
  my $ct=shift;
  my %o;
  my @ll;
  my $ph=1;

Arbitrarily, I'll try up to twice as many primes as I actually want answers. (This should be enough, but I can't prove it.)

  foreach my $p (@{genprimes(nthprimelimit($ct*2))}) {

If we have enough solutions, and the new one would be higher than any solution we have, leave.

    if (scalar keys %o >= $ct) {
      if ($p >= max(keys %o)) {
        last;
      }
    }

Get the next primorial.

    $ph *= $p;

Iterate through candidate numbers, testing for primality.

    my $l=$p+1;
    while (!isprime($l+$ph)) {
      $l++;
    }
    $o{$l}=1;

If we've got more than enough answers, sort them and prune the highest.

    if (scalar keys %o > $ct) {
      @ll=sort {$a <=> $b} keys %o;
      splice @ll,$ct;
      %o=map {$_ => 1} @ll;
    }
  }
  return \@ll;
}

This is all quite inefficient; I suspect a reverse BTree, with quick access to the maximum value, would do better.

Task 2: Pisano Period

Write a script to find the period of the 3rd Pisano Period.

where π(3) is the period of the Fibonacci sequence modulo 3. The whole set of periods is at the OEIS.

How do we find a period? The sequence must start with (0,1) and continue with (1) (as the Fibonacci sequence itself does), and eventually will come back to (0,1), so we can run through looking for the first recurrence of that pattern.

We can optimise* slightly by taking the prime factors of the period length. If we are looking for π(36), this can be regarded as two overlapping sequences, 2² and 3², and the result is the lcm of π(2²) and π(3²). (Not necessary for or even relevant to this specific problem of finding π(3), but I solved for the more interesting general case.)

* (For values of "optimise" which involve writing a lot more code: I imported my existing prime-generation routines, and wrote a prime-factoriser too, except in Ruby where it's all included. Probably many of these languages have existing libraries that would do the job too.)

This is where I got frustrated with Raku; it turns out that it interprets if ($n!=0) { as something other than if ($n != 0) { and then, having changed the value of $n, throws an error message from an entirely different part of the code. I don't regard this as reasonable behaviour for a modern language.

In Rust, some utility functions:

fn primefactor(n: usize) -> HashMap<usize, usize> {
    let mut f: HashMap<usize, usize> = HashMap::new();
    let mut m = n;

Start out using easy small primes.

    for p in [2, 3, 5, 7] {
        while m % p == 0 {
            m /= p;
            let en = f.entry(p).or_insert(0);
            *en += 1;
        }
    }

Because if I've reduced it, I don't need to generate such large primes for the rest of the job.

A further optimisation would generate primes only up to sqrt(m), and any remaining value of m at the end is itself a prime factor. Maybe next time there's a prime factorisation problem. (Hint from Future Roger: I'll do this in a few weeks' time.)

    if m > 1 {
        for p in genprimes(m).iter().skip(4) {
            while m % p == 0 {
                m /= p;
                let en = f.entry(*p).or_insert(0);
                *en += 1;
            }
            if m == 1 {
                break;
            }
        }
    }
    f

One of the standard GCD algorithms, and LCM as a corollary of it.

fn gcd(m0: usize, n0: usize) -> usize {
    let mut m = m0;
    let mut n = n0;
    while n != 0 {
        let t = n;
        n = m % n;
        m = t;
    }
    m
}

fn lcm(m: usize, n: usize) -> usize {
    m / gcd(m, n) * n
}

Rust doesn't have integer exponentiation; not sure why, but this is my solution. (Break down the power into a binary sequence, so that x⁶ would be calculated as x² × x⁴.)

fn pow(x0: usize, pow0: usize) -> usize {
    let mut x = x0;
    let mut pow = pow0;
    let mut ret = 1;
    while pow > 0 {
        if (pow & 1) == 1 {
            ret *= x;
        }
        x *= x;
        pow >>= 1;
    }
    ret
}

So to do the actual job we factorise the number, then for each distinct prime in the factorisation work out the length of that sequence; then combine them with lcm.

fn pisano(n: usize) -> usize {
    let mut a = 1;
    for (k, v) in primefactor(n).iter() {
        let nn = pow(*k, *v);
        let mut t = 1;
        let mut x = [1, 1];
        while x != [0, 1] {
            t += 1;
            x = [x[1], (x[0] + x[1]) % nn];
        }
        a = lcm(a, t);
    }
    a
}

Full code on github.