TWC 148: Numbrs, and Cubic Cards

In which we see that you can spell "numbers" without an e, and realize that first implies an ordering.

TWC Task #1 - Eban Numbers

Eban: A number that has no letter e when spelled in English.

Observations:

Since the Lingua::EN::Numbers exists for both Perl and Raku, the solutions are just expanded one-liners. Both output:

2 4 6 30 32 34 36 40 42 44 46 50 52 54 56 60 62 64 66

Raku

use Lingua::EN::Numbers;
put grep *.&cardinal.contains('e').not, 0..100;

Perl

use Modern::Perl;
use Lingua::EN::Numbers qw<num2en>;

say join ' ', grep { !(num2en($_) =~ /e/) } 0..100;

Note: I missed using Perl's !~ operator, which would have been clearer than negating the whole expression: grep { num2en($_) !~ /e/ }

TWC Task #2 - Cardano Triplets

  ________     ________
ยณโˆš ๐‘Ž + ๐‘โˆš๐‘  + ยณโˆš ๐‘Ž - ๐‘โˆš๐‘  == 1

Observations:

  • Most of my calculations are in the README.

  • Exact wording of the task is "Write a script to generate first 5 Cardano Triplets".

  • The order that triplets occur depends on the method you use to search for them. Therefore, if more than one method exists to search for Cardano triplets, the output of all the participants might not match each other.

  • Several methods exists to search for Cardano triplets.

  • a,b,c are independent and not interchangeable; they each play a different role in the equation. So, we cannot use the common technique of "for a = 1..Inf", then have inner loops that bound b and c to <= a, or we risk missing some triplets.

  • Whenever searching for N-tuples (in this case, triples) that satisfy some condition, a good approach is to see if one piece of the tuple can be calculated from the other pieces. If so, you can remove a loop!

  • c can be calculated from a,b.
    WolframAlpha says that
    Solve[CubeRoot[a + b Sqrt[c]] + CubeRoot[a - b Sqrt[c]] == 1, c] has the solution c = (8a - 1)(a + 1)ยฒ / 27bยฒ.
    So, if (8a - 1)(a + 1)ยฒ is evenly divisible by 27bยฒ, then we have found a triplet and computed c.

  • All triplets must have ๐‘Ž โ‰ก 2 ๐‘š๐‘œ๐‘‘ 3, which is the same as saying a = 3k+2 where k is a positive integer.
    (Stated halfway through Jean Marie's StackExchange answer to Parametrization of Cardano triplet)

  • Combining the last two points, we can replace a with 3k+2 in (8a - 1)(a + 1)ยฒ, and simplify it to 27(8k + 5)(k + 1)ยฒ.
    Since that has a factor of 27, and every triple must be evenly divisible by 27bยฒ, then every valid a (that is, ๐‘Ž โ‰ก 2 ๐‘š๐‘œ๐‘‘ 3) has a least one triple, where b=1. If a has square divisors, then there are multiple triples for that a.

  • A straightforward coding of is_Cardano_triplet() uses floating-point math then checks for an equality, which is the classic FP bug. We must compare to some tolerance, so that we do not exclude a valid triple.

    say "Bad!" if sqrt(3)**2 != 3' # Output: Bad!

Raku

  • Since the math involves more than just division, even Raku's invisible Rat(s) will not save us; ($r - 1).abs < $some_tolerance, or maybe $r =~= 1.

  • Raku's is-almost-equal-to operator, =~= uses the current $*TOLERANCE setting, whose default of 1e-15 is too strict to find some triplets.
    (Aside: if you do something like my $*TOLERANCE = 1e-14; and omit the minus sign, then suddenly every triple is a Cardano triple, because now one is almost equal to a trillion. Blew my mind until I saw my error.)

  • I found it concise to use two uncommon techniques: find_Cardano_triplet() and is_Cardano_triplet() both take a single argument of a triplet (notice the doubled parenthesis in the sub declarations), and return either a complete triplet or Nil. You might expect find_Cardano_triplet() to return c, is_Cardano_triplet() to return a Bool and both to take 3 arguments. This allowed a slick combining of candidate sources and filters/generators via map and grep.

  • While I used $m % $n in the Perl code to check for divisibility, Raku's Rat types always reduce to GCD, so I can just divide and then examine .denominator == 1. The .narrow method is needed so that I don't return a Rat where one would expect an Int.

  • Similar to using Unicode superscript numerals to express exponents ($nยฒ), you can use Unicode "vulgar fractions" in Raku. This was convenient for cube roots, although special handling was needed for negatives.

    say [15 * โ…–, 27 ** โ…“]; # Output: [6 3]

    sub cbrt (\n) { n.sign * n.abs ** โ…“ } # Handles positives and negatives uniformly.

  • While Raku provides the Cross operator (and meta-operator), it does not work usefully with more than one infinite list, so we can only use it (or the equivalent three nested for loops) if we already know the answer to the problem. I have always disliked this limitation (across many problems), but I did use the pre-knowledge this time for comparison, as constants @fixed_X_triplets and @fixed_X_doublets with 1..21 ranges.

  • I invented (probably rediscovered) an algorithm to iterate over N-tuple cross-products of N infinite lazy lists, never generating (or even having to check for) duplicates. It works well in any language with generators, and is a good fit here, but the ordering that it produces the tuples looks irregular to humans.

I expanded the output count from 5 triplets to 6, to better show the differences in the orderings. Here are the 4 lines that my program produces (plus one from elsewhere).

(2 1 5) (5 2 13) (17 18   5) (17 9  20) ( 8  3  21) (11 4   29) # 1
(2 1 5) (5 2 13) ( 8  3  21) (17 9  20) (17 18   5)             # 2
(2 1 5) (5 1 52) ( 5  2  13) ( 8 1 189) ( 8  3  21) (11 1  464) # 3
(2 1 5) (5 1 52) ( 8  1 189) (11 1 464) (14  1 925) (17 1 1620) # 4
(2 1 5) (5 2 13) ( 8  3  21) (11 4  29) (14  5  37) (17 6   45) # 5
  1. triplet_generator().grep( &is_Cardano_triplet)
    "smallest" for minimizing max(a,b,c)

  2. @fixed_X_triplets.grep( &is_Cardano_triplet)
    Notice that it omits the 6th triplet, because the (1..21) that was just enough for 5; it "ran out" of triplets before finding the 6th!

  3. @fixed_X_doublets.map(&find_Cardano_triplet)
    a,b can be smaller because calculated c is not restricted by a range. This one will "run out" too.

  4. (^Inf).map({find_Cardano_triplet( (3 * $_ + 2, 1), )})
    Locks b=1, and walks through all the a=3k+2 since we know they all work with b=1, finding the correct c. This is reasonable; you would get the same answer from a big enough triple-loop if b was the outer variable, find_Cardano_triplet would return Nil for all the a that were not ๐‘Ž โ‰ก 2 ๐‘š๐‘œ๐‘‘ 3 and therefore be suppressed.

  5. This line does not come from my challenge solution; Flavio Poletti adapted Jean Marie's complete answer (that I only managed half of, above) to derive a formula that, given b, produces the smallest-working a and its corresponding c. Splendid! And also, yet another ordering, with a reasonable take on "first 5", that gives different output from any of my four. The cost in finding (5 2 13) so easily, is never seeing other b=2 like (11 2 116).

    say (1..Inf).map({ 3 * $_ - 1, $_, 8 * $_ - 3 }).head(6)
    
  6. Inspired by Flavio, I wrote this just before posting. It was experimentally arrived-at (doublet_generator+find_Cardano_triplet|sort, OEIS), not algebraically derived, but it produces correct output (with tolerance=1e-12) for the first 1000 triplets, locking c=5.

    say (1..Inf).map(-> \n { (15 * (2 * n - 1)ยฒ + 1)/8, (2 * n - 1)*(5 * nยฒ - 5 * n + 2)/2, 5 }).head(6)

    Output: ((2 1 5) (17 18 5) (47 80 5) (92 217 5) (152 459 5) (227 836 5))

Six first 5's, all for valid interpretations of "first". Surprising, and Fun!

Perl

I simply translated find_Cardano_triplet from Raku, and called it with a=0..4, b=1 for the most humorous definition of "first 5".

sub find_Cardano_triplet ( $x, $y ) {
    my $m = ($x + 1)**2 * (8*$x - 1);
    my $n = 27 * $y * $y;

    return if $m % $n;
    return [ $x, $y, $m / $n ];
}
say sprintf('%3d %3d %3d', @{$_}) for map { find_Cardano_triplet( (3 * $_ + 2, 1), ) } 0..4;

Output is the same triples as #4 in the Raku section.

Leave a comment

About Bruce Gray

user-pic "Util" on IRC and PerlMonks. Frequent speaker on Perl and Raku, but infrequent blogger.