wilson: C implementation of depth-first search

Given a depth limit, this will search for best solutions.

wilson/Makefile

@@ -1,2 +1,7 @@
+all: search-breadth search-depth
+
 search-breadth: search-breadth.c breadth.c breadth.h mygmp.c mygmp.h
 	gcc -g -O6 -Wall -o search-breadth search-breadth.c breadth.c mygmp.c -lgmp
+
+search-depth: search-depth.c depth.c depth.h mygmp.c mygmp.h
+	gcc -g -O6 -Wall -o search-depth search-depth.c depth.c mygmp.c -lgmp

wilson/depth.c

@@ -0,0 +1,151 @@
+#include "depth.h"
+
+mpq_t limit_calc;
+mpz_t limit_div;
+
+extern mpq_t r;
+extern mpq_t rone;
+extern mpq_t limit;
+
+typedef struct {
+    mpq_t q;
+    ulong count;
+    ulong actual;
+    ulong best_bits_num;
+    ulong best_bits_den;
+} frame_t;
+
+frame_t *stack;
+ulong max_depth;
+bool solved;
+ulong count;
+
+void init_depth(ulong depth) {
+    ulong i;
+    max_depth = depth;
+    /* we count from 0 to n */
+    stack = (frame_t *)malloc((depth + 1) * sizeof(frame_t));
+    for (i = 0; i <= max_depth; ++i) {
+        QINIT(&stack[i].q, "stack[%lu].q", i);
+        stack[i].actual = (ulong)0;
+        stack[i].best_bits_num = (ulong)0;
+        stack[i].best_bits_den = (ulong)0;
+    }
+    solved = 0;
+    QINIT(&limit_calc, "depth limit_calc");
+    ZINIT(&limit_div, "depth limit_div");
+}
+
+void finish_depth(void) {
+    ulong i;
+    ZCLEAR(&limit_div, "depth limit_div");
+    QCLEAR(&limit_calc, "depth limit_calc");
+    for (i = 0; i <= max_depth; ++i)
+        QCLEAR(&stack[i].q, "stack[%lu].q", i);
+    free(stack);
+}
+
+void report_depth(ulong depth) {
+    ulong i;
+    gmp_printf("%Qd is solved in %lu steps by [ ", r, depth);
+    for (i = depth; i > 0; --i)
+        printf("%lu ", stack[i].count);
+    printf("]\n");
+}
+
+extern double timing(void);
+void report_progress(ulong depth) {
+    printf("tried %lu values (%.2fs)\n", count, timing());
+}
+
+void report_final(void) {
+    ulong i;
+    for (i = 0; i < max_depth; ++i) {
+        frame_t *f = &stack[i];
+        gmp_printf("%Qd - g > %lu: %lu, best_bits %lu/%lu\n",
+                r, i + 1, f->actual, f->best_bits_num, f->best_bits_den);
+    }
+}
+
+bool try_depth(ulong depth) {
+    frame_t *cur = &(stack[depth]);
+    frame_t *next = &(stack[depth + 1]);
+    int final = (depth + 1 >= max_depth) ? 1 : 0;
+    bool locally_solved = 0;
+    ulong bits_num, bits_den;
+
+    ++count;
+    if ((count & 0x3ffffff) == 0)
+        report_progress(depth);
+
+    mpq_inv(next->q, cur->q);
+
+    /* if final, we only care whether we can reach 1 */
+    if (final) {
+        mpq_sub(next->q, next->q, rone);
+        if (mpq_sgn(next->q) < 0)
+            return 0;
+        mpq_div(next->q, next->q, r);
+        if (mpz_cmp_ui(mpq_denref(next->q), (ulong)1) != 0)
+            return 0;
+        next->count = mpz_strict_get_ui(mpq_numref(next->q));
+        report_depth(depth + 1);
+        return 1;
+    }
+
+    if (mpq_cmp(next->q, limit) <= 0) {
+        next->count = (ulong)0;
+    } else {
+        /* count = floor((q - limit) / r) */
+        mpq_sub(limit_calc, next->q, limit);
+        mpq_div(limit_calc, limit_calc, r);
+        mpz_fdiv_q(limit_div, mpq_numref(limit_calc), mpq_denref(limit_calc));
+        next->count = mpz_strict_get_ui(limit_div);
+
+        /* q -= count * r */
+        mpq_set(limit_calc, r);
+        mpz_mul_ui(mpq_numref(limit_calc), mpq_numref(limit_calc), next->count);
+        mpq_canonicalize(limit_calc);
+        mpq_sub(next->q, next->q, limit_calc);
+    }
+
+    /* now we've guaranteed curq - r <= limit, so all values we find are
+     * useful.
+     */
+    /* while ((q -= r) > 0) { ... } */
+    while (1) {
+        mpq_sub(next->q, next->q, r);
+        if (mpq_sgn(next->q) <= 0)
+            return locally_solved;
+        ++next->count;
+        if (mpq_equal(next->q, rone)) {
+            report_depth(depth + 1);
+            max_depth = depth + 1;
+            solved = 1;
+            return 1;
+        }
+        ++next->actual;
+        bits_num = mpz_bitsize(mpq_numref(next->q));
+        if (!next->best_bits_num || next->best_bits_num >= bits_num) {
+            next->best_bits_num = bits_num;
+            bits_den = mpz_bitsize(mpq_denref(next->q));
+            if (!next->best_bits_den || next->best_bits_den > bits_den)
+                next->best_bits_den = bits_den;
+        }
+
+        /* and recurse */
+        if (try_depth(depth + 1))
+            locally_solved = 1;
+    }
+    /* not reached */
+}
+
+bool search_depth(void) {
+    bool result;
+    mpq_set(stack[0].q, rone);
+    stack[0].count = 0;
+    stack[0].actual = 1;
+    result = try_depth(0);
+    report_final();
+    return result;
+}

