## Advent Calendar - December 18, 2022

Sunday, Dec 18, 2022| Tags: Perl, Raku

### |   Day 17   |   Day 18   |   Day 19   |

The gift is presented by `Ryan Thompson`. Today he is talking about his solution to “The Weekly Challenge - 171”. This is re-produced for Advent Calendar 2022 from the original post by him.

## PWC 171 › Odd Abundant Numbers

This post is part of a series on Mohammad Anwar’s excellent Weekly Challenge, where hackers submit solutions in Perl, Raku, or any other language, to two different challenges every week. (It’s a lot of fun, if you’re into that sort of thing.)

This week’s tasks include a simple number theory calculation, and a language feature. Here’s a look at `task 1`.

## Odd Abundant Numbers

`Abundant numbers` are numbers where the `sum of the proper divisors is greater than the number`. The first odd abundant number is `945`. 945’s proper divisors are `1, 3, 5, 7, 9, 15, 21, 27, 35, 45, 63, 105, 135, 189`, and `315`. (Recall that we exclude `945` itself as a proper divisor.)

The sum of those divisors is `975`, so therefore `945` is an abundant number.

Equivalently, a number is abundant if the sum of all divisors is greater than twice the number. We’ll see both variations in the solutions below.

I’m going to take you through a few different approaches. Since I don’t like repeating myself, let’s get some foreshadowing for `task #2` going, by using a `first class function`:

``````sub n_odd_abundant(&\$) {
my (\$is_abundant, \$N) = @_;

my @r;
for (my \$n = 3; \$N; \$n += 2) {
if (\$is_abundant->(\$n)) {
push @r, \$n;
\$N--
}
}

@r;
}
``````

The above code takes in a code ref (`\$is_abundant`) and a limit (`\$N`). We loop over all odd numbers, pushing any that pass the `\$is_abundant check` to our result.

Although the above version suited my purposes better, it’s also possible to do this with an iterator, to avoid having to store the intermediate result:

``````sub odd_abundant_iterator(&) {
my \$is_abundant = shift;
my \$n = 1;

sub {
do { \$n += 2 } until \$is_abundant->(\$n);

\$n;
}
}
``````

Now that we have a framework for gathering abundant numbers, let’s try it a few different ways.

## Brute force

The first way you might think to try is simply to brute force your way through every divisor of every number. This is `O(n)` for each number, and takes over a second to find the first `20 numbers`:

``````sub n_abundant_naive {
n_odd_abundant {
my \$n = shift;
\$n < sum grep { \$n % \$_ == 0 } 1..\$n/2;
} \$_[0];
}
``````

## Using `sqrt`

Stopping at `√n` perhaps unsurprisingly brings the asymptotic time down to `O(√n)`. Since divisors come in pairs, we can simply calculate the other divisor and avoid looping through most of the numbers:

``````sub n_abundant_sqrt {
n_odd_abundant {
my \$n    = shift;
my \$sqrt = sqrt(\$n);

my \$sum  = sum map { \$_,  \$n / \$_ }
grep { \$n % \$_ == 0 } 1..\$sqrt;

\$sum -= \$sqrt if \$sqrt == int \$sqrt;

2*\$n < \$sum;
} \$_[0];
}
``````

It might not seem like a huge change, but the above code runs about `27 times` faster than the naïve version, when asked to find the `first 20 odd abundant numbers`.

## Using `Math::Prime::Util` divisor_sum

Just for fun, our old friend, `Math::Prime::Util` has a function that seems perfect for our needs: divisor_sum. It does what it says on the tin: it calculates the sum of the divisors of whatever number we give it.

``````sub n_abundant_mpu {
n_odd_abundant {
my \$n = shift;
my \$sum = divisor_sum(\$n);

2*\$n < \$sum;
} \$_[0];
}
``````

This one is another `12 times` faster than the `sqrt` solution, and `356 times` faster than the `naïve method`. Great, right? Well, not so fast. Under the hood, the `divisor_sum` function is still finding all divisors for every number, so we’re still at `O(n√n)` time. It’s only faster because of the very tightly optimized `C code` under the hood.

## Sieve

We can still do quite a bit better by realizing that since we’re looking for a whole bunch of abundant numbers, we’re repeating the same calculations over and over again, for every multiple of a number we’ve seen already. So as long as we’re willing to tweak the requirement slightly to find all abundant numbers below a given limit (although we’ll see how we can still accommodate the old calling syntax if we really want to), we can do much better:

``````sub n_abundant_sieve {
my \$lim = shift;
my @r;

my @div_sum; # Sum of divisors for each number
for my \$n (1..\$lim) {
\$div_sum[\$n*\$_] += \$n for 1..\$lim/\$n+1;
push @r, \$n if \$n % 2 and 2*\$n <= \$div_sum[\$n];
}

@r;
}
``````

When called as `n_odd_abundant_sieve(10000)`, this returns the first `23` odd abundant numbers, about `twice` as fast as the `sqrt` version returns `20`. That’s because this algorithm runs in `O(n logn)` time, which is strictly less than `O(n√n)` time for all `n>0`. It grows significantly slower.

## Sneaking up on the limit

In number theory circles, finding all numbers below a limit is usually just fine, but if we wanted to be a stickler for the parameters of the challenge and return only the `first 20`, we could simply call `n_odd_abundant_sieve` multiple times, doubling the limit every time, until we have at least `20` results. I wouldn’t bother, though.

## You might be wondering, `“how does this perform compared to the MPU version?”`

For small values (the first 23 numbers in the sequence), the `MPU` version is `faster` by about `200%`. However, when given more realistic limits, our better algorithm pulls ahead. It breaks even at about `100` numbers on my machine. By the time we ask for `150` numbers, the `sieve` is already `60%` faster.

Benchmarks are an incredibly useful—and essential—tool when you need to write high performance code. Knowing that one algorithm is faster than another doesn’t always translate to real code. When dealing with implementation details, those hidden constants can really make a difference. What’s also really important, however, is knowing your requirements. If we’re only ever going to need to find less than `100 abundant primes`, the `MPU` version might be a worthy choice. However, if we want to push it beyond that, then the better algorithm wins.

Of course, translating the `sieve algorithm` to tight `C code` would blow both of them out of the water at any value of `n`, since the `sieve` version grows strictly slower than the `MPU` version.

## Raku

Here’s a quick port to `Raku`. If you’re only concerned with outputting the results, there’s no point in storing them first:

``````sub MAIN(Int \$lim = 10000) {
my @div_sum; # Sum of divisors for each number

for 1...\$lim -> \$n {
@div_sum[\$n*\$_] += \$n for 1..\$lim/\$n+1;
\$n.say if \$n % 2 and 2*\$n <= @div_sum[\$n];
}

}
``````

If you have any suggestion then please do share with us perlweeklychallenge@yahoo.com.

## SO WHAT DO YOU THINK ?

If you have any suggestions or ideas then please do share with us.