#include <stdio.h>
#include <math.h>
#include "mpfi.h"
#include "mpfi_io.h"

/* include sin and cos, which aren't yet implemented in mpfi */
/* these could be much more efficient... oh well             */
int mpfi_cos(mpfi_t rop,mpfi_srcptr op) {
	mpfi_t temp;
	mpz_t  left_z,right_z;
	mpfr_t temp_fr;
	mpfi_t result;

	mpfi_init2(result,mpfi_get_prec(rop));
	mpfi_init2(temp,mpfi_get_prec(rop));
	mpfr_init2(temp_fr,mpfi_get_prec(rop));
	mpz_init(left_z); mpz_init(right_z);

	mpfr_cos(temp_fr,&op->left,GMP_RNDD);
	mpfi_set_fr(result,temp_fr);
	mpfr_cos(temp_fr,&op->left,GMP_RNDD);
	mpfi_put_fr(result,temp_fr);
	mpfr_cos(temp_fr,&op->right,GMP_RNDU);
	mpfi_put_fr(result,temp_fr);
	mpfr_cos(temp_fr,&op->right,GMP_RNDU);
	mpfi_put_fr(result,temp_fr);

	mpfi_const_pi(temp);
	mpfi_div(temp,op,temp);
	mpfr_get_z(left_z,&temp->left,GMP_RNDD);
	mpfr_get_z(right_z,&temp->right,GMP_RNDU);

	while (mpz_cmp(left_z,right_z) <= 0) {
		if (mpfi_is_inside_z(left_z,temp))
			mpfi_put_si(result,mpz_fdiv_ui(left_z,2) ? -1 : 1);
		mpz_add_ui(left_z,left_z,1);
	}

	mpfi_set(rop,result);
	mpfi_clear(result);
	mpfi_clear(temp); mpfr_clear(temp_fr);
	mpz_clear(left_z); mpz_clear(right_z);
	return 3; /* both endpoints inexact */
}

int mpfi_sin(mpfi_t rop,mpfi_srcptr op) {
	mpfi_t temp;

	mpfi_init2(temp,mpfi_get_prec(op));
	mpfi_const_pi(temp);
	mpfi_div_ui(temp,temp,2);
	mpfi_sub(temp,temp,op);
	mpfi_cos(rop,temp);
	mpfi_clear(temp);
	return 3;
}

/* the function mpfi_inp_str seems to be broken */
void myread(mpfi_t rop,FILE *fp) {
	mpfr_t temp_fr;
	mpfr_init2(temp_fr,mpfi_get_prec(rop)+128);
	mpfr_inp_str(temp_fr,fp,10,GMP_RNDN);
	mpfi_set_fr(rop,temp_fr);
	mpfr_clear(temp_fr);
}

/* we'll have d < dmax, M < Mmax, N < Nmax */
#define dmax 64
#define Mmax 64
#define Nmax 256
int d,M,N;

/* sample points: x,y and [xmin,xmax], [ymin,ymax] */
mpfi_t sample[4][Nmax];

/* rectangular derivative variables */
mpfi_t fderiv[Nmax][dmax][dmax],fderiv_nbhd[Nmax][dmax][dmax];
mpfi_t s[Mmax][dmax],w[Mmax][dmax];

/* the precision to use for reading constants,
   must be sufficiently high for all bessel function computations! */
#define const_prec 2048

/* polar conversion variables */
mpz_t P[dmax][dmax][dmax],Q[dmax][dmax][dmax],Qnext[dmax][dmax][dmax];
mpfi_t xp[Nmax][dmax],yp[Nmax][dmax],nxp[Nmax][dmax],nyp[Nmax][dmax];

/* variables depending only on the eigenfunction */
mpfi_t R,R2,lambda,sqrt_lambda,bessel_const_mag,bessel_const_phase,A[Mmax];
int parity;

/* Bessel function extrema */
struct {
	mpfi_t u,w;
} extrema0[16],extrema1[16];
int ecount0,ecount1;

/* this must be sufficiently large to catch all extrema */
#define bessel_points 63

mpz_t  c[dmax][dmax];     /* recurrence constants c_{lj} */
mpfi_t twopi;             /* 2*Pi with default precision */

