Attractive 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, Mohammad asks us to output a list of all Attractive Numbers between 1 and 50. Attractive numbers, as described by the Online Encyclopedia of Integer Sequences (OEIS) are:

Numbers with a prime number of prime divisors (counted with multiplicity)

OEIS Sequence A063989

The first numbers are 4, 6, 8, 9, 10, 12, 14, 15, 18, 20, 21, 22, 25, 26, 27, 28, 30, 32, 33, 34, 35, 38, 39, 42, 44, 45, 46, 48, 49, 50, …

This is a straightforward problem, but I’m not going to let that stop me from finding an interesting way to solve it. If I were looking for a sensible way to solve it, I’d do something like this:

use Math::Prime::Util ':all';
say for grep { is_prime( factor($_) ) } 1..50;

With only core Perl modules at our disposal, we’re on our own when it comes to prime number testing. I’ve already implemented more than my share of prime sieves, Miller-Rabin, and other tests for primality, so I opted for something that maybe no one else in the Challenge this week will be foolish enough to implement:

Wilson’s Theorem

One of the most basic theorems of prime numbers is Wilson’s Theorem, which states, simply, that a number n is prime iff \((n-1)! \equiv -1 \pmod n\).

Or, in more familiar Perl 5 syntax, $n is prime iff:

(factorial( $n - 1 ) + 1) % $n == 0

Wilson’s Theorem is almost never used as a test for prime numbers, since it scales very quickly thanks to the factorial, and much faster methods exist. However, it is an amazing bit of number theory with many applications, and that’s why I chose to showcase this theorem today.

In our case, one little tweak can remove an entire term from the complexity for free. Since we are checking consecutive primes between 1 and 50, we can avoid a whole lot of multiplying by building each n! up from the previous results. Otherwise, my code is about as literal of an implementation of Wilson’s as you can get.

I hereby present possibly the least efficient prime number checking function I have ever written:

# Wilson's theorem, generates primes < $N
# (or melts your CPU trying)
sub primes_to {
    use bigint;
    my $N = shift;
    my $fac = 1;
    my @r;
    for my $n (2..$N) {
        $fac *= $n - 1;

        push @r, $n unless ($fac + 1) % $n;
    }
    @r;
}

Complexity

(It isn’t good.)

It would be tempting to think that the algorithm has a big-O complexity of O(n), since for n primes we do n multiplications (and n modulo divisions). Unfortunately, the multiplications themselves are a factor here, since the length of the numbers increases with n, thanks to the factorial. After the numbers get too large for the CPU’s multiply instruction, we can no longer count on them being multiplied within one clock cycle. This happens at just 13! on 32-bit systems, and 21! on 64-bit systems. Most people could easily list primes below 21 without a computer.

bigint picks up the slack for us, so we need to figure out how bigint multiplies large numbers. After scouring the C and Perl sources, I believe that bigint multiplies large numbers in O(n²) time, or O(n^log₂3) ≈ O(n^1.585) time with GMP installed. I only had time for a cursory analysis, so this claim should have a big asterisk. I’d welcome some verification. Still, it’s a pretty good bet that GMP isn’t going to do any better than O(n^log₂3), and doing worse than O(n²) would be a challenge in itself.

The whole Wilson’s algorithm to check the numbers 1..n for primality is therefore between O(n^2.585) and O(n³). It does grow quickly, but at n = 50, we are still at some constant multiple of only 24,651 or 125,000 operations, which is nothing to a modern computer. Indeed, checking all numbers from 1..50 for primality using this method takes about 3ms on my VM.

Of course Math::Prime::Util‘s primes() runs about 1000x faster, even on that small input!

Full Solution

Once we can get a list of primes < 50, we can make a function to find prime factors. We need these in multiplicity. In other words, given 48, which is 2 × 2 × 2 × 2 × 3, we would expect to get a list containing the 5 elements, (2, 2, 2, 2, 3). That’s easy enough to implement:

@primes50 = primes_to(50);

sub prime_div_mult {
    my $n = shift;
    my @div;
    for my $div (@primes50) {
        last if $div > $n;
        while ($n % $div == 0) {
            $n /= $div;
            push @div, $div;
        }
    }

    @div;
}

Now all that’s left is to grep through the numbers 1..50 for the ones where prime_div_mult($_) (in scalar context) is itself a prime number.

my @primes50   = primes_to_N(50);
my %primes50   = map  { $_ => 1 } @primes50;

say for grep { $primes50{ prime_div_mult($_) } } 1..50;

Feel free to check out my full Perl 5 solution (GitHub).

Raku Version

Raku has a built-in is-prime() function that does what it says on the tin, so we’ll gleefully use that. There is no divisors/factors function, though, so we’ll have to build our own, just as we did with the Perl 5 solution.

We could do this as an ordinary sub, and the code would be almost the same. However, I’d like to show off an interesting but lesser-known Raku feature, MONKEY-TYPING. With this pragma, you can add methods to existing classes and grammars via the augment declarator.

Here is the divisors method:

use MONKEY-TYPING;

augment class Int {
    #| Return an Array of prime divisors (with multiplicity)
    method divisors returns Array {
        my @div;
        my $n = self;
        for (2..self.sqrt).grep: *.is-prime -> $div {
            while $n % $div == 0 {
                $n /= $div;
                @div.push: $div;
            }
        }
        @div;
    }
}

Warning: The docs are very vocal about augment being a bad idea, since augment allows you to change class behavior globally, in lexical scope. And I’d agree, usually. However, I’m mainly using this pragma to show off a unique language feature, and it would also be trivial to write divisors() as an ordinary function.

To find Attractive Numbers, we now need only grep: *.divisors.is-prime on whatever range we want! Our MAIN() sub looks like this:

sub MAIN( Int :$max = 50 ) {
    .say for (1..$max).grep: *.divisors.is-prime;
}

And, you can have a look at my full Raku solution on GitHub as well.

That’s it!

I hope you enjoyed this long read on a short problem. It was certainly a deep dive on a couple of tangential topics, like our analysis of the complexity of Wilson’s theorem, and the MONKEY-TYPING pragma in Raku. While my Perl and Raku solutions were both more or less useless for any sane practical application, my intent was to show some interesting concepts you might be able to apply to other problems.