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.
Comments on this post are now closed. If you have particular grounds for adding a late comment, comment on a more recent post quoting the URL of this one.