getmaninsymbol(coefftable,N,t,u) = 
{
	local(temp,temp2,i,r,s,M1,t1,M,h,u1,s1,Mh,U,SUM,sumlength,q1,q2,gammaM,gammaN,gammaMMh,T1,T2,n);

	/* set up parameters in notation of Goldfeld */
	temp = bezout(t,u); r = temp[1]; s = -temp[2];
	M1 = gcd(t,N); t1 = t/M1; M = N/M1;
	temp = bezout(t,M); h = temp[1]*u/temp[3]; u1 = temp[2]*u/temp[3];
	s1 = s - r*h;
	Mh = gcd(M,h);
	if( M/Mh == 1 , l = 1 , l = M/(h*N) % (M/Mh) );
	U = sqrt(M*Mh);

	/* compute the gammas */
	temp = factor(M); temp2 = matsize(temp);
	gammaM = prod(i=1,temp2[1],(-coefftable[temp[i,1]])^temp[i,2]);
	temp = factor(N); temp2 = matsize(temp);
	gammaN = prod(i=1,temp2[1],(-coefftable[temp[i,1]])^temp[i,2]);
	temp = factor(M/Mh); temp2 = matsize(temp);
	gammaMMh = prod(i=1,temp2[1],(-coefftable[temp[i,1]])^temp[i,2]);

	/* compute the sum */
	SUM = 0; sumlength = length(coefftable);
	q1 = exp( -(2*Pi/M) * (U/sqrt(N) - I*h) );
	q2 = exp( -(2*Pi*Mh) * ( 1/(U*sqrt(N)) + I*l/M ) );
	T1 = gammaM; T2 = -gammaN*gammaM*gammaMMh;
	for(n=1,sumlength,
		T1 *= q1;
		T2 *= q2;
		SUM += (coefftable[n]/n)*(T1+T2);
	);
	
	return(SUM);
}

ellmaninraw0(E) =
{
	local(N,sumlength,an,i,j,ik,jk,k,temp,gcdij,outmat,lcmdenom1,lcmdenom2);

	N = ellglobalred(E)[1];
	sumlength = 10*N;
	an = ellan(E,sumlength);
	outmat = matrix(N,N);
	lcmdenom1 = 1;
	lcmdenom2 = 1;

	/* fill in table based on omega+ and omega- determined by elliptic curve */
	for(i=1,N, for(j=1,N, 
		gcdij = gcd(i,j);
		if(gcd(gcdij,N) == 1,
		   if(outmat[i,j] == 0,
		      temp = lindep([E.omega[1],imag(E.omega[2])*I,getmaninsymbol(an,N,i/gcdij % N,j/gcdij % N)]);
		      outmat[i,j] = [-temp[1]/temp[3],-temp[2]/temp[3]];
		      lcmdenom1 = lcm(lcmdenom1,denominator(outmat[i,j][1]));
		      lcmdenom2 = lcm(lcmdenom2,denominator(outmat[i,j][2]));
		      for(k=1,N, if(gcd(k,N) == 1,
		         ik = i*k % N; if(ik == 0,ik = N);
		         jk = j*k % N; if(jk == 0,jk = N);
		         outmat[ik,jk] = outmat[i,j];
		      ));
		    );
		,
		   outmat[i,j] = [1];
		);
	));

	/* renormalize the table for omega+ and omega- determined by modular form */
	print([lcmdenom1,lcmdenom2]);
	for(i=1,N, for(j=1,N,
		if(length(outmat[i,j]) == 2, outmat[i,j] = [outmat[i,j][1]*lcmdenom1,outmat[i,j][2]*lcmdenom2]);
	));

	return(outmat);
}

ellmanincomb(matr,a,b) =
{
	local(i,j);
	for(i=1,matsize(matr)[1], for(j=1,matsize(matr)[2],
		if(length(matr[i,j]) == 2, matr[i,j] = [a*matr[i,j][1] + b*matr[i,j][2],0]);
	));

	return(matr);
}


ellcuspint0(table,a,b) =
{
	return(ellcuspint1(table,b) - ellcuspint1(table,a));
}

ellcuspint1(table, x) =
{
	local(N,aa,kappa,pp,qq,k,msym,SUM);

	N = matsize(table)[1];

	if(x == [1],
	   SUM = table[N,1];
	, if(denominator(x) == 1,
	   SUM = [0,0];
	,

	   /* compute Farey sequence from 0 to x */
	   /* aa[k+1] = a[k], pp[k+3] = p[k], qq[k+3] = q[k] */
	   aa = contfrac(x);
	   kappa = length(aa) - 1;
	   qq = vector(kappa+3);
	   qq[1] = 1; qq[2] = 0;
	   for(k=0,kappa,
		qq[k+3] = aa[k+1]*qq[k+2] + qq[k+1];
	   );
	
	   /* add up corresponding Manin symbols */
	   SUM = [0,0];
	   for(k=1,kappa,
		msym = [(-1)^(k-1)*qq[k+3] % N, qq[k+2] % N];
		if(msym[1] == 0, msym[1] = N);
		if(msym[2] == 0, msym[2] = N);
		SUM += table[msym[1],msym[2]];
	   );

	));

	return(SUM);
}

ellmanindbadd0(cremonaentry) =
{
	local(rawtable,outfilename);
	cload(cremona); /* WARNING! */
	rawtable = ellmaninraw(ellini(cremonaentry));
	outfilename = concat([manindbfolder,cremonaentry,".raw"]);
	write(outfilename,rawtable);
}

ellmanindbget0(cremonaentry) =
{
	local(infilename);
	infilename = concat([manindbfolder,cremonaentry,".raw"]);
	read(infilename);
}