Primorials and Kronecker products

Primorials

I am always amazed at the way mathematicians find a meaning or a use for obscure sequences and mathematical operations. And here's another one: primorial P(n) is defined as the product of the first n prime numbers. That is true whether or not you accept 1 as the first prime, but apparently P(0) is 1 which neatly sidesteps that argument.

It's easy to see that successive primorials increase fast. We are asked to generate the first 10, and P(9) is already over 223 million, and after that we are quickly into BigInt territory.

As this is a relatively simple task I chose to generate the primes rather than using a module and my algorithm looks like this:

k = 1
Loop j from 2 until we've found P(9)
    Next j if j is divisible by any prime we've found so far
    So j is prime, save it
    k = k * j, and k is the next primorial

I'm sure someone will come with a one-liner for this, but as always I tens to prioritise lucidity over brevity.

Kronecker products

You come across Leopold Kronecker's musings all over the place, and I see that at least a dozen mathematical 'things' are named after him.  Fortunately, this one isn't too hard to follow, thanks to someone giving the key to solving it in Wikipedia, albeit in very small print.

As with some other recent challenges, it takes more lines code to display the data neatly  than to perform the required calculation.

I used two helper functions. First is quorem which takes two arguments, a dividend and a divisor, and returns the quotient and the remainder as per primary school division.  The second is show which simply displays a matrix with the columns neatly lined up. I use that to display the two input matrices (A and B) and the resulting product matrix (K).

Wikipedia kindly gives us the formula for each element K(i, j) as a function of two elements of A and B. The calculation of which elements of A and B are needed involves some modular arithmetic, which is where quorem comes in.

The show function makes two passes over the matrix to be displayed. In the first pass it determines which element is the widest, for which Perl's length function is useful as it intrinsically copes with minus signs. The cell width is then made one space wider than the widest number for a pleasing result.

I chose to store the 3 matrices as 2-dimensional arrays - actually arrays of references, so that element X(i, j) is $X->[$i]->[$j]. This nicely parallels the mathematical representation.

Brilliant and Achilles

Brilliant numbers

A brilliant number has exactly two prime factors, and they both have the same number of digits.

That sounds like an easily met pair of criteria.  For example, there are 21 2-digit prime numbers, yielding 231 brilliant numbers ranging from 441 to 9409, so they're not rare.

The ever-helpful Math::Prime::Util has a function factor which returns a list of prime factors of its argument. so the algorithm looks like this:

loop over test from 1 to something big until we've found 20 brilliant numbers
    get the prime factors (pf) of test
    skip to the next test unless pf contains exactly 2 same-length members 
    test is brilliant!
end

That generates the desired 20 numbers in microseconds, 2000 in under a second and 20 000 in under 15 seconds on my lowly Raspberry Pi.  The 20 000th brilliant number is 2 733 299, which is 1117 x 2447.

Achilles numbers

Write a script to generate first 20 Achilles Numbers: numbers that are powerful but imperfect. A positive integer n is a powerful number if, for every prime factor p of n, p^2 is also a divisor. A number is imperfect if it has at least two distinct prime factors.

Mohammad doesn't mention the restriction that an Achilles number also cannot be a perfect power (ie a square, cube etc), although the Wikipedia article he refers to does mention that.

So, my algorithm looks like this:

loop over test from 1 to something big until we've found 20 Achilles numbers

    # check for powerful
    for each prime factor (p) of test
        if test isn't a multiple of p^2 -- no good
        while we're here, count how many times (seen(f)) p is a factor f

    # check for imperfection
    if any seen(f) < 2 -- no good

    # check for not being a perfect power
    find the least frequently occurring prime factor: say it occurs m times
    check that all the values in seen are *not* multiples of m

    ... and if that is the case, test is an Achilles number

The 'not a perfect power' criterion is quite tricky, although for the current exercise it could be done simply by calculating all the squares, cubes etc less than the number under test.  To see why my way works, consider a number n with prime factors:

n = a^z * b^y * c^x * d^w

where a, b, c etc are primes, and they are ordered such that z <= y <= x <= w.  Let's suppose that y, x and w are all multiples of z, so that:

n = a^z * b^2z * c^6z * d^8z 

We have already checked n for imperfection, so z must be at least 2. But if z is 2, then:

n = (a * b^z * c^3z * d^4z) ^2

and n is therefore a perfect square.  Similarly, if z = 3 and y, x, w etc are all multiples of 3, then n is a perfect cube.  And so on.

So the 'not a perfect power' rule boils down to any of the coefficients - ie the number of times each prime factor occurs - not being a multiple of the number of times the least frequent factor occurs.

Having worked all that out, the algorithm is quite fast, yielding 1000 Achilles numbers in under 10 seconds. The 1000th is 1 151 064, which is 2^3 * 3^3 * 73^2.

More funny numbers ... and a very big one

Perrin primes

The Perrin sequence is defined to start with 3, 0, 2 and after that, term n is the sum of terms n - 2 and n - 3. So far as I can see this sequence has no obvious practical use, but that's what pure mathematicians like.  Perrin primes are then those terms in the sequence which are prime numbers.

Looking at the supplied value of f(13) it seems wise to invoke our old friend Math::BigInt, and the only other quirk is that some early numbers - such as 5 - appear more than once so duplicates need to be avoided.

f(23) is (allegedly):
3631640163992448158050321979101634523424467070532989940589376200895999542521324121865744873084026078365592113103829057044319371093458267314150632770247926037880226504980936257910602481948018841362454143562440537190514898173776176693598426395086189616722660098879586330664613823090197360409779437591689520837492830513163054777061491401259817572546197753109857199993236881971656255401039799820579630315398215866349742611432329007353997352494443986017317922833363523351835711663212252398827126207580953779469798265218514506497114067477064259789799733135324524166520280952689291443318735365943242441087374207019201381566622887361047383284786893087439845660097204995566088460835424395601898782600822606786314286293

... and beyond that is_prime begins to struggle.

Home primes

The home prime HP(n) of an integer n is obtained by concatenating its prime factors (in ascending order) to form a new number and then repeating this process until the result is prime.  Perl is good at switching between the binary and character representation of a number as the context demands, so if we are allowed factor and is_prime from Math::Util this task is fairly trivial:

sub home {

    my @prime_factors = factor(shift);
    my $concat .= $_ for @prime_factors;
    return $concat if is_prime($concat);
    home($concat);
}

This works in microseconds up to the 48th home prime, but the 49th is apparently yet to be solved.