##################################################################### # standard `algcurves` macros # read("../edgardo/AbelMap_3Aug.mpl"); macro( submatrix = linalg['submatrix'] , linsolve = linalg['linsolve'] , coldim = linalg['coldim'] , Whereis = `algcurves/Whereis` , homology = `algcurves/homology` , Propagate = `algcurves/Propagate` , Partintegral = `algcurves/Partintegral` , Pathintegral = `algcurves/Pathintegral` , AC = `algcurves/Acontinuation` , puiseux = algcurves['puiseux'] , genus = algcurves['genus'] , monodromy = algcurves['monodromy'] , differentials = algcurves['differentials'] , periodmatrix = algcurves['periodmatrix'] ): # rcv ------ RiemannConstantVector := module() #options package; export CanonicalDivisor , AllBranches , SmallCD ; local ModuleApply , input_check , kappa_base , MakeH , RefineViaThetaDivisor , FastTheta , Abel_diff , In_SGamma , Mylimit , Valuation , Next_P , nearest , ProblemPoints , StandardCurve ; # rcv ------------------------------------------------------------------------------------------ # the traffic director for the Riemann Constant Vector ModuleApply := proc(curve, x :: name, y :: name, P0, accuracy) local F, kappa, g, abel, i, K; F := StandardCurve(curve, x, y); # Check that the input is formatted as needed and that `P0` can be # determined input_check(F, x, y, P0); g := genus(F, x, y); # Abel transform from P0 to base # 'ZERO' is a flag to cause the use of basepoint chosen by monodromy abel := AbelMap(F, x, y, P0, 'ZERO', accuracy); # Hand off to computational engine for point-independent bits # Any other RCV is an Abel map away kappa := kappa_base(F, x, y, accuracy); # Put it all together K := [seq(evalf((kappa[1][i] - (g - 1)*abel[i])), i = 1..g)]; # Reduce modulo the lattice AbelMap:-ModPeriodLattice(K, periodmatrix(F, x, y, 'Riemann'), `if`(has([args], 'fraction'), fraction, NULL)); end: ##################################################################### # @input_check # Cheeck that the input is in correct form, on curve # and one of `algcurves[puiseux]` expansions can be chosen. input_check := proc(curve, x::name, y::name, inP) local temp, defaultMessage, P, Px, L, Y, test, match, r, res, i, j, k, low_y, tau, t; # supress mint complaint t := t; # If `inP` has the index 'valuation' then it does not need checked. if has(inP, 'valuation') then return inP; fi; # Avoid mint noise, `tau` and `Px` are both used as name/flag Px, tau := Px, tau; defaultMessage := "`algcurves[puiseux]` may be helpful."; # Basic type-checking first. # Check that inP is either a list or the flag 'ZERO' ( = basepoint) if evalb(inP = 'ZERO') then return 'ZERO'; elif not type(inP, list) and not type(inP, table) then error("The starting point of the Abel map should be input as a list or as flag 'ZERO'\ , but received ", inP, defaultMessage); fi; if hastype(inP, float) then error("Floats not supported by RiemannConstantVector.", defaultMessage); fi; P := `if`(hastype(inP, 'table'), copy(inP), convert(inP, 'table')); # Convert to RootOf P[x] := convert(P[x], RootOf); P[y] := convert(P[y], RootOf); # check that there is only one name in P[x] if nops(indets(P[x])) > 1 then error("Please enter %1-component as either numeric value or series in one parameter.", x); fi; # check if the user gave a parameter, introduce one if not P['parameter'] := `if`(nops(indets(P[x])) = 1, op(indets(P[x])), t); if not(has(P[x], P['parameter'])) then P[x] := P[x] + P['parameter'] fi; # make sure no extra parameters floating about if not({op(map(op, [indices(P)]))} minus {x, y, 'valuation', 'dfdy_degree', 'made_float', 'parameter'} = {}) or not(indets(P[x]) = {P['parameter']}) or not(indets(P[y]) minus {P['parameter']} = {}) then error("Please input point as a list of form [x = x_0+b*t^r, y = c*t^s+c'*t^(s')...],\ or [x = x_0, y = y_0]", defaultMessage); fi; temp := coeffs(P[x], P['parameter'], Px); temp := {op(map(degree, [Px], P['parameter']))}; if evalb(nops(temp) > 2 or nops(temp minus {0}) <> 1 or not(type(op(temp minus {0}), integer))) then error("Expansion for %1 should be of form a+b*%2^r with {a, b} algebraic and r an integer." , x, P['parameter'], defaultMessage); fi; # passed basic type-checking. # Make sure user's expansion can be matched to algcurves[puiseux]. # If user's point is finite and far enough from discriminant points, # then it is easy to determine which point is meant. L := LIMIT(P, x, y, P['parameter']); r := degree(P[x], P['parameter']); # The lowest `t`-degree in `y` needed for comparison low_y := `if`(P[y] = 0, 0, ldegree(P[y], P['parameter'])); if not(has(L[x], infinity)) and not(has(L[y], infinity)) then # !!!!!!!!!! MAYBE REPLACE THE CONDITION BELOW WITH ONE THAT CECKS THAT, # GIVEN X_VALUE, Y_VALUE CAN BE PICKED WITHIN DIGITS!!!!!!!!!!!!!!!!!!!!! if not evala(eval(curve, {x = L[x], y = L[y]})) = 0 then error("The given point %1 is not on %2", inP, curve, defaultMessage); fi; Y := ListTools:-MakeUnique([fsolve(eval(curve, x = L[x]), y, complex)]); if nops(Y) = degree(curve, y) and is_resolved_abel(Y, L[y]) then # This is a regular point, make sure branching number is unity. if r = 1 then # If `y`-value is resolved and there are `N` fibre-values just ask Next_P(uiseux). return Next_P(curve, x, y, P['parameter'], [x = P[x], y = L[y]]); else error("Exponent in %1 expansion should be\ unity around a regular point. %2", x, defaultMessage); fi; else # Finite and not resolved. Implies "near" discriminant point. # Re-use name `Y`. Y := AllBranches(curve, x, L[x], y, 0, tau); match := [seq(true, i = 1..nops(Y))]; for i to nops(Y) do test := convert(Y[i], table); if degree(test[x], tau) <> r or (ldegree(test[y], tau) <> low_y and low_y <> 0) then # The exponent of `t` in `x` expansion must be branching number. # Thus the low degrees in `y` must also match (or be zero if under-resolved). match[i] := false; else # Compute the `r` possible `res`iduals of P[y]-test[y] res := [seq(eval(test[y] , tau = P['parameter']*RootOf(_Z^r-coeff(P[x], P['parameter']^r) /coeff(test[x], tau^r) , index = k)) - P[y], k = 1..r)]; # Now test each index to see if it matches to order. for j to r do # Test order by order. for k from ldegree(P[y], P['parameter']) to min(low_y, degree(test[y], tau)) do # This coeff is not zero, break as it is not a match. if not evalc(evala(coeff(res[j], P['parameter'], k))) = 0 then res[j] := false; break; fi: od; od: # Now test if any `res` was zero to order; set match[i] to true if it is. match[i] := evalb(op(ListTools:-MakeUnique(res)) <> false); fi; od: fi; else # Infinity. Y := AllBranches(curve, x, L[x], y, 0, tau); match := [seq(true, i = 1..nops(Y))]; for i to nops(Y) do test := convert(Y[i], table); if degree(test[x], tau) <> r or (ldegree(test[y], tau) <> low_y and low_y <> 0) then # The exponent of `t` in `x` expansion must be branching number. # Thus the low degrees in `y` must also match (or be zero if under-resolved). match[i] := false; else # Compute the `r` possible `res`iduals of P[y]-test[y] res := [seq( eval(test[y], tau = P['parameter']*RootOf(_Z^abs(r) - coeff(test[x], tau^r) /coeff(P[x], P['parameter']^r), index = k)) - P[y] , k = 1..abs(r))]; # Now test each index to see if it matches to order. for j to abs(r) do # If it is zero it is a match. if res[j] = 0 then match[i] = true; break; fi; # Test order by order to avoid over-evala-ing. for k from ldegree(P[y], P['parameter']) to min(low_y, degree(test[y], tau)) do # This coeff is not zero, break as it is not a match. if not evala(coeff(res[j], P['parameter']^k)) = 0 then res[j] := false; break; fi: od; od: # Now test if any `res` was zero to order; set match[i] to true if it is. match[i] := evalb(op(ListTools:-MakeUnique(res)) <> false); fi; od; fi; # Check that there is but one match. temp := nops(map(z -> `if`(z, true, NULL), match)); if temp = 1 then # There is one match, return it. eval(Y[seq(`if`(match[i], i, NULL), i = 1..nops(match))], tau = P['parameter']); elif temp > 1 then error("The expansion %1 matches more than one branch. %2", inP, defaultMessage); else error("The expansion %1 matches no branch. %2", inP, defaultMessage); fi; end: # kappa_base --------------------------------------------------------------------------------- # Compute the kappa-bits for P0=basepoint (x=monodromy[1],y=monodromy[2][1]) kappa_base := proc(curve, x :: name, y :: name, accuracy) local CD, B, g, K_fog, A, H, abs_K, Index, i, FiniteSum, cold_filtered, dK , FS_args, rho, Rb, Ib, IbInv, IbChol, j, k, eps1, eps2, R1, R2, dA; options cache; g := genus(curve, x, y); # Digits must be increased if genus is high enough Digits := max(Digits, round(log10(2^(2*g + 1)))); B := periodmatrix(curve, x, y, 'Riemann'); userinfo(3, 'FastTheta', `Finished computing Riemann Matrix`); # Compute an efficient cannonical divisor (Maybe use SMALLCD) CD := AbelMap:-CanonicalDivisor(curve, x, y, 't'); # Compute the Abel transform of the divisor. A := AbelMap:-DivisorAbelMap(curve, x, y, 'ZERO', CD, Digits); # "Divide by -2," at least in fundamental parallelogram K_fog := AbelMap:-ModPeriodLattice(-A, B)/2; # Compute decompositions needed to create FiniteSum and # compute shortest lattice vector # Make `B` symmetric in order to get Cholesky decomp of Im(B). for i from 1 to g do for j from i to g do if i <> j then B[i, j] := (B[i, j] +B[j, i])/2; B[j, i] := B[i, j]; fi; od: od: # The compiled version must have Re/Im perts separate. Rb, Ib := map(Re, B), map(Im, B); # Need inverse and Cholesky decomposition of Im(Riemann matrix). IbInv, IbChol := linalg:-linsolve(Ib, array('identity', 1..g, 1..g)), linalg:-cholesky(Ib); # Compute the shortest lattice vector in Lambda (via LLL) rho := sqrt(evalf(Pi)*add(j^2, j in sort(IntegerRelations:-LLL([seq([seq(IbChol[i, j], j = 1..g)] , i = 1..g)]) , (a,b) -> evalb(add(i^2, i in a) < add(i^2, i in b) ))[1])); # Compute epsilon and R (the "radius of the bounding ellipsoid")for final pass as per # disertation (MP) eps2 := evalf(1/2^(2*g + 1)); R2 := max((sqrt(2*g) + rho)/2 , fsolve(eps2 = g*2^(g - 1)/rho^g*GAMMA(g/2, (_Z - rho/2)^2), _Z = rho, 'fulldigits')); # Compute epsilon_1 for the initial filtering, first need R1 R1 := max((sqrt(2*g) + rho)/2 , fsolve(g*2^(g - 1)/rho^g*GAMMA(g/2, (_Z - rho/2)^2)*R2^g - _Z^g, _Z = R2, 'fulldigits')); eps1 := g*2^(g - 1)/rho^g*GAMMA(g/2, (R1 - rho/2)^2); # Not much point in having an eps1 > 0.1 (??) eps1 := min(eps1, 0.1); # Use the Riemann theta function to check which half-lattice vector gives correct K FiniteSum, FS_args := FastTheta(K_fog, eps1, B, Rb, Ib, IbInv, IbChol); # Test for zeros: compute them all on first pass, as per MSP_DSRT abs_K := [seq(FiniteSum(i, op(FS_args)), i = 1..2^(2*g))]; # Test that AT LEAST one is still a candidate if min(op(abs_K)) < eps1 then # Test if the lowest value is MUCH smaller than the rest. cold_filtered := [seq(`if`(abs_K[i] < eps1, i, NULL), i = 1..2^(2*g))]; if nops(cold_filtered) > 1 then # Compute a perturbation to K as per MSP_DSRT, # `In_SGamma` returns (g - 1) Abel maps. dA := In_SGamma(curve, x, y, max(accuracy, -ceil(log10(eps2)))); # Perturb the candidates by dK, of the many possible choices... #dK := [seq(0, i = 1..g)]; #dK := [seq(add((g - 1)*dA[j][i], j = 1..1), i = 1..g)]; dK := [seq(add(dA[j][i], j = 1..(g - 1)), i = 1..g)]; # Refinement step # Note the "perturbation" by a vector in A(SGamma^{g - 1}) FiniteSum, FS_args := FastTheta([seq(K_fog[i] + dK[i], i = 1..g)] , eps2, B, Rb, Ib, IbInv, IbChol); abs_K := convert(abs_K, table); for i in cold_filtered do abs_K[i] := FiniteSum(i, op(FS_args)); od: # Make sure AT LEAST one is still valid if min(op(map(op, [entries(abs_K)]))) < eps2 then cold_filtered := [seq(`if`(abs_K[i] < eps2, i, NULL), i in cold_filtered)]; # Test if nops(cold_filtered) > 1 then # Final stage. Try to refine by using Theta-Divisor # Ref...isor returns appropriate error # messages if it fails. Sucsess is assumed here Index := RefineViaThetaDivisor(K_fog, cold_filtered, dA, B, eps2, Rb, Ib, IbInv, IbChol); H := MakeH(Index, g); return K_fog - [seq(add(B[j, k]*H[k], k = 1..g), j = 1..g)] - H[g + 1..2*g], H; else # If there is only one, then return the correctly refined `kappa_base`. Index := cold_filtered[1]; H := [seq((1/2)*`if`(j = 1,`mod`(Index, 2), `mod`(trunc(Index/(2^(j - 1))), 2)) , j = 1..2*g)]: return K_fog - [seq(add(B[j, k]*H[k], k = 1..g), j = 1..g)] - H[g + 1..2*g], H; fi; else error "%1 failed to find a zero of the Riemann theta-function.",procname; fi; else # There was only one that survived the first filter. Index := Whereis(abs_K, min(op(abs_K))); H := [seq((1/2)*`if`(j = 1,`mod`(Index, 2), `mod`(trunc(Index/(2^(j - 1))), 2)) , j = 1..2*g)]: FiniteSum, FS_args := FastTheta([seq(K_fog[i], i = 1..g)] , eps2, B, Rb, Ib, IbInv, IbChol); return K_fog - [seq(add(B[j, k]*H[k], k = 1..g), j = 1..g)] - H[g + 1..2*g], H; fi; else error "%1 failed to find a zero of the Riemann theta-function.",procname; fi; end: # end kappa_base ----------------------------------------------------------------------------- # MakeH # Construct the `Index`-th half-lattice vector # Note: Be carfull about [I|B] vs. [B|I] MakeH := proc(Index, g) local h, j; h := seq((1/2)*`if`(j = 1,`mod`(Index, 2), `mod`(trunc(Index/(2^(j - 1))), 2)) , j = 1..2*g): [h[g + 1..2*g], h[1..g]] end: # RefineViaThetaDivisor ---------------------------------------------------------------------- RefineViaThetaDivisor := proc(K, cold_filtered :: list, A, B, epsilon, Rb, Ib, IbInv, IbChol) local i, j, g, CF, FiniteSum, FS_args, dK, test; # There are (g - 1) vectors in A. The later ones in the list arise # as the Abel map to places that are farther from the left place, # one would suppose that these would tend to be "more different." # As such, use these first. Give up after trying twice for each dude. CF := cold_filtered; g := nops([entries(A)]) + 1; for i from g - 1 to 1 by -1 do # Create the perturbation. loop goes backwards as larger Abel # maps tend to be at the "end" of A dK := [seq((g - 1)*A[i][j], j = 1..g)]; # Create a compiled version FiniteSum, FS_args := FastTheta([seq(dK[j] + K[j], j = 1..g)], epsilon, B, Rb, Ib, IbInv, IbChol); for j in CF do test := FiniteSum(j, op(FS_args)); if test < epsilon then # Leave it in else # Take it out CF := remove(z -> evalb(z = j), CF); fi; od: od: # Exception handling if nops(CF) = 1 then return CF[1]; elif nops(CF) = 0 then error("Computing the Riemann constant vector of the left place failed."); else error("Computing vector of Riemann constants of the left place resulted in non-unique vector."); fi; end: # end RefineViaThetaDivisor ------------------------------------------------------------------ # FastTheta ---------------------------------------------------------------------------------- # Compute the kappa-bits for P0=basepoint (monodromy[1]), attempt to compile FastTheta := proc(K, epsilon, B, InRb, InIb, InIbInv, InIbChol) local temp, g, i, FiniteSum, F, carol, ReK, ImK, Z, ans, Rb, Ib, IbInv, IbChol, datatype; g := linalg:-rowdim(B); # Mint supression datatype := datatype; Rb, Ib := convert(InRb, vector), convert(InIb, vector); IbInv, IbChol := convert(linalg:-transpose(InIbInv), vector), convert(linalg:-transpose(InIbChol), vector); # Vectors to arrays. Rb := Array(1..g^2, Rb, `datatype` = float[8]); Ib := Array(1..g^2, Ib, `datatype` = float[8]); IbInv := Array(1..g^2, IbInv, `datatype` = float[8]); IbChol := Array(1..g^2, IbChol, `datatype` = float[8]); userinfo(5, 'FastTheta', `Finished conversions; computed radius costing.`, time() - temp); # The Riemann theta function to compile. Z := [seq('zeta'[i], i = 1..g)]; temp := time(); F := op(RiemannTheta(Z, B, [], epsilon, 'output' = list)[2]); temp := time() - temp; userinfo(3, 'FastTheta', `Computing the needed lattice-vectors took.`, temp); # F[4]=Lattice vectors needed for accurate approximation userinfo(5, 'FastTheta', `Number of lattice-vectors.`, nops(F[4])); temp := time(); # Build a local copy of Bernard's `RiemannTheta/finitesum` to compile. # Why does it take longer to compile as genus goes up?? FiniteSum:=proc(Index :: posint, gee :: posint, NumIVs :: posint , ReK :: Array(`datatype` = float[8]), ImK :: Array(`datatype` = float[8]) , ReB :: Array(`datatype` = float[8]), ImB :: Array(`datatype` = float[8]) , InvImB :: Array(`datatype` = float[8]), CholImB :: Array(`datatype` = float[8]) , H :: Array(`datatype` = float[8]) , ReZ :: Array(`datatype` = float[8]), ImZ :: Array(`datatype` = float[8]) , IVs :: Array(`datatype` = float[8]) , intshift :: Array(`datatype` = float[8]) , damp :: Array(`datatype` = float[8]), frag :: Array(`datatype` = float[8]) ) :: float; local i :: integer, j :: integer, k :: integer, temp :: float[8], ReT :: float[8] , ImT :: float[8], Retemp :: float[8], Imtemp :: float[8]; # First make the half-lattice vector. for i from 1 to 2*gee do if i=1 then H[i]:=-0.5*`mod`(Index, 2); else temp:=evalhf(Index)/2.^(i - 1); H[i]:=-0.5*`mod`(trunc(temp), 2); end if; end do: # Make the candidate vector for i from 1 to gee do ReZ[i]:=ReK[i]+H[gee+i]; ImZ[i]:=ImK[i]; for j from 1 to gee do ReZ[i]:=ReZ[i]+ReB[(i-1)*gee+j]*H[j]; ImZ[i]:=ImZ[i]+ImB[(i-1)*gee+j]*H[j]; end do: end do: # Compute the vector "damping" bits in exponetial (See DvH2, p.3) for i from 1 to gee do damp[i]:=0.0; for j from 1 to gee do damp[i]:=damp[i]+InvImB[(i-1)*gee+j]*ImZ[j]; end do; intshift[i]:=round(damp[i]); frag[i]:=damp[i]-intshift[i]; end do: ReT:=0; ImT:=0; for i from 1 to NumIVs do # Look into excluding certain lattice-vectors(U-wise vs P-wise). Retemp:=0.0; Imtemp:=0.0; for j from 1 to gee do # Add the bits from and 0.0. Asymmetry due to factoring out growth as per DvH2. No shift in yet. temp:=0.0; Imtemp:=Imtemp+ReZ[j]*(IVs[(i-1)*gee+j]-intshift[j]); for k from 1 to gee do temp:=temp+CholImB[(j-1)*gee+k]*(IVs[(i-1)*gee+k]+frag[k]); if k=j then Imtemp:=Imtemp+ 0.5*ReB[(j-1)*gee+k]*(IVs[(i-1)*gee+j]-intshift[j])*(IVs[(i-1)*gee+k]-intshift[k]); elif k>j then Imtemp:=Imtemp+ ReB[(j-1)*gee+k]*(IVs[(i-1)*gee+j]-intshift[j])*(IVs[(i-1)*gee+k]-intshift[k]); fi; #Imtemp:=Imtemp+ 0.5*ReB[(j-1)*gee+k]*(IVs[(i-1)*gee+j]-intshift[j])*(IVs[(i-1)*gee+k]-intshift[k]); end do: Retemp:=Retemp+temp^2; end do: ReT:=ReT+exp(-Pi*Retemp)*cos(2*Pi*Imtemp); ImT:=ImT+exp(-Pi*Retemp)*sin(2*Pi*Imtemp); od: return (ReT^2+ImT^2)^(0.5); end proc: # Make the arguments that go into compiled function ReK := Array(1..g, map(Re, K), `datatype` = float[8]); ImK := Array(1..g, map(Im, K), `datatype` = float[8]); carol:=[g, nops(F[4]) # The vector(Re/Im) around which we check the half-latti. , ReK, ImK # The real, imaginary, inverse imaginary and cholesky imaginary periodmatrix , Rb, Ib, IbInv, IbChol # A spot for the half-lattice offsets ,Array(1..2*g, [seq(0.0, i = 1..2*g)], `datatype` = float[8]) # Spots for the Re&Im bits of the vector argument. ,Array(1..g, [seq(0.0,i = 1..g)], `datatype` = float[8]) ,Array(1..g, [seq(0.0,i = 1..g)], `datatype` = float[8]) # The lattice vectors for `accuracy` ,Array(1..g*nops(F[4]), ListTools:-Flatten(convert(F[4], list)), `datatype` = float[8]) # Spots for the shift and fractional part. ,Array(1..g, [seq(0.0, i = 1..g)], `datatype` = float[8]) ,Array(1..g, [seq(0.0, i = 1..g)], `datatype` = float[8]) # Spot for the damping bits of the exponential argument. ,Array(1..g, [seq(0.0, i = 1..g)], `datatype` = float[8]) ]; temp := time() - temp; ans := Compiler:-Compile(FiniteSum, 'optimize' = true), carol; userinfo(3, 'FastTheta', `Time to compile`, temp); return ans; end: # end FastTheta ------------------------------------------------------------------------------ # @SmallCD ----------------------------------------------------------------------------------- # An attempt to create a small CannonicalDivisor at infinity SmallCD:=proc(curve, x::name, y::name, t::name) local temp, g, Omega, omega, disc, XINF, target, DSL, ans, J, dfdy, d, Om, P, i, j; options remember; # Get a basis of (`g`=genus) differentials without denomiantor. dfdy := diff(curve, y); Omega := Abel_diff(curve, x, y, false); g := nops(Omega); # Puiseux expansions of the points at (infinity,infinity) XINF := [seq(convert(i, table), i in AllBranches(curve, x, infinity, y, 0, t))]; # A canonical divisor should be degree 2g - 2 target := 2*g - 2; # sequence of [differential numerator, degree of differenial] temp := [seq([Omega[i], add( degree(diff(XINF[j][x], t), t) + degree(eval(Omega[i], {x = XINF[j][x], y = XINF[j][y]}), t) - XINF[j]['dfdy_degree'] , j = 1..nops(XINF))], i = 1..g)]; # Select the diff with largest degree less than the target (remove # call SHOULD be superfluous) omega := op(select(z -> evalb(z[2] = max(op(remove(z -> evalb(z > target), map2(op, 2, temp))))), temp)); # The divisor based just on places at infinity ans := [seq([degree(diff(XINF[j][x], t), t) + degree(eval(omega[1], {x = XINF[j][x], y = XINF[j][y]}), t) - XINF[j]['dfdy_degree'] # the rest of this line is formatting as downstream procs need , [x = XINF[j][x], y = XINF[j][y], 'valuation' = XINF[j]['valuation'] , 'dfdy_degree' = XINF[j]['dfdy_degree'], 'parameter' = XINF[j]['parameter']]] , j = 1..nops(XINF))]; if omega[2] = target then # There are the requisite number of zeros, so return divisor. return ans; else # Check the discriminant points. The new target is old minus # the zeros already found target := target-omega[2]; # Calculate the discriminant points with respect to `y` and # roots of leading polynomial coefficient _EnvExplicit := false; disc := ProblemPoints(curve, x, y); # Make sure eavla thinks these numbers are in a "simple form". disc := map(evala, disc); # Group "like" discriminant points together and arrange by the size of these groups. # D(iscriminantpoint)S(tring)L(ength): length of the dpoint as a string. # This is purely to avoid complicated (ie, computationally # involved )RootOf expresions DSL := [seq([nops(convert(convert(i, string), list)), i], i in disc)]; # How many of each "type" are there? DSL := [seq([nops(select(z -> evalb(z[1] = DSL[i][1]), DSL)), DSL[i][2]], i = 1..nops(DSL))]; # Sort them by how many there are. If alpha_1 and alpha_2 are # "similar" algebraic numbers, then BOTH would be likely to be # involved in a particular divisor. If LOTS are simliar, then # a divisor would perhaps contain all of them, and require a # lot of computation. disc := sort(DSL, (a,b) -> evalb( a[1] < b[1] )); for d in disc do # Compute the lifts over discpoint... P := [seq(convert(j, table), j in AllBranches(curve, x, d[2], y, 0, t))]; # ...and for each lift for J in P do # Differential with the denominator Om := eval(omega[1]/dfdy, {x = J[x], y = J[y]})*diff(J[x], t); # Determine if place(or lift) is a zero if limit(Om, t = 0) = 0 then # If place is a zero, update the divisor ans := [op(ans), [ldegree(eval(omega[1]*diff(J[x], t), {x = J[x], y = J[y]}), t) - J['dfdy_degree'] # rest of line is formatting (downstream) , [x = J[x], y = J[y], 'valuation' = J['valuation'] ,'dfdy_degree' = J['dfdy_degree']]]]; target := target - ans[-1][1]; if target = 0 then return ans fi; fi; od: od: error "The differential with adjoint polynomial %1 appears to have %2 too few zeros", omega[1], 2*g - 2 - target; fi; end: # end SmallCD -------------------------------------------------------------------------------- # Abel_diff ---------------------------------------------------------------------------------- # Compute differentials "in order" # That is, sorted on {x, y} degree, then on {x} degree Abel_diff := proc(curve, x::name, y::name, flag::boolean) local dFdy, diffs; options remember; diffs := sort(differentials(curve, x, y, 'skip_dx') , (a, b) -> evalb((degree(a, {x, y}) < degree(b, {x, y})) or ((degree(a, {x, y}) = degree(b, {x, y})) and (degree(a, x) > degree(b, x))))); if flag = false then # Return without denominator return diffs else dFdy := diff(curve, y); # Return with denominator return map(z -> z/dFdy, diffs); fi; end proc: # end Abel_diff ------------------------------------------------------------------------------ # @In_SGamma ---------------------------------------------------------------------------------- # Compute a Abel map (g - 1) places (with 'ZERO' as initial place) # The sum of any effective combo of degree = (g-1) is in the (g - 1)th # symetric power of Gamma In_SGamma := proc(curve, x :: name, y :: name, acc) local g, M, base, X, i, tau, Y, A; g := genus(curve, x, y); M := monodromy(curve, x, y, `give paths`); _Env_Explicit := false; # Generate (g - 1) perturbed x-values to the left of left-point base := convert(M[1], fraction); X := [seq(base - i, i = 1..(g - 1))]; # Compute the coresponding Ys (all on sheet 1) for i from 1 to nops(X) do # Analytically continue to stay on sheet 1 Y := convert(AC(curve, x, y, (1 - 'tau')*base + 'tau'*X[i], 'tau')[-1][2][1], fraction); # Back to algebraic numbers Y := nearest([solve(eval(curve, x = X[i]), y)], Y)[1]; A[i] := AbelMap(curve, x, y, 'ZERO', [x = X[i], y = Y], acc); od: A; end proc: # end In_SGamma ------------------------------------------------------------------------------- # @CanonicalDivisor -------------------------------------------------------------------------- # Compute `g` Canonical Divisors. CanonicalDivisor:=proc(curve, x::name, y::name, t::name, accuracy) local P, ans, ANS, F, g, omega, disc, zeros, test_places, count, i, j, k, Omega, found; # Set the environment to not give radicals _EnvExplicit := false; # Curve in preferred form (As per DvH code). F := StandardCurve(curve, x, y); # Genus g := genus(F, x, y); # Sort the differentials with denoms, higher degree may get 2g - 2 zeros in # one evaluation omega := ListTools:-Reverse(Abel_diff(F, x, y, true)); # Calculate the discriminant points with respect to `y` disc := ProblemPoints(F, x, y); # Make sure Mark thinks these numbers are reduced. disc:=map(evala,disc); for i from 1 to nops(omega) do # Initialize some stuff. ans := NULL; found := 0; # Calculate the resultant of the numerator of the differential # and the curve; here and the discpoints are where this differential may vanish. zeros := ListTools:-MakeUnique([op(map(allvalues , [solve(resultant(F, numer(omega[i]), y), x)], 'implicit')), op(disc)]); for j from 1 to nops(zeros) do test_places := AllBranches(F, x, zeros[j], y, 0, t); for k from 1 to nops(test_places) do # Evaluate the numerator of differential_i at each of these test places. P := convert(test_places[k], table); Omega := eval(omega[i], {x = P[x], y = P[y]})*diff(P[x], t); # THIS IS THE LINE THAT IS (OR WAS) TAKING ALLLLLL THE TIME. # Maybe judicious use of testeq can help. if Mylimit(Omega, t = 0) = 0 then # THIS IS THE LINE THAT IS TAKING ALLLLLL THE TIME. found := found + 1; ans := ans,[ldegree(numer(Omega), t) - ldegree(denom(Omega), t) , [x = P[x], y = P[y], 'valuation' = P['valuation'], 'dfdy_degree' = P['dfdy_degree']]]; if found=2*(g - 1) then break fi; fi; od: od: # Check that the degree of the divisor is 2g-2 if nops([ans]) = 0 or not( add([ans][count][1], count = 1..nops([ans])) = 2*g - 2 ) then # There is an emergency going on. error("There is an emergency going on."); else ANS[i] := [ans]; if not(has({args}, 'all')) then break; fi; fi; od: forget('AllBranches'); ANS end: # end CannonicalDivisor ----------------------------------------------------------------------- # @Mylimit ------------------------------------------------------------------------------------ # Just used for profiling as of DEC6 Mylimit:=proc(expr, point) limit(expr,point) end: # end Mylimit --------- ----------------------------------------------------------------------- # @Allbranches -------------------------------------------------------------------------------- # get all the representations over x=a AllBranches:=proc(curve, x::name, A, y::name, N::integer, t::name) options cache; local alpha, places, place, Place, ans, a, L; # Supress mint alpha := alpha; a := `if`(hastype(A, radical), convert(A, RootOf), A); # Get the Puiseux expansions from van Hoeij in form I want. places := map(collect, map(expand, convert(map( op, {puiseux(curve, x = a, y, N, t)}), list)), t); # Conjugate by `allvalues` places := ListTools:-MakeUnique(map(allvalues, places, 'implicit')); # Initialize some stuff: `L` counts branches L:=0; ans:=NULL; # Go. for place in places do Place := convert(place, table); Place['parameter'] := t; # Attach the `valuation`s for use in finding subsequent series. Place['valuation'] := Valuation(eval(curve, {x = Place[x], y = Place[y] + alpha*t^(`if`(Place[y] = 0, 0, degree(Place[y], t)) + 1)}), t, alpha) - 1; # Calculate the degree in `t` like which dfdy vanishes for when computing integrals Place['dfdy_degree'] := Valuation(eval(diff(curve, y), {x = Place[x], y = Place[y]}), t, alpha); # Update branching number L := L + abs(degree(Place[x],t)); # Downstream methods want ONE `y`-expansion (rest obtained by conjugation) Place[y] := [allvalues(Place[y], 'implicit')][1]; # If the expansion for `y` is zero, then replace with `next_p` if not(has(Place[y], t)) then Place := Next_P(curve, x, y, t, Place); fi; ans := [x = Place[x], y = Place[y], 'valuation' = Place['valuation'], 'dfdy_degree' = Place['dfdy_degree'] , 'parameter' = Place['parameter']], ans; od: # Check to see if there are the correct number of branches if L = degree(curve, y) then return [ans] else error(`AllBranches expects to find`, degree(curve,y), `branches, but found`, L); fi: end: # end AllBranches ----------------------------------------------------------------------------- # valuation ----------------------------------------------------------------------------------- # find the valuation (quick,safe?) of an expresion in `t` # try float, then `flag` and check with `evala` if needed Valuation:=proc(expr,t::name) local i; if expr=0 then # zero is as zero does return -infinity; fi; # from the low degree to the high degree for i from ldegree(expr,t) to degree(expr,t) do if not evala(coeff(expr,t,i))=0 then # it is not zero, return the degree return i; fi; od: # if there is not one non-zero, then well, it is zero return -infinity; end proc: # end Valuation ------------------------------------------------------------------------------- # Next_P -------------------------------------------------------------------------------------- # get next Puiseux series Next_P := proc(curve, x::name, y::name, t::name, p::{list,table}) local j, alpha, P, balance, y_deg, found, Alpha; # Mint supression alpha := alpha; # Convert to table P := `if`(hastype(p, table), copy(p), convert(p, table)); # Check for a valuation, attatch if missing if not has({indices(P)}, 'valuation') then P['valuation'] := Valuation(eval(curve ,{x = P[x] , y = P[y] + alpha*t^(1 + `if`(P[y] = 0, 0, degree(P[y], t)))}), t, alpha) - 1; fi; # find degree of `y(t)` y_deg := `if`(P[y] = 0, 0, degree(P[y], t)); # initialize loop: found sets to TRUE when a new expansion is found found,j:=false,1; while not found do # Compute the coefficient of the next power of 't' balance := coeff(eval(curve,{x = P[x], y = P[y] + alpha*t^(j + y_deg)}), t, P['valuation'] + j); if testeq(balance) or balance = 0 or degree(balance, alpha) = 0 then # No term here, go on. elif degree(balance,alpha)=1 then # Proceed. _EnvExplicit:=false added Dec 15 _EnvExplicit := false; Alpha := evala(solve(balance, alpha)); if Alpha <> 0 then if not hastype(p, table) then return [x = P[x], y = P[y] + collect(Alpha*t^(j + y_deg), t) ,'valuation'=P['valuation'] + j, 'dfdy_degree' = P['dfdy_degree'] , 'parameter' = P['parameter']]; elif hastype(p, table) then P[y] := collect(P[y] + Alpha*t^(j + y_deg), t); P['valuation'] := P['valuation'] + j; return P; fi; fi; else # complain about non-linear alpha error("next_p found more than one branch."); fi; # increment on the power of t j := j + 1; od: end proc: # end Next_P ---------------------------------------------------------------------------------- # @nearest ------------------------------------------------------------------------------------ # procedure to find (finite) point nearest to (finite) z0 in a lisjt = [a, b, c, ...] nearest := proc(lisjt::list, z0) local i, dumb, temp, dist; # Mint suppresion i, temp := i, temp; dist := seq(`if`(has(i, infinity), infinity, abs(evalf(i-z0))), i in lisjt); dumb := member(min(dist), [dist], `temp`); if dumb then lisjt[temp], temp else error("nearest failed."); fi; end: ##################################################################### # @ProblemPoints # Compute the "problem points" of the `curve` ProblemPoints := proc(curve, x :: name, y :: name) local F, PP; _EnvExplicit := false; F := StandardCurve(curve, x, y); # Problem points, (discriminant points and leading coefficient zeros) PP := map(allvalues, {solve(resultant(F, diff(F, y), y), x)}, 'implicit'); PP := ListTools:-MakeUnique(convert(`union`(PP,{solve(coeff(F, y, degree(F, y)), 'implicit')}), list)); end: ##################################################################### # @StandardCurve StandardCurve := proc(curve, x::name, y::name) # The curve in standard algcurves form. if indets(curve) <> {x, y} then error "Polynomial defining algebraic curve should be in %1 and %2.\ The object(s) %3 are scanning as variable", x, y, indets(curve) minus {x, y}; fi; return collect(primpart(curve, {x, y}), {x, y}, 'distributed', Normalizer); end: end module: