level = 37; md(a,b,res)=\ res=a%b;\ if(2*res>abs(b),res=res-abs(b),);\ res symbol(a,b,w,ww,aa,bb,res) =\ until( ((2*abs(a) <= abs(b)*level) && (2*abs(b)<= abs(a))) || (b==0),\ a=md(a,b*level);\ b=md(b,a) );\ if(b<0,b=-b;a=-a,);\ if(b==0,res=0,\ if(b==1,res=table[a%level],\ w=round(-a/b);\ ww = (1/w)% level;\ bb= b*w + a;\ aa= (a* ww + b*(w*ww-1));\ res=symbol(aa,bb)-table[ww] ));\ res \\ Sets up the computation with disc = d. setup(d0,f,d,u)=\ d=d0*f^2;\ gd0=d0;gc=f;\ gs=sqrt(d);\ gsp=factmod(x^2-d,level);\ gsp=compo(gsp,1);\ gsp=compo(gsp,1);\ gsp=compo(gsp,1);\ gsp=compo(gsp,2);\ u=unit(d)^2;\ gu=compo(u,2);gv=compo(u,3);\ if(d%4==1,gu=gu+gv/2; gv=gv/2;gid=try(1,1),\ gv=gv/2;gid=try(0,1));\ print(gu^2-d*gv^2);\ adjust(f,a,b,c,res)=\ a=f[1];b=f[2];c=f[3];\ if (a%level==0,if( (b+gsp) % level ==0,res=f,res=adjust([c,-b,a])),\ k=( ( -b+gsp)/(2*a)) % level;\ res=shift(f,k);\ res=[res[3],-res[2],res[1]]);\ res geodesic(f,ff,a,b)=\ ff=adjust(f);\ a=ff[1];b=ff[2];\ symbol(gu+gv*b,-2*a*gv/level) computetheta(f,n, th,dummy,ff,j)=\ th = vector(n,dummy,0);\ ff=f;\ for(j=1,n,\ th[j]=geodesic(ff);\ print1(j);print1(" ");print(th[j]);\ ff=cmp(f,ff));\ th computetheta2(f1,n1,f2,n2, th,dummy1,dummy2,ff,j)=\ th = matrix(n1,n2,dummy1,dummy2,0);\ ff1=f1;\ for(j1=1,n1,\ ff=ff1;\ for(j2=1,n2,\ th[j1,j2]=geodesic(ff);\ print1(j1);print1(" ");print1(j2);\ print1(" ");print(th[j1,j2]);\ ff=cmp(f2,ff));\ ff1=cmp(ff1,f1) );\ th complete(u,v,x,y) =\ if(v==0,x=1/u;y=0,x=(1/u)%v;y=(u*x-1)/v);\ [u,y;v,x] nextN(f)=\ a=level*f[3];b=-f[2];c=f[1]/level;\ while(c%level==0,c=c/level;a=a*level);\ k=(-b/c)%level;\ if(k>5,k=k-level,);\ print1(k);print1(" ");\ aa=c*k^2+b*k+a;\ if((aa%121)==0,if(k<0,k=k+level,k=k-level);aa=c*k^2+b*k+a,);\ bb=b+2*k*c;\ reduceN([aa,bb,c]))