/**
 * Author:      Jesse Short-Gershman    * Tight design parameter searcher
 * Email:       jesseasg@yahoo.ca       * File: tightsearch.c
 * Date:        May 24, 2011            * Module: Executable
 *
 * Recommended compilation: gcc -Wall -O3 tightsearch.c -o ts -lm.
 *
 * Usage: ts [s >= 2] [beta_0 > 0] [derivative threshold >= 0] [v]
 * The first two parameters must be included in the order specified. Each of the
 * last two parameters is optional, and if both are included then they can be in
 * either order. Including the "v" flag causes verbose output; each of the pairs
 * (v,k) such that a_2/a_0 is an integer is printed along with the corresponding beta.
 * Examples of arguments: [5 40], [5 40 v], [5 40 .25], [5 40 .25 v], [5 40 v .25].
 *
 * Purpose: Check all of the values of v and k such that v >= 2k+1 and beta <= beta_0
 * to determine if there could possibly be a tight 2s-design with such parameters.
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#ifdef _WIN32
# define PRI64 "I64"
#else
# define PRI64 "ll"
#endif

#define DEFAULT_THRESHOLD(b) (.15+.002*(b))

typedef struct {
	long long num, den;
} Fraction;

long long ftoi(Fraction a) {
	return a.num/a.den;
}

double ftod(Fraction a) {
	return (double)a.num/a.den;
}

long long gcd(long long a, long long b) {
	while (b != 0) {
		long long temp = b;
		b = a%b;
		a = temp;
	}
	return a;
}

void reduce(Fraction *a) {
	long long d = gcd(a->num, a->den);
	a->num /= d;
	a->den /= d;
}

Fraction newfrac(long long a, long long b) {
	Fraction ans = {a,b};
	reduce(&ans);
	return ans;
}

Fraction multbyint(Fraction a, int n) {
	return newfrac(a.num*n, a.den);
}

Fraction mult(Fraction a, Fraction b) {
	return newfrac(a.num*b.num, a.den*b.den);
}

Fraction dividellby(long long n, Fraction a) {
	return newfrac(a.den*n, a.num);
}

void addto(Fraction *a, Fraction *b) {
	a->num = a->num*b->den+a->den*b->num;
	a->den = a->den*b->den;
	reduce(a);
}

void pf(Fraction a) {
	printf("%"PRI64"d/%"PRI64"d",a.num,a.den);
}

void hyphens(int n) {
	int i;
	for (i = 0; i < n; i++)
		printf("-");
}

void print_runtime(double seconds, char *alg_name) {
	char temp[30];
	sprintf(temp,"%.3f",seconds);
	int num_hyphens = strlen(alg_name)+strlen(temp)+17;
	printf("+"); hyphens(num_hyphens);
	printf("+\n| %s took %.3f seconds. |\n+",alg_name,seconds);
	hyphens(num_hyphens); printf("+\n");
}

unsigned long long min(unsigned long long a, unsigned long long b) {
	return a <= b ? a : b;
}

unsigned long long max(unsigned long long a, unsigned long long b) {
	return a >= b ? a : b;
}

double round(double x) {
	return floor(x+0.5);
}

unsigned long long factorial(int n) {
	unsigned long long ans = 1;
	int i;
	for (i = 1; i <= n; i++)
		ans *= i;
	return ans;
}

unsigned long long choose(unsigned long long n, int k) {
	if (k < 0 || k > n)
		return 0;
	unsigned long long ans = 1;
	int t;
	for (t = 0; t < min(k,n-k); t++)
		ans = ans*(n-t)/(t+1);
	return ans;
}

int s, sc2, verbose, any_possible, vktot;
unsigned long long attot;

Fraction psi(unsigned long long v, unsigned long long k, unsigned long long x) {
	Fraction ans = {0,1}, temp;
	int i;
	for (i = 0; i <= s; i++) {
		temp = newfrac(((s-i)&1 ? -1 : 1)*choose(v-s,i)*choose(k-i,s-i)*choose(k-1-i,s-i)*choose(x,i),choose(s,i));
		addto(&ans,&temp);
	}
	return ans;
}

double g(double n, double a) {
	return sc2*a*(n+1)*(n+2)/((n*n+n)/a+1);
}

double ginverse(unsigned long long c, double a) {
	return (3*sc2*a-c/a+sqrt(pow(3*sc2*a-c/a,2)-4*(c/a-sc2*a)*(c-2*sc2*a)))/(2*(c/a-sc2*a));
}

double dgdn(double n, double a) {
	return 2*sc2*a/((n*n+n)/a+1)-sc2*(2*n-a+2)*(2*n+1)/pow((n*n+n)/a+1,2);
}

void check_tight_design(unsigned long long n, Fraction a) {
	if (dividellby(n*n+n,a).den == 1) {
		unsigned long long k = n+s;
		unsigned long long v = ftoi(dividellby(n*n+n,a))+2*s-1;
		if (verbose) {
			char temp[20];
			sprintf(temp,"%f",(1-ftod(a)/n)*sqrt(ftod(a)));
			printf("| %-10"PRI64"u%-8"PRI64"u%.8s |",v,k,temp);
		}
		vktot++;
		if ((choose(v,s)*choose(k,2*s)) % choose(v,2*s) == 0) {
			if (!verbose)
				printf("v = %"PRI64"u, k = %"PRI64"u:",v,k);
			printf(" The lambda works");
			unsigned long long l = (choose(v,s)*choose(k,2*s))/choose(v,2*s);
			int check_poly = 1;
			int u;
			for (u = 1; u < 2*s; u++) {
				if ((l*choose(v-u,2*s-u)) % choose(k-u,2*s-u) != 0) {
					check_poly = 0;
					break;
				}
			}
			if (check_poly) {
				printf(", checked the polynomial");
				long long a_s = (s&1 ? -1 : 1)*choose(k,s)*choose(k-1,s);
				Fraction a_0 = newfrac(choose(v-s,s),factorial(s));
				if (dividellby(a_s,a_0).den == 1) {
					a_s = ftoi(dividellby(a_s,a_0));
					int count = 0;
					unsigned long long x;
					for (x = 1; x < k; x++) {
						if (a_s % x == 0 && psi(v,k,x).num == 0)
							count++;
					}
					if (count == s) {
						printf(", a tight %d-(%"PRI64"u,%"PRI64"u,%"PRI64"u) design might exist",2*s,v,k,l);
						any_possible = 1;
					}
				}
			}
			printf(".\n");
		} else if (verbose)
			printf("\n");
	}
}

int main(int argc, char *argv[]) {
	s = atoi(argv[1]);
	if (s < 2) {
		fprintf(stderr,"%s: s must be an integer greater than or equal to 2\n",argv[0]); exit(1);
	}
	double b = atof(argv[2]);
	if (b <= 0) {
		fprintf(stderr,"%s: beta_0 must be a positive real number\n",argv[0]); exit(1);
	}
	double threshold;
	if (argc >=4) {
		if ((threshold = atof(argv[3]))) {
			if (argc >= 5) {
				if (strcmp(argv[4], "v") == 0) {
					verbose = 1;
					if (argc >= 6) {
						fprintf(stderr,"%s: invalid argument: %s\n",argv[0],argv[5]); exit(1);
					}
				} else {
					fprintf(stderr,"%s: invalid argument: %s\n",argv[0],argv[4]); exit(1);
				}
			} else
				verbose = 0;
		} else {
			if (strcmp(argv[3], "v") == 0) {
				verbose = 1;
				if (argc >= 5) {
					if ((threshold = atof(argv[4]))) {
						if (argc >= 6) {
							fprintf(stderr,"%s: invalid argument: %s\n",argv[0],argv[5]); exit(1);
						}
					} else {
						fprintf(stderr,"%s: invalid argument: %s\n",argv[0],argv[4]); exit(1);
					}
				} else
					threshold = DEFAULT_THRESHOLD(b);
			} else {
				fprintf(stderr,"%s: invalid argument: %s\n",argv[0],argv[3]); exit(1);
			}
		}
	} else {
		verbose = 0;
		threshold = DEFAULT_THRESHOLD(b);
	}
	if (verbose) {
		printf("+"); hyphens(28);
		printf("+\n| %-10s%-8s%-8s |\n+", "v", "k", "beta");
		hyphens(28); printf("+\n");
	}
	clock_t start_time = clock();
	sc2 = choose(s,2);
	vktot = attot = any_possible = 0;
	Fraction a = {1,s}, one_over_s = {1,s};
	int i;
	unsigned long long n, n_min, c, c_max, m, d;
	double n_max, ad = ftod(a);
	for (; ad < 4*b*b; addto(&a,&one_over_s), ad = ftod(a)) {
		n_min = max(s,ftoi(multbyint(a,2))+1);
		c = ftoi(mult(a,multbyint(a,sc2)))+2;
		n_max = ginverse(c-1,ad);
		c_max = (unsigned long long)g((double)n_min,ad);
		if (ad > b*b && ad/(1-b/sqrt(ad)) < n_max) {
			n_max = ad/(1-b/sqrt(ad));
			c = (unsigned long long)g(n_max,ad)+1;
			if ((fabs(n_max-round(n_max)) >= 1e-7 || fabs(g(n_max,ad)-round(g(n_max,ad))) >= 1e-7) && c <= c_max)
				n_max = ginverse(c++,ad);
		}
		while (-dgdn(n_max,ad) < threshold && c <= c_max && n_max >= n_min) {
			for (i = 0; i < 50 && c <= c_max && n_max >= n_min; i++, attot++) {
				if (fabs(n_max-round(n_max)) < 1e-7)
					check_tight_design((unsigned long long)round(n_max),a);
				n_max = ginverse(c++,ad);
			}
		}
		m = a.num;
		d = a.den;
		for (n = n_min; n <= n_max; n++, attot++) {
			if ((s*(s-1)*m*m*(n*n+3*n+2)) % (2*d*(n*n*d+n*d+m)) == 0)
				check_tight_design(n,a);
		}
	}
	if (verbose) {
		printf("+"); hyphens(28); printf("+\n");
	}
	printf("Total pairs (v,k) checked: %d.\n",vktot);
	printf("Total pairs (alpha,t) enumerated: %"PRI64"u.\n",attot);
	char alg_name[70];
	sprintf(alg_name, "tightsearch with s = %d, beta_0 = %g and threshold = %g",s,b,threshold);
	print_runtime((clock()-(double)start_time)/CLOCKS_PER_SEC, alg_name);
	if (!any_possible)
		printf("There are no nontrivial tight %d-designs with beta <= %g.\n",2*s,b);
	return 0;
}