/* read parameters and compute recurrence constants */
void initialize(FILE *fp) {
	int j,l,n;
	mpfi_t temp;

	/* initialize global constants */
	/* Bessel function constants must be treated as highly precise,
	   even though they may not be accurate to this precision.
		 This is so that the power series is not overrun with error. */
	mpfi_init2(R,const_prec);
	mpfi_init2(R2,const_prec);
	mpfi_init2(lambda,const_prec);
	mpfi_init2(sqrt_lambda,const_prec);
	mpfi_init2(bessel_const_mag,const_prec);
	mpfi_init2(bessel_const_phase,const_prec);
	mpfi_init(twopi); /* 2Pi is computed in default precision */

	fscanf(fp,"%d",&parity);
	printf("PARITY: %d\n",parity);
	myread(R,fp);
	myread(bessel_const_mag,fp);
	myread(bessel_const_phase,fp);
	mpfi_init(temp);
	for (n=1;n<=M;n++) {
		mpfi_init(A[n]);
		myread(A[n],fp);
		mpfi_set_ui(temp,n);
		mpfi_sqrt(temp,temp);
		mpfi_div(A[n],A[n],temp);
	}
	mpfi_clear(temp);

	mpfi_sqr(R2,R);
	mpfi_add_d(lambda,R2,0.25);
	mpfi_sqrt(sqrt_lambda,lambda);

	mpfi_const_pi(twopi);
	mpfi_mul_ui(twopi,twopi,2);

	/* compute c_{lj} */
	mpz_init_set_si(c[2][0],-1);
	for (l=2;l<d;l++) {
		mpz_init(c[l+1][0]);
		mpz_mul_si(c[l+1][0],c[l][0],-l);
		for (j=1;j<=l-2;j++) {
			mpz_init(c[l+1][j]);
			mpz_mul_si(c[l+1][j],c[l][j],j-l);
			mpz_add(c[l+1][j],c[l+1][j],c[l][j-1]);
		}
		mpz_init_set_si(c[l+1][l-1],-1);
	}

	/* initialize s and w arrays */
	for (n=1;n<=M;n++)
		for (l=0;l<=d;l++) {
			mpfi_init(s[n][l]);
			mpfi_init(w[n][l]);
		}
	}

/* compute sample points: x,y,xmin,ymin,xmax,ymax */
void samples(void) {
	mpfi_t th,h,emh;
	int i,j;

	mpfi_init(th); mpfi_init(h); mpfi_init(emh);
	mpfi_div_ui(h,twopi,24*N);
	mpfi_neg(emh,h);
	mpfi_put(h,emh);
	mpfi_exp(emh,h);

	for (j=0;j<=N;j++) {
		/* theta = Pi/3*(1+j/2N) */
		mpfi_mul_ui(th,twopi,2*N+j);
		mpfi_div_ui(th,th,12*N);
		for (i=0;i<4;i++)
			mpfi_init(sample[i][j]);
		mpfi_cos(sample[0][j],th);
		mpfi_sin(sample[1][j],th);

		/* over a neighborhood, theta varies by +/- h
		   and the radius between exp(-h) and exp(h) */
		mpfi_add(th,th,h);
		mpfi_cos(sample[2][j],th);
		mpfi_mul(sample[2][j],sample[2][j],emh);
		mpfi_sin(sample[3][j],th);
		mpfi_mul(sample[3][j],sample[3][j],emh);
	}
	mpfi_clear(th); mpfi_clear(h); mpfi_clear(emh);
}

