part: add find2, find2c.c

Find partitions of an n-cube of side 2 into parts of size at most 2.
This finds A045310.

Author: hvds

part/find2

@@ -0,0 +1,48 @@
+#!/opt/maths/bin/perl -w
+use strict;
+@ARGV == 1 or die "Usage: $0 <dimension=0..6>\n";
+
+my($dim) = @ARGV;
+my @fullmask = (1);
+for (1 .. 6) {
+    my $prev = $fullmask[$_ - 1];
+    $fullmask[$_] = $prev | ($prev << (1 << ($_ - 1)));
+}
+
+my $count = count_d($dim, $fullmask[$dim]);
+printf "200 a(2, %s) = %s (%.2f)\n", $dim, $count, scalar times();
+exit 0;
+
+=head2 count_d ( d, mask )
+
+Calculates and returns the number of ways the specified shape can be
+partitioned into polyominoes of size at most 2. The shape is some or all
+of a I<d>-dimensional hypercube of side 2; I<mask> is a vector of C<2 ** d>
+bits specifying the unit I<d>-cubes that are included in the shape.
+
+The approach is recursive: we pick a dimension, split the shape in 2 over
+that dimension giving two new C<d-1>-dimensional shapes, and then try
+assigning size 2 polyominoes that straddle the break in each way possible.
+
+Calculation uses perl's inbuilt types without regard for overflow or loss
+of precision; for n=6 you'll need perl built with 64-bit integers.
+
+=cut
+
+sub count_d {
+    my($d, $mask) = @_;
+    return 1 if $d == 0;
+    my $hd = $d - 1;        # dimension of half shape
+    my $hp = 1 << $hd;      # number of unit I<hd>-cubes in half shape
+    my $halfmask = $fullmask[$hd];  # mask of all bits in half shape
+    my $left = $mask & $halfmask;   # mask representing one of the half shapes
+    my $right = $mask >> $hp;
+    my $shared = $left & $right;    # mask of cube-pairs straddling the break
+    my $x = 0;
+    my $sum = 0;
+    do {
+        $sum += count_d($hd, $left ^ $x) * count_d($hd, $right ^ $x);
+        $x = ($x - $shared) & $shared;  # cycle over the bits of $shared
+    } while $x != 0;
+    return $sum;
+}

part/find2c.c

