Perl Weekly Challenge 148: Cardano Triplets

(2,1,5) is the first Cardano Triplets.

First point of business is, what defines the first 5 triplets? Do we order them on smallest a, then smallest b? Do we order them on the minimum of a, b, and c (breaking ties with the second smallest)? Do we order them on the sum of a, b and c? Something else?

We decide to order them on the sum of a, b, and c.

Cordano triplets are the subject of task 251 of Project Euler.

We start off by manipulating the formula above to eliminate the cube and square roots ([U3]):

\[ \begin{array}{lcl} \sqrt[3]{a + b \sqrt{c}} + \sqrt[3]{a - b \sqrt{c}} = 1 & \implies & \\ \sqrt[3]{a + b \sqrt{c}} = 1 - \sqrt[3]{a - b \sqrt{c}} & \implies & \text{(cubing both sides)} \\ a + b \sqrt{c} = 1 - 3 \left(\sqrt[3]{a - b \sqrt{c}}\right) \left(1 - \sqrt[3]{a - b \sqrt{c}}\right) - a + b \sqrt{c} & \implies & \\ 2 a - 1 = -3 \left(\sqrt[3]{a - b \sqrt{c}}\right) \left(1 - \sqrt[3]{a - b \sqrt{c}}\right) & \implies & \text{(identity on second line)} \\ 2 a - 1 = -3 \left(\sqrt[3]{a - b \sqrt{c}}\right) \left(\sqrt[3]{a + b \sqrt{c}}\right) & \implies & \text{(cubing both sides)} \\ 8 a^3 - 12 a^2 + 6 a - 1 = -27 \left(a^2 - b^2 c\right) & \implies & \\ 8 a^3 + 15 a^2 + 6 a - 27 b^2 c = 1 \end{array} \tag{1} \]

Given the formula \(8 a^3 + 15 a^2 + 6 a - 27 b^2 c = 1\), we can think of a parametrization.

Let's take this modulo \(3\) [JM]:

\[ \begin{array}{lcl} 8 a^3 + 15 a^2 + 6 a - 27 b^2 c \equiv 1 \text{ mod } 3 & \implies & \text{(}15 a^2 \text{, } 6 a \text{, and } 27 b^2 c \text{ are all divisible by } 3 \text{)} \\ 2 a^3 \equiv 1 \text{ mod } 3 & \implies & \text{(multiply both sides by 2)} \\ 4 a^3 \equiv 2 \text{ mod } 3 & \implies & \\ a^3 \equiv 2 \text{ mod } 3 & \implies & \\ a \equiv 2 \text{ mod } 3 & \implies & \\ a = 3 k + 2 \text{ for some integer } k \geq 0 & & \end{array} \tag{2} \]

Now, if we plug in the result of (2) into the result of (1), we get:

\[ \begin{array}{lcl} 8 a^3 + 15 a^2 + 6 a - 27 b^2 c = 1 & \implies & \\ 27 b^2 c = 8 a^3 + 15 a^2 + 6 a - 1 & \implies & \text{(result of (2))} \\ 27 b^2 c = 8 (3 k + 2)^3 + 15 (3k + 2)^2 + 6 (3k + 2) - 1 & \implies & \text{(Mathematica)} \\ 27 b^2 c = 27 (k + 1)^2 (8k + 5) & \implies & \\ b^2 c = (k + 1)^2 (8k + 5) & & \end{array} \]

So, for each \(k \geq 0\) we find an \(a\), and a \(b^2 c\). But there can be multiple solutions for \(b\) and \(c\):

1. Each divisor \(d_1\) of \(k + 1\) will be a solution for \(b\).
2. For each divisor \(d_2\) of \(8k + 5\), if \(d_2\) is a perfect square, \(\sqrt{d_2}\) will be a solution for \(b\).
3. For each divisor \(d_1\) found in step 1, and each divisor \(\sqrt{d_2}\) found in step 2, \(d_1 \sqrt{d_2}\) will be a solution for \(b\).

Once we have \(k\) and \(b\), we can calculate \[ c = \frac{(k + 1)^2 (8 k + 5)}{b^2} \]

How do we know we can stop generating triangles? We will keep an array of the best 5 triangles found so far. If \(3 k + 2 = a\) is larger than the largest of the sums of the best 5 triangles found so far, we know we cannot do better with larger \(k\), so we can stop.