/* z^{1/2}K(z) and its first derivative */
void mybessel_fr(mpfi_t result1,mpfi_t result2,mpfr_t zin) {
	mp_prec_t pr;
	mp_exp_t realexp,imagexp;
	mpfi_t freal,fimag,temp,temp2,z,z2,s1,s2;
	mpfr_t temp_fr;
	double zd;
	int k,kmin,thresh;

	/* add extra precision to the computation b/c of cancellation--
	   the terms in the series get as large as about exp(z), the answer is
	   of size exp(-z), and we want to be accurate to the precision of z */
	pr = mpfr_get_prec(zin);
	zd = mpfr_get_d(zin,GMP_RNDU);
	thresh = -(pr+(int)ceil(zd*M_LOG2E)+32);
	pr += (int)ceil(zd*2*M_LOG2E)+32;

	/* abort if we need more precision than we have for R,lambda,etc. */
	if (pr > const_prec) {
		fprintf(stderr,"constant precision too low\n");
		exit(1);
	}
	kmin = (int)ceil(zd);

	/* initialize variables */
	mpfi_init2(freal,pr); mpfi_init2(fimag,pr);
	mpfi_init2(z,pr);     mpfi_init2(z2,pr);
	mpfi_init2(s1,pr);    mpfi_init2(s2,pr);
	mpfi_init2(temp,pr);  mpfi_init2(temp2,pr);
	mpfr_init2(temp_fr,pr);

	/* upgrade the precision of zin */
	mpfr_set(temp_fr,zin,GMP_RNDN);
	mpfi_set_fr(z,temp_fr);
	mpfi_div_ui(z,z,2);
	mpfi_sqr(z2,z);

	/* set f = C*z^(1/2-iR), where
	   C = -i*sqrt(2)*Pi/(sinh(Pi*R)*gamma(1-iR)) is precomputed
		 and given in polar coordinates by bessel_const_{mag,phase} */
	mpfi_log(temp,z);
	mpfi_mul(temp,temp,R);
	mpfi_sub(temp,bessel_const_phase,temp);
	mpfi_cos(freal,temp);
	mpfi_sin(fimag,temp);
	mpfi_sqrt(temp,z);
	mpfi_mul(temp,temp,bessel_const_mag);
	mpfi_mul(freal,freal,temp);
	mpfi_mul(fimag,fimag,temp);

	/* set initial terms: s1 = real(f), s2 = real((1/2-iR)f); */
	mpfi_set(s1,freal);
	mpfi_div_ui(s2,freal,2);
	mpfi_mul(temp,fimag,R);
	mpfi_add(s2,s2,temp);

	k = 1;
	do {
		/* f *= (1+iR/k); */
		mpfi_div_ui(temp,R,k);
		mpfi_mul(temp2,temp,fimag);
		mpfi_mul(temp,temp,freal);
		mpfi_sub(freal,freal,temp2);
		mpfi_add(fimag,fimag,temp);

		/* f *= z^2/(r^2+k^2); */
		mpfi_add_ui(temp,R2,k*k);
		mpfi_div(temp,z2,temp);
		mpfi_mul(freal,freal,temp);
		mpfi_mul(fimag,fimag,temp);

		/* s1 += real(f), s2 += real((2k+1/2-iR)f); */
		/* note: doubles are treated as exact!      */
		mpfi_add(s1,s1,freal);
		mpfi_mul_d(temp,freal,2.0*k+0.5);
		mpfi_add(s2,s2,temp);
		mpfi_mul(temp,fimag,R);
		mpfi_add(s2,s2,temp);

		mpfi_mag(temp_fr,freal);
		realexp = mpfr_get_exp(temp_fr);
		mpfi_mag(temp_fr,fimag);
		imagexp = mpfr_get_exp(temp_fr);
		k++;
	} while (k <= kmin || realexp > thresh || imagexp > thresh);

	/* temp2 = |f|*k^2/(k^2-z^2) */
	mpfi_sqr(temp,freal);
	mpfi_sqr(temp2,fimag);
	mpfi_add(temp2,temp,temp2);
	mpfi_sqrt(temp2,temp2);
	mpfi_mul_ui(temp2,temp2,k*k);
	mpfi_ui_sub(temp,k*k,z2);
	mpfi_div(temp2,temp2,temp);

	/* W series tail is bounded by |f|*k^2/(k^2-z^2) */
	mpfi_neg(temp,temp2);
	mpfi_put(temp,temp2);
	mpfi_add(result1,s1,temp);

	/* W' series tail is bounded by |f|*k^2/(k^2-z^2)*(2k+sqrt(lambda)) */
	mpfi_add_ui(temp,sqrt_lambda,2*k);
	mpfi_mul(temp2,temp2,temp);
	mpfi_neg(temp,temp2);
	mpfi_put(temp,temp2);
	mpfi_add(temp,s2,temp);
	mpfi_div_ui(temp,temp,2);
	mpfi_div(result2,temp,z);

	/* free variables */
	mpfi_clear(freal); mpfi_clear(fimag);
	mpfi_clear(z);     mpfi_clear(z2);
	mpfi_clear(s1);    mpfi_clear(s2);
	mpfi_clear(temp);  mpfi_clear(temp2);
	mpfr_clear(temp_fr);
}