wilson/depth.h

@@ -0,0 +1,10 @@
+#ifndef DEPTH_H
+#define DEPTH_H
+
+#include "mygmp.h"
+
+extern void init_depth(ulong depth);
+extern void finish_depth(void);
+extern bool search_depth(void);
+
+#endif

wilson/search-depth.c

@@ -0,0 +1,61 @@
+#include <unistd.h>
+#include <sys/times.h>
+#include "mygmp.h"
+#include "depth.h"
+
+long clock_tick;
+
+mpq_t r;        /* the rational we're testing */
+mpq_t rone;     /* handy 1 */
+mpq_t limit;    /* r + 1/r */
+
+double timing(void) {
+    struct tms ttd;
+    times(&ttd);
+    return ((double)ttd.tms_utime) / clock_tick;
+}
+
+void init(int p, int q) {
+    QINIT(&r, "r");
+    mpq_set_ui(r, (ulong)p, (ulong)q);
+    QINIT(&rone, "rone");
+    mpq_set_ui(rone, (ulong)1, (ulong)1);
+    QINIT(&limit, "limit");
+    mpq_inv(limit, r);
+    mpq_add(limit, limit, r);
+}
+
+void finish(void) {
+    QCLEAR(&limit, "limit");
+    QCLEAR(&rone, "rone");
+    QCLEAR(&r, "r");
+}
+
+int main(int argc, char** argv) {
+    int p, q, d;
+
+    if (argc != 4) {
+        fprintf(stderr, "Usage: search-depth <p> <q> <depth>\n");
+        return 1;
+    }
+    p = atoi(argv[1]);
+    q = atoi(argv[2]);
+    d = atoi(argv[3]);
+    if (!(p > 0 && q > p)) {
+        fprintf(stderr, "Value error: need 0 < p/q < 1\n");
+        return 1;
+    }
+    if (!(d > 1)) {
+        fprintf(stderr, "Value error: need depth > 1\n");
+        return 1;
+    }
+
+    clock_tick = sysconf(_SC_CLK_TCK);
+    init(p, q);
+    init_depth((ulong)d);
+    if (!search_depth())
+        gmp_printf("%Qd - g > %d (%.2fs)\n", r, d, timing());
+    finish_depth();
+    finish();
+    return 0;
+}