\\ Package to compute the "derived period integrals" attached to a
\\  (modular) elliptic curve E/Q and a real quadratic field K of narrow 
\\  class number one for which E has sign (-1,-1).
\\ 

\\ The function integral(r,t1,t2) returns a numerical approximation to 
\\  the integral, from r+it1 to r+it2, of the function e^z/z.
\\  It is assumed that r is strictly negative.
\\

intezoverz(r,t1,t2,\
 er,rp,ip,t) =\
 er = exp(r);\
 rp = -er*intnum(t=t1,t2, (r*sin(t)-t*cos(t))/(r^2+t^2) );\
 ip =  er*intnum(t=t1,t2, (r*cos(t)+t*sin(t))/(r^2+t^2) );\
 return(rp + I*ip);


\\ The function period(gammapsi) returns the derived BD period 
\\ attached to gammapsi, correspoding to a fundamental unit of
\\ the corresponding real quadratic field.
\\

period(E,gammapsi,\
 a,b,c,d,z1,gz1, xpsi,ypsi, lim1,lim2,lim3,lim4, res, r1, r2, t1,t2,\
 digits, M1, M2, q, n, fc, trc) =\
 digits = 5;\
 a = gammapsi[1,1]; b = gammapsi[1,2];\
 c = gammapsi[2,1]; d=gammapsi[2,2];\
 trc = a+d;\
 print([a,b,c,d]);\
 z1 = -d/c+ I/abs(c);\
 gz1 = a/c +I/abs(c);\
 epsilon = (trc+sqrt(trc^2-4))/2;\
 print("unit = ", epsilon);\
 xpsi = (epsilon-d)/c;\
 ypsi = (1/epsilon -d)/c;\
 lim1 = 2*Pi*I*(z1-xpsi) ;\
 lim2 = 2*Pi*I*(gz1-xpsi) ;\
 lim3 = 2*Pi*I*(z1-ypsi) ;\
 lim4 = 2*Pi*I*(gz1-ypsi) ;\
 r1 = real(lim1);\
 if (r1<> real(lim2), print("ERROR!!!"),print("ok"));\
 r2 = real(lim3);\
 if (r2<> real(lim4), print("ERROR!!!"),print("ok"));\
 print("r2 =", r2);\
 M1 = floor(digits/abs(r1)*log(10)+ 1/abs(r1)*abs(log(abs(r1)/(2*Pi*(a+d))))) + 1;\
 print("M1 = ", M1);\
 M2 = floor(digits/abs(r2)*log(10)+ 1/abs(r2)*abs(log(abs(r2)/(2*Pi*(a+d))))) + 1;\
 M=max(M1,M2);\
 fc = ellan(E,M);\
 gq1 = exp(2*Pi*I*gz1); q1 = exp(2*Pi*I*z1);\
 per= sum(n=1,M,fc[n]/n*(gq1^n-q1^n));\
 print("period = ", per);\
 res=0;\
 t1 = imag(lim1); t2 = imag(lim2);\
 print( "Computing sum with ", M1, " terms");\
 q = exp(2*Pi*I*xpsi);\
 for(n=1,M1, \
   res = res - fc[n]*q^n* intezoverz(n*r1,n*t1,n*t2) ;\
   print1([n]); );\
 t1 = imag(lim3); t2 = imag(lim4);\
 print( "Computing sum with ", M2, " terms");\
 q = exp(2*Pi*I*ypsi);\
 for(n=1,M2, \
   res = res + fc[n]*q^n* intezoverz(n*r2,n*t1,n*t2) ;\
   print1([n]); );\
 return(res);
 
 
   
E = ellinit([0,0,-1,-1,0]);

ellglobalred(E);


/* 
This is some calculation for Q(sqrt(3)) 
gammapsi = [17,6; -37, -13];
period(E,gammapsi)

gammapsi = [17,-6; 37, -13];
period(E,gammapsi)  

gammapsi = [-20,-13; 37, 24];
period(E,gammapsi)  
*/


/*  Calculation for Q(sqrt(53))    */

gammapsi = [-90,-49; 259, 141 ];
per = 2*Pi*I*period(E,gammapsi)
print(gammapsi)