# $Source: /u/maple/research/lib/algcurves/src/RCS/ratpar,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ # #--> ratpar: Compute a rational bijective morphism from a conic or a line # to a given curve f with genus 0. Here "rational" means: using no # algebraic extensions. # # # Input: f - A polynomial in 2 variables x and y, describing # an irreducible algebraic curve C. # x,y,s,t - variables # - optional argument `parametrize by line` to tell # ratpar to compute a bijection with a line, even if # that involves algebraic extensions # # Output: If ratpar is able to compute a rational parametrization (i.e. # a rational bijection with a line) then this will be the output. # The syntax is [X(s),Y(s)] where the variable s is given in the # input and X(s) and Y(s) are rational functions. # # If ratpar can not compute a rational parametrization, then # - if the option `parametrize by line` is given then a # parametrization over an algebraic extension of degree 2 is # given # - Otherwise the output is of the form [X(s,t),Y(s,t),F2(s,t)] # where [X(s,t),Y(s,t)] is a rational bijective morphism from a # conic C2 (defined by F2) to the curve C defined by f. # # Calling sequences: # ratpar(f,x,y,s,t); # ratpar(f,x,y,s,t,`parametrize by line`); # # Mark van Hoeij # The code will search for an odd place or point. By "odd place" or # "odd point" is meant: defined over an algebraic extension of odd degree. If # the code finds an odd place then it can compute a rational (i.e. using no # algebraic extensions) parametrization. # If no odd place is found, then (depending on whether the extra option is # given or not) either a parametrization (defined over a quadratic extension) # or a rational bijection with a conic is returned. # # Approach: # First a bijection to a conic C2 is computed. Then we try to find # a rational point on C2. This is done by searching an odd point on C2, # and for this we search for odd regular points on C. If this fails # then we compute the inverse morphism, i.e. we compute a bijection from C2 # to C. If now we have a place (not a regular point) of odd degree on C, then # we use this bijection to construct an odd point on C2 (note that this can be # done differently, one could do this without computing the inverse morphism). # If it succeeds we can still compute a rational parametrization, if it fails # then we'll either quit if `parametrize by line` is not given, and introduce # and algebraic extension of degree 2 if this option is given. # # Procedures: # odd_root: give (if possible) and odd root of a polynomial # odd_regpoint_C: try to find an odd regular point on C # odd_point_on_C2: try to find an odd point on C2 # rat_point_on_C2: try to find a rational point on C2 # search_rat_param: given C2, try to find a rational parametrization of C. # Each of these procedures uses the previous one. The result of each # of these procedures can be "failed". If search_rat_param fails then # ratpar has to figure out something else. In that it tries the procedure # odd_singularity_on_C: try to find a line on which C has an odd place # if this works, then compute the inverse morphism from C2 to C, and use # that to find an odd point on C2, that results in a rational point and # that gives a rational parametrization. If odd_singularity_on_C fails then use # an algebraic extension of degree 2. # # Other procedures: # inverse_of_g(f,f2,s,t,u,G,x,y): Compute the inverse of # the map [G[1],G[2],G[3]]: C -> C2 # genus1_charpol: this is section 4.1 in my ISSAC'95 paper. # express_in_p(f,x,y,s,v1,v2): if v1/v2 generates the function field # then this procedure computes a parametrization # L_inf, frac_integral_a: see "Rational Parametrizations of Algebraic # Curves", submitted to JSC. macro( ratpar=`algcurves/ratpar`, pyth2=`algcurves/pyth2`, odd_point_on_a_old=`algcurves/oddp_a`, inverse_of_g=`algcurves/inv_g`, genus1_charpol=`algcurves/g1_charpol`, odd_root=`algcurves/odd_root`, regpoint_from_sing=`algcurves/rp_from_s`, odd_regpoint_C=`algcurves/oddrp_C`, odd_point_on_C2=`algcurves/oddp_C`, rat_point_on_C2=`algcurves/ratp_C`, search_rat_param=`algcurves/s_param`, odd_singularity_on_C=`algcurves/odds_C`, parametrize_cube=`algcurves/param_cube`, parametrize_conic=`algcurves/param_conic`, express_in_p=`algcurves/expr_in_p`, express_x_in_p=`algcurves/expr_x_in_p`, compute_x_old=`algcurves/comp_x`, frac_integral_a=`algcurves/frac_int_a`, FFDIV=`algcurves/FFDIV`, L_inf=`algcurves/L_inf`, g_gcd=`algcurves/g_gcd` ): # This code uses the procedures in `algcurves/g_expand`, # `algcurves/puiseux`, `algcurves/integral_basis` and # `algcurves/singularities` macro( solve=`solve/linear`, homogeneous=`algcurves/homogeneous`, singularities=`algcurves/singularities`, find_points=`algcurves/find_points`, degree_ext=`algcurves/degree_ext`, integral_basis=`algcurves/integral_basis`, ext_to_coeffs=`algcurves/e_to_coeff`, `puiseux/technical_answer`=`algcurves/puiseux_te`, g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, g_factors=`algcurves/g_factors`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, z=`algcurves/z` ): # f is a polynomial in x and y # f must be irreducible over the algebraic closure of the constants # Output: a parametrization [X(s,t),Y(s,t),F(s,t)] where X(s,t),Y(s,t) are # the points on the curve f and where s,t are the points on the curve F. Here # F is a curve of degree 2. Or Output: a parametrization with by a line. ratpar:=proc(ff,x,y,s,t) global z; local u,j,co,f,n,i,res,b,v,a,den,d,d1,i0,bb,f2,vars,w,y_powers ,zero,rf; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; f:=collect(ff,{x,y},'distributed',Normalizer); n:=degree(f,{x,y}); if n=0 then error "degree should be at least 1" elif degree(f,x)*degree(f,y)=0 and n>1 then error "wrong input" elif degree(f,y)=1 then return subs(x=s,[x,subs(solve({f},{y}),y)]) elif degree(f,x)=1 then return subs(y=s,[subs(solve({f},{x}),x),y]) elif coeff(f,y,n)=0 then # Use one of these 2 linear transformations to avoid this case. if coeff(f,x,n)<>0 then res:=procname(f,y,x,args[4..nargs]); return [res[2],res[1],op(3..nops(res),res)] elif coeff(coeff(f,x,0),y,0)<>0 and member(`parametrize by line`,{args}) then # switch z and y f:=subs(y=1,z=y,homogeneous(f,x,y,z,polynom)); res:=procname(f,x,y,args[4..nargs]); return map(Normalizer,[res[1]/res[2],1/res[2]]) fi; res:=procname(subs(x=x+y,f),args[2..nargs]); return [Normalizer(res[1]+res[2]),op(2..nops(res),res)] fi; # Now coeff(f,y,n)<>0, so [0,1,0] is not a point on the curve if Normalizer(discrim(subs(z=0,x=1,homogeneous(f,x,y,z,polynom) ),y))=0 then res:=procname(numer(subs({x=x/(x+1),y=y/(x+1)},f)),args[2..nargs]); # Apply the same transformation on the result and normalize # the result in Q(s)[t]/(res[3]) if nops(res)=2 then return [Normalizer(res[1]/(res[1]+1)), Normalizer(res[2]/(res[1]+1))] else return [FFDIV(res[3],s,t,res[1]/(res[1]+1)), FFDIV(res[3],s,t,res[2]/(res[1]+1)),res[3]] fi fi; # now there are no ramification points nor singularities at infinity b:=integral_basis(f,x,y); if add(degree(denom(lcoeff(i,y)),x),i=b)<>(n-1)*(n-2)/2 then error "genus is not zero" fi; if n=2 then if member(`parametrize by line`,{args}) then return parametrize_conic(f,x,y,s) else return subs(x=s,y=t,[x,y,f]) fi elif n=3 then return parametrize_cube(f,x,y,s) fi; # elimate denominators: den:=denom(b[n]); b:=[seq(g_normal(den*i),i=b)]; w:=L_inf(b,x,y,den,n); v:=expand([1*den,x*den,x^2*den,op(w),seq(i*x,i=w)]); # Now v is a basis for L(2*(line at infinity)) # L(-canonical divisor) = L(2*(line at infinity)) intersected # with the set {b | b*diff(all integral elements,x) is integral} # First treat the ramification points at the multiplicity 1 factor of # the discriminant. We take those elements of the vector space spanned # by v that have roots in this points. To do this take a function a # that has poles in these points and then multiply with v d:=numer(discrim(f,y)); d1:=evala(Sqrfree(d,x),'expanded')[2][1]; if d1[2]=1 then # the discriminant does have a multiplicity 1 factor d1:=expand(d1[1]); # (*) Commented the following out, because this case is so rare so it's # not worth the trouble to make a special case for it, and caused # a bug (namely the 2*n+1 a little further should be nops(v)). # if member(n-degree(d1,x),{1,2}) then # # then we'll have 2 or 3 functions left afterwards # v:=[1,x,op(w)] # fi; a:=numer(FFDIV(f,x,y,diff(f,x),diff(f,y))); # Now a/d1 has poles of order 1 at the ramification points of d1 # Take a generic number i0: i0:=0; while has(g_gcd(subs(y=i0,f),d1),x) do i0:=i0+1 od; # a*y^i mod f, with y=i0 substituted. So these are polynomials in x. y_powers:=[seq(subs(y=i0,rem(expand(y^i*a),f,y)),i=0..n-1)]; zero:=add(co[i]*evala(Rem(collect( add(coeff(v[i],y,j)*y_powers[j+1],j=0..n-1) ,x),d1,x)),i=1..2*n+1); vars:={seq(co[i],i=1..2*n+1)}; v:=subs(solve({coeffs(collect(zero,x),x)},vars) ,add(co[i]*v[i],i=1..2*n+1)); vars:=indets(v) intersect vars; v:=[seq(subs(solve(vars,vars),subs(i=1,v)),i=vars)]; # The dimension of v should have dropped degree(d1,x) because we # have specified roots for the functions in v in degree(d1,x) points if nops(v) <> 2*n+1-degree(d1,x) then error "found wrong number of solutions" # = bug fi else d1:=1 fi; # switch the order in the integral basis so the basis elements are # used in a different order in the following loop. bb:=[b[1],b[2],seq(b[n-i],i=0..n-3)]; if nops(v)>3 then # integral basis is not recomputed because of the # options remember rf:=add(denom(i),i=integral_basis(f,x,y)); # The following contains the relevant factors of the # discriminant for frac_integral_a rf:=numer(Normalizer(d/rf^2/d1)) fi; i:=1; while nops(v)>3 do i:=i+1; # a:=diff(generic integral element,x) # (each integral basis element need not be generic, but we # treat them one by one, and together they are generic). # The denominator d1 has already been taken care of and # can be disregarded a:=FFDIV(f,x,y,subs(RootOf(f,y)=y,diff(subs(y=RootOf(f,y), bb[i]/den),x)*d1)); # Compute a basis for those elements b in the vectorspace # generated by v which satisfy a*b = integral v:=frac_integral_a(f,x,y,v,a,b,den,n,rf) # For a=diff(generic integral element,x) the dimension of the # L(D) we end up with must be 3. od; v:=[seq(evala(Expand(i/icontent(i))),i=v)]; # The following can't occur now anymore because of (*). # if nops(v)=2 then # return express_in_p(f,x,y,s,v[1],v[2]) # fi; # Write v[1]=s v[2]=t v[3]=u and search for a degree 2 relation: i0:=-1; f2:=add(add(co[i,j]*s^i*t^j*u^(2-i-j),i=0..2-j),j=0..2); vars:=indets(f2,name) minus {s,t,u}; while nops(vars)>1 do i0:=i0+1; f2:=subs(solve({coeffs(collect(evala(Rem(expand(subs(x=i0, subs(s=v[1],t=v[2],u=v[3],f2))),subs(x=i0,f),y)),y),y)}, vars),f2); vars:=indets(f2) intersect vars od; f2:=subs(op(vars)=1,f2); vars:=[s,t,u]; for i from 1 to 3 do if degree(f2,vars[i])=1 then # The quotient of the other two generates the function field return express_in_p(f,x,y,s,op({op(v)} minus {v[i]})) fi od; # the degree is 2 in all three variables s,t,u w:=search_rat_param(f,x,y,z,f2,s,t,u,v,n,d, op({`parametrize by line`} intersect {args})); if w<>'failed' then return w fi; f2:=subs(u=1,f2); w:=inverse_of_g(f,f2,s,t,u,v,x,y); j:=odd_singularity_on_C(f,x,y,d); if not member(`parametrize by line`,{args}) and j='failed' then # Then parametrization by a conic is fine return w fi; # Now we have the bijection to C2 in two directions # Lets see if we can find an odd degree place: if j='failed' then userinfo(2,'algcurves',`parametrizing the conic`,f2); j:=parametrize_conic(f2,s,t,x); userinfo(2,'algcurves', `composing this parametrization with a map from the conic to`,f); # Compose this with the morphism from C2 to C: map(evala@Normal,subs(s=j[1],t=j[2],x=s,[w[1],w[2]])) else # solve w[1]=j[1] a:=evala(Resultant(w[1]-j[1],f2,t)); # The following speeds things up a little, but is not strictly # necessary: a:=evala(Norm(a,{},g_ext_r(f))); # Now a is a polynomial in s, on which there is # an odd degree algebraic point. search_rat_param(f,x,y,z,f2,s,t,u,v,n,d, `point is given`,odd_point_on_a_old(f2,s,t,a,g_ext_r([f,a]))) fi end: odd_point_on_a_old:=proc(f,x,y,a,ext) local r,v; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; userinfo(4,'algcurves',`searching for a point, algebraic of odd-degree`); r:=odd_root(a,x,ext); v:=g_factors(subs(x=r,f),y,g_ext_r([ext,r])); if degree(v[1][1],y)<>1 then return procname(f,x,y, numer(g_normal(a/evala(Norm(x-r)))),ext) fi; [r,subs(solve({v[1][1]},{y}),y),1] end: # Let g be a bijective morphism from C to C2. This procedure computes # the inverse morphism. inverse_of_g:=proc(f,f2,s,t,u,v,x,y) # local res,resp,i0,p,i; local i; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; # res:= [seq(FFDIV(f2,s,t,subs(solve( {genus1_charpol(v[1]/v[3],op(i),f,s)-genus1_charpol(v[2]/v[3], op(i),f,t)},{i[1]}),i[1])),i=[[x,y],[y,x]]),f2]; # ### This correctness test is not necessary because `algcurves/parametrization` ### has a builtin test as well. # #if not has([res,f],RootOf) then # # Check the result (not in the RootOf case, then mod p is problematic) # p:=nextprime(10000); # i0:=0; # while Expand(subs(s=i0,denom(res[1])*denom(res[2])*numer(res[3]))) # mod p=0 do # i0:=i0+1; # p:=nextprime(p) # od; # resp:=subs(s=i0,[Expand(res[1]),Expand(res[2]) # ,Expand(numer(res[3]))]) mod p; # resp:=[Rem(resp[1],resp[3],t),Rem(resp[2],resp[3],t),resp[3]] mod p; # if Rem(Expand(subs(x=resp[1],y=resp[2],Expand(subs(s=i0,numer(f))) # mod p)) mod p,resp[3],t) mod p <> 0 then error "failed check" # fi; # userinfo(3,'algcurves',`Checked correctness, OK`) #fi; # res end: # Compute the characteristic polynomial of x over Qbar(x0v) # f is a polynomial in x and y giving the algebraic relation between x and y # x0v is an element of Qbar(x,y) # We use the variable t for x0v. genus1_charpol:=proc(x0v,x,y,f,t) local r,ansatz,a,n,j,k,i,aa; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; n:=degree(f,y); # Not necessarily true: a[n,2]:=1; ansatz:=add(t^k*(a[k,2]*x^2+a[k,1]*x+a[k,0]),k=0..n); aa:=indets(ansatz,name); j:=-1; i:=0; while i<5 do j:=j+1; r:=map(Normalizer,subs(x=j,[denom(x0v),f])); if r[1]=0 or not has(r[2],y) or has(g_gcd(op(r)),y) then next fi; r:=numer(t-subs(x=j,x0v)); r:=collect(evala(Resultant(r,subs(x=j,f),y)),t,g_normal); r:=collect(r/lcoeff(r,t),t,g_normal); if degree(r,t)<>n then next fi; ansatz:=subs(solve({seq(subs(x=j,coeff(ansatz,t,k)- coeff(r,t,k)*coeff(ansatz,t,n)),k=0..n-1)},aa),ansatz); i:=i+1 od; r:=numer(g_normal(ansatz/lcoeff(ansatz,t))); collect(r/lcoeff(r,x),x,g_normal) end: #-------------------------------------------------------------- # ratpar_odd_points # Input: a polynomial # Output: a root of the smallest factor of odd degree # ext: give the base field. # RootOf syntax odd_root:=proc(f,x,ext) local v,m,i,w; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; v:=g_factors(args); if v=[] then return 'failed' fi; w:={}; for i in v do if irem(degree(i[1],x),2)=1 then w:={op(w),i[1]} fi od; if w={} then return 'failed' fi; m:=min(seq(degree(i,x),i=w)); v:={}; for i in w do if degree(i,x)=m then v:={op(v),i} fi od; m:=min(seq(length(i),i=v)); for i in v do if length(i)=m then return RootOf(i,x) fi od; error "should not reach this point" end: # Construct an odd regular point, from an odd singular point # with an odd multiplicity. # RootOf syntax, except P is in rootof syntax. regpoint_from_sing:=proc(f,x,y,z,n,disc,ex,P) local ext,d,v,w,a,i,j; global g_conversion2; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; ext:=g_ext([args]); d:=0; while (d+2)*(d+1)/2 <= degree_ext(ext,g_ext(f))+1 do d:=d+1 od; userinfo(4,'algcurves',`intersecting the curve with a curve of degree`, d,`through the singularity`,subs(g_conversion2,P)); v:=add(add(a[i,j]*x^i*y^j,i=0..d-j),j=0..d); w:=g_expand(subs(x=P[1][1],y=P[1][2],v),ext); w:=solve(subs(g_conversion2,{ext_to_coeffs(w,g_ext(f))}), indets(v,name) minus {x,y}); v:=subs(w,v); w:=[op(indets(v) intersect {seq(seq(a[i,j],i=0..d-j),j=0..d)})]; if member(`second try`,{args}) then w:=subs(seq(w[i]=1,i=1..nops(w)),v) else w:=subs(w[1]=1,seq(w[i]=0,i=2..nops(w)),v) fi; a:=denom(FFDIV(f,x,y,1,w)); a:=g_normal(a/g_gcd(a,disc)); # to avoid roots of the discriminant do v:=odd_root(a,x,ex); if v='failed' then if member(`second try`,{args}) then userinfo(5,'algcurves', `failed to find a point which is algebraic of odd-degree`); # This case is in principle very unlikely, if it occurs it is # probably a bug. return 'failed' else # Now try with degree d one more. return procname(args,`second try`) fi fi; w:=odd_root(subs(x=v,f),y,g_ext_r([v,f])); if w='failed' then # To make sure that odd_root(a,x,ex) takes a different root # next time: a:=g_normal(a/evala(Norm(x-v,{},ex))) else userinfo(4,'algcurves',`Found a point`,[v,w,1], `which is algebraic of odd-degree on`,f); # Furthermore v is not a root of the discriminant, which # makes it easier to compute the image of this point under # the map from f to the conic. return [v,w,1] fi od end: # Try to find a regular point, avoiding disc=0, on C, algebraic of odd degree. # RootOf syntax odd_regpoint_C:=proc(f,x,y,z,n,disc,ext) local v,i,m,w; global g_conversion1; # First look at infinity option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; v:=odd_root(subs(z=0,x=1,homogeneous(f,x,y,z,polynom)),y,ext); if v<>'failed' then return [1,v,0] fi; # First look on x=0, x=1, x=-1 and x=infinity if irem(n,2)=1 then m:=100 else m:=1 fi; for i from -1 to m do if not testeq(subs(x=i,disc)) then v:=odd_root(subs(x=i,f),y,ext); if v<>'failed' then return [i,v,1] fi fi od; if irem(degree(f,x),2)=1 then m:=100 else m:=1 fi; for i from -1 to m do v:=odd_root(subs(y=i,f),x,ext); if v<>'failed' and not testeq(subs(x=v,disc)) then return [v,i,1] fi od; v:=find_points(subs(g_conversion1,f),x,y,g_ext(f)); w:={}; for i in v do if irem(i[2]*degree_ext(i,g_ext(f)),2)=1 then # Both the multiplicity i[2], and the degree of the # field extension of this point is odd. w:=w union {i} fi od; if w<>{} then m:=min(seq(degree_ext(i,g_ext(f)),i=w)); # Construct an odd regular point from this: for i in w do if degree_ext(i,g_ext(f))=m then return regpoint_from_sing(args,i) fi od fi; 'failed' end: # Assume here that the degree of f is more than 2. # Try to find a point on C2 which is algebraic of odd degree # RootOf syntax odd_point_on_C2:=proc(f,x,y,z,f2,s,t,u,g,n,d) local ext,G,i,m,P; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; ext:=g_ext_r(f); G:=g_factors(subs(u=0,s=1,f2),t,ext); if degree(expand(G[1][1]),t)=1 then # Found a rational point by luck return [1,subs(solve({G[1][1]},{t}),t),0] fi; P:=odd_regpoint_C(f,x,y,z,n,d,ext); if P='failed' then return P fi; m:=max(seq(degree(i,{x,y}),i=g)); G:=[seq(z^(m-degree(i,{x,y}))*homogeneous(i,x,y,z,polynom),i=g)]; G:=evala(subs({x=P[1],y=P[2],z=P[3]},G)); if G=[0,0,0] then # odd_point_on_C shouldn't give a point that has # this effect. userinfo(3,'algcurves', `Should have (but did not) find an odd degree point.`); return 'failed' fi; G end: # Search for a rational (i.e. defined over the same field as f) point # on f2. # RootOf syntax rat_point_on_C2:=proc(f,x,y,z,f2,s,t,u,g,n,d) local P,d0,co,vars,i,j,G,ext; global g_conversion1,g_conversion2; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; if member(`point is given`,{args}) then P:=args[nargs] else P:=odd_point_on_C2(f,x,y,z,f2,s,t,u,g,n,d) fi; if P='failed' then return P fi; d0:=degree_ext(P,f); if d0=1 then return P fi; userinfo(3,'algcurves',`Computing a rational point`); G:=add(add(co[i,j]*s^i*t^j*u^((d0+1)/2-i-j) ,i=0..(d0+1)/2-j),j=0..1); vars:=indets(G,name) minus {s,t,u}; ext:=g_ext([f,P]); G:=subs(solve(subs(g_conversion2,{ext_to_coeffs(g_expand( subs(s=P[1],t=P[2],u=P[3],g_conversion1,G),ext),g_ext(f))}) ,vars),G); vars:=[op(vars intersect indets(G))]; G:=g_normal(subs(vars[1]=0,vars[2]=1,G)/subs(vars[2]=0,vars[1]=1,G)); G:=FFDIV(op(subs(u=1,[f2,s,t,G]))); P:=subs(solve({denom(G)},{s}),s); [P,subs(solve({subs(s=P,numer(G))},{t}),t),1] end: # Try to find a rational point on C2. Then from it construct a # parameter, then express x and y in terms of it. # RootOf syntax search_rat_param:=proc(f,x,y,z,f2,s,t,u,g,n,d) local P,P1,P2; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; P:=rat_point_on_C2(args); if P='failed' then # This part is obsolete. # if member(`parametrize by line`,{args}) and # `dont compute map from C2 to C`=true and # odd_singularity_on_C(f,x,y,d)='failed' then # P:=[1,RootOf(subs(u=0,s=1,f2),t),0] # else # Exit this method and we'll find a rational # parametrization using the result of odd_singularity_on_C # or if the option `parametrize by line` is not given # we can also compute a rational bijection with a # conic. return P # fi fi; userinfo(4,'algcurves',`constructing a parameter from the point`,P); if P[3]<>0 then P1:=P[3]*g[1]-P[1]*g[3]; P2:=P[3]*g[2]-P[2]*g[3] else # switch 1 and 3 P1:=P[1]*g[3]-P[3]*g[1]; P2:=P[1]*g[2]-P[2]*g[1] fi; # P1/P2 is a parameter, i.e. generates the function field express_in_p(f,x,y,s,P1/icontent(P1),P2/icontent(P2)); end: # Give a singular point which contains a place that is algebraic of odd degree. # Only to be used if no odd degree regular point could be found. If an # odd degree odd multiplicity point can be found then we'll use # regpoint_from_sing first, so this procedure will be only useful in rare # cases (i.e. the case where no odd degree odd multiplicity point was found, # but still an odd degree place was found). # RootOf syntax odd_singularity_on_C:=proc(F,x,y,disc) local v,ext,res,i,m,p,extp,f,f_translated,pl_transl; global g_conversion1,g_conversion2; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; ext:=g_ext(F); f:=subs(g_conversion1,F); v:=find_points(f,x,y,ext,[`search points on`,{1,3}]); res:={}; for i in v do if nops(i)=4 and i[4]=3 then if i[1][3]=0 then error "should not reach this point" fi; # case not occur because the # existence of odd places at infinity has already been checked # in odd_regpoint_C. # Now there is an odd place on the line x=i[1][1] res:=res union {i} fi od; if res={} then return 'failed' fi; # Take the one with the smallest algebraic degree: m:=min(seq(degree_ext(i,ext),i=res)); # Now determine which singularity it was (in fact this has already # been determined inside find_points, but is not given in the output # of that procedure. Recomputing won't take time because of the # options remember in the Puiseux algorithm). for i in res do if degree_ext(i,ext)=m then p:=i[1][1]; extp:=g_ext([p,ext]); f_translated:=g_expand(subs(x=x+p,f),extp); pl_transl:=[op(`puiseux/technical_answer`(f_translated,x,y ,0,extp))]; res:=NULL; for i in pl_transl do if irem(i[6],2)=1 then res:=res,i fi od; res:=[res]; m:=min(seq(i[6],i=res)); for i in res do if i[6]=m then return subs(g_conversion2,[p,subs(x=0,i[1]),1]) fi od fi od end: # Input: f of degree 3, degree(f,y) should be 3 (necessary for express_in_p). # Output: a parametrization # RootOf syntax parametrize_cube:=proc(f,x,y,t) local s; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; s:=singularities(f,x,y); if s={} then error "genus is not 0" fi; s:=s[1][1]; if s[3]=0 then # This case shouldn't occur anyway, it is caught in ratpar express_in_p(f,x,y,t,y*s[1]-x*s[2],1) else express_in_p(f,x,y,t,x-s[1],y-s[2]) fi end: # Input: f of degree 2. degree(f,y) should be 2. # Output: a parametrization # RootOf syntax parametrize_conic:=proc(f,x,y,t) local s,P; # Take a point at infinity and compute a parametrization option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; s:=RootOf(g_factors(subs(t=0,x=1,homogeneous(f,x,y,t,polynom)) ,y,g_ext_r(f))[1][1],y); if indets(f,{'name','RootOf'}) = {x,y} and has(s,RootOf) then # try to find a rational point on the curve using `isolve/pyth` P:=pyth2(f,x,y); if P<>FAIL then userinfo(2,'algcurves',`Found a rational point on the conic`); return `algcurves/iss94`(f,x,y,t,[P[1],P[2],1]) fi fi; express_in_p(f,x,y,t,y-x*s,1) end: # find a rational finite point on a conic, that does # not have a rational point at infinity. pyth2:=proc(f,x,y) local P,c,i; option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; if indets(f,{'name','RootOf'})<>{x,y} or {degree(f,y),degree(f,{x,y})}<>{2} then error "wrong input" elif degree(f,x)=1 then P:=[subs(solve({subs(y=0,f)},{x}),x),0] elif coeff(f,y,1)<>0 then c:=-coeff(f,y,1)/coeff(f,y,2)/2; P:=procname(expand(subs(y=y+c,f)),x,y); if P<>FAIL then P:=[P[1],P[2]+subs(x=P[1],c)] fi elif coeff(f,x,1)<>0 then c:=-coeff(f,x,1)/coeff(f,x,2)/2; P:=procname(expand(subs(x=x+c,f)),x,y); if P<>FAIL then P:=[P[1]+c,P[2]] fi else _Env_ifactor_easy:=true; # avoid getting stuck in ifactor, note # that this can cause `isolve/pyth` to generate errors. P:=traperror(`isolve/pyth`( primpart(homogeneous(f,x,y,c)),x,y,c,[])); if P=lasterror or subs(P,[x,y,c])=[0,0,0] then P:=FAIL else P:=subs(P,[x,y,c]); while indets(P)<>{} do c:=indets(P,'name')[1]; i:=0; while normal(subs(c=i,numer(P[3])))=0 do i:=i+1 od; P:=map(normal,subs(c=i,P)) od; P:=[P[1]/P[3],P[2]/P[3]]; if subs(x=P[1],y=P[2],f)<>0 then P:=FAIL; userinfo(1,'algcurves', `Error in finding point on the conic`) # I don't know if this can occur fi fi fi; P end: #-------------------------------------------------------------- # ratpar_various # read(ratpar_OK); # Assumption: (0,1,0) is not a point on the curve # v1/v2 generates the function field. # Compute X(s) and Y(s) such that x=X(v1/v2) and y=Y(v1/v2) # RootOf syntax express_in_p:=proc(f,x,y,s,v1,v2) global z; local R,z1,z2,zero; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; zero:=expand(v1-s*v2); z1:=subs(x=1,homogeneous(zero,x,y,z,polynom)); # x=1 can be assumed at z=0 because (0,1,0) is not on the curve z2:=subs(z=0,expand(z1/z^ldegree(z1,z))); # This is zero at the denominator of the parametrization. Now we # try to find the denominator by solving z2=subs(z=0,f)=0 R:=subs(z=0,x=1,homogeneous(f,x,y,z,polynom)); R:=evala(Resultant(z2,R,y)); if R<>0 then R:=R/icontent(R) fi; [seq(express_x_in_p(f,op(z1),s,zero,R),z1=[[x,y],[y,x]])] end: # d is a polynomial in s which vanishes if x=infty. The denominators # of X(s) and Y(s) divide d. express_x_in_p:=proc(f,x,y,s,zero,d) local X,cX,vars,i,R; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; if d=0 # Denominator is not yet known, so I'll have to use the old algorithm or testeq(discrim(d,s)) then return compute_x_old(args) fi; X:=add(cX[i]*s^i,i=0..degree(f,{x,y})); vars:=indets(X,name) minus {s}; i:=-1; while has(X,vars) do i:=i+1; if i=6 then userinfo(5,'algcurves', `giving up, trying a different method`); return compute_x_old(args) fi; R:=subs(x=i,[zero,f]); if degree(R[1],s)*degree(R[2],y)*degree(R[1],y)=0 or has(g_gcd(op(R)),y) then next fi; R:=evala(Resultant(op(R),y)); X:=evala(Expand(subs(solve({coeffs(collect(evala(Rem(X-i*d, expand(R),s)),s),s)},vars),X))) od; g_normal(X/d) end: # f is a genus 0 curve in x and y. p is the parameter, i.e. p generates # the function field # Output: a rational function in t that gives x as a function of p. compute_x_old:=proc(f,x,y,t,zero) local i,X,np,dp,R,ansatz,C,ar; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; # The following case shouldn't occur anyway: # if not member(y,indets(f)) then return solve(f,x) fi; dp:=-coeff(zero,t,1); np:=coeff(zero,t,0); X[1]:=-1; for i from 1 to 3 do R[i]:=0; while degree(R[i],t) degree(collect(g_gcd(dena,i^(z+1)),x),x) do z:=z+1 od; dena:=g_gcd(dena,i^z); numa:=evala(Rem(expand(numer(a)),dena,x)); ansatz:=add(co[i]*v[i],i=1..nops(v)); z:=0; for i from 1 to nops(v) do z:=z+co[i]*evala(Rem(expand(numa*v[i]),f,y)) od; # now z/dena/den should be integral so:=solve(evala({coeffs(expand(evala(Rem(expand(z) ,dena,x))),[x,y])}),vars); ansatz:=subs(so,ansatz); z:=subs(so,z); vars:=vars intersect indets(ansatz); userinfo(5,'algcurves',`Done solving linear equations.`, nops(vars),`variables remaining`); if nops(vars)>3 then # Now z is divisible by dena z:=subs(g_conversion1,g_normal(z/dena)); for i from d-1 by -1 to 0 do userinfo(5,'algcurves',`step number`,i); while degree(coeff(z,y,i),x) >=degree(lcoeff(b[i+1],y),x) and coeff(z,y,i)<>0 do z:=g_expand(z-b[i+1]* g_normal(lcoeff(coeff(z,y,i),x)/lcoeff(lcoeff(b[i+1],y),x)) *x^(degree(coeff(z,y,i),x)-degree(lcoeff(b[i+1],y),x)),ext); # Without normalizing you never know if coefficients are # going to vanish, so just in case: if indets(f) minus {x,y} <> {} then z:=collect(z,{x,y},'distributed',g_normal) fi od od; equations:=subs(g_conversion2,{coeffs(z,[x,y])}); userinfo(4,'algcurves',`solving linear equations...`); z:=subs(solve(equations,vars),ansatz); vars:=indets(z) intersect vars; userinfo(4,'algcurves',`Done solving.`,nops(vars),`variables remaining`) fi; [seq(g_normal(subs(solve(vars,vars),subs(i=1,z))),i=vars)] end: # Function Field DIVision, compute a/b in Qbar(x)[y]/(f) # Algebraic numbers in RootOf syntax FFDIV:=proc(f,x,y,a,b) local d,q,i,qc; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; if nargs=4 then d:=normal(a); return FFDIV(f,x,y,numer(d),denom(d)) fi; userinfo(4,'algcurves',`doing a division in the function field...`); d:=degree(f,y); q:=add(qc[i]*y^i,i=0..d-1); q:=evala(Normal(subs(solve(evala({coeffs(collect(rem(collect(q*b-a,y), f,y),y,numer),y)}),{seq(qc[i],i=0..d-1)}),q))); userinfo(4,'algcurves',`Done division.`); q end: # Output: L(line at infinity), i.e. a basis for all integral functions # with a pole of order <= 1 at infinity (we skip the elements 1 and x). # Assumptions: no singularities nor ramification points at infinity # b/den=integral basis, the genus is 0, the curve is irreducible over the # algebraic closure of the constants and y is integral over C[x]. # Note that the output should still be divided by den to get the correct # answer. # Algebraic numbers in RootOf syntax L_inf:=proc(b,x,y,den,n) local R,vars,i,j,d,Cr,Cb,di,ai,S,co,C; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; userinfo(4,'algcurves',`Computing L(divisor of line at infinity)...`); R:=0; # variables which must not be eliminated: vars:=indets([args,seq(C[i],i=0..n)],name); # The basis for L(1*line) will be obtained # by substituting (0,0,..,1,0,0,..) in the C[i] in R. # The coefficient of y^i*(suitable power of x)/den must be C[i]. # At infinity this looks like: y^i/z + ... for i from n-1 by -1 to 0 do userinfo(5,'algcurves',`step number`,i); d:=1+degree(den,x)-i; if d<0 then error "wrong degree estimate" # = bug fi; Cr:=coeff(R,y,i)-C[i]*x^d; Cb:=lcoeff(b[i+1],y); di:=degree(Cr,x)-degree(Cb,x); ai:=add(co[i,j]*x^j,j=0..di); S:=collect(Cr-ai*Cb,x); # from x^d and up this should be zero S:={co[0,0],seq(coeff(S,x,j),j=d..degree(S,x))}; R:=collect(subs(solve(evala(S),indets(S,name) minus vars), R-ai*b[i+1]),{x,y},'distributed',g_normal) od; vars:={seq(C[i]=0,i=0..n)}; userinfo(4,'algcurves',`Done computing L(divisor).`); [seq(subs(vars,subs(C[i]=1,R)),i=1..n-1)] end: # Compute the gcd # RootOf syntax g_gcd:=proc() local a,b; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; a:=numer(args[1]); # Otherwise evala/Gcd could return an error message b:=numer(args[2]); if has([a,b],RootOf) then primpart(evala(Gcd(a,b))) else primpart(gcd(a,b)) fi end: #savelib ('ratpar','odd_point_on_a_old','inverse_of_g','genus1_charpol',\ 'odd_root','regpoint_from_sing','odd_regpoint_C','odd_point_on_C2',\ 'rat_point_on_C2','search_rat_param','odd_singularity_on_C',\ 'parametrize_cube','parametrize_conic','express_in_p',\ 'express_x_in_p','compute_x_old','frac_integral_a','FFDIV',\ 'pyth2',\ 'L_inf','g_gcd'):