# $Source: /u/maple/research/lib/algcurves/src/RCS/periodmatrix,v $ # $Notify: deconinc@newton.colorado.edu $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ macro( TRY_CCQUAD = 11, # higher accuracy in Int than this -> use _CCquad ERR1 = `Newton method failed, using more Digits...`, periodmatrix=`algcurves/periodmatrix`, Periodmatrix=`algcurves/Periodmatrix`, Contourintegrate=`algcurves/Contourintegrate`, Pathintegral=`algcurves/Pathintegral`, Partintegral=`algcurves/Partintegral`, Acy=`algcurves/Acy`, integrand2=`algcurves/integrand2`, integrate=`algcurves/integrate`, A1=`algcurves/A1`, # global variables. A2=`algcurves/A2`, Analyticcontinuation=`algcurves/Acontinuation`, Permuted=`algcurves/Permuted`, genus=algcurves['genus'], differentials=algcurves['differentials'], homology=algcurves['homology'], Whereis=`algcurves/Whereis`, Putfirst=`algcurves/Putfirst` ): #R6:=evalb(traperror(type(0,nonreal))<>lasterror): #if R6 then macro( # linsolve=LinearAlgebra['LinearSolve'], # delcols=LinearAlgebra['DeleteColumn'], # coldim=LinearAlgebra['ColumnDimension'] # ) #else macro( linsolve=LinearAlgebra:-LA_Main:-LinearSolve, delcols=LinearAlgebra:-LA_Main:-DeleteColumn, coldim=LinearAlgebra:-LA_Main:-ColumnDimension ): #fi: # This procedure computes a Riemann matrix for a Riemann surface, # given as a plane algebraic curve. For a Riemann surface of genus # g, the Riemann matrix is dependent on the choice of the homology, # but all Riemann matrices are equivalent upto a class of symplectic # transformations. These transformations then correspond to changes # in the choice of the homology. periodmatrix:=proc(curve,x::name,y::name) local pm,alpha,beta,g,F,A,i,j,m; # make sure that the input is of an allowed form option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved. Authors: B. Deconinck and M. van Hoeij.`; if indets(curve,'name') <> {x,y} or not type(curve,polynom(anything, [x, y])) then error "Curve should be a polynomial in both variables %1 and no other variables", {x,y} elif indets(curve,float)<>{} then error "No floating point coefficients allowed" fi; # make sure degree gives the correct result F:=collect(primpart(curve,{x,y}),{x,y},'distributed',Normalizer); g:=genus(F,x,y); if g<=0 then error "Input must have positive genus" fi; A:=[args[4..-1]]; # pass options to monodromy: _Env_algcurves_opt:=op(A); pm:=Periodmatrix(Digits,F,x,y); if member('Riemann',A) or member('normalized',A) then alpha:=delcols(pm,g+1..2*g,outputoptions=[]); beta:=delcols(pm,1..g,outputoptions=[]); pm:=linsolve(alpha,beta,method=none,free=_t0,conjugate=true, inplace=false,outputoptions=[]); # works better than evalm(alpha^(-1)) if g>1 then m:=max(seq(seq(abs(pm[i,j]-pm[j,i]),i=j+1..g),j=1..g-1)); i:=cat("About: ",convert(evalf(-log10(m),3),string), " accurate Digits"); if m>10^(3-Digits) then WARNING( i ) # Shouldn't happen. else userinfo(2,'algcurves',i) fi fi; # Test if Imaginary part is positive definite. m:=Matrix(g,g,[seq([seq(Im(pm[i,j]) ,j=i..g)],i=1..g)],scan=[triangular[upper]],shape=symmetric); if not LinearAlgebra:-LA_Main:-IsDefinite(m,query='positive_def') then error "Imaginary part not positive definite" fi; fi; # if member('Riemann',A) then # pm:=evalm(2*evalf(Pi)*I*pm) # fi; pm end: Periodmatrix:=proc(DI,curve,x,y) local h,g,pm,kappa,i,j,k,contours,diffs; options remember, `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; h:=homology(curve,x,y,`give paths`); kappa:=h[1]['linearcombination']; contours:=h[1]['cycles']; diffs:=differentials(curve,x,y,'skip_dx'); g:=nops(diffs); pm:=Matrix(g,2*g); for i to g do userinfo(1,'algcurves',`Integrating differential:`,diffs[i]/diff(curve,y)); for j to 2*g do pm[i,j]:=0; for k to coldim(kappa) do if kappa[j,k]<>0 then pm[i,j]:=pm[i,j]+kappa[j,k]*Contourintegrate( curve,diffs[i],x,y,contours[k], h[1]['basepoint'],h[1]['sheets'],h[3],h[3]['parameter'],h[4]) fi od od od; forget(Contourintegrate); forget(Pathintegral); forget(Partintegral); forget(integrand2); pm end: Contourintegrate:=proc(curve,differential,x,y,contour,x0,sheets,paths,t,radii) local n,m,i,integ,bpoint,firstsheet,lastsheet,perm,cycle; options remember, `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; integ:=0; n:=nops(contour)/2; cycle:=[op(contour),contour[1]]; for i from 1 to n do bpoint:=cycle[2*i][1]; perm:=cycle[2*i][2]; firstsheet:=cycle[2*i-1]; lastsheet:=cycle[2*i+1]; perm:=Putfirst(perm,firstsheet); for m from 1 to Whereis(perm,lastsheet)-1 do # Go around the branchpoint, and add the integral # over that path to integ: integ:=integ+ Pathintegral(curve,differential,x,y,paths,t,x0 ,sheets[perm[m]],sheets[perm[m+1]],bpoint ,`if`(bpoint=infinity, map(op,{indices(paths)}) minus {infinity,'parameter'}, radii[bpoint]) ) od od; integ end: Pathintegral:=proc(curve,differential,x,y,paths,t,x0,y0,y1,bpoint,radius) local i,n,integ,yi,yj,parint,rpoints,ipoints, rmax,imin,imax,v,path; options remember, `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; integ, yi, yj := 0, y0, y1; if bpoint=infinity then # integrate around the branchpoint rpoints:=[seq(Re(i),i=radius)]; ipoints:=[seq(Im(i),i=radius)]; imin, imax, rmax := min(op(ipoints))-1, max(op(ipoints))+1, max(op(rpoints))+1; v:=[x0, Re(x0)+I*imax, rmax+I*imax, rmax+I*imin, Re(x0)+I*imin, x0]; for i to 5 do parint:=Partintegral(curve,differential,x,y, v[i]*(1-t)+v[i+1]*t ,t,yi,Digits-2); integ:=integ+parint[1]; yi:=parint[2] od elif `algcurves/realcurve`(curve,x,y) and Im(bpoint)>0 then # Integrate around conjugate(bpoint) instead, so that the options # remember is more effective. n:=procname(args[1..7],seq(conjugate(i),i=[y1,y0,bpoint]),radius); return -conjugate(n) else # go to branchpoint path:=paths[bpoint]; n:=nops(path); for i to n do parint:=Partintegral(curve,differential,x,y,path[i],t,yi,Digits-2); yi:=parint[2]; integ:=integ+parint[1]; od; # integrate around the branchpoint parint:=Partintegral(curve,differential,x,y,bpoint-radius*exp(2*Pi*I*t),t,yi,Digits-2); yi:=parint[2]; integ:=integ+parint[1]; # integrate following the path in the reverse direction for i to n do parint:=Partintegral(curve,differential,x,y,path[i],t,yj,Digits-2); yj:=parint[2]; integ:=integ-parint[1]; od fi; if abs(yi-yj)>(1+abs(yi))*10^(-iquo(Digits,2)) then error "Path was not followed correctly during the integration" fi; integ end: Partintegral:=proc(curve,differential,x,y,path,t,y0,acc) local k,integrand,N; global A1,A2; options remember, `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; # userinfo(4,'algcurves',`Integrating over path`,path); Digits:=max(10,Digits); A1:=curve,x,y,path,t,Digits; A1:=Analyticcontinuation(A1),A1; N:=op(Permuted(A1[1][1][2],[y0])); # Number of current sheet. A1:=[N,A1]; A2:=subs(y=_X,differential); integrand:=proc(S) global A1,A2; local v; option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; if indets(S,name)={} then Digits:=A1[-1]; v, Digits:=integrand2(A1,evalf(S)); evalf(v[1]*subs(A1[4]=v[2],_X=v[3],A2)) else 'procname(args)' fi end: k:=traperror(integrate(integrand,acc)); if k=lasterror then Digits:=Digits+5; procname(args) else [k,A1[2][-1][2][N]] fi end: # Compute evalf(Int(integrand,0..1)) with acc accurate digits. integrate:=proc(integrand,acc) local k,method; global A1; option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; method:=`if`(nargs<4 and acc>TRY_CCQUAD,_CCquad,_NCrule); userinfo(10,'algcurves',cat("integration method ",method,", accuracy " ,acc,", Digits ",A1[-1]," (",Digits,"), nested ",nargs-2)); k:=traperror(evalf(Int(integrand,0..1,acc,method))); if k=lasterror then userinfo(1,'algcurves',"evalf/Int returned :",k); # Increase Digits by 3. A1:=subsop(-1=A1[-1]+3,A1); if k=ERR1 then error k else # Divide the interval, makes the integrand smoother. k:=args[2..-1],0; (procname(proc(s) integrand(s/2) end,k) +procname(proc(s) integrand(s/2+1/2) end,k))/2 fi else k fi end: integrand2:=proc(A,s) options remember, `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; local x0,y0,F,dFdy; x0:=evalf(subs(A[7]=s,A[6])); F:=evalf(subs(A[4]=x0,A[5]=_X,A[3])); y0, dFdy, Digits:=Acy(A[1],A[2],s,F); [evalf(subs(A[7]=s,diff(A[6],A[7])) / dFdy), x0,y0], Digits end: # s = number between 0 and 1 giving the point on the path. # x0 = subs(t=s,path) is the is the current point on the path. # F = subs(x=x0,y=_X,curve) in C[_X] # G = diff(F,_X) # Ac = Analytic continuation on that path. # N = the number of the sheet. # # Output: The y-value at x0 on sheet N (label numbers correspond to Ac). # Method: Use Ac to get an approximation for this y-value, and then use # Newton's method to improve the accuracy of this y-value as a # root of F. # Acy:=proc(N,Ac,s,f) local i,d,l,r,L,R,F,G; option `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved.`; i:=1; while s>Ac[i][1] do i:=i+1 od; r,R:=Ac[i][1],Ac[i][2][N]; if s-r=0 then d:=R; else l,L:=Ac[i-1][1],Ac[i-1][2][N]; d:=(s-l)/(r-l); d:=evalf(R*d+L*(1-d)); fi; F:=convert(f,horner); G:=convert(diff(f,_X),horner); r:=evalf(subs(_X=d,G)); if abs(r)<1/32 then Digits:=Digits-round(.435*ln(abs(r))); G:=convert(diff(f,_X),horner) fi; r:=evalf(subs(_X=d,F))/r; d:=d-r; i,l:=1,Digits; Digits:=max(5,l-2); R:=0.1^l; while i<4 or abs(r)>R do i:=i+1; if i>10 and i>evalf(2*ln(l)) then # Digits:=l; # return `fsolve/refine2`(d,F,G,l) error ERR1 fi; L:=evalf(subs(_X=d,G)); r:=evalf(subs(_X=d,F))/L; d:=d-r; Digits:=Digits+2 od; d, L, l end: #savelib('periodmatrix','Periodmatrix','Contourintegrate',\ 'Pathintegral','Partintegral','Acy','integrand2','integrate'):