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

void generate_primes(unsigned int primes[], int n);
void recurse(mpf_t result, unsigned int primes[], int pos);

int main(int argc, char *argv[])
{
    if (argc == 2) {
	mpf_t result, partial;
	int i;
	int prime_count = atoi(argv[1]);
	if (prime_count < 1)
	    prime_count = 1;
	unsigned int primes[prime_count];

	generate_primes(primes, prime_count);
	mpf_set_default_prec(64);
	mpf_init(result);
	mpf_init(partial);

	// Sum up percentage of used numbers
	// Time is O(2^(prime_count+1))
	for (i = 0; i < prime_count; i++) {
	    recurse(partial, primes, i);
	    mpf_add(result, result, partial);
	}

	mpf_ui_sub(result, 1, result);
	mpf_mul_ui(result, result, 100);

	gmp_printf("Best guess: %Fg%%\n", result);

	mpf_clear(partial);
	mpf_clear(result);

	return 0;
    }

    fprintf(stderr, "Usage: %s primes_to_use\n", argv[0]);
    return 1;
}

void recurse(mpf_t result, unsigned int primes[], int pos)
{
    int sign, i;
    mpz_t mask, total;
    mpf_t accum, num, signf;

    mpf_init(accum);
    mpf_init(num);
    mpf_init(signf);
    mpz_init(mask);
    mpz_init(total);

    mpz_ui_pow_ui(total, 2, pos);

    mpf_set_ui(num, 1);
    mpf_div_ui(num, num, primes[pos]);

    mpz_init_set_ui(mask, 0);

    mpf_set_ui(result, 0);
    for (mpz_init_set_ui(mask, 0);
	 mpz_cmp(mask, total) < 0; mpz_add_ui(mask, mask, 1)) {
	mpf_set_ui(accum, 1);
	for (sign = 1, i = 0; i < pos; i++) {
	    if (mpz_tstbit(mask, i)) {
		mpf_div_ui(accum, accum, primes[i]);
		sign = -sign;
	    }
	}
	mpf_set_si(signf, sign);
	mpf_mul(accum, accum, signf);
	mpf_add(result, result, accum);
    }
    mpf_mul(result, result, num);

    mpf_clear(signf);
    mpf_clear(num);
    mpf_clear(accum);
    mpz_clear(total);
    mpz_clear(mask);
}


void generate_primes(unsigned int primes[], int n)
{
    int l = 3, c = 1, i, prime;

    if (n > 0) {
	primes[0] = 2;
	if (n > 1) {
	    for (c = 1; c < n; l += 2) {
		prime = 1;
		for (i = 0; i < c; i++) {
		    if (l % primes[i] == 0) {
			prime = 0;
			break;
		    }
		}
		if (prime) {
		    primes[c] = l;
		    c++;
		}
	    }
	}
    }
}