@@ -0,0 +1,176 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <gmp.h>
+#include "clock.h"
+
+/* we rely on mask_t to be an integer type
+ * with at least 2 ** MAX_DIMENSION bits
+ */
+#define MAX_DIMENSION (6)
+typedef unsigned long long mask_t;
+
+mask_t fullmask[MAX_DIMENSION + 1] = {
+    0x0000000000000001,
+    0x0000000000000003,
+    0x000000000000000f,
+    0x00000000000000ff,
+    0x000000000000ffff,
+    0x00000000ffffffff,
+    0xffffffffffffffff
+};
+unsigned int step = 0;
+clock_t t0;
+
+mpz_t sum[MAX_DIMENSION + 1];
+mpz_t scratch[MAX_DIMENSION + 1];
+
+void init_gmp(void) {
+    int i;
+    for (i = 0; i <= MAX_DIMENSION; ++i) {
+        mpz_init(sum[i]);
+        mpz_init(scratch[i]);
+    }
+}
+
+void clear_gmp(void) {
+    int i;
+    for (i = 0; i <= MAX_DIMENSION; ++i) {
+        mpz_clear(sum[i]);
+        mpz_clear(scratch[i]);
+    }
+}
+
+/*
+ * count_part(d, mask): set sum[d] to the number of ways the specified
+ * shape can be partitioned into polyominoes of size at most 2.
+ * The shape is some or all of a I<d>-dimensional hypercube of side 2;
+ * I<mask> is a vector of C<2 ** d> bits specifying the unit I<d>-cubes
+ * that are included in the shape.
+ *
+ * Overwrites scratch[0..d] and sum[0..d] during calculation.
+ */
+void count_part(int d, mask_t mask) {
+    mask_t left, right, shared, u;
+    if (d == 0) {
+        mpz_set_ui(sum[d], 1);
+        return;
+    }
+    left = mask & fullmask[d - 1];
+    right = mask >> (1 << (d - 1));
+    shared = left & right;
+
+    u = 0;
+    mpz_set_ui(sum[d], 0);
+    while (1) {
+        count_part(d - 1, left ^ u);
+        mpz_set(scratch[d], sum[d - 1]);
+        count_part(d - 1, right ^ u);
+        mpz_mul(scratch[d], scratch[d], sum[d - 1]);
+        mpz_add(sum[d], sum[d], scratch[d]);
+        u = (u - shared) & shared;
+        if (u == 0)
+            break;
+    }
+}
+
+/*
+ * count_full(d): same as count_part(d, fullmask[d])
+ */
+void count_full(int d) {
+    mask_t lim, u;
+    if (d == 0) {
+        mpz_set_ui(sum[d], 1);
+        return;
+    }
+    lim = fullmask[d - 1];
+    mpz_set_ui(sum[d], 0);
+    for (u = 0; u <= lim; ++u) {
+        count_part(d - 1, u);
+        mpz_mul(scratch[d], sum[d - 1], sum[d - 1]);
+        mpz_add(sum[d], sum[d], scratch[d]);
+        if ((++step & 0xfffff) == 0) {
+            gmp_printf("300 at (%d, %llx) sum is %Zu (%.2f)\n",
+                    d, u, sum[d], difftime(t0, curtime()));
+        }
+    }
+}
+
+#define LIM4 (1 << (1 << 4))
+unsigned int cache4[LIM4];  /* the sum of these is 41025, int is ample room */
+void init_cache4(void) {
+    int u;
+    for (u = 0; u < LIM4; ++u) {
+        count_part(4, (mask_t)u);
+        cache4[u] = mpz_get_ui(sum[4]);
+    }
+}
+
+/*
+ * count_5_part(mask): same as count_part(5, mask)
+ *
+ * Uses a prepared cache of results for d=4, to minimize recursion.
+ */
+void count_5_part(mask_t mask) {
+    mask_t left, right, shared, u;
+    left = mask & fullmask[4];
+    right = mask >> (1 << 4);
+    shared = left & right;
+    u = 0;
+    mpz_set_ui(sum[5], 0);
+    while (1) {
+        mpz_add_ui(sum[5], sum[5], cache4[left ^ u] * cache4[right ^ u]);
+        u = (u - shared) & shared;
+        if (u == 0)
+            break;
+    }
+}
+
+/*
+ * count_6_full(): same as count_full(6)
+ *
+ * Uses a prepared cache of results for d=4 (via count_5_part()) to reduce
+ * recursion at a reasonable memory cost.
+ */
+void count_6_full(void) {
+    mask_t lim, u;
+    init_cache4();
+    lim = fullmask[5];
+    mpz_set_ui(sum[6], 0);
+    for (u = 0; u <= lim; ++u) {
+        count_5_part(u);
+        mpz_mul(scratch[6], sum[5], sum[5]);
+        mpz_add(sum[6], sum[6], scratch[6]);
+        if ((++step & 0x3ffffff) == 0) {
+            gmp_printf("300 at (%d, %llx) sum is %Zu (%.2f)\n",
+                    6, u, sum[6], difftime(t0, curtime()));
+        }
+    }
+}
+
+/*
+ * Calculates the number of ways a I<d>-dimensional hypercube of side 2
+ * can be partitioned into polyominoes of size at most 2.
+ */
+int main(int argc, char** argv) {
+    int d;
+
+    if (argc != 2) {
+        fprintf(stderr, "usage: %s <dimension>\n", argv[0]);
+        return 1;
+    }
+    d = atoi(argv[1]);
+    if (d < 0 || d > MAX_DIMENSION) {
+        fprintf(stderr, "dimension must be in the range 0 to %u\n",
+                MAX_DIMENSION);
+        return 1;
+    }
+
+    setup_clock();
+    t0 = curtime();
+    init_gmp();
+    (d == 6) ? count_6_full() : count_full(d);
+    gmp_printf("200 a(2, %d) = %Zu (%.2f)\n",
+            d, sum[d], difftime(t0, curtime()));
+    clear_gmp();
+    return 0;
+}