First, we start off with some initialization. We keep the best triples found so far in an array @out — for each triple, we also store its sum. We initialize @out with high numbers, much larger than any triple we will report. We also keep track of the index of the largest triple in @out:

use Math::Prime::Util qw [divisors];
use List::Util        qw [sum];    

my $COUNT = 5;   
my $A     = 0;
my $B     = 1;   
my $C     = 2;
my $SUM   = 3;
my @out;
foreach my $i (0 .. $COUNT - 1) {
    $out [$i]        = [(999999) x 3];  
    $out [$i] [$SUM] = sum @{$out [$i]};
}

my $max_index = 0;

We then start a loop, incrementing k. Given a k, we find a, and store the value k + 1 and 8 * k + 5.

for (my $k = 0; 3 * $k + 2 <= $out [$max_index] [$SUM]; $k ++) {
    my $a  = 3 * $k + 2;
    my $f1 =     $k + 1;
    my $f2 = 8 * $k + 5;

We then find the divisors of f1, and the square roots of the divisors of f2, discarding the ones which are integers:

my @d1 = divisors ($f1);
my @d2 = grep {$_ == int ($_)} map {sqrt $_} divisors ($f2);

We know multiply each element of @d1 and @d2, assign this to b, and calculate c. If the sum of a, b and c is better than the fifth best triple found so far, we discard the fifth best triple, replacing it with a, b and c. But we only do this if a, b and c are not yet in @out: we don't need duplicates:

foreach my $d1 (@d1) {
  D2:
    foreach my $d2 (@d2) {
        my $b = $d1 * $d2;
        my $c = $f1 * $f1 * $f2 / ($b * $b);
        if ($a + $b + $c < $out [$max_index] [$SUM]) {
            #
            # Avoid duplicate entries
            #
            foreach my $info (@out) {
                next D2 if $$info [$A] == $a && $$info [$B] == $b;
            }

            #
            # Put triple in the output structure
            #
            $out [$max_index] = [$a, $b, $c, $a + $b + $c];

We now have to find the index of the newest fifth best triple:

#
# Find the index of the new highest value
#
$max_index = 0;
my $max_sum = $out [$max_index] [$SUM];
for (my $i = 1; $i < $COUNT; $i ++) {
    if ($out [$i] [$SUM] > $max_sum) {
        $max_index = $i;
        $max_sum   = $out [$i] [$SUM];
    }
}

Once we have found all triples, we can print them:

say "@$_[$A, $B, $C]" for @out;

Find the full program on GitHub.

Most languages don't have a divisor method (or they may have in some module, but I generally don't go searching for them), so we have to find the divisors ourselves. Python for instance.

We'll skip the initialization, and dive directly into the body of the main loop, where we have to find the divisors of f1 and f2:

a  = 3 * k + 2
f1 =     k + 1
f2 = 8 * k + 5

d1 = []
for i in range (1, f1 + 1):
    if i * i > f1:
        break
    if f1 % i == 0:
        d1 . append (i)
        if i != f1 // i: 
            d1 . append (f1 // i)


d2 = []
for i in range (1, f2 + 1):
    if i * i > f2:
        break
    if f2 % i == 0:
        s1 = math . isqrt (i)
        s2 = math . isqrt (f2 // i)
        if s1 * s1 == i:
            d2 . append (s1)
        if s2 * s2 == f2 // i and s1 != s2:
            d2 . append (s2)

Finding the possibilities for b and c isn't very different from the Perl solution:

for d1v in d1:
    for d2v in d2:
        b = d1v * d2v
        c = f1 * f1 * f2 // (b * b)

        if a + b + c < out [max_index] [SUM]:
            seen = False
            for i in range (COUNT):
                if out [i] [A] == a and out [i] [B] == b:
                    seen = true
                    break
            if seen:
                break

            out [max_index] = [a, b, c, a + b + c]

            max_index = 0
            max_sum   = out [max_index] [SUM]
            for i in range (1, COUNT):
                if max_sum < out [i] [SUM]:
                    max_sum   = out [i] [SUM]
                    max_index = i

Find the full program on GitHub.

We also have implementation in some other languages, all using a very similar algorithm:

AWK, bc, C, Go, Java, Lua, Node.js, Pascal, R, Ruby, and Tcl.