hpxhint0(tau1,tau2,z1,z2,table,prec) =
{
	local(balls,s,j,k,PROD,numballs,p,temp,tnumer,tdenom,realnumer,imagnumer,denomin,val,t,radius,mu,musum,a1,a2,b1,b2,gammaz1,gammaz2);
	
	/* collect info on tau1 and tau2 */
	if(type(tau1) != "t_QUAD" || type(tau2) != "t_QUAD", error("tau1 and tau2 must be in H_p"));
	a1 = real(tau1); b1 = imag(tau1);
	a2 = real(tau2); b2 = imag(tau2);
	if(type(b1) != "t_PADIC" || type(b2) != "t_PADIC", error("tau1 and tau2 must be in H_p with explicit precision"));
	p = component(b1,1);
	if(component(b2,1) != p, error("tau1 and tau2 must be in H_p for a fixed p"));
	s = cpnonsquare(p);

	/* prepare a list of balls */
	balls = vector(p^prec+p^(prec-1));
	k = 1;
	for(j=0,p^prec-1,
		balls[k] = [j,prec]; k++;
	);
	balls[k] = [1/p^(prec+2),prec]; k++;
	for(j=1,p^(prec-1)-1,
		balls[k] = [1/(p*j),prec]; k++;
	);

	/* main loop */
	PROD = 1;
	musum = 0;
	numballs = length(balls);
	for(j=1,numballs,

		t = balls[j][1]; radius = balls[j][2];

		/* evaluate mu{z1,z2}(U_j) */
		if(valuation(t,p) >= 0,
			gammaz1 = (z1 - t)/p^radius;
			gammaz2 = (z2 - t)/p^radius;
		,
			if(z1==0, gammaz1 = [1], gammaz1 = (-1/z1 + 1/t)/p^radius);
			if(z2==0, gammaz2 = [1], gammaz2 = (-1/z2 + 1/t)/p^radius);
		);
		/* absgamma = (-1)^(radius % 2); */
		absgamma = 1;
		temp = absgamma*ellcuspint(table,gammaz1,gammaz2);
		mu = temp[1] + temp[2];

		musum += mu;
		/*
		print("["j" t="t", gz1="gammaz1", gz2 = "gammaz2", radius="radius", mu="mu", musum="musum"]");
		*/

		if(mu == 0, next());

		/* evaluate (t_j - tau1)/(t_j - tau2) or inverse according to sign of mu */
		realnumer = b1*b2*s - t^2 + t*a2 + a1*t - a1*a2;
		if(mu>0,
			imagnumer = -t*b2 + a1*b2 + b1*t - b1*a2;
			denomin = -t^2 + 2*t*a2 - a2^2 + b2^2*s;
		,
			imagnumer = -t*b1 + a2*b1 + b2*t - b2*a1;
			denomin = -t^2 + 2*t*a1 - a1^2 + b1^2*s;
			mu = -mu;
		);
		val = cptoquad([realnumer/denomin,imagnumer/denomin,s]);


		PROD *= val^mu;

	);

	print("musum: " musum);
	return(PROD);
}

hpxhgammatau(disc) =
{
        local(w,u,gammatau);
        w = quadgen(disc);
        u = quadunit(disc);
        if(norm(u) == -1, u = u^2);
        gammatau = real(u)*[1,0;0,1] + imag(u)*[trace(w),-norm(w);1,0];
        return(gammatau);
}

hpxhlfx(gamm,tau) =
{
	local(numer,denom);
	numer = cpadd(cpmul(gamm[1,1],tau),gamm[1,2]);
	denom = cpadd(cpmul(gamm[2,1],tau),gamm[2,2]);
	return(cpdiv(numer,denom));
}

hpxhmodsymb8(tau,fareyseq,manintable,prec) =
{
	local(lgth,eta01,j,a,b,c,d,xformj,INTEGRAL,tau2);
	lgth = length(fareyseq) - 1;
	eta01 = hpxhint(hpxhlfx([1,-1;1,0],tau),hpxhlfx([1,0;1,1],tau),0,1,manintable,prec);
	PROD = 1;
	for(j=1,lgth,
		a = numerator(fareyseq[j]); b = denominator(fareyseq[j]);
		c = numerator(fareyseq[j+1]); d = denominator(fareyseq[j+1]);
		if(a*d-b*c == -1,
			c = -c;
			d = -d;
		);
		if(a*d - b*c != 1, error("ad-bc != 1"));
		xformj = [-b,a;-(d+b),c+a];
		tau2 = hpxhlfx(xformj,tau);
		if(tau == tau2,
			INTEGRAL = 1;
		,
			INTEGRAL = hpxhint(tau,tau2,0,1,manintable,prec);
		);
		PROD = cpmul(cpmul(INTEGRAL,PROD),eta01);
	);
	return(PROD);
}

