Collatz experiment

collatz.c

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <gmp.h>

#define S(n) mpz_get_str(NULL, 10, n)
#define EQ(z, n) (mpz_cmp_ui(z, n) == 0)

/* Print out collatz path for given number. */
int path(char *str)
{
    int count = 0;
    mpz_t n;
    mpz_init(n);
    mpz_set_str(n, str, 10);
    while (mpz_cmp_ui(n, 1) != 0) {
	count++;
	if (mpz_divisible_ui_p(n, 2)) {
	    mpz_cdiv_q_ui(n, n, 2);
	    printf("DIV (%s)\n", S(n));
	} else {
	    mpz_mul_ui(n, n, 3);
	    mpz_add_ui(n, n, 1);
	    printf("MUL (%s)\n", S(n));
	}
    }
    mpz_clear(n);
    return count;
}

/* Find the lowest number with collatz length len, starting at res. */
void find(mpz_t res, int len)
{
    if (len == 0) return;

    mpz_t n, m;
    mpz_init(n);
    mpz_init(m);

    /* Try n/2. */
    mpz_set(n, res);
    mpz_mul_ui(n, n, 2);
    find(n, len - 1);

    /* Try 3n+1. */
    int fm = 0;
    mpz_set(m, res);
    mpz_sub_ui(m, m, 1);
    if (mpz_cmp_ui(m, 3) > 0 && mpz_divisible_ui_p(m, 3)) {
	mpz_cdiv_q_ui(m, m, 3);
	if (!mpz_divisible_ui_p(m, 2)) {
	    fm = 1;
	    find(m, len - 1);
	}
    }

    /* Return best result. */
    if (fm && mpz_cmp(m, n) < 0) {
	mpz_set(res, m);
    } else {
	mpz_set(res, n);
    }
    mpz_clear(m);
    mpz_clear(n);
}

/* Return current time as a double. */
double now()
{
    struct timeval t;
    gettimeofday(&t, NULL);
    return t.tv_sec + (t.tv_usec / 1000000.0);
}

/* Perform a find() from 1, and print out timing information. */
void time_find(int d)
{
    mpz_t n;
    mpz_init(n);
    mpz_set_ui(n, 1);

    double t = now();
    find(n, d);
    t = now() - t;

    char *r = S(n);
    printf("%d\t%f\t%s\n", d, t, r);
    fflush(stdout);
    free(r);
    mpz_clear(n);
}

int main(int argc, char **argv)
{
    //printf("%d steps\n", path("2367363789863971985761"));

    int i = 0;
    while (++i) time_find(i);

    if (argc < 2) {
	fprintf(stderr, "error: need a depth argument\n");
	exit(1);
    }
    time_find(atoi(argv[1]));

    return 0;
}