##################################################################### # standard `algcurves` macros 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'] ): ##################################################################### # Compute the Abel transform of divisor given as [multiplicity, point] # Abel map for points AbelMap := module() #options package; export ModPeriodLattice # ,Divisor ,CanonicalDivisor ,DivisorAbelMap ; local ModuleApply ,input_check ,departure ,IntegratePuiseuxVector ,Refine_diff ,ResolveBranch ,Resolve ,AllBranches ,Valuation ,Next_P ,trip ,abelpart ,FNext_p ,evalf_sum ,how_big_can_it_get ,bound_rootof ,swirl ,normalizer ,navigate ,pointBypoint ,short_path ,short ,connect_sheet ,leftmost ,whereis ,LIMIT ,is_resolved_abel ,nearest ,Sort_arg ,fix_rootof_table ,termlist ,on_0_2Pi ,size_log10 ,powers_raised ,powers_multiplied ,lowest_new_term ,FRAC ,AbelDiff ,AllCD ,SmallCD ,De_Minting ,plus # ,CommonDenominator # ,DivisorPlace # ,ProblemPoints ; ##################################################################### # @De_Minting # proc to stop mint complaint De_Minting:= proc() De_Minting(); ModPeriodLattice(); # Divisor(); CanonicalDivisor(); AllCD(); SmallCD(); end: ##################################################################### # @ ModuleApply ModuleApply := proc(curve, x :: name, y :: name, inP, inQ, Acc :: integer) local i, P, Q, g, escape_P, escape_Q, pilgrimage_P, pilgrimage_Q, twist, F, acc, V; options cache; # Quick check that `input_check` does not cover. if evalb(indets(curve) <> {x, y}) then error("polynomial should be in two variables."); fi; acc := `if`(nargs = 6, Acc, Digits); F := collect(primpart(curve, {x, y}), {x, y}, 'distributed', Normalizer); g := genus(F, x, y); Digits := max(Digits, acc) + 1; # Check that P <> Q, which will happen sometimes when calling from `DivisorAbelMap` if evalb(inP = inQ) then return [seq(0.0, i = 1..g)]; fi; # Check that input is in required form, is on curve etc. P, Q := input_check(curve, x, y, inP), input_check(curve, x, y, inQ); # `departure` returns integral, y-value, discpoint, x-value, t-value; escape_P := departure(F, x, y, P, acc); escape_Q := departure(F, x, y, Q, acc); # `trip` is the trip on one sheet. pilgrimage_P := trip(F, x, y, escape_P[2..5], acc); pilgrimage_Q := trip(F, x, y, escape_Q[2..5], acc); # `swirl` computes the integral for travel from sheet(P) to sheet(Q) twist := swirl(F, x, y, pilgrimage_P[2], pilgrimage_Q[2], acc); V := [seq(escape_P[1][i] - escape_Q[1][i] + pilgrimage_P[1][i] - pilgrimage_Q[1][i] + twist[i] , i = 1..g)]; if has({args}, {'reduce', 'Reduce', 'Jacobian'}) then return ModPeriodLattice(V, periodmatrix(F, x, y, 'Riemann')); elif has({args}, {'fraction', 'Fraction'}) then return ModPeriodLattice(V, periodmatrix(F, x, y, 'Riemann'), 'fraction'); else return V; fi; end: ##################################################################### # @DivisorAbelMap # Compute the Abel transform of divisor given as [multiplicity,point] DivisorAbelMap:=proc(F, x :: name, y :: name, P0, Divisor, Acc :: integer) local i, g, A, spot, acc; acc := `if`(nargs = 6, Acc, Digits); g := genus(F, x, y); A := [seq(0, i = 1..g)]; for spot in Divisor do A:=plus(A, spot[1]*AbelMap(F, x, y, P0, spot[2], acc)); od: A 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, default, 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; default := "`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, default); fi; if hastype(inP, float) then error("Floats not supported by AbelMap.", default); fi; P := `if`(hastype(inP, 'table'), copy(inP), convert(inP, 'table')); # 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]", default); 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'], default); 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, default); 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, default); 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 # if evala(res[j]) = 0 then # match[i] = true; # break; # fi; # 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, default); else error("The expansion %1 matches no branch. %2", inP, default); fi; end: ##################################################################### # @departure # A function of (curve, x::name, y::name, point::list, acc) # point is a list[x = f(t), y = g(t)], 'valuation' need not be present. # Returns: # V = the evaluated integral # planet = the nearest discpoint # intecept = the `x`-value where intersection with DvH path happens # y_star = the `y`-value over intercept; F(intercept, y_star) = 0 # theta = the argument on_0_2Pi of the intersection WRT planet (unnecessary) departure := proc(curve, x::name, y::name, point, acc) local N, P, Planet, Roots, V, c, diffs, discpoints, g, h, intercept, n, offset, t ,ord, pad, planet, t1, t2, t_max, theta, x0, y0, y_star, i, j, tau, Dig, df; options cache; tau := tau; # avoid mint noise # 'ZERO' is the flag for basepoint. if has(point, 'ZERO') then return [seq(0.0, i = 1..genus(curve, x, y))], seq(false, i = 1..4); fi; # Pick up the homology, differential normalizer, genus of `curve` and a conversion # Ask for extra digits when normalizing, can get costly in a hurry. # Save current value of Digits for `nearness` tests. Dig := Digits; Digits := max(Digits, acc) + 2; h, g := homology(curve, x, y, `give paths`), genus(curve, x, y); c, P := normalizer(curve, x, y, g), convert(point, table); t := P['parameter']; # Return zero if escaping the 'basepoint' (Flaged as a call from `ZERO`) # in Riemann constant vector calculation if has(point, 'ZERO') then return [seq(0.0, i = 1..genus(curve, x, y))], seq(false, i = 1..4); fi; # Obtain the disc...points and their radii from `homology` discpoints := map(op, [indices(h[4])]); Roots := fix_rootof_table(h[5], discpoints); # Find the point we are leaving and particulars about it. x0, n := evalf(limit(P[x], t = 0)), degree(P[x], t); y0 := evalf(limit(P[y], t = 0)); if not(has(P[y], t)) then # Constant `y`-value induces problems, so fix. P := Next_P(curve, x, y, t, P); fi; if n > 0 then # User point is finite. Find nearest (finite) discriminant point planet := nearest(map(op, [discpoints]), x0)[1]; # Path starts at `t` = 0, now find end. theta := on_0_2Pi(argument( `if`( # If the desired x-value is very close abs(x0 - planet) < 10^(1 - Dig), # Then use the end of the path (left/right) `if`( # If the path is not-empty nops(h[3][planet])>0, # Use left/right end of path eval(h[3][planet][-1], h[3]['parameter'] = 1), # Or if path empty then this is near basepoint and path ends on left planet - 10^(2 - Dig) ), # Otherwise use the x-value itself x0 ) - planet)); offset := h[4][planet]*exp(I*theta); intercept := planet + offset; # Find the `y0` that AC will expect. Idiot-proofing should check. y0, N := nearest([`algcurves/fsolve`(eval(curve, x = x0), y, complex)], y0); if abs(planet-x0) >= h[4][planet] then # The point is outside of safety circle so the point is outside the # safety circle of the nearest discriminant point then the sheets are # separarted and `y0` should be easy to pick up pad := x0*(1 - tau) + intercept*tau; # `differentials` without denominators diffs := differentials(curve, x, y, 'skip_dx'); # Compute integral! V := [seq(Partintegral(curve, diffs[i], x, y, pad, tau, y0, acc)[1], i = 1..g)]; # Catch that last `y`-value. y_star := AC(curve, x, y, pad, tau, acc)[-1][2][N]; # Normalize. V := [seq(add(c[i, j]*V[j], j = 1..g), i = 1..g)]; return V, planet, intercept, y_star, theta; else # Integrate out from the nearby discriminant point. # Get branches `resolved` over `x0`, and refined until resolved over intercept if not(abs(planet - x0) < 10^(1 - Dig)) then # `x0` is close but not VERY CLOSE to `planet.` Planet := Resolve(curve, x, y, t, Roots[planet], x0, y0, intercept); else # `x0` is VERY CLOSE to `planet. Use the provided PE, making sure resolved over intercept. Planet := ResolveBranch(curve, x, y, t, P, intercept); fi; # Have the `Planet`, now get the `t`-values. # The `evalf`s in the next line are needed to keep Digits same in both # `Planet[x]` and `x0` t1 := Sort_arg([`algcurves/fsolve` (evalf[Digits - 1](Planet[x]) - evalf[Digits - 1](x0), t, complex)])[1]; t2 := Sort_arg([`algcurves/fsolve` (evalf[Digits - 1](Planet[x]) - evalf[Digits - 1](intercept), t, complex)])[1]; # Get the differentials. diffs := differentials(curve, x, y, 'skip_dx'); df := diff(curve, y); # Check if the error (at least for `y(t2)`) is acceptable, refine if not. Planet, ord := Refine_diff(curve, x, y, t, Planet, t2, diffs, acc + 1); # Compute. V := IntegratePuiseuxVector (x, y, t, Planet, t1, t2, [seq(i/diff(curve, y), i in diffs)], ord, acc); V := [seq(add(c[i, j]*V[j], j = 1..g), i = 1..g)]; y_star := eval(Planet[y], t = t2); return V, planet, intercept, y_star, theta; fi; elif n < 0 then # Go to leftpoint. Most likely will include an analytic continuation (done by `trip`). planet := h[1]['basepoint']; t_max := abs([`algcurves/fsolve` (evalf[Digits - 1](P[x]) - evalf[Digits - 1] (max(seq(abs(i), i in discpoints))), t, complex)][1]); intercept := -abs(eval(P[x], t = t_max/2)); theta := on_0_2Pi(argument(intercept - planet)); # If user handed off `infinity` then it is not VERY CLOSE (it is ON) Planet := ResolveBranch(curve, x, y, t, P, intercept); # Have the `Planet`, now get the `t`-values. t1, t2 := 0, Sort_arg ([`algcurves/fsolve`(evalf[Digits - 1](Planet[x]) - evalf[Digits - 1](intercept), t, complex)])[1]; diffs := differentials(curve, x, y, 'skip_dx'); df := diff(curve, y); # Check if the error is acceptable, refine if not. Planet, ord := Refine_diff(curve, x, y, t, Planet, t2, diffs, acc + 1); # Compute. V := IntegratePuiseuxVector(x, y, t, Planet, t1, t2, [seq(i/df, i in diffs)], ord, acc); V := [seq(add(c[i, j]*V[j], j = 1..g), i = 1..g)]; y_star := eval(Planet[y], t = t2); return V, planet, intercept, y_star, theta; else error("the expansion about %1 is constant in %2", x0, t); fi; end: ##################################################################### # Ah...well, integrate then I suppose. IntegratePuiseuxVector := proc(x, y, t, planet, t1::constant, t2::constant , diffs::list, ord::integer, acc::integer) local diff_temp, dx, van, i, down, P, extra, up, answer, temp_ord; P := `if`(hastype(planet, 'table'), copy(planet), convert(planet, 'table')); van := P['dfdy_degree']; dx := diff(P[x], t); # If the denominator has 1/t terms it is possible to get ord = 0. No # good for `series` extra := `if`(van < 0, abs(van), 0) + 1; # Make the functions, all denominators are the same, take the first. down := convert(series(eval(evalf[acc+1](denom(diffs[1])), {y = P[y], x = P[x]}), t = 0, ord + extra), polynom); # Any term in the denominator (D_y F) with t-degree less than `van` # must be removed. down := select(z -> evalb(degree(z, t) >= van), down); # Turn into series for denominator down := series(1/down, t = 0, ord + extra); # Compute the numerator up := map(z -> convert(series(eval(z*dx,{y = P[y], x = P[x]}), t = 0, ord + extra) , polynom), map(numer, diffs)); # CHANGE TEH FOLLOWINHG LINES TO... # deal with the fact that is the series for down starts at 'ord', a # bunch of zeroes are returned # The order `ord` may be too low to convince Maple to compute any terms # hence "max(ord, ldegree(...))" temp_ord := [seq(max(ord, -van, ldegree(convert(i, polynom), t )) + 1, i in up)]; # Multiply through by numerators and dx diff_temp := [seq(convert(series(up[i]*down, t, temp_ord[i]), polynom), i = 1..nops(up))]; # Check that there is a series there if has(diff_temp, 0) then error ("Series for differential is zero to order of %1-series near %2=%3",y,x,limit(P[x],t=0)); fi; # Make sure series has no negative powers. if has (diff_temp,t^(-anything)) then error ("Series for differential is not holomorphic near %1=%2",x,limit(P[x],t=0)); fi; answer := [seq(int(diff_temp[i], t = t1..t2) , i = 1..nops(diff_temp))]; answer; end: ##################################################################### # Compute ever more terms in resolved `y`-expansion until last term in # worst differential is tolerable # @Refine_diff Refine_diff := proc(curve, x, y, t, planet, t_max, diffs, acc::integer) local dumb, P, df, dx, g, ord, counter, i, new_diffs, old_diffs, last_one_small, CM, Ex, Why, sigma, denom_ord, temp; P := planet; dx := evalf[acc + 1](eval(diff(P[x], t), t = t_max)); df := diff(curve, y); # `last_y` only used in here, so here is where it gets updated #if not has({indices(P)}, 'last_y') then # P['last_y'] := `if`(P[y] = 0, 0, coeff(P[y], t, degree(P[y]) ) * t ^ degree(P[y],t)); #fi; g := nops(diffs); # Iterate. 100 is MAGIC. # !!! These are convergent series. The refinement should be fast. I warn # the user if the refinement is not fast. I have seen the refinement # go past ten. # Leave it as 100 for now and I will think about this more. # compute the powers of the terms dumb,CM,Ex,Why := termlist(df,x,y,t,P); # mint noise avoidence dumb := dumb; # I want small difference twice in a row for counter to 100 do # $$$$$$$$$$$$$$$$$$ TIME SINK $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ temp := evalf[acc + 2](eval(eval(df, {x = P[x], y = P[y]}), t = t_max)); new_diffs := [seq(evalf[acc + 2] (eval(eval(diffs[i] * dx / temp, {x = P[x], y = P[y]}), t = t_max)), i=1..g)]; # $$$$$$$$$$$$$$$$$$ TIME SINK $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ # Compare with the old version of the differential at the end of the path if counter > 2 and degree(P[y],t) > 1 and max(seq( abs(new_diffs[i] - old_diffs[i]), i = 1..g)) < 10^(-acc) then # Assume the series for differential is valid to same degree(order) as y-series # Make sure `small` at least twice in a row. if last_one_small then # Compute the lowest new term added to denominator by refinement in y denom_ord := lowest_new_term(CM, Ex, Why, degree(df,x), degree(df,y)); # Find the the most restrictive numerator, and use that to determine the order sigma := min(seq(Ex[1]*degree(i,x) + Why[1]*degree(i,y) , i in map(numer,diffs))) +degree(dx,t) - P['dfdy_degree']; ord := denom_ord - P['dfdy_degree'] + sigma; # `ord` is now the valid order of the series for most restrictive # differential. This must be modified (in IntegratePuiseuxVector) to # accomadate for negative powers . if ord < 1 then error("The order of a differential is zero in, %1", t); fi; return P, ord; else # It is good, but get the next Puisex expansion to get two small # differences in a row. Update the powers lisjt. last_one_small := true; P := FNext_p(curve, x, y, t, P, acc); Why := [op(Why),degree(P[y],t)]; fi; else # Not within tolerence. Get new expansion. Update the powers lisjt last_one_small := false; P := FNext_p(curve, x, y, t, P, acc); Why := [op(Why),degree(P[y],t)]; fi; old_diffs:=new_diffs; od: error("Did not make it to requested accuracy when refining Puiseux expansion."); end: ##################################################################### # @ResolveBranch # Refine until resolved over `intercept` ResolveBranch := proc(curve, x, y, t, P, intercept) local Planet, Y, i, initial_dist, t_star, j; Planet := `if`(hastype(P, 'table'), copy(P), convert(P, 'table')); # Bump up `Digits` locally. Digits := Digits+2; # Find the values of `t` such that P[x](t) = x1, sort by arg(on_0_2Pi) t_star := Sort_arg([`algcurves/fsolve`(evalf[Digits-1](Planet[x])-evalf[Digits-1](intercept), t, complex)])[1]; # What is the actual floating point fibre over `x1`? Y := [`algcurves/fsolve`(eval(curve, x = intercept), y, complex)]; initial_dist := min(seq(abs(Y[i]-evalf(eval(Planet[y], t = t_star))), i = 1..nops(Y))); for i to 1000 do if is_resolved_abel(Y, eval(Planet[y], t = t_star)) then # It is resolved, away we go. return Planet else # Refine Planet := FNext_p(curve, x, y, t, Planet, Digits); fi; # Warn if over-iterating. There should be a solution suggested (will increasing Digits help?) if `mod`(i, 20) = 0 then WARNING("%1 is working quite hard to find a match at %2", procname, intercept); elif `mod`(i, 100) = 0 then WARNING("%1 is working VERY hard to find a match at %2. Suggest canceling computation", procname, intercept); elif i>3 and min(seq(abs(Y[j]-eval(Planet[y], t = t_star)), j = 1..nops(Y))) > initial_dist then # The expansion calculation is blowing up. error("The calculation of the expansion around %1 has failed.", eval(P[x], t = 0)); fi; od: end: ##################################################################### # @Resolve # Compute the branches from (algebraic)`planet` that are resolved over # `x1` and matches (algebraic, float)`y`-value `y1` # then refine until resolved over `x2` (optional) Resolve := proc(curve, x, y, t, planet, x1, y1, x2) local Branch, Planet, Y, branch_number_plus, fibre, fibre_index, i, j, k, t_star, y2; Planet := map(convert, AllBranches(curve, x, planet, y, 0, t), table); # Find the values of `t` such that P[x](t) = x1, sort by arg(on_0_2Pi) t_star := [seq(Sort_arg([ `algcurves/fsolve`(evalf[2*Digits](Planet[i][x])-evalf[2*Digits](x1), t, complex)])[1] , i = 1..nops(Planet))]; # What is the actual floating point fibre over `x1`? Y := [`algcurves/fsolve`(eval(curve, x = x1), y, complex)]; # From this list pick the match to `y1` (assumed to be resolved). # Calculate the `y`-fibre over `x1` according to the PE from `Planet` fibre := [seq(`if`(degree(Planet[i][x], t) = 1, evalf(eval(Planet[i][y], t = t_star[i])), seq(evalf(eval(Planet[i][y], evalf(t = t_star[i]*exp(2*Pi*I*(j-1)/degree(Planet[i][x], t))))) , j = 1..degree(Planet[i][x], t))), i = 1..nops(t_star))]; # Check that the `fibre` according to PE is resolved, 1000 is MAGIC for i to 1000 do if is_resolved_abel(fibre, y1) then # It is resolved, away we go. Get the right branch (careful of conjugation.) fibre_index := [seq(`if`(degree(Planet[k][x], t) = 1 , [k, 0] , seq([k, (j-1)], j = 1..degree(Planet[k][x], t))) , k = 1..nops(t_star))]; # Get the branch number and angle offset and assign. branch_number_plus := fibre_index[nearest(fibre, y1)[2]]; Branch := Planet[branch_number_plus[1]]; Branch[y] := eval(Branch[y], t = t*exp(2*Pi*I*branch_number_plus[2]/degree(Branch[x], t))); if nargs = 8 then # Iterate until resolved over `x2`, first get the true-float-fibre Y := [`algcurves/fsolve`(eval(curve, x = x2), y, 'complex')]; t_star := Sort_arg([`algcurves/fsolve` (evalf[Digits-1](Branch[x]) - evalf[Digits-1](x2), t, 'complex')])[1]; for j to 1000 do y2 := eval(Branch[y], t = t_star); if is_resolved_abel(Y, y2) then return Branch else # Refine Branch := FNext_p(curve, x, y, t, Branch, Digits); fi; # Warn if over-iterating. There should be a solution suggested if `mod`(i, 20) = 0 then WARNING("%1 is working quite hard to find a match at %2 from %3" , procname, x2, Planet); elif `mod`(i, 100) = 0 then WARNING("%1 is working VERY hard to find a match at %2 from %3.\ Suggest canceling computation", procname, x2, Planet); fi; od: else return Branch fi; else # Iterate until resolved; get next expansion and calculate fibre Planet := [seq(FNext_p(curve, x, y, t, Planet[j], Digits), j = 1..nops(Planet))]; fibre := [seq(`if`(degree(Planet[k][x], t) = 1, eval(Planet[k][y], t = t_star[k]), seq(eval(Planet[k][y], evalf(t = t_star[k]*exp(2*Pi*I*(j - 1) /degree(Planet[k][x], t)))) , j = 1..degree(Planet[k][x], t))), k = 1..nops(t_star))]; fi; # Warn if over-iterating. There should be a solution suggested if `mod`(i, 20) = 0 then WARNING("%1 is working quite hard to find a match at %2 from %3" , procname, x1, Planet); elif `mod`(i, 100) = 0 then WARNING("%1 is working VERY hard to find a match at %2 from %3.\ Suggest canceling computation", procname, x1, Planet); fi; od: end: ##################################################################### # get all the representations over x = a # @AllBranches AllBranches := proc(curve, x::name, A, y::name, N::integer, t::name) local L, Place, Places, a, ans, Ans, places, alpha, flag, extra # ,VCheck, i ; options cache; # `alpha` is local. it is used as dummy, solved for and tossed. alpha := alpha; # `extra` is incremented by 2. Used to make sure series all contain # `t` when `A` =0 extra := 0; # `a` is the users point in my format a := `if`(hastype(A, 'radical'), convert(A, 'RootOf'), A); # Get the Puiseux expansions from van Hoeij in form I want. # Make sure `t` shows up in every series by iteration. flag := true; while flag do extra := extra + 1; places := map(collect, map(expand, convert(map(op, {puiseux(curve, x = a, y, N + extra, t)}), list)) , t); flag := has( map(z->has(rhs(z),t), ListTools:-Flatten(places)), false); od: # Conjugate by `allvalues` places := ListTools:-MakeUnique(map(allvalues, places, 'implicit')); Places := map(convert, places, 'table'); # Initialize some stuff L := 0; ans := NULL; Ans := []; # Go. for Place in Places do # Attach the 'valuation's for use in finding subsequent series. Place['parameter'] := t; 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; # Need the degree of vanishing of 'dfdy' to avoid singularities from small error at leading order in denom Place['dfdy_degree'] := Valuation(eval(diff(curve, y), {x = Place[x], y = Place[y]}), t, alpha); L := L + abs(degree(Place[x], t)); # The methods that follow want ONE PARTICULAR `y`-expansion Place[y] := [allvalues(Place[y], implicit)][1]; # an attempt to compute the valuations faster (as per vH suggestion) # VCheck := add(`if`( degree(Place[x], t) = 1, 1, abs(degree(Place[x], t) - 1))*Valuation(Place[y] - Places[i][y] # + alpha*t^(`if`(Place[y] = 0, 0, degree(Place[y], t)) + 1), t) # , i = 1..nops(Places)) # + Valuation(eval(coeff(curve, y, degree(curve, y)), x = Place[x]), t); # VCheck2 := add(`if`(Places[i] = Place # , `if`(Place[y] = 0, 0, degree(Place[y], t)) + (abs(degree(Place[x], t)) - 1)*ldegree(Place[y], t) # , abs(degree(Place[x], t))*Valuation(Place[y] - Places[i][y], t) ) # , i = 1..nops(Places)) + Valuation(eval(coeff(curve, y, degree(curve, y)), x = Place[x]), t); # If the expansion for `y` is constant, then replace with `Next_p` if not(has(Place[y], t)) then Place := Next_P(curve, x, y, t, Place); fi; if not has({args}, 'table') then ans := [x = Place[x], y = Place[y], 'valuation' = Place['valuation'] , 'dfdy_degree' = Place['dfdy_degree'], 'parameter' = t] , ans; else Ans := [copy(Place), op(Ans)]; fi; od: # Check to see if there are the correct number of branches if L = degree(curve, y) then return `if`(not has({args}, 'table'), [ans], Ans) else error("AllBranches expects to find %1 branches, but found %2", degree(curve, y), L); fi: end: ##################################################################### # @Valuation # find the valuation of an expresion in `t` # (try float, then `flag` and check with `evala` if needed?) Valuation := proc(expr, t::name) local i; if not has(indets(expr), t) or expr = 0 then return 0; 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, and the degree is... -infinity; end: ##################################################################### # @Next_P # get next Puiseux series Next_P := proc(curve, x::name, y::name, t::name, place::{list, table}) local Alpha, P, balance, found, j, y_deg, alpha; # Avoid mint noise; alpha used as name/flag alpha := alpha; # avoid mint noise # convert to table P := `if`(hastype(place, 'table'), copy(place), convert(place, '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; # Need the degree of vanishing of `dfdy` to avoid singularities from small error at leading order in denom if not has({indices(P)}, 'dfdy_degree') then P['dfdy_degree'] := Valuation(eval(diff(curve, y), {x = P[x], y = P[y]}), t, alpha); fi; if not has({indices(P)}, 'made_float') then P['made_float'] := false; fi; # Find degree of `y(t)` y_deg := `if`(P[y] = 0, 0, degree(P[y], t)); # Initialize loop found, j := false, 1; while not found do balance := coeff(eval(curve, {x = P[x], y = P[y] + alpha*t^(j + y_deg)}) , t, P['valuation'] + j); if 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(place, 'table') then return [x = P[x], y = P[y] + Alpha*t^(j + y_deg) , 'valuation' = P['valuation'] + j , 'dfdy_degree' = P['dfdy_degree'], 'parameter' = t]; else 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. Are your expansions separated?"); fi; # Increment j := j + 1; od: end: ##################################################################### # @trip # this does the "walk"(integration) on one sheet from `x_star` to 'basepoint' # `x_star` = the place in C where puiseux meets the path # `y_star` = the correct cover over `x_star` # `t_star` = the parameter value (vis-a-vis the path, an angle in radians) # `planet` = the (not necsesarily) branchpoint about which the meeting takes place trip := proc(curve, x, y, Planet, x_star, y_star, t_star, acc) local planet, base, c, diffs, fibre, g, j, m, pad, paths, r, radii, reentry , sheets, side, sweep, t, temp, y_base, i; options cache; if Planet = false then # This means escaping from base itself, so sheet is always number one. return [seq(0.0, i = 1..genus(curve, x, y))] , monodromy(curve, x, y, `give paths`)[2][1]; fi; m := monodromy(curve, x, y, `give paths`); base, sheets, paths, radii := m[1], m[2], m[5], m[6]; # Select nearest discpoint to planet to avoid accuracy mismatch problems planet := nearest(map(op, [indices(homology(curve, x, y, `give paths`)[4])]), Planet)[1]; # Pluck out the parameter name and the genus. t, g := paths['parameter'], genus(curve, x, y); # Construct the path in the complex x plane from basepoint to x_star if planet = base then # go nowhere to get to `base.` The path is from `base` to the intersection. pad := [base*(1-t)+t*x_star]; else r := radii[planet]; # Check if path from homology ends to left or right of `planet`. temp := paths[planet]; if nops(temp) = 0 or abs(eval(temp[-1], t = 1) - (planet - r))<10^(3 - Digits) then # left point of `planet` side := 1; elif abs(eval(temp[-1], t = 1) - (planet + r))<10^(3 - Digits) then # right point of `planet`' side := 0; else error("the path to %1 given by algcurves[homology] is not standard.", planet); fi; # `sweep` is the path from the end of `monodromy` path to the # intersection of avoidence circle and ray leaving discriminant # point to `x_star, ` where the Puiseux expansion quit sweep := evalf(planet+r*exp(I*(side*Pi*(1-t)+t*t_star))); # `reentry` is path inwards to `x_star` if necesary reentry := evalf((planet+r*exp(I*t_star))*(1-t)+t*x_star); pad := [op(paths[planet]), `if`(abs(eval(sweep, t = 0)-eval(sweep, t = 1))>10^(1-Digits), sweep, NULL), `if`(abs(eval(reentry, t = 0)-eval(reentry, t = 1))>10^(1-Digits), reentry, NULL) ]; fi; # find the starting sheet fibre := sheets; for j in pad do fibre := Propagate(AC(curve, x, y, j, t, Digits), fibre); od; y_base := sheets[whereis(fibre, y_star)]; # normalizing coefficients for the differentials c := normalizer(curve, x, y, g); diffs := differentials(curve, x, y, 'skip_dx'); # The (-1) below is due to fact that `abelpart` integrates "out" Digits := max(Digits, acc) + 1; [seq(add((-1)*c[i, j]*abelpart(curve, diffs[j], x, y, pad, y_base, t), j = 1..g), i = 1..g)], y_base end: ##################################################################### # Integrate a non-normalized differential along a path that starts on a given sheet s # used as integration engine by `trip` abelpart := proc(curve, diff, x, y, pad, s, t) local integ, y0, i, parint; integ := 0; y0 := s; for i to nops(pad) do parint := Partintegral(curve, diff, x, y, pad[i], t, y0, Digits-2); integ := integ+parint[1]; y0 := parint[2]; od; integ end: ##################################################################### # Get next Puiseux series floating point # @FNext_p FNext_p := proc(curve, x::name, y::name, t::name, place::{list, table}, acc) local Alpha, P, balance, found, j, y_deg, alpha; # alpha is used as a name/flag; supress mint complaint alpha := alpha; P := `if`(type(place, 'table'), copy(place), convert(place, 'table')); # Check for valuation, error if absent if not has({indices(P)}, 'valuation') then error("Finding valuation in floats is unstable.") fi; # If 'made_float' is not assigned then this point has not been made float if not has({indices(P)}, 'made_float') then P['made_float'] := false; fi; # See if there are non-float coeffs, and float them if there are # Find degree of `y(t)` y_deg := `if`(P[y] = 0, 0, degree(P[y], t)); # Initialize loop found, j := false, 1; while not found do # $$$$$$$$$$$$$$$$$$$$$$$$$ TIME SINK $$$$$$$$$$$$$$$$$$$$$$$$$ balance := evalf_sum(coeff(eval(curve, {x = P[x], y = P[y]+alpha*t^(j+y_deg)}), t, P['valuation']+j), alpha, acc + 1); # $$$$$$$$$$$$$$$$$$$$$$$$$ TIME SINK $$$$$$$$$$$$$$$$$$$$$$$$$ Alpha := `if`(balance = 0 or not has(balance, alpha), 0 , `algcurves/fsolve`(balance, alpha)); if nops([Alpha])>1 then error("The branches should be separated."); elif Alpha <> 0 then # Return in proper format if not hastype(place, 'table') then return [x = `if`(P['made_float'], P[x], evalf_sum(P[x], t , acc + 1)) , y = `if`(P['made_float'], P[y], evalf_sum(P[y], t , acc + 1)) + Alpha*t^(j+y_deg) , 'valuation' = P['valuation']+j, 'parameter' = P['parameter'] ,'dfdy_degree' = P['dfdy_degree'], 'made_float'=true]; else P[x] := `if`(P['made_float'], P[x], evalf_sum(P[x], t , acc + 1)); P[y] := `if`(P['made_float'], P[y], evalf_sum(P[y], t , acc + 1)) + Alpha*t^(j+y_deg); P['made_float'] := true; P['valuation'] := P['valuation'] + j; return P; fi ; fi; # increment j := j + 1; od: end: ##################################################################### # @evalf_sum # `evalf` an `expression` in `t` (involving `RootOf`s etc.) to `accuracy` # must determine size of (un-evaluated) numbers being subtracted evalf_sum := proc(expression, t, accuracy) local i; if expression = 0 then return 0 fi; #if not type(ldegree(expression, t),integer) then # return 0 #fi; # The code commented out below is directly translated to one call to # `add` afterward. Left in for clarity. #l, h := ldegree(expression, t), degree(expression, t); #for i from l to h do # if coeff(expression, t, i) = 0 then # a[i] := 0.; # else # # guess how big coefficient "could" be and compute ceil(log10(of that)) # # log10_size := size_log10(how_big_can_it_get(coeff(expression, t, i))); # a[i] := evalf[accuracy+max(log10_size, 0)](coeff(expression, t, i)); # fi; #od: # # `if`(Digits Insheet2 then m := monodromy(curve, x, y, `give paths`); whereto := navigate(m); # Monodromy may have been computed in different places, # make sure floats match when they should sheet1 := nearest(m[2], Insheet1)[1]; sheet2 := nearest(m[2], Insheet2)[1]; # Pick up the directions from sheet to sheet gate := [whereto[sheet1, sheet2]]; brpoints := map(z -> op(1, z), m[3]); basepoint, sheets, permutations, discpaths, discradii := m[1], m[2], m[3], m[5], m[6]; t := discpaths['parameter']; integ := [seq(0, i = 1..g)]; # Normalizing coefficients for the differentials, ask for extra # accuracy as ana...c...tion asks for Digits - 3 Digits := max(Digits, acc); c := normalizer(curve, x, y, g); # This is the Digits setting that Partintegral wants Digits := max(10, Digits); # Compute the differentials and start/end sheet diffs := differentials(curve, x, y, 'skip_dx'); paris := sheet1; roubaix := sheet2; for k in gate do # Find the beginning and final values for y. indexper := whereis(brpoints, k[1]); i := 1; while not member(whereis(sheets, paris), permutations[indexper][2][i], 'h') do i := i+1 od; perpart := permutations[indexper][2][i]; if k[2] = true then if h = nops(perpart) then westvleteren := sheets[perpart[1]] else westvleteren := sheets[perpart[h+1]] fi; else if h = 1 then westvleteren := sheets[perpart[-1]] else westvleteren := sheets[perpart[h-1]] fi; fi; # And... integrate! if k[2] = true then integ := integ + [seq(add(c[kay, j]*Pathintegral(curve, diffs[j], x, y, discpaths, t, basepoint, paris, westvleteren, k[1], discradii[k[1]]) , j = 1..g), kay = 1..g)]; else integ := integ - [seq(add(c[kay, j]*Pathintegral(curve, diffs[j], x, y, discpaths, t, basepoint, westvleteren, paris, k[1], discradii[k[1]]) , j = 1..g), kay = 1..g)]; fi; paris := westvleteren; od; if roubaix <> westvleteren then error("Path was not constructed correctly during numerical integration of swirl.") fi; integ else [seq(0, i = 1..g)]; fi; end: ##################################################################### # @normalizer # weights for deferential differentials normalizer := proc(curve, x, y, g) options cache; linsolve(submatrix(periodmatrix(curve, x, y), 1..g, 1..g), array('identity', 1..g, 1..g)) end: ##################################################################### # @navigate # takes the return of Algcurves\monodromy and returns an array of directions to "navigate" a # R.S. the [i, j] component of the returned array is a list of points to encircle and in which sense # e.g. [a+I*b, true] means encircle the point a+I*b counterclockwise, # [c, false] means encircle the point c-clockwise navigate := proc(monod) local A, branchpoint, bpoint, perm, perms, permutations, sheets, numsheets, i, j, directions; sheets := monod[2]; numsheets := nops(sheets); permutations := monod[3]; A := array(1..numsheets, 1..numsheets, [seq([seq({}, j = 1..numsheets)], i = 1..numsheets)]); for branchpoint in permutations do bpoint := branchpoint[1]; perms := branchpoint[2]; for perm in perms do for i to nops(perm) do j := i+1; if evalb(j>nops(perm)) then j := 1 fi; A[perm[i], perm[j]] := {op(A[perm[i], perm[j]]), [bpoint, true ]}; A[perm[j], perm[i]] := {op(A[perm[j], perm[i]]), [bpoint, false]}; od: od: od: # Take the left-most point to minimize integration A := map(leftmost, A); # Now turn it into `directions` indexed by branchpoints for i to numsheets do for j to numsheets do if i <> j then directions[sheets[i], sheets[j]] := pointBypoint(short_path(short(i, j, A)), A); else directions[sheets[i], sheets[j]] := []; fi; od: od: directions; end: ##################################################################### # @pointBypoint # converts from list of sheets to point by point directions # where `A` is a table of "shortest" ways between sheets (given as complex numbers) pointBypoint := proc(path, A) local i; if op(path) = NULL then return NULL fi; seq(op(A[path[i], (path[i+1])]), i = 1..nops(path)-1) end: ##################################################################### # @short_path # takes a list of pairs of directly connected sheets, and # returns an ordered list of sheets short_path := proc(chain) local i, link, path; # next line deals with directly connected sheets if not hastype(op(chain)[1], list) then return chain fi; path := chain[-1][-1], chain[-1][-2]; for i from (nops(chain)-1) to 1 by -1 do for link in chain[i] do if link[2] = path[-1] then path := path, link[1] fi; od: od: ListTools:-Reverse([path]) end: ##################################################################### # @short # is a depth-first search, returns a "chain" of direct connections to get from start to finish # stops when `finish` is found, needs `short_path` to return actual `chain` short := proc(start, finish, A) local i, link, connect, found, chain; if start = finish then return [] fi; found := {start, op(connect_sheet(start, A))}; chain[1] := [seq([start, connect_sheet(start, A)[i]], i = 1..nops(connect_sheet(start, A)))]; for link in chain[1] do if link[2] = finish then return link fi; od: for i to coldim(A) do chain[i+1] := []; for link in chain[i] do for connect in connect_sheet(link[2], A) do if connect = finish then chain[i+1] := [link[2], finish];return convert(chain, list) fi; if not member(connect, found) then found := {connect} union found; chain[i+1] := [op(chain[i+1]), [link[2], connect]]; fi; od: od: od: error("Riemann Surface does not seem to be connected from sheet %1 to %2.", start, finish); end: ##################################################################### # @connect_sheet # returns a list of sheets directly(only one "elevator") connected to start connect_sheet := proc(start, A) local N, connected, LC, i; N := coldim(A); connected := {seq([i, leftmost(A[start, i])], i = 1..N)}; N := nops(connected); for i to N do LC[i] := `if`(nops(connected[i]) <> 1, connected[i][1], NULL); od: [seq(LC[i], i=1..N)]; end: ##################################################################### # @leftmost # takes a list of sense-consious points and returns the point with least, but positive, # imaginary part, of those points with the least real part # preference is to clockwise, or "true" # e.g. leftmost([[1, true], [-1+I, true], [-1-I, true], [-1+I, false]]); gives [-1+I, true] # doesn't work for non-sensed(ha!) leftmost := proc(points) local left, point; if nargs = 0 or nops(points) = 0 then return NULL fi; left := points[1]; for point in points do if point[1] = infinity then left := left elif Re(point[1]) < Re(left[1]) # it is to the left or Re(point[1]) = Re(left[1]) and abs(Im(point[1]))sign(Im(left[1])) # it has abs positive or Re(point[1]) = Re(left[1]) and Im(point[1]) = Im(left[1]) and point[2] = true # clockwise then left := point fi; od: [left] end: ##################################################################### # @whereis whereis := proc(lijst, element) local i, j, m; m := evalf(abs(lijst[1]-element)); if nops(lijst) = 1 then return(1); else j := 1; for i from 2 to nops(lijst) do if evalb(evalf(m>abs(lijst[i]-element))) then j := i; m := evalf(abs(lijst[i]-element)); fi; od; fi; if m <> 0 then # warn the user if element not close `enough` if not is_resolved_abel(lijst, element) then WARNING("Analytic continuation not done properly. More digits recommended") fi; fi; j; end: ##################################################################### # @LIMIT # , a procedure to determine behavior of point P as `t` -> 0 # point comes in as list or table LIMIT := proc(point::{list, table}, x::name, y::name, t::name) local P, Q; if not type(P, table) then P := convert(point, table); elif type(P, table) then P := point; else return FAIL; fi; Q[x], Q[y] := limit(P[x], t = 0), limit(P[y], t = 0); table([ x = `if`(type(evalf(Q[x]), finite), Q[x], infinity) , y = `if`(type(evalf(Q[y]), finite), Q[y], infinity)]) end: ##################################################################### # @is_resolved_abel # determine if cover `Y` can by picked from `Y_list` is_resolved_abel := proc(Y_list::list, Y::constant) local a1, d1, a2, d2; a1 := nearest(Y_list, Y)[1]; d1 := abs(Y-a1); a2 := nearest(convert(convert(Y_list, set) minus {a1}, list), Y)[1]; d2 := abs(Y-a2); evalb(evalf(d1 <= d2/10)) end: ##################################################################### # @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; i, temp := i, temp; # avoid mint noise # if has(z0, infinity) or not member(false, {seq(has(i, infinity), i in lisjt)}) then # error "can not determine distance to infinity."; # fi; 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 is confused."); fi; end: ##################################################################### # @ModPeriodLattice # a function to reduce vector to fundamental parralelogram ModPeriodLattice := proc(V, B) # !!! local g, EyE, Omega, L_V, v, i, j; g := LinearAlgebra:-RowDimension(B); EyE := 1.*Matrix(g, g, 'shape' = 'identity'); Omega := <, >; L_V := LinearAlgebra:-LinearSolve(Omega, <, >); v := map(FRAC, L_V); L_V := map(FRAC, L_V-v); if has([args], fraction) then return [seq(i, i in v)]; else return [seq(add([j, i]*v[i, 1], i = 1..2*g), j = 1..g)]; end if; end: ##################################################################### # @Sort_arg # Sort a list of complexs on argument (on_0_2Pi) Sort_arg := proc(lijst::list) sort(lijst, (a, b) -> evalb(on_0_2Pi(argument(a)) < on_0_2Pi(argument(b)))); end: ##################################################################### # @fix_rootof_table # RootOf table returned by `homology` does not index roots, I need indices # if a `disc_point` is not in table, assign it a converted fraction fix_rootof_table := proc(R::table, dpoints::list) options cache; local answer, all, float; # Compute allvalues of the `rootof`s all := convert(map(allvalues, map(op, {entries(R)}), 'implicit'), 'list'); # Get the floats according to the table for float in map(op, [indices(R)]) do answer[float] := nearest(all, float)[1]; od: # Infinity may be a discriminant point, but will not be in `rootof` table if member(infinity, dpoints) then answer[infinity] := infinity; fi; answer end: ##################################################################### # @termlist # list the terms in the series for dFdy termlist := proc(dF,x::name,y::name,t::name,Point) local dumb, CM, M, N, i, j, r, sigma, rho, lisjt; M,N := degree(dF,x),degree(dF,y); r := degree(Point[x],t); if ldegree(Point[x],t)<>r then rho := [ldegree(Point[x],t),r]; else rho := [r]; fi; dumb := coeffs(Point[y],t,`sigma`); # mint noise avoidence dumb := dumb; sigma := sort([seq(degree(i,t),i in [sigma])]); CM:= (i,j)->`if`( coeff(dF,x,i-1)<>0, coeff(coeff(dF,x,i-1),y,j-1), `if`( coeff(dF,y,j-1)<>0, coeff(coeff(dF,y,j-1),x,i-1), 0 ) ); lisjt :=sort (ListTools:-MakeUnique([seq(seq(`if`(CM(i,j)=0 , NULL , op(powers_multiplied(powers_raised(rho, i - 1, rho) , powers_raised(sigma, j - 1, sigma)))) ,i = 1..M + 1),j = 1..N + 1)])); return lisjt,CM,rho,sigma; end: ##################################################################### # @lowest_new_term # The lowest degree new term in the series for dFdy lowest_new_term := proc(CM,Ex::list,Why::list,M::integer,N::integer) local i, j, X, Y; # Double check, does not work if not sorted. X,Y := sort(Ex),sort(Why); min(seq(seq( `if`(CM(i,j)<>0, Y[-1] + Y[1]*(j-2) + X[1]*(i-1), NULL) ,i = 1..M + 1),j = 2..N + 1)) end: ##################################################################### # @powers_raised # raise power_list to exponent, # start with arguments (power_list,exponent,power_list) powers_raised := proc(power_list,exponent,storage) local i, j; # mint appeasal i, j := i, j; if exponent = 0 then return [0]; elif exponent = 1 then return sort(storage,`>`); else return powers_raised(power_list,exponent - 1 , ListTools:-MakeUnique([seq (seq(i + j, j in storage) , i in power_list)])); end: end: ##################################################################### # @powers_multiplied # Multiply two `powers` lists (a and b). # That is, the output is a list of the exponents that show up in # add(t^i, i in a) * add(t^j, j in b) powers_multiplied := proc(a :: list, b :: list) local i,j; # mint appeasal i, j := i, j; # Perform needed function, could be Quickie. Writing it this way # ensures that proc is only used on lists. ListTools:-MakeUnique([seq(seq(i + j, i in a), j in b)]); end: ##################################################################### # @CanonicalDivisor # Front-end to compute Canonical Divisors CanonicalDivisor:=proc(curve, x :: name, y :: name, t :: name) local ans, F; options cache; # The curve in standard algcurves form. if evalb(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; F := collect(primpart(curve, {x, y}), {x, y}, 'distributed', Normalizer); # Two choices, (i) small: (one Divisor, attempted first at infinity) # (ii) all: keep computing until 'g' Divisors are found if has([args], 'All') then ans := AllCD(F, x, y, t); if ans[1] = FAIL then error "Computing the divisor implied by differential with numerator %1 did not return %2 zeros.\ It is possible that %3 will return a valid cannonical divisor.", ans[2], 2*ans[3] -2 , procname(args[1..4]); fi; else ans := SmallCD(F, x, y, t); if ans[1] = FAIL then error "Attempting to compute a small cannonical divisor at infinity failed.\ The differential with numerator %1 gave %2 zeros, a deficiency of %3.\ Try %4.", ans[2], ans[3], ans[4], procname(args, 'All'); fi; fi; ans end: ##################################################################### # @SmallCD # An attempt to create a small CanonicalDivisor at infinity SmallCD:=proc(F, x :: name, y :: name, t :: name) local temp, g, Omega, omega, disc, XINF, target, DSL, ans, J, dfdy, d, Om, P, i, j; options cache; # Get a basis of (`g` = genus) differentials without denomiantor. dfdy := diff(F, y); Omega := AbelDiff(F, x, y, false); g := nops(Omega); # The points at (infinity,infinity) XINF := [seq(convert(i, table), i in AllBranches(F, x, infinity, y, 0, t))]; # There are 2g - 2 zeros of a Canonical Divisor target := 2*g - 2; # Creat a list of differentials, and the sum of the orders of their zeros at infinity 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 differential that gives the highest order divisor # (that is not too high [remove(w if w is greater than target)]) omega := op(select(z -> evalb(z[2] = max(op(remove(w -> evalb(w > target), map2(op, 2, temp))))) ,temp)); # May be more than one, choose one with low degree (first, as using AbelDiff) if nops([omega]) > 1 then omega := omega[1]; fi; # Create a divisor (possibly not complete) 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'] , [x = XINF[j][x], y = XINF[j][y], 'valuation' = XINF[j]['valuation'] , 'parameter' = XINF[j]['parameter'], 'dfdy_degree' = XINF[j]['dfdy_degree']]] , j = 1..nops(XINF))]; # Remove any for which the valuation is zero! ans := remove(z -> evalb(z[1] = 0), ans); if omega[2] = target then # There are the requisite number of zeros, so return it. return ans; else # Better check the discriminant points The new target is... target := target - omega[2]; # Calculate the discriminant points with respect to `y` _EnvExplicit := false; disc := map(allvalues, [solve(resultant(F, diff(F, y), y), x)], 'implicit'); # Add in the the roots of y^N. These are finite x-points with y = infinity disc := [op(disc), solve(coeff(F, y^degree(F, y)), x)]; # And add in zeros of the numerator disc := [op(disc), op(map(allvalues, [solve(resultant(F, omega[1], y), x)], 'implicit'))]; # Make sure Mark thinks these numbers are reduced, and Maple thinks # list is unique. disc := ListTools:-MakeUnique(map(evala, disc)); # Group "like" discriminant points together and arrange by the size of these groups. # D(isc)S(tring)L(ength): length of the dpoint as a string. # this is done as it is faster to compute with an integer than a # complicated RootOf 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. disc := sort(DSL, (a, b) -> evalb( a[1] < b[1] )); for d in disc do # Compute all the branches over this discriminant point P := [seq(convert(j, table), j in AllBranches(F, x, d[2], y, 0, t))]; # For each place, check for zeros/orders for J in P do Om := eval(omega[1]/dfdy, {x = J[x], y = J[y]})*diff(J[x], t); if limit(Om,t = 0) = 0 then # If you are here we have found a zero, calculate the order. ans := [op(ans) , [ldegree(eval(omega[1]*diff(J[x], t), {x = J[x], y = J[y]}), t) - J['dfdy_degree'] , [x = J[x], y = J[y] , 'valuation' = J['valuation'] , 'parameter' = J['parameter'] , 'dfdy_degree' = J['dfdy_degree']]]]; # Update the target, then check. target := target - ans[-1][1]; if target = 0 then return ans fi; fi; od: od: # Did not find what we wanted. return FAIL, omega[1], 2*g - 2 - target, target; fi; end: ##################################################################### # @AllCD # Compute `g` Canonical Divisors. AllCD:=proc(F, x :: name, y :: name, t :: name) local P, ans, ANS, g, omega, disc, zeros, test_places, count, i, j, k, Omega, found, deg_t # ,ti ; # Set the environment to not give radicals _EnvExplicit := false; # Genus g := genus(F, x, y); # Sort the differentials. No denominators, use dfdy_degree omega := ListTools:-Reverse(AbelDiff(F, x, y, false)); # Calculate the discriminant points with respect to `y` disc := ListTools:-MakeUnique(map(allvalues,[solve(resultant(F, diff(F, y), y), x)], 'implicit')); # Add in the point x=infinity, and all the roots of y^N disc := ListTools:-MakeUnique([op(disc), solve(coeff(F, y^degree(F,y)), x), infinity]); # Make sure Mark thinks these numbers are reduced. disc := map(evala, disc); for i from 1 to nops(omega) do # Initialize some stuff prior to search. 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, omega[i], y), x)], 'implicit')), op(disc)]); for j from 1 to nops(zeros) do # Must compute the branches over each of these points. 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. # ti := time(); deg_t := ldegree(evala(Omega), t) - P['dfdy_degree']; if deg_t > 0 then # THIS IS THE LINE THAT IS TAKING ALLLLLL THE TIME. found := found + 1; ans := ans, [deg_t, [x = P[x], y = P[y], 'valuation' = P['valuation'] , 'dfdy_degree' = P['dfdy_degree'], 'parameter' = t]]; # This is a somewhat silly check. How to fix? if found = 2*(g - 1) then break fi; elif deg_t < 0 then error("Holomorphic differential %1 appears to have a pole at %2", omega[i], evalf(P[x]) ); 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 # error("The degree of a computed divisor is not # 2*genus(curve) - 2."); return FAIL, omega[i], g; else ANS[i] := [ans]; fi; od: # forget(AllBranches); ANS end: ##################################################################### # @AbelDiff # Compute differentials "in order" AbelDiff:=proc(curve, x :: name, y :: name, flag :: boolean) local dFdy, diffs; options cache; 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 # Do nothing return diffs else dFdy:=diff(curve, y); return map(z -> z/dFdy, diffs); fi; end proc: ##################################################################### # Quickies on_0_2Pi := x -> evalf(`if`(evalf(x) >= 0, frac(x/2/Pi)*2*Pi, (1+frac(x/2/Pi))*2*Pi)): #--------------------------------- FRAC := x -> `if`(x >= 0, frac(x), frac(1+frac(x))): #--------------------------------- size_log10 := (x) -> ceil(log10(max(abs(Re(evalf(x))), abs(Im(evalf(x)))))): #--------------------------------- plus := proc(x, y) local i; `if`(nops(x) = nops(y) , [seq(x[i] + y[i], i = 1..nops(x))], FAIL) end: end module: # end of AbelMap # Some commands to set up 'AbelMap'in a repository. This allows for # profiling # march('create', "/home/mpatters/maple10/lib/AbelMap.mla"); # libname := "/home/mpatters/maple10/lib/AbelMap.mla", libname; # savelibname := "/home/mpatters/maple10/lib/AbelMap.mla": # savelib('AbleMap'):