hpxhmodsymb9(tau,fareyseq,manintable,prec) =
{
	local(lgth,eta01,j,a,b,c,d,xformj,INTEGRAL,tau2);
	lgth = length(fareyseq) - 1;
	eta01 = hpxhint9(hpxhlfx([1,-1;1,0],tau),hpxhlfx([1,0;1,1],tau),0,1,manintable,prec);
	PROD = 1;
	for(j=1,lgth,
	        a = numerator(fareyseq[j]); b = denominator(fareyseq[j]);
                c = numerator(fareyseq[j+1]); d = denominator(fareyseq[j+1]);
		if(a*d-b*c == -1,
			c = -c;
                        d = -d;
                );
                if(a*d - b*c != 1, error("ad-bc != 1"));
                xformj = [-b,a;-(d+b),c+a];
                tau2 = hpxhlfx(xformj,tau);
                if(tau == tau2,
                        INTEGRAL = 1;
                ,
                        INTEGRAL = hpxhint9(tau,tau2,0,1,manintable,prec);
                );
                PROD = cpmul(cpmul(INTEGRAL,PROD),eta01);
        );
        return(PROD);
}


hpxhappr(x,N=0)=
 {  local(prec,res,denom, num);
    prec=component(x,2); res=[];
   if(N==0, N=sqrtint(prec),);
   for(denom=1,N,
     num = (denom*x) % prec;
     if(gcd(num,denom)==1,
     if(num<=N,res=concat(res,num/denom),),);
     num = (-denom*x) % prec;
     if(gcd(num,denom)==1,
     if(num<=N,res=concat(res,-num/denom),),)  );
    return(res)  }

hpxhgammaF(F) =
{
	local(A,B,C,delta,u,u1,u2,M);
	A = component(F,1);
	B = component(F,2);
	C = component(F,3);
	delta = B^2 - 4*A*C;
	u = quadunit(delta);
	if(norm(u) == -1, u = u^2);
	u1 = real(u); u2 = imag(u);
	if(delta % 4 == 0,
		M = (u1 + u2*B/2) + u2*[-B,-C;A,0];
	,
		M = (u1 + u2*(B+1)/2) + u2*[-B,-C;A,0];
	);
	return(M);
}

cfcgs(L) =
{
	local(a,p,q,cgs,j);
	a = L;
	p = vector(length(L));
	q = vector(length(L));
	cgs = vector(length(L));
	p[1] = a[1]; p[2] = a[2]*a[1] + 1;
	q[1] = 1; q[2] = a[2];
	cgs = vector(length(L));
	cgs[1] = p[1]/q[1]; cgs[2] = p[2]/q[2];
	for(j=3,length(L),
		p[j] = a[j]*p[j-1] + p[j-2];
		q[j] = a[j]*q[j-1] + q[j-2];
		cgs[j] = p[j]/q[j];
	);
	return(cgs);
}

hpxhtauF(F,Oprec) =
{
	local(A,B,C,coreF,D,f,outs);
	A = component(F,1);
	B = component(F,2);
	C = component(F,3);
	coreF = core(B^2-4*A*C,1);
	D = coreF[1];
	f = coreF[2];
	outs = [-B/(2*A),f/(2*A) + Oprec,D];
	return(outs);
}

hpxhmodsymb17(F) =
{
	local(z1,z2,tauF,gammaF,fareyseq0,fareyseq,fareyintro,lgth,intlist,powlist,tau,tau2,j,k,a,b,c,d,xformj);

	/* get tau and Farey sequence from z1 = 0 to z2 */
	/* SHOULD GET BETTER FAREY SEQUENCE GENERATOR!!! */
	tauF = hpxhtauF(F);
	gammaF = hpxhgammaF(F);
	z1 = 0;
	z2 = hpxhlfx(gammaF,z1);
	/* */
	fareyseq0 = cfcgs(contfrac(z2));
	k = length(fareyseq0);
	while(denominator(fareyseq0[k]) != 1, k--);
	fareyseq = vector(length(fareyseq0) - k + 1, j, fareyseq0[j + k - 1]);
	if(fareyseq[1] == 1 || fareyseq[1] == -1, fareyseq = concat([0],fareyseq));
	if(fareyseq[1] >= 1 || fareyseq[1] <= -1, fareyseq = concat([0,[1]],fareyseq));
	/* while(fareyseq[1] != z1, fareyseq = concat(fareyseq[1] + sign(z1-fareyseq[1]),fareyseq)); */
	/* */
	if(fareyseq[1] != z1 || fareyseq[length(fareyseq)] != z2, error("Fatal error in strategy: hpxhmodsymb17"));
	tau = tauF;

	print(fareyseq);

	/* pass through Farey sequence to get list of integrals to what exponents needed */
	lgth = length(fareyseq) - 1;
	intlist = [[hpxhlfx([1,-1;1,0],tau),hpxhlfx([1,0;1,1],tau),0,1]];
	powlist = [0];
	for(j=1,lgth,
		if(fareyseq[j] == [1], 
			a = 1; b = 0;
		,
		        a = numerator(fareyseq[j]); b = denominator(fareyseq[j]);
		);
		if(fareyseq[j+1] == [1],
			c = 1; d = 0;
		,
	                c = numerator(fareyseq[j+1]); d = denominator(fareyseq[j+1]);
		);
		if(a*d-b*c == -1,
			c = -c;
                        d = -d;
                );
                if(a*d - b*c != 1, error("ad-bc != 1"));
                xformj = [-b,a;-(d+b),c+a];
                tau2 = hpxhlfx(xformj,tau);
                if(tau != tau2,
                        intlist = concat(intlist,[[tau,tau2,0,1]]);
			powlist = concat(powlist,1);
                );
		powlist[1] = powlist[1] + 1;
        );

        return([intlist,powlist]);
}

