Merge pull request #5589 from Util/branch-for-challenge-149

Add Raku, Perl, and C solutions, blog URL, for TWC 149 by Bruce Gray

8 files changed, 303 insertions, 76 deletions

diff --git a/challenge-149/bruce-gray/README b/challenge-149/bruce-gray/README
index 19a8ae9e9b..cec9d84ed8 100644
--- a/challenge-149/bruce-gray/README
+++ b/challenge-149/bruce-gray/README
@@ -1,79 +1,72 @@
-Solutions by Bruce Gray for https://theweeklychallenge.org/blog/perl-weekly-challenge-148/
+Solutions by Bruce Gray for https://theweeklychallenge.org/blog/perl-weekly-challenge-149/
+Languages: Raku, Perl, C
 
-The Raku solution to Task#2 shows four different results for "first 5",
-to show the different orderings produced by different algorithms.
+Task 1: Generate https://oeis.org/A028840
+Task 2: Generate https://oeis.org/A287298
 
-Output:
-$ perl perl/ch-1.pl
-    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 4 6 30 32 34 36 40 42 44 46 50 52 54 56 60 62 64 66
-
-$ perl perl/ch-2.pl 
-     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))
+Sample runs:
+$ perl perl/ch-1.pl 30
+$ raku raku/ch-1.raku 30
+$ gcc -Wall c/ch-1.c && ./a.out 30
+    All have the output:
+    0, 1, 2, 3, 5, 8, 10, 11, 12, 14, 17, 20, 21, 23, 26, 30, 32, 35, 41, 44, 49, 50, 53, 58, 62, 67, 71, 76, 80, 85
     
