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

/*
check.c:  Computes the sum S_Y as described in "A Test for
Identifying Fourier Coefficients of Automorphic Forms and Application
to Kloosterman Sums."

Andrew Booker (arbooker@math.princeton.edu)
Run the program with no arguments for details.
*/

#define F(t) ((t)*(t)*exp(-(t)))
struct plist { int prime; double a; } *l;
struct ylist { double sum,yinv; } *ystart, *yend;
int n;

/* Perform the sum recursively.  The Hecke relations are done as we go. */
static void sum(int m,struct plist *p,double co,double a1,double a2) {
	int end = n/m;
	double temp = co * a1;
	struct ylist *y;

	for (y=ystart;y<yend;y++)
		y->sum += temp * F(y->yinv*m);
	if (p->prime <= end) {
		sum(m * p->prime,p,co,p->a*a1-a2,a1);
		while ((++p)->prime <= end)
			sum(m * p->prime,p,temp,p->a,1.0);
	}
}

int main(int argc,char *argv[]) {
	int i,j,k,level=1,q=1,sqrtn,n2;
	struct plist *p;
	unsigned int *sieve;
	double ymin,ymax,t;
	int sign = 0;

	/* read parameters from command line */
	if (argc < 5) {
		printf("Usage: %s [-n] <Ymin> <Ymax> <k> <n> [level] < inputfile\n"
"Description:\n"
"Computes S_Y at k geometrically spaced values of Y between Ymin and Ymax,\n"
"including all terms a(i) with i <= n and i relatively prime to the level.\n"
"The input file should contain a(p) in IEEE double precision format for\n"
"each prime p (including primes dividing the level, although those values\n"
"are not used).  Specify -n to multiply all prime coefficients by -1.\n",
			argv[0]);
		return 0;
	}
	if (!strcmp(argv[1],"-n")) {
		sign = 1;
		--argc, ++argv;
	}
	sscanf(argv[1],"%lf",&ymin);
	sscanf(argv[2],"%lf",&ymax);
	sscanf(argv[3],"%d",&k);
	sscanf(argv[4],"%d",&n);
	if (argc > 5)
		sscanf(argv[5],"%d",&level);

	/* compute k geometrically spaced Y values from ymin to ymax */
	ystart = malloc(sizeof(*ystart)*k);
	yend = &ystart[k];
	if (k == 1)
		ystart->sum = 0.0, ystart->yinv = 1.0/ymin;
	else {
		ymin = log(ymin), ymax = log(ymax), t = (ymax-ymin)/(k-1);
		for (i=k-1;i>=0;i--)
			ystart[i].sum = 0.0, ystart[i].yinv = exp(t*i-ymax);
	}

	/* sieve to find all odd primes <= n */
	sieve = calloc(1+(n>>6),sizeof(*sieve));
	sqrtn = (int)sqrt((double)n), n2 = (n-1)>>1;
	for (i=3;i<=sqrtn;i+=2)
		if ( !(sieve[i>>6] & (1 << ((i & 63) >> 1))) )
			for (j=i*i>>1;j<=n2;j+=i)
				sieve[j>>5] |= 1 << (j&31);

	/* allocate space for >= pi(n) primes and a(p)'s */
	t = 1.0/log((double)n);
	p = l = malloc(sizeof(*l)*(int)(n*t*(1.0+t+2.6*t*t)));

	/* read a(2) */
	fread(&p->a,sizeof(double),1,stdin);
	if (sign) p->a = -p->a;
	if (level & 1)         /* use only primes not dividing level */
		(p++)->prime = 2;
	else if (level & 2)
		q = 3;

	/* read a(p), p > 2 */
	for (i=3;i<n;i+=2)
		if ( !(sieve[i>>6] & (1 << ((i & 63) >> 1))) ) {
			fread(&p->a,sizeof(double),1,stdin);
			if (sign) p->a = -p->a;
			if (level % i)
				(p++)->prime = i;
			else if ((level/i) % i)
				q *= (1+i);
		}
	p->prime = n+1;

	/* do the sum */
	sum(1,l,1.0,1.0,0.0);

	/* print the results */
	printf("        Y           normalized sum     eigenvalue >=\n");
	for (i=k-1;i>=0;i--) {
		ystart[i].sum *= sqrt(ystart[i].yinv);
		t = 42.88 * sqrt(fabs(ystart[i].sum)) / (level * ystart[i].yinv)-3.0;
		if (t < 0.25) t = 0.25;
		printf("%14.5f%20.12f%18.2f\n",1.0/ystart[i].yinv,ystart[i].sum,t);
	}
	return 0;
}