RogerBW's Blog

 The Weekly Challenge 155: Pisano's Fortune 10 March 2022 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 { let mut f: HashMap = 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. Search 10 per page 20 per page 50 per page 100 per page Archive Tags 1920s 1930s 1940s 1950s 1960s 1970s 1980s 1990s 2000s 2010s 3d printing action advent of code aeronautics aikakirja anecdote animation anime army astronomy audio audio tech aviation base commerce battletech beer boardgaming book of the week bookmonth chain of command children chris chronicle church of no redeeming virtues cold war comedy computing contemporary cornish smuggler cosmic encounter coup covid-19 crime crystal cthulhu eternal cycling dead of winter doctor who documentary drama driving drone ecchi economics en garde espionage essen 2015 essen 2016 essen 2017 essen 2018 essen 2019 essen 2022 essen 2023 existential risk falklands war fandom fanfic fantasy feminism film firefly first world war flash point flight simulation food garmin drive gazebo genesys geocaching geodata gin gkp gurps gurps 101 gus harpoon historical history horror hugo 2014 hugo 2015 hugo 2016 hugo 2017 hugo 2018 hugo 2019 hugo 2020 hugo 2021 hugo 2022 hugo 2023 hugo 2024 hugo-nebula reread in brief avoid instrumented life javascript julian simpson julie enfield kickstarter kotlin learn to play leaving earth linux liquor lovecraftiana lua mecha men with beards mpd museum music mystery naval noir non-fiction one for the brow opera parody paul temple perl perl weekly challenge photography podcast politics postscript powers prediction privacy project woolsack pyracantha python quantum rail raku ranting raspberry pi reading reading boardgames social real life restaurant reviews romance rpg a day rpgs ruby rust scala science fiction scythe second world war security shipwreck simutrans smartphone south atlantic war squaddies stationery steampunk stuarts suburbia superheroes suspense television the resistance the weekly challenge thirsty meeples thriller tin soldier torg toys trailers travel type 26 type 31 type 45 vietnam war war wargaming weather wives and sweethearts writing about writing x-wing young adult Special All book reviews, All film reviews
Produced by aikakirja v0.1