/*
 *  fractran.c
 *  Copyright (C) 2002  Frank Buá (fb@frank-buss.de)
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You can get the GNU General Public License at
 *  http://www.gnu.org/licenses/gpl.html
 */

#include <stdio.h>
#include <time.h>
#include <gmp.h>

int main()
{
	unsigned long numerators[] =   { 17, 78, 19, 23, 29, 77, 95, 77,  1, 11, 13, 15,  1, 55 };
	unsigned long denominators[] = { 91, 85, 51, 38, 33, 29, 23, 19, 17, 13, 11,  2,  7,  1 };
	char* start = "2";
	unsigned long max_count = 100000000;

	clock_t start_time, stop_time;
	mpz_t accu, quot;
	unsigned long count = 0;
	mpz_init(quot);
	mpz_init_set_str(accu, start, 10);
	start_time = clock();
	while (count < max_count) {
		unsigned long remainder, pos;
		int i;
		count++;
		for (i = 0; i < sizeof(numerators); i++) {
			remainder = mpz_tdiv_q_ui(quot, accu, denominators[i]); 
			if (!remainder) {
				mpz_mul_ui(accu, quot, numerators[i]);
				break;
			}
		}
		if (remainder) {
			printf("remainder > 0, stopping\n");
			break;
		}
		pos = mpz_scan1(accu, 0);
		if (mpz_sizeinbase(accu, 2) == pos + 1) {
			printf("steps = %i, accumulator = 2^%i, fraction = %i (%i/%i)\n", count, pos, i + 1, numerators[i], denominators[i]);
		}
	}
	stop_time = clock();
	printf("running %i steps in %i seconds.\n", count, (stop_time - start_time) / CLOCKS_PER_SEC);
	fflush(stdout);
	return 0;
}