/* interval Bessel function */
void mybessel(mpfi_t rop1,mpfi_t rop2,mpfi_t zin) {
	mpfi_t l0,l1,r0,r1;
	int i;

	mpfi_init(l0); mpfi_init(l1);
	mpfi_init(r0); mpfi_init(r1);
	mybessel_fr(l0,l1,&zin->left);
	mybessel_fr(r0,r1,&zin->right);
	mpfi_put(r0,l0); mpfi_put(r1,l1);
	/* check for extrema inside the interval */
	for (i=0;i<ecount0;i++)
		if (!mpfi_cmp(zin,extrema0[i].u))
			mpfi_put(r0,extrema0[i].w);
	for (i=0;i<ecount1;i++)
		if (!mpfi_cmp(zin,extrema1[i].u))
			mpfi_put(r1,extrema1[i].w);
	mpfi_set(rop1,r0); mpfi_set(rop2,r1);
	mpfi_clear(l0); mpfi_clear(l1);
	mpfi_clear(r0); mpfi_clear(r1);
}

/* compute Bessel function extrema */
void bessel_extrema(void) {
	int i,pr;
	mpfr_t u,umin,step,u1,u2,usamp[bessel_points+1];
	mpfi_t w0,w1,wsamp0[bessel_points+1],wsamp1[bessel_points+1];
	mpfi_t temp;

	mpfr_init(u);  mpfr_init(umin);
	mpfr_init(u1); mpfr_init(u2);
	mpfi_init(w0); mpfi_init(w1);
	mpfr_init(step); mpfi_init(temp);

	/* lowest value we compute is 2*Pi*(ymin of first sample point) */
	mpfr_mul(umin,&twopi->left,&sample[3][0]->left,GMP_RNDD);

	/* oscillations stop at sqrt(lambda) */
	mpfr_sub(step,&sqrt_lambda->right,umin,GMP_RNDU);
	mpfr_div_ui(step,step,bessel_points,GMP_RNDU);

	/* compute initial samples */
	for (i=0;i<=bessel_points;i++) {
		mpfr_mul_ui(u,step,i,GMP_RNDN);
		mpfr_add(u,u,umin,GMP_RNDN);
		mpfr_init(usamp[i]);
		mpfr_set(usamp[i],u,GMP_RNDN);
		mpfi_init(wsamp0[i]); mpfi_init(wsamp1[i]);
		mybessel_fr(wsamp0[i],wsamp1[i],usamp[i]);
		/* This is pedantic, but we need to be sure that the sample points
			 separate the zeros and don't straddle them */
		if (mpfi_is_inside_ui(0,wsamp0[i]) || mpfi_is_inside_ui(0,wsamp1[i])) {
			fprintf(stderr,"Bessel function sample contains 0!!\n");
			exit(1);
		}
	}

	ecount0 = 0, ecount1 = 0;
	for (i=0;i<bessel_points;i++) {
		/* check for extremum of W (zero of W') */
		mpfi_union(w1,wsamp1[i],wsamp1[i+1]);
		if (mpfi_is_inside_ui(0,w1)) {
			/* We know the computed sample points don't contain 0,
			   so we can be sure that W'(u1) >= 0, W'(u2) <= 0 */
			if (mpfi_is_pos(wsamp1[i])) {
				mpfr_set(u1,usamp[i],GMP_RNDN);
				mpfr_set(u2,usamp[i+1],GMP_RNDN);
			} else {
				mpfr_set(u1,usamp[i+1],GMP_RNDN);
				mpfr_set(u2,usamp[i],GMP_RNDN);
			}
			pr = mpfr_get_prec(u);
			do {
				mpfr_add(u,u1,u2,GMP_RNDN);
				mpfr_div_ui(u,u,2,GMP_RNDN);
				mybessel_fr(w0,w1,u);
				if (mpfi_is_pos(w1)) mpfr_set(u1,u,GMP_RNDN);
				if (mpfi_is_neg(w1)) mpfr_set(u2,u,GMP_RNDN);
			} while (--pr && !mpfi_is_inside_ui(0,w1));

			mpfi_init(extrema0[ecount0].u);
			mpfi_init(extrema0[ecount0].w);

			/* we know the zero is contained between u1 and u2 */
			mpfi_interv_fr(extrema0[ecount0].u,u1,u2);

			/* since the function is concave or convex here, we can
			   bound the extremum by a linear approximation */
			mpfi_sub_fr(temp,extrema0[ecount0].u,u);
			mpfi_mul(temp,temp,w1);
			mpfi_add(extrema0[ecount0].w,temp,w0);
			ecount0++;
		}

		/* check for extremum of W' (zero of W) */
		mpfi_union(w0,wsamp0[i],wsamp0[i+1]);
		if (mpfi_is_inside_ui(0,w0)) {
			/* We know the computed sample points don't contain 0,
			   so we can be sure that W(u1) >= 0, W(u2) <= 0 */
			if (mpfi_is_pos(wsamp0[i])) {
				mpfr_set(u1,usamp[i],GMP_RNDN);
				mpfr_set(u2,usamp[i+1],GMP_RNDN);
			} else {
				mpfr_set(u1,usamp[i+1],GMP_RNDN);
				mpfr_set(u2,usamp[i],GMP_RNDN);
			}
			pr = mpfr_get_prec(u);
			do {
				mpfr_add(u,u1,u2,GMP_RNDN);
				mpfr_div_ui(u,u,2,GMP_RNDN);
				mybessel_fr(w0,w1,u);
				if (mpfi_is_pos(w0)) mpfr_set(u1,u,GMP_RNDN);
				if (mpfi_is_neg(w0)) mpfr_set(u2,u,GMP_RNDN);
			} while (--pr && !mpfi_is_inside_ui(0,w0));

			mpfi_init(extrema1[ecount1].u);
			mpfi_init(extrema1[ecount1].w);

			/* we know the zero is contained between u1 and u2 */
			mpfi_interv_fr(extrema1[ecount1].u,u1,u2);

			/* replace w0 by W'' */
			mpfi_set_fr(temp,u);
			mpfi_sqr(temp,temp);
			mpfi_div(temp,lambda,temp);
			mpfi_ui_sub(temp,1,temp);
			mpfi_mul(w0,w0,temp);

			/* since the function is convex here, we can
			   bound the extremum by a linear approximation */
			mpfi_sub_fr(temp,extrema1[ecount1].u,u);
			mpfi_mul(temp,temp,w0);
			mpfi_add(extrema1[ecount1].w,temp,w1);
			ecount1++;
		}
	}

	/* add one more extremum of W' at sqrt(lambda) */
	mpfi_init(extrema1[ecount1].u);
	mpfi_init(extrema1[ecount1].w);
	mpfi_set(extrema1[ecount1].u,sqrt_lambda);
	mpfi_get_fr(u,extrema1[ecount1].u);
	mybessel_fr(w0,w1,u);

	/* replace w0 by W'' */
	mpfi_set_fr(temp,u);
	mpfi_sqr(temp,temp);
	mpfi_div(temp,lambda,temp);
	mpfi_ui_sub(temp,1,temp);
	mpfi_mul(w0,w0,temp);

	/* since the function is convex here, we can
	   bound the extremum by a linear approximation */
	mpfi_sub_fr(temp,extrema1[ecount1].u,u);
	mpfi_mul(temp,temp,w0);
	mpfi_add(extrema1[ecount1].w,temp,w1);
	ecount1++;

	mpfr_clear(u);  mpfr_clear(umin);
	mpfr_clear(u1); mpfr_clear(u2);
	mpfi_clear(w0); mpfi_clear(w1);
	mpfr_clear(step); mpfi_clear(temp);
	for (i=0;i<=bessel_points;i++) {
		mpfr_clear(usamp[i]);
		mpfi_clear(wsamp0[i]);
		mpfi_clear(wsamp1[i]);
	}
}

