I’ve been doing the Weekly
Challenges. The
latest
involved a Fibonacci-like word construction and a search for
square-free numbers. (Note that this is open until 6 February 2022.)
Task 1: Fibonacci Words
You are given two strings having same number of digits, $a
and $b
.
Write a script to generate Fibonacci Words
by concatenation of the
previous two strings. Finally print 51st digit of the first term
having at least 51 digits.
I wondered whether there might be a way of working out the relevant
character without actually assembling the strings, but the problem
requires the strings to be produced as well, so that's the way I did
it. I see no advantage to recursion here. Raku:
sub fibstr($aa,$bb,$limit) {
my $a=$aa;
my $b=$bb;
loop {
my $c=$a ~ $b;
say $c;
if (chars($c) >= $limit) {
return substr($c,$limit-1,1);
}
($a,$b)=($b,$c);
}
}
Some of the other languages I'm using distinguish between a 1-length
string and a single char
value.
Task 2: Square-free Integer
Write a script to generate all square-free integers <= 500.
A square-free integer (or squarefree integer) is defined as an
integer which is divisible by no perfect square other than 1. That
is, its prime factorization has exactly one factor for each prime
that appears in it. For example, 10 = 2 ⋅ 5
is square-free, but 18 = 2 ⋅ 3 ⋅ 3
is not, because 18 is divisible by 9 = 3².
This sequence is in the OEIS, though only
up to 113.
The approach I thought of first (there's another one later) is to look
at each number and factorise it. In Ruby:
def squarefree(max)
Initialise the output list.
out=[]
We'll prime-factor each candidate number, with short-cuts. Say I'm
considering candidate 500. Since I'm only interested in factors that
are square or higher powers of primes, the highest prime that
(squared) could be a factor is 19; 23² is 529. So I'll never need a
prime larger than the square root of the highest candidate.
primes=genprimes(Integer.sqrt(max))
And for each individual candidate, I don't need to look at primes higher
than its own square root.
plimit=1
Iterate over each candidate.
1.upto(max) do |candidate|
A number is square-free if I find no squared prime factors.
squarefree=true
Increase the limit until the square of the limit is no lower than the
candidate.
while plimit * plimit < candidate do
plimit += 1
end
Start factoring here.
cc=candidate
Check each prime in turn, and break out if we reach the end of the
list or exceed the limit.
primes.each do |pr|
if pr > plimit then
break
end
Count the number of times the number is divisible by this prime.
n=0
while cc % pr == 0 do
n += 1
If there's more than one, it's not square-free, so bail out.
if n > 1 then
squarefree=false
break
end
Otherwise, reduce the candidate and continue. We don't need to store
the prime factors, just confirm that none of them is a square.
cc=(cc/pr).to_i
end
If this isn't a square-free number, we don't need to test any more primes.
if !squarefree then
break
end
end
If we got through all the primes and found nothing to disqualify it,
add it to the output list.
if squarefree then
out.push(candidate)
end
end
Having tested all candidates, return the list of ones that qualified.
return out
end
Full code on
github.
As a bonus for blog readers, I also came up with another algorithm,
which starts with an output list containing just 1, then for each
prime up to max multiplies each entry in the list by that prime, and
stores that in the list if it's not too high. So the list starts as
(1)
; multiplying by 2 makes it (1, 2)
; multiplying by 3 makes it
(1, 2, 3, 6)
, etc. Because we never multiply by the same prime
twice, there can be no non-squarefree numbers in the result.
This doesn't form part of my solution submission, for reasons I'll
mention at the end, but it might be useful somewhere. Here's the Rust
version, relying on genprimes
as above.
fn squarefree(max: u32) -> Vec<u32> {
This time we take primes all the way up to max
(because the largest
prime within the limit is itself a squarefree number).
let primes = genprimes(max);
Initialise the output list. The great thing about a BTree is that,
while storing to it may be fractionally more expensive than simply
slapping a value onto the end of a list, I can cheaply get a list of
its members in ascending order. (In fact, I'm already using it in the
Rust version of the prime generator.)
let mut sf = BTreeSet::new();
sf.insert(1u32);
Loop over each prime, and over each entry in the output list. Note
that the list could be changed inside the loop, but in some languages
including Rust that will cause this algorithm to go awry; so instead I
build a list of values to be inserted, then insert them after the
inner loop is complete.
for pr in primes {
let mut ii=Vec::new();
for k in sf.iter() {
Generate a new solution number. If it's low enough, store it (and
we'll multiply it in future along with the other entries); if not,
bail out, because all the rest of the products in this inner loop will
also be too high. (Which is why we need the members of sf
in order,
because especially with later primes that'll be the majority of the
list.)
let y=k*pr;
if y <= max {
ii.push(y);
} else {
break;
}
}
Insert all the new values into the BTree. (Not that it matters, but
none of them will be there already as they're all products of a new
prime.)
if ii.len() > 0 {
for y in ii {
sf.insert(y);
}
}
}
Finally, return the BTree as a plain vector of numbers.
return sf.iter().map(|i| *i).collect::<Vec<u32>>();
}
This gets relatively faster with higher numbers; break-even time seems
to be with max
somewhere in the range 2^17-2^18, and by about 2^21
on my test machine it's reaching an answer twice as fast as the Rust
version of my first solution. However, it does rely for performance on
having a working BTree implementation; without that, I'd have to
re-sort the keys in sf
for each pass through the prime loop, and I
don't feel like writing a BTree for the languages that don't have it
readily available. (Though I'll probably do it eventually.)
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.