-
-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².
-
-
-Jean Marie goes further, in https://math.stackexchange.com/questions/1885095/parametrization-of-cardano-triplet ,
-showing (halfway through) that 𝑎 ≡ 2 𝑚𝑜𝑑 3. 
-
-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.
+$ perl perl/ch-2.pl 2-12 14-16 18-19
+     2           1                     1         1                 1
+     3           1                     1         1                 1
+     4          15                   225        33              3201
+     5          24                   576        44              4301
+     6         195                 38025       523            452013
+     7         867                751689      2346           6250341
+     8        3213              10323369      6215          47302651
+     9       18858             355624164     27773         823146570
+    10       99066            9814072356     99066        9814072356
+    11      528905          279740499025    331413       A8701245369
+    12     2950717         8706730814089    BA3711      B8750A649321
+    14   105011842     11027486960232964   DD3789C    DC71B30685A924
+    15   659854601    435408094460869201  3CDE271B   EDAC93B24658701
+    16  4285181505  18362780530794065025  FF6AAE41  FED5B39A42706C81
+    18  198009443151  3.92077395769691e+22  HH7CF68B9  HGF80ADC537126GBH2
+    19  1404390324525  1.97231218361943e+24  46D29B1F53  IHGFD3408C68ID1IBG7
+(21m38s runtime)
+$ raku raku/ch-2.raku 2-12 14-16 18-19
+     2           1                     1         1                 1
+     3           1                     1         1                 1
+     4          15                   225        33              3201
+     5          24                   576        44              4301
+     6         195                 38025       523            452013
+     7         867                751689      2346           6250341
+     8        3213              10323369      6215          47302651
+     9       18858             355624164     27773         823146570
+    10       99066            9814072356     99066        9814072356
+    11      528905          279740499025    331413       A8701245369
+    12     2950717         8706730814089    BA3711      B8750A649321
+    14   105011842     11027486960232964   DD3789C    DC71B30685A924
+    15   659854601    435408094460869201  3CDE271B   EDAC93B24658701
+    16  4285181505  18362780530794065025  FF6AAE41  FED5B39A42706C81
+    18  198009443151  39207739576969100808801  HH7CF68B9  HGF80ADC53712EB649
+    19  1404390324525  1972312183619434816475625  46D29B1F53  IHGFD3408C6E715A2B9
+(4m51s runtime)
+$ gcc -Wall -lgmp c/ch-2.c && ./a.out 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16    18 19 20   22 23 24 
+f(2)='1'
+f(3)='1'
+f(4)='3201'
+f(5)='4301'
+f(6)='452013'
+f(7)='6250341'
+f(8)='47302651'
+f(9)='823146570'
+f(10)='9814072356'
+f(11)='a8701245369'
+f(12)='b8750a649321'
+f(13)='cba504216873'
+f(14)='dc71b30685a924'
+f(15)='edac93b24658701'
+f(16)='fed5b39a42706c81'
+f(18)='hgf80adc53712eb649'
+f(19)='ihgfd3408c6e715a2b9'
+f(20)='jihg03dac457bfe96281'
+f(22)='lkjig5d14b9032fhac867e'
+f(23)='mlkjefg5ic1d9h8042ab376'
+f(24)='nmlkjbgc6a0d579482i3efh1'
+(13m20s runtime)
diff --git a/challenge-149/bruce-gray/blog.txt b/challenge-149/bruce-gray/blog.txt
new file mode 100644
index 0000000000..9d904325e7
--- /dev/null
+++ b/challenge-149/bruce-gray/blog.txt
@@ -0,0 +1 @@
+http://blogs.perl.org/users/bruce_gray/2022/01/twc-149-limited-fibs-and-bases-of-unusual-size.html
\ No newline at end of file
diff --git a/challenge-149/bruce-gray/c/ch-1.c b/challenge-149/bruce-gray/c/ch-1.c
new file mode 100644
index 0000000000..ca9a492ea5
--- /dev/null
+++ b/challenge-149/bruce-gray/c/ch-1.c
@@ -0,0 +1,61 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+
+int
+is_fibonacci(const int a) {
+    // The maximum digit sum of a 64-bit int is 171,
+    // so Fibonacci 233 and higher are not needed.
+    return (a ==  0 || a ==  1 || a ==  2 || a ==  3 || a ==  5 || a ==  8
+         || a == 13 || a == 21 || a == 34 || a == 55 || a == 89 || a ==144)
+        ? 1 : 0;
+}
+
+int
+sum_of_digits(int a) {
+    int sum = 0;
+    int digit;
+
+    while (a) {
+        digit = a % 10;
+        sum  += digit;
+        a    -= digit;
+        a    /= 10;
+    }
+
+    return sum;
+}
+
+int
+main(int argc, char **argv) {
+    assert(argc == 2);
+    int count;
+    int i;
+
+    if (argv[1][0] == 'b') {
+        // Exactly reproduce https://oeis.org/A028840/b028840.txt
+        for (i = 0, count = 1 ; count <= 10000 ; i++) {
+            if ( is_fibonacci(sum_of_digits(i)) ) {
+                printf("%d %d\n", count, i);
+                count++;
+            }
+        }
+        for (count = 1 ; count < 64 ; count++) {
+            printf("\n");
+        }
+    }
+    else {
+        // Produce first N values, formated as per Task example.
+        for (i = 0, count = atoi(argv[1]) ; count ; i++) {
+            if ( is_fibonacci(sum_of_digits(i)) ) {
+                count--;
+                printf("%d", i);
+                if (count)
+                    printf(", ");
+            }
+        }
+        printf("\n");
+    }
+    return 0;
+}
diff --git a/challenge-149/bruce-gray/c/ch-2.c b/challenge-149/bruce-gray/c/ch-2.c
new file mode 100644
index 0000000000..59ff9f7d78
--- /dev/null
+++ b/challenge-149/bruce-gray/c/ch-2.c
@@ -0,0 +1,78 @@
+// Adapted from Python+gmpy2 code at https://oeis.org/A287298
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include "gmp.h"
+
+#define MAX_BASE 62
+
+char digits[MAX_BASE+2];        // Never using more
+char repetitions_array[128];    // max ord of '9','Z','z' is 122
+mpz_t m1;
+mpz_t mpz_square;
+mpz_t initial_all_digits;
+
+int
+has_no_repeating_digits(char *s) {
+    memset(repetitions_array, 0, sizeof(repetitions_array));
+    unsigned char *c;
+    for (c = (unsigned char *) s ; c[0] != '\0' ; c++) {
+        if (repetitions_array[ c[0] ])
+            return 0;
+        repetitions_array[ c[0] ] = 1;
+    }
+    return 1;
+}
+
+int
+square_has_no_repeating_digits(const mpz_t i, int base) {
+    mpz_mul(mpz_square, i, i);
+    mpz_get_str(digits, base, mpz_square);
+
+    return has_no_repeating_digits(digits);
+}
+
+void
+set_all_digits_descending(mpz_t m, int base) {
+    int digit;
+    // e.g. For base=5, initial_all_digits will be 43210
+    mpz_set_ui(m, 0);
+    for (digit = base - 1 ; digit >= 0 ; digit--) {
+        mpz_add_ui(m, m, digit);
+        mpz_mul_ui(m, m, base);
+    }
+    mpz_divexact_ui(m, m, base);
+}
+
+mpz_t*
+A287298(int base) {
+    assert(base <= 62);
+    assert(base >=  2);
+
+    set_all_digits_descending(m1, base);
+    for ( mpz_sqrt(m1, m1) ; mpz_cmp_si(m1, 0) > 0 ; mpz_sub_ui(m1, m1, 1) ) {
+        if (square_has_no_repeating_digits(m1, base))
+            return &m1;
+    }
+    assert(1 == 0);  // Should not ever hit this point!
+}
+
+int
+main(int argc, char **argv) {
+    int i;
+    assert(argc > 1);
+
+    mpz_init(m1);
+    mpz_init(mpz_square);
+    mpz_init(initial_all_digits);
+
+    for (i = 1 ; i < argc ; i++) {
+        int base = atoi(argv[i]);
+        mpz_t *r = A287298(base);
+        mpz_mul(mpz_square, *r, *r);
+        printf("f(%d)='%s'\n", base, mpz_get_str(NULL, base, mpz_square));
+    }
+
+    return 0;
+}
diff --git a/challenge-149/bruce-gray/perl/ch-1.pl b/challenge-149/bruce-gray/perl/ch-1.pl
new file mode 100644
index 0000000000..f4923e0c76
--- /dev/null
+++ b/challenge-149/bruce-gray/perl/ch-1.pl
@@ -0,0 +1,14 @@
+use Modern::Perl;
+use List::Util      qw<sum>;
+use Math::Fibonacci qw<isfibonacci>;
+
+my $count = abs( shift // 20 );
+
+my @out;
+for ( my $i = 0 ; $count ; ++$i ) {
+    next unless isfibonacci( sum split '', $i );
+    $count--;
+    push @out, $i;
+}
+
+say join ', ', @out;
diff --git a/challenge-149/bruce-gray/perl/ch-2.pl b/challenge-149/bruce-gray/perl/ch-2.pl
new file mode 100755
index 0000000000..3ad4fedb7b
--- /dev/null
+++ b/challenge-149/bruce-gray/perl/ch-2.pl
@@ -0,0 +1,37 @@
+#!/usr/bin/env perl
+use Modern::Perl;
+use experimental qw<signatures>;
+use List::Util   qw<uniq>;
+use Math::BigInt;
+# use Math::BigInt only => 'GMP,Pari';
+
+sub largest_rep_free_square_in_base ($base) {
+    # 'BA9876543210' for $base==12
+    my $highest = join '', map { Math::BigInt->new($_)->to_base($base) }
+                  reverse 0 .. $base-1;
+
+    # $square will be a _single_ Math::BigInt. 
+    # $root will be plain perl IV for speed. (Will overflow past base==27)
+    my $root   = Math::BigInt->from_base($highest, $base)->bsqrt->numify();
+    my $square = Math::BigInt->new($root) ** 2;
+
+    while ($root) {
+        my @digits = split '', $square->to_base($base);
+        return $root if @digits == uniq(@digits);
+
+        $square->bsub( $root + $root - 1 ); # (a - 1)² == a² - 2a + 1
+        $root--;
+    }
+    return;
+}
+
+my @bases = map { /(\d+)-(\d+)/ ? ($1..$2) : $_ } @ARGV;
+@bases = (2 .. 12, 14 .. 16, 18) if not @bases;
+
+for my $base (@bases) {
+    my $root = largest_rep_free_square_in_base($base);
+    my @rs = ( $root, $root**2 );
+
+    printf "%2s  %10s  %20s  %8s  %16s\n",
+        $base, @rs, map { Math::BigInt->new($_)->to_base($base) } @rs;
+}
diff --git a/challenge-149/bruce-gray/raku/ch-1.raku b/challenge-149/bruce-gray/raku/ch-1.raku
new file mode 100644
index 0000000000..63c2a373e4
--- /dev/null
+++ b/challenge-149/bruce-gray/raku/ch-1.raku
@@ -0,0 +1,20 @@
+multi sub is-fibonacci (       0 --> True ) { }
+multi sub is-fibonacci ( UInt $n --> Bool ) {
+    constant φ = ( 1 + 5.sqrt ) / 2;
+    constant @fib = 0, 1, { $^a + $^b } … Inf;
+
+    # my $k = ( ($n * 5.sqrt).log / φ.log ).round;    #  20% Faster
+    my $k = @fib.first: :k, * >= $n;              # 100% Clearer
+
+    return @fib[$k] == $n;
+}
+
+constant @s = grep *.comb.sum.&is-fibonacci, ^Inf;
+
+sub MAIN ( UInt $count = 20 ) {
+    say @s.head($count).join(", ");
+}
+
+# Exactly reproduce https://oeis.org/A028840/b028840.txt
+# say "{.key + 1} {.value}" for @s.head(10_000).pairs;
+# say '' for ^63;
diff --git a/challenge-149/bruce-gray/raku/ch-2.raku b/challenge-149/bruce-gray/raku/ch-2.raku
new file mode 100644
index 0000000000..6cc25f4897
--- /dev/null
+++ b/challenge-149/bruce-gray/raku/ch-2.raku
@@ -0,0 +1,23 @@
+subset Base of Int where 2..36;
+
+sub no-reps ( Str $s --> Bool ) { $s.comb.repeated.not }
+
+sub largest_rep-free_square_in_base ( Base $base ) {
+    my Str $highest = (^$base)».base($base).join.flip; # 'BA9876543210' for $base==12
+
+    my UInt $root = $highest.parse-base($base).sqrt.round;
+
+    $root-- until ($root²).base($base).&no-reps;
+
+    return $root;
+}
+
+my @bases = flat map { /(\d+)\-(\d+)/ ?? (+$0 .. +$1) !! $_ },
+            (@*ARGS ?? @*ARGS !! <2-12 14-16 18>);
+
+for @bases -> Base $base {
+    my $root = largest_rep-free_square_in_base($base);
+
+    printf "%2s  %10s  %20s  %8s  %16s\n",
+        $base, $root, $root², $root.base($base), ($root²).base($base);
+}