/* derivatives 0 through d of z^{1/2}*K(z) at z=2*Pi*n*y for n=1..M */
/* return results in global variable w */
void bessel_derivatives(mpfi_t y) {
	int j,l,n;
	mpfi_t twopiy,zinv,zp,temp;

	mpfi_init(twopiy); mpfi_init(temp);
	mpfi_init(zinv);   mpfi_init(zp);
	mpfi_mul(twopiy,twopi,y);

	for (n=1;n<=M;n++) {
		/* compute first two derivatives */
		mpfi_mul_ui(temp,twopiy,n);
		mybessel(w[n][0],w[n][1],temp);

		/* compute higher derivatives by recurrence */
		mpfi_ui_div(zinv,1,temp);
		for (l=2;l<=d;l++) {
			mpfi_sqr(zp,zinv);
			mpfi_mul(zp,zp,lambda);
			mpfi_set(w[n][l],w[n][l-2]);
			for (j=l-2;j>=0;j--) {
				mpfi_mul(temp,zp,w[n][j]);
				mpfi_mul_z(temp,temp,c[l][j]);
				mpfi_add(w[n][l],w[n][l],temp);
				mpfi_mul(zp,zp,zinv);
			}
		}
	}
	mpfi_clear(twopiy); mpfi_clear(temp);
	mpfi_clear(zinv);   mpfi_clear(zp);
}