hpxhmodsymb18(rawmaninid,functional,F,Oprec) =
{
	local(rawtable,mtable,temp,intlist,powlist,numints,j,intouts,currint,tau1,tau2,z1,z2,xprec,outval);

	rawtable = ellmanindbget(rawmaninid);
	mtable = ellmanincomb(rawtable,functional[1],functional[2]);

	temp = hpxhmodsymb17(F);
	intlist = temp[1]; powlist = temp[2];
	numints = length(intlist);
	intouts = vector(numints);

	print("Computing " numints " integrals...");
	for(j=1,numints,
		currint = intlist[j];
		tau1 = currint[1] + [0,Oprec^3,0];
		tau2 = currint[2] + [0,Oprec^3,0];
		z1 = currint[3];
		z2 = currint[4];
		xprec = padicprec(Oprec,component(Oprec,1));
		print("Computing integral " j "...  ");
		/* intouts[j] = hpxhint(tau1,tau2,z1,z2,mtable,xprec); */
                intouts[j] = hpxhint9(tau1,tau2,z1,z2,mtable,xprec);
	);

	outval = 1;
	for(j=1,numints,
		/* print([outval,intouts[j],powlist[j]]);
		if(j=numints,print(intlist)); */
		outval = cpmul(outval,cppow(intouts[j],powlist[j]));
	);
	
	print("Modular symbol computation complete!");

	return(outval);
}

hpxhmodsymb19(rawmaninid,functional,F,Oprec) =
{
	local(modsymbval,A,B,C,D,outfile,outstuff,prevvec,testline);

	modsymbval = hpxhmodsymb18(rawmaninid,functional,F,Oprec);

	A = component(F,1);
	B = component(F,2);
	C = component(F,3);
	D = B^2 - 4*A*C;
	outfile = concat([componentfolder,"db/J/",rawmaninid,"-",D]);
	testline = concat(["if test ! -e ",outfile,"; then echo [] > ",outfile,"; fi"]);
	system(testline);
	prevvec = read(outfile);
	outstuff = concat(["[Qfb(",A,",",B,",",C,"),",functional,",",Oprec]);
	outstuff = concat([outstuff,",",Str(modsymbval)]);
	outstuff = concat(outstuff,"]");
	prevvec = concat(prevvec,[outstuff]);
	system(concat("rm ",outfile));
	write(outfile,prevvec);
}

hpxhshowJ(Jid) =
{
	local(infile,invec,curr,F,V,tau,fnl,prec);
	infile = concat([componentfolder,"db/J/",Jid]);
	invec = read(infile);
	print("Record for J database " Jid);
	print("\tQfb\t\ttau\t\tfnl\t\tprec");
	for(j=1,length(invec),
		curr = eval(invec[j]);
		F = curr[1]; 
		V = [component(F,1),component(F,2),component(F,3)];
		tau = hpxhtauF(F);
		fnl = curr[2];
		prec = curr[3];
		print(j "\t" V "\t" tau "\t" fnl "\t" prec);
	);
}

hpxhgetJ(Jid,j,fl) =
{
	local(infile,invec,curr);
	infile = concat([componentfolder,"db/J/",Jid]);
	invec = read(infile);
	curr = eval(invec[j]);
	if(fl==0,
		print([curr[1],curr[2],curr[3]]);
		return(curr[4]);
	,
		return([curr[1],curr[2],curr[3],curr[4]]);
	);
}

hpxhFtau(tau) = 
{
	local(a,b,c,A,B,C,d,V);
	a = tau[1]; b = tau[2]; c = tau[3];
	if(b == 0, error("tau should be in upper half plane!"));
	if(b < 0, A = -1, A = 1);
	B = -2*A*a;
	C = A*(a^2-b^2*c);
	d = lcm(denominator(B),denominator(C));
	V = d*[A,B,C];
	return(Qfb(V[1],V[2],V[3]));
}


/* load additional routines */
read(concat(componentfolder,"src/hpxhgetpoint.pari"));
read(concat(componentfolder,"src/twisttools.pari"));