Merge pull request #7104 from jo-37/contrib

Solutions to challenge 191

3 files changed, 247 insertions, 0 deletions

diff --git a/challenge-191/jo-37/maxima/ch-2.wxm b/challenge-191/jo-37/maxima/ch-2.wxm
new file mode 100644
index 0000000000..d7e0cbef2d
--- /dev/null
+++ b/challenge-191/jo-37/maxima/ch-2.wxm
@@ -0,0 +1,17 @@
+/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
+/* [ Created with wxMaxima version 13.04.2 ] */
+
+/* [wxMaxima: input   start ] */
+n: 15;
+/* [wxMaxima: input   end   ] */
+
+/* [wxMaxima: input   start ] */
+a: genmatrix(lambda([i, j], if mod(i, j) = 0 or mod(j, i) = 0 then 1 else 0), n, n);
+/* [wxMaxima: input   end   ] */
+
+/* [wxMaxima: input   start ] */
+permanent(a);
+/* [wxMaxima: input   end   ] */
+
+/* Maxima can't load/batch files which end with a comment! */
+"Created with wxMaxima"$
diff --git a/challenge-191/jo-37/perl/ch-1.pl b/challenge-191/jo-37/perl/ch-1.pl
new file mode 100755
index 0000000000..5b446c688e
--- /dev/null
+++ b/challenge-191/jo-37/perl/ch-1.pl
@@ -0,0 +1,75 @@
+#!/usr/bin/perl -s
+
+use v5.16;
+use Test2::V0 '!float';
+use PDL;
+
+our ($tests, $examples);
+
+run_tests() if $tests || $examples;	# does not return
+
+die <<EOS unless @ARGV;
+usage: $0 [-examples] [-tests] [--] [N]
+
+-examples
+    run the examples from the challenge
+ 
+-tests
+    run some tests
+
+N...
+    List of integers
+
+EOS
+
+
+### Input and Output
+
+say twice_largest(@ARGV);
+
+
+### Implementation
+#
+# The task requires careful reading: Operating on *items*, not on
+# *values*.  This differentiation becomes crucial if the maximum is not
+# taken on a unique item.  Interpreting this task as: Take the largest
+# two items from the list and check if the first is at least twice the
+# second.
+
+sub twice_largest {
+    state $w = long(1, -2)->dummy(0);
+
+    my $n = long(@_);
+    # Find the indices of the two largest items.
+    maximum_n_ind($n, my $max2 = zeroes(indx, 2));
+
+    # Check if the maximum item is at least twice the next.  As items
+    # are integers, we may compare the scalar product of the maximum two
+    # items with the weight vector against -0.5 to gain the requested
+    # -1/1 result.
+    ($n->index($max2) x $w)->sclr <=> -0.5;
+}
+
+
+### Examples and tests
+
+sub run_tests {
+    SKIP: {
+        skip "examples" unless $examples;
+
+        is twice_largest(1, 2, 3, 4), -1, 'Example 1';
+        is twice_largest(1, 2, 0, 5),  1, 'Example 2';
+        is twice_largest(2, 6, 3, 1),  1, 'Example 3';
+        is twice_largest(4, 5, 2, 3), -1, 'Example 4';
+    }
+
+    SKIP: {
+        skip "tests" unless $tests;
+
+        is twice_largest(1, 3, 3), -1, 'non-unique max';
+        is twice_largest(-10, 1, -20), 1, 'negative items';
+	}
+
+    done_testing;
+    exit;
+}
diff --git a/challenge-191/jo-37/perl/ch-2.pl b/challenge-191/jo-37/perl/ch-2.pl
new file mode 100755
index 0000000000..aca77d20f0
--- /dev/null
+++ b/challenge-191/jo-37/perl/ch-2.pl
@@ -0,0 +1,155 @@
+#!/usr/bin/perl -s
+
+use v5.16;
+use Test2::V0 '!float';
+use PDL;
+use Math::Prime::Util qw(forsetproduct forcomb);
+use experimental 'signatures';
+
+our ($tests, $examples, $ryser);
+
+run_tests() if $tests || $examples;	# does not return
+
+my $permanent = $ryser ? \&ryser : \&bbfg;
+
+die <<EOS unless @ARGV;
+usage: $0 [-examples] [-tests] [-ryser] [N]
+
+-examples
+    run the examples from the challenge
+ 
+-tests
+    run some tests
+
+-ryser
+    Use an implementation of the Ryser formula instead of
+    Balasubramanian-Bax-Franklin-Glynn.
+
+N
+    length of cute lists to be counted.
+
+EOS
+
+
+### Input and Output
+
+say count_cute(shift, $permanent);
+
+
+### Implementation
+
+# First approach:
+# Loop over permutations, check for cute lists and count these.  This
+# looks kind of insane as 15! = 1307674368000.  Not even tried it.
+# 
+# Second approach:
+# Loop over the possible values at a given position and then recurse to
+# the next position with a reduced list of values.  This works, but
+# isn't really satisfying.  Dropped it.
+#
+# Third approach:
+# Do not count the number of cute lists but calculate it instead.  A
+# cute list may be regarded as a complete matching between the numbers
+# from 1 to N and the N slots in the list.  A number is permitted in a
+# slot, if the corresponding element in the adjacency matrix is one.
+# The number of perfect matchings in an unweighted bipartite graph
+# defined by its adjacency matrix A equals the permanent of A. See [2].
+# For a cute list the adjacency matrix A is defined by:
+#   a(i,j) = 1 if i is divisible by j or j is divisible by i and
+#   a(i,j) = 0 otherwise.
+# Calculating the permanent of a matrix is rather expensive.  However,
+# utilizing the Balasubramanian-Bax-Franklin-Glynn formula (see [3]) the
+# computation for N>15 becomes possible.  The running time roughly
+# doubles when N is incremented by 1.  Here N=20 takes about 10 s and
+# N=28 takes about 40 min.  Beyond that the given PDL implementation
+# runs into longlong overflows on the "cute" adjacency matrix and will
+# produce false results at some point.
+#
+# NB: Tried some different ways of implementing the better known Ryser
+# formula.  None was able to keep up with B-B-F-G.  See below for the
+# best of these attempts (that has more than twice the running time
+# compared to the B-B-F-G implementation).
+#
+# NB2: What a shame! wxMaxima finds the solution for N=28 in 5 min.
+#
+# References:
+# [1] http://oeis.org/A320843
+# [2] https://en.wikipedia.org/wiki/Permanent_(mathematics)
+# [3] https://en.wikipedia.org/wiki/Computing_the_permanent
+
+sub count_cute ($n, $permanent) {
+    # Build the upper right triangle of the adjacency matrix A for a
+    # "cute list" of size N.
+    my $a = !((xvals(long, $n) + 1) % (yvals(long, 1, $n) + 1));
+    # Build the full adjacency matrix utilizing its symmetry.
+    $a->inplace->or2($a->xchg(0, 1), 0);
+
+    # Find the number of cute lists.
+    $permanent->($a);
+}
+
+# Calculate the permanent of a square matrix using the
+# Balasubramanian-Bax-Franklin-Glynn formula.
+sub bbfg ($a) {
+    my $n = $a->dim(1);
+    # Allocate the vector δ.  As the summands will become very large,
+    # calculations for N > 14 need to be done on longlong types coming
+    # in here.
+    my $d = zeroes longlong, $n;
+    my $dr = $d->get_dataref;
+    my $perm = longlong 0;
+    # Straightforward implementation of the
+    # Balasubramanian-Bax-Franklin-Glynn formula, see [3]:
+    forsetproduct {
+        # Set d to a single "1" followed by all of {-1, 1}**(N-1)
+        # Manipulating d's internal data to gain maximum speed.
+        $$dr = pack 'q*', 1, @_;
+        $d->upd_data;
+        # The outer sum:
+        $perm = $perm + prodover($d) * prodover($d x $a);
+    } ([-1, 1]) x ($n - 1);
+
+    # The permanent of A:
+    $perm->sclr / 2 ** ($n - 1);
+}
+
+# Calculate the permanent of a square matrix using the Ryser formula.
+sub ryser ($a) {
+    # Operate on longlong type.
+    $a = convert $a, longlong;
+    my $n = $a->dim(0);
+    my $minus_one = long -1;
+    my $perm = longlong 0;
+    # Straightforward implementation of the Ryser formula, see [3].
+    # Loop over all subsets of {0,...,N-1}
+    forcomb {
+        $perm = $perm + $minus_one**@_ * prodover sumover $a->dice(\@_);
+    } $n;
+
+    # The permanent of A:
+    (-1)**$n * $perm->sclr;
+}
+
+
+### Examples and tests
+
+sub run_tests {
+    SKIP: {
+        skip "examples" unless $examples;
+
+        is count_cute(2, \&bbfg), 2, 'Example';
+    }
+
+    SKIP: {
+        skip "tests" unless $tests;
+
+        # See [1].
+        is count_cute(10, \&bbfg),  700,   'N=10 bbfg';
+        is count_cute(10, \&ryser), 700,   'N=10 ryser';
+        is count_cute(15, \&bbfg),  24679, 'Last from task';
+        is count_cute(16, \&bbfg),  87328, 'From A320843';
+	}
+
+    done_testing;
+    exit;
+}