/* derivatives 0 through d of cos^{(epsilon)}(z) at z=2*Pi*n*x for n=1..M */
/* return results in global variable s */
void sincos_derivatives(mpfi_t x) {
	mpfi_t twopix,z,cs[4];
	int k,n;

	mpfi_init(twopix); mpfi_init(z);
	mpfi_mul(twopix,twopi,x);
	for (k=0;k<4;k++) mpfi_init(cs[k]);
	for (n=1;n<=M;n++) {
		mpfi_mul_ui(z,twopix,n);
		mpfi_cos(cs[0],z);     mpfi_sin(cs[3],z);
		mpfi_neg(cs[2],cs[0]); mpfi_neg(cs[1],cs[3]);
		for (k=0;k<=d;k++)
			mpfi_set(s[n][k],cs[(k+parity) & 3]);
	}
	for (k=0;k<4;k++) mpfi_clear(cs[k]);
	mpfi_clear(twopix); mpfi_clear(z);
}

/* compute powers of x up to the nth */
void powers(mpfi_t p[],mpfi_t x,int n) {
	int i;

	mpfi_set_ui(p[0],1);
	mpfi_set(p[1],x);
	for (i=2;i<=n;i++)
		mpfi_mul(p[i],p[i-1],x);
}

/* Compute derivatives in rectangular coordinates:
   for each sample point j and each k and l, give
   the value of the (k,l) derivative at point j.
	 Results returned in f, which is initialized as we go.
*/
void rectangular(mpfi_t f[][dmax][dmax],mpfi_t x[],mpfi_t y[]) {
	int j,k,l,n;
	mpfi_t temp,power[dmax];

	mpfi_init(temp);
	for (k=0;k<=d;k++)
		mpfi_init(power[k]);

	for (j=0;j<=N;j++) {
		sincos_derivatives(x[j]);
		bessel_derivatives(y[j]);
		for (k=0;k<=d;k++)
			for (l=0;l<=d-k;l++)
				mpfi_init_set_ui(f[j][k][l],0);
		for (n=1;n<=M;n++) {
			mpfi_mul_ui(temp,twopi,n);
			powers(power,temp,d);
			for (k=0;k<=d;k++)
				for (l=0;l<=d-k;l++) {
					mpfi_mul(temp,A[n],power[k+l]);
					mpfi_mul(temp,temp,w[n][l]);
					mpfi_mul(temp,temp,s[n][k]);
					mpfi_add(f[j][k][l],f[j][k][l],temp);
				}
		}
	}
	mpfi_clear(temp);
	for (k=0;k<=d;k++)
		mpfi_clear(power[k]);
}

