Merge pull request #5557 from Util/branch-for-challenge-148

Add Raku and Perl solutions for TWC 148 by Bruce Gray.

5 files changed, 142 insertions, 17 deletions

diff --git a/challenge-148/bruce-gray/README b/challenge-148/bruce-gray/README
index 7acad0424d..19a8ae9e9b 100644
--- a/challenge-148/bruce-gray/README
+++ b/challenge-148/bruce-gray/README
@@ -1,25 +1,79 @@
-Solutions by Bruce Gray for https://theweeklychallenge.org/blog/perl-weekly-challenge-147/
+Solutions by Bruce Gray for https://theweeklychallenge.org/blog/perl-weekly-challenge-148/
 
-NOTE: Both Task#2 solutions deliberately avoid the need for is_pentagon_number(),
-[as per quadratic transformation of `n(3n-1)/2 = P`
-to `sqrt(24P + 1) must be 5(mod 6)` and `24P + 1 must be a perfect square`]
-, just for fun.
+The Raku solution to Task#2 shows four different results for "first 5",
+to show the different orderings produced by different algorithms.
 
-Sample runs:
+Output:
 $ perl perl/ch-1.pl
-    2 3 5 7 13 17 23 37 43 47 53 67 73 83 97 103 107 113 137 167
-    2 3 5 7 13 17 23 37 43 47 53 67 73 83 97 113 137 167 173 197
-
+    2 4 6 30 32 34 36 40 42 44 46 50 52 54 56 60 62 64 66
 $ raku raku/ch-1.raku
-    2 3 5 7 13 17 23 37 43 47 53 67 73 83 97 103 107 113 137 167
-    2 3 5 7 13 17 23 37 43 47 53 67 73 83 97 113 137 167 173 197
+    2 4 6 30 32 34 36 40 42 44 46 50 52 54 56 60 62 64 66
 
 $ perl perl/ch-2.pl 
-    7042750 5482660
+     2   1   5
+     5   1  52
+     8   1 189
+    11   1 464
+    14   1 925
+$ raku raku/ch-2.raku
+    ((  2   1   5) (  5   2  13) ( 17  18   5) ( 17   9  20) (  8   3  21))
+    ((  2   1   5) (  5   2  13) (  8   3  21) ( 17   9  20) ( 17  18   5))
+    ((  2   1   5) (  5   1  52) (  5   2  13) (  8   1 189) (  8   3  21))
+    ((  2   1   5) (  5   1  52) (  8   1 189) ( 11   1 464) ( 14   1 925))
+    
+
+Analysis of Task#2:
+
+If I get a blogpost written, I plan to delve into how `=~=` is insufficient for this task, as the default $*TOLERANCE misses some cases.
+
+Original equation: (a + b√c)^⅓ + (a - b√c)^⅓ = 1
+When solved for `c` via https://www.wolframalpha.com/ :
+  Solve[ Cbrt[a + bSqrt[c]] + Cbrt[a - bSqrt[c]] = 1, c ]
+      c = (a + 1)² * (8a - 1) / 27b²
+  Useful!
+  Also means that:
+      (a + 1)² * (8a - 1) / 27b²c = 1
+
+# Full derivation:
+# https://math.stackexchange.com/questions/2160805/cardano-triplet-transformation
+Original equation:
+    (a + b√c)^⅓ + (a - b√c)^⅓ = 1
+Move 2nd term across:
+    (a + b√c)^⅓  = 1 - (a - b√c)^⅓
+Cubing just removes the cube-root on the left, and expands on the right to:
+    Expand[ (1 - Cbrt[a - bSqrt[c]])³ ]
+                    3 ((a - b√c)^⅓)² - 3 (a - b√c)^⅓ - a + b√c + 1
+     a + b√c == 1 + 3 ((a - b√c)^⅓)² - 3 (a - b√c)^⅓ - a + b√c
+     a       == 1 + 3 ((a - b√c)^⅓)² - 3 (a - b√c)^⅓ - a      
+    2a       == 1 + 3 ((a - b√c)^⅓)² - 3 (a - b√c)^⅓          
+    2a - 1   ==     3 ((a - b√c)^⅓)² - 3 (a - b√c)^⅓
+    2a - 1   ==    -3 ((a - b√c)^⅓) (1 - (a - b√c)^⅓)
+Use original equality to substitute the last part from (a-...) to (a+...) :
+    2a - 1   ==    -3 ((a - b√c)^⅓) ((a + b√c)^⅓)
+Cube both sides again:
+    Expand[ (2a - 1)³ ]
+        8a³ - 12a² + 6a - 1
+    Expand[ (-3 ((a - b√c)^⅓) ((a + b√c)^⅓))³ ]
+        27b²c - 27a²
+    8a³ - 12a² + 6a - 1 == 27b²c - 27a²
+    8a³ + 15a² + 6a - 1 == 27b²c
+Factor[8a³ + 15a² + 6a - 1]
+    (a + 1)² (8a - 1) == 27b²c
+Now solving for `c` is easy:
+    (a + 1)² (8a - 1) / 27b² == c
+`c` can only be a whole number if (a + 1)² (8a - 1) can be evenly divided by 27b².
+
 