void polar(mpfi_t result) {
	mpfi_t rssum[Nmax],klsum,hf[dmax],term,temp;
	mpz_t term_z,temp_z;
	int j,k,l,m,r,s;

	mpfi_init(term);  mpfi_init(temp);
	mpz_init(term_z); mpz_init(temp_z);
	mpfi_div_ui(temp,twopi,24*N);
	mpfi_neg(term,temp);
	mpfi_put(temp,term);

	mpfi_init_set_ui(hf[0],1);
	for (k=1;k<=d;k++) {
		mpfi_init(hf[k]);
		mpfi_mul(hf[k],hf[k-1],temp);
		mpfi_div_ui(hf[k],hf[k],k);
	}
	for (j=0;j<=N;j++) {
		for (k=0;k<=d;k++) {
			mpfi_init(xp[j][k]);
			mpfi_init(yp[j][k]);
			mpfi_init(nxp[j][k]);
			mpfi_init(nyp[j][k]);
		}
		powers(xp[j],sample[0][j],d);
		powers(yp[j],sample[1][j],d);
		powers(nxp[j],sample[2][j],d);
		powers(nyp[j],sample[3][j],d);
	}

	for (k=0;k<=d;k++)
		for (l=0;l<=d-k;l++)
			for (m=0;m<=k+l;m++) {
				mpz_init(P[k][l][m]);
				mpz_init(Q[k][l][m]);
				mpz_init(Qnext[k][l][m]);
			}
	mpz_set_ui(P[0][0][0],1);
	mpz_set_ui(Q[0][0][0],1);

	mpfi_init(klsum);
	for (j=0;j<=N;j++)
		mpfi_init_set_ui(rssum[j],0);
	for (r=0;r<=d;r++) {
		if (r > 0) {
			/* diff wrt t */
			for (k=0;k<=r;k++)
			for (l=0;l<=r-k;l++)
			for (m=0;m<=k+l;m++) {
				mpz_set_ui(term_z,0);
				if (k+l < r) {
					mpz_mul_ui(temp_z,P[k][l][m],k+l);
					mpz_add(term_z,term_z,temp_z);
				}
				if (k > 0 && m > 0)
					mpz_add(term_z,term_z,P[k-1][l][m-1]);
				if (l > 0 && m < k+l)
					mpz_add(term_z,term_z,P[k][l-1][m]);
				mpz_set(Q[k][l][m],term_z);
			}
			for (k=0;k<=r;k++)
			for (l=0;l<=r-k;l++)
			for (m=0;m<=k+l;m++)
				mpz_set(P[k][l][m],Q[k][l][m]);
		}
		for (s=0;s<=d-r;s++) {
			if (s > 0) {
				/* diff wrt theta */
				for (k=0;k<=r+s;k++)
				for (l=0;l<=r+s-k;l++)
				for (m=0;m<=k+l;m++) {
					mpz_set_ui(term_z,0);
					if (m > 0) {
						if (k+l < r+s) {
							mpz_mul_ui(temp_z,Q[k][l][m-1],k+l-m+1);
							mpz_add(term_z,term_z,temp_z);
						}
						if (l > 0)
							mpz_add(term_z,term_z,Q[k][l-1][m-1]);
					}
					if (m < k+l) {
						if (k+l < r+s) {
							mpz_mul_ui(temp_z,Q[k][l][m+1],m+1);
							mpz_sub(term_z,term_z,temp_z);
						}
						if (k > 0)
							mpz_sub(term_z,term_z,Q[k-1][l][m]);
					}
					mpz_set(Qnext[k][l][m],term_z);
				}
				for (k=0;k<=r+s;k++)
				for (l=0;l<=r+s-k;l++)
				for (m=0;m<=k+l;m++)
					mpz_set(Q[k][l][m],Qnext[k][l][m]);
			}
			if (r+s == d)
				for (j=0;j<=N;j++) {
					mpfi_set_ui(klsum,0);
					for (k=0;k<=r+s;k++)
					for (l=0;l<=r+s-k;l++) {
						mpfi_set_ui(term,0);
						for (m=0;m<=k+l;m++) {
							mpfi_mul(temp,nxp[j][m],nyp[j][k+l-m]);
							mpfi_mul_z(temp,temp,Q[k][l][m]);
							mpfi_add(term,term,temp);
						}
						mpfi_mul(term,term,fderiv_nbhd[j][k][l]);
						mpfi_add(klsum,klsum,term);
					}
					mpfi_mul(temp,hf[r],hf[s]);
					mpfi_mul(temp,temp,klsum);
					mpfi_add(rssum[j],rssum[j],temp);
				}
			else if ((r+parity) & 1)
				for (j=0;j<=N;j++) {
					mpfi_set_ui(klsum,0);
					for (k=0;k<=r+s;k++)
					for (l=0;l<=r+s-k;l++) {
						mpfi_set_ui(term,0);
						for (m=0;m<=k+l;m++) {
							mpfi_mul(temp,xp[j][m],yp[j][k+l-m]);
							mpfi_mul_z(temp,temp,Q[k][l][m]);
							mpfi_add(term,term,temp);
						}
						mpfi_mul(term,term,fderiv[j][k][l]);
						mpfi_add(klsum,klsum,term);
					}
					mpfi_mul(temp,hf[r],hf[s]);
					mpfi_mul(temp,temp,klsum);
					mpfi_add(rssum[j],rssum[j],temp);
				}
		}
	}

	mpfi_set_ui(result,0);
	for (j=0;j<=N;j++) {
		mpfi_abs(temp,rssum[j]);
		mpfi_put(result,temp);
	}
	mpfi_mul_ui(result,result,2);  /* the factor of 2 in equation (49) */
	mpfi_clear(term);  mpfi_clear(temp);
	mpz_clear(term_z); mpz_clear(temp_z);
	for (k=0;k<=d;k++) mpfi_clear(hf[k]);
	mpfi_clear(klsum);
	for (j=0;j<=N;j++) mpfi_clear(rssum[j]);
}

/* Lower Riemann sum for the integral from a to b of W(2*Pi*y)^2/y^2 */
void riemann_sum(mpfi_t answer,double a,double b) {
#define riemann_steps 256
	mpfi_t u,sum,w0,w1,temp;
	double h,u1,u2;
	int i;

	mpfi_init2(u,53);   mpfi_init2(temp,53);
	mpfi_init2(w0,53);  mpfi_init2(w1,53);
	mpfi_init2(sum,53);	mpfi_set_ui(sum,0);

	/* really want an integral out to infinity, so ok if not exact here */
	u1 = a, h = (b-a)/riemann_steps;

	for (i=1;i<=riemann_steps;i++) {
		u2 = a+i*h;
		mpfi_interv_d(u,u1,u2);

		mpfi_mul(temp,u,twopi);
		mybessel(w0,w1,temp);
		mpfi_div(w0,w0,u);
		mpfi_sqr(w0,w0);

		mpfi_set_fr(temp,&u->right);
		mpfi_sub_fr(temp,temp,&u->left);
		mpfi_mul(w0,w0,temp);

		mpfi_add(sum,sum,w0);
		u1 = u2;
	}
	mpfi_set(answer,sum);
	mpfi_clear(u);   mpfi_clear(temp);
	mpfi_clear(w0);  mpfi_clear(w1);
	mpfi_clear(sum);
}

int main(int argc,char *argv[]) {
	mpfi_t bound,sum,temp;
	FILE *fp=stdin;
	int pr=75;
	d=10, M=10, N=50;

	if (argc > 1) fp=fopen(argv[1],"r");
	if (argc > 2) sscanf(argv[2],"%d",&d);
	if (argc > 3) sscanf(argv[3],"%d",&M);
	if (argc > 4) sscanf(argv[4],"%d",&N);
	if (argc > 5) sscanf(argv[5],"%d",&pr);

	mpfr_set_default_prec(pr);
	mpfi_init(bound); mpfi_init(sum); mpfi_init(temp);

	initialize(fp);
	samples();
	bessel_extrema();
	rectangular(fderiv,sample[0],sample[1]);
	rectangular(fderiv_nbhd,sample[2],sample[3]);
	polar(bound);

	/* Below we use the fact that the hyperbolic delta-nbhd fits within the
	Euclidean nbhd of size 1-exp(-delta), even though in the paper we just
	say delta.  In other words, it's fine to put simply delta = Pi/12N */

	/* temp = delta^(-1) */
	mpfi_ui_div(temp,24*N,twopi);

	/* bound *= 40*delta^(-3/2) */
	mpfi_mul(bound,bound,temp);
	mpfi_sqrt(temp,temp);
	mpfi_mul(bound,bound,temp);
	mpfi_mul_ui(bound,bound,40);

	/* multiply by C_{f,2} */
	if (!parity) {
		/* temp = 2*cos(R*log(2)) */
		mpfi_set_ui(temp,2);
		mpfi_log(temp,temp);
		mpfi_mul(temp,temp,R);
		mpfi_cos(temp,temp);
		mpfi_mul_ui(temp,temp,2);

		/* bound *= 2*(sqrt(2)+1/sqrt(2)) */
		mpfi_set_ui(sum,2);
		mpfi_sqrt(sum,sum);
		mpfi_mul(bound,bound,sum);
		mpfi_mul_ui(bound,bound,3);

		/* sum = T_2 eigenvalue = sqrt(2)*A[2] */
		mpfi_mul(sum,sum,A[2]);

		mpfi_sub(temp,temp,sum);
		mpfi_abs(temp,temp);
		mpfi_div(bound,bound,temp);
	}

	mpfi_div_ui(temp,twopi,24*N);
	mpfi_exp(temp,temp);
	if (!parity)
		mpfi_mul_ui(temp,temp,2);

	riemann_sum(sum,mpfr_get_d(&temp->right,GMP_RNDU),5.0);
	mpfi_sqrt(sum,sum);
	mpfi_div(bound,bound,sum);

	printf("%14g\n",mpfr_get_d(&bound->right,GMP_RNDU));
	return 0;
}