-$ time raku raku/ch-2.raku
-    7042750 5482660
+Jean Marie goes further, in https://math.stackexchange.com/questions/1885095/parametrization-of-cardano-triplet ,
+showing (halfway through) that 𝑎 ≡ 2 𝑚𝑜𝑑 3. 
 
-    real	0m2.432s
-    user	0m2.534s
-    sys 	0m0.067s
+Humor:
+Easy to prove that, if we lock `b` to always be 1, and `𝑎 ≡ 2 𝑚𝑜𝑑 3` , then `c` will be integer for all `a` generated from 3k+2.
+    (a + 1)² (8a - 1) / 27 == c
+    Expand[((3k + 2) + 1)² * (8(3k + 2) - 1)]
+        216k² + 567k² + 486k + 135
+    Factor[Expand[((3k + 2) + 1)² * (8(3k + 2) - 1)]]
+        27 (k + 1)² (8k + 5)
+    Aha! Always divisible by 27!
+So, a cheap way to generate Cardano triplets is (3k+2, 1, (k + 1)² (8k + 5)) for k=0..Inf.
+Since 1 in the lowest possible `b`, using k=0..4 would give a reasonable (although unexpected) answer for the task.
diff --git a/challenge-148/bruce-gray/perl/ch-1.pl b/challenge-148/bruce-gray/perl/ch-1.pl
new file mode 100644
index 0000000000..200522d8f4
--- /dev/null
+++ b/challenge-148/bruce-gray/perl/ch-1.pl
@@ -0,0 +1,5 @@
+use Modern::Perl;
+use Lingua::EN::Numbers qw<num2en>;
+
+say join ' ', grep { !(num2en($_) =~ /e/) } 0..100;
+
diff --git a/challenge-148/bruce-gray/perl/ch-2.pl b/challenge-148/bruce-gray/perl/ch-2.pl
new file mode 100644
index 0000000000..725cafb3bc
--- /dev/null
+++ b/challenge-148/bruce-gray/perl/ch-2.pl
@@ -0,0 +1,17 @@
+use Modern::Perl;
+use experimental qw<signatures>;
+
+# See README for analysis common to Perl and Raku.
+
+# Uses WolframAlpha solution for `c`.
+# Empty return if c would be non-integer, otherwise aref of all three.
+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 ];
+}
+
+# Locking b=1, to use a humorous definition of "first 5".
+say sprintf('%3d %3d %3d', @{$_}) for map { find_Cardano_triplet( (3 * $_ + 2, 1), ) } 0..4;
diff --git a/challenge-148/bruce-gray/raku/ch-1.raku b/challenge-148/bruce-gray/raku/ch-1.raku
new file mode 100644
index 0000000000..973d83788c
--- /dev/null
+++ b/challenge-148/bruce-gray/raku/ch-1.raku
@@ -0,0 +1,3 @@
+use Lingua::EN::Numbers;
+
+put grep *.&cardinal.contains('e').not, 0..100;
diff --git a/challenge-148/bruce-gray/raku/ch-2.raku b/challenge-148/bruce-gray/raku/ch-2.raku
new file mode 100644
index 0000000000..b8cb232279
--- /dev/null
+++ b/challenge-148/bruce-gray/raku/ch-2.raku
@@ -0,0 +1,46 @@
+# See README for analysis common to Perl and Raku.
+
+# Technical definition from the task.
+sub is_Cardano_triplet ( (\a,\b,\c) --> Bool ) {
+    sub cbrt (\n) { n.sign  *  n.abs ** ⅓ }
+
+    my \bsc = b * sqrt(c);
+    return 1e-14 > abs( cbrt(a + bsc) + cbrt(a - bsc) - 1 );
+}
+
+# Given (A,B), solve for C, and return all three if C is integer.
+sub find_Cardano_triplet ( (\a, \b) ) {
+    my Rat \c = (a + 1)² * (8*a - 1) / (27 * b²);
+    return (a, b, c.narrow) if c.denominator == 1;
+}
+
+# My own algorithm for efficient lazy infinite N-tuples.
+# May be a module someday, or at least a blog post.
+sub triplet_generator ( ) {
+    return lazy gather for 1..* -> $limit {
+        my @inner = 1 ..^ $limit;
+
+        # 3 outer planes minus the edges, then the 3 edges, then the corner.
+        .take for @inner X @inner X $limit;
+        .take for @inner X $limit X @inner;
+        .take for $limit X @inner X @inner;
+
+        .take for $limit X $limit X @inner;
+        .take for $limit X @inner X $limit;
+        .take for @inner X $limit X $limit;
+
+        .take for $limit X $limit X $limit;
+    }
+}
+
+constant @fixed_X_triplets = [X] (1..21) xx 3;
+constant @fixed_X_doublets = [X] (1..21) xx 2;
+
+my @answers =
+    triplet_generator().grep( &is_Cardano_triplet), # Small 3-tuples, no bias
+      @fixed_X_triplets.grep( &is_Cardano_triplet), # Small 3-tuples, left-biased, c-range-limit forces higher a&b
+      @fixed_X_doublets.map(&find_Cardano_triplet), # Small 2-tuples, left-biased, unlimited range for `c`
+                 (^Inf).map({find_Cardano_triplet( (3 * $_ + 2,  1), )}), # Force b=1
+;
+
+say .head(5)».fmt('%3d') for @answers;