# This file is just the concatenation of the source files for the # algcurves (used to be called IntBasis) package. These source files # themselves can be downloaded from # http://www-math.sci.kun.nl/math/compalg/IntBasis/ # Simply type read(this file); in Maple 5.3 (should work under 5.4 as well) # If you find bugs or other problems please send them to hoeij@sci.kun.nl # macro( readlib=proc() args end ); # This macro is un-done at the end of the file. #------------------------ filename: algcurves/g_expand --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/g_expand,v $ # $Notify: hoeij@sci.kun.nl $ # Purpose of this file: (was first called: ground_field) # Procedures to allow code that works over different kinds of ground # fields. So for example g_expand looks what kind of expression is given, # and then decides what to do, expand, or evala(Expand()), or expanding # with normalization of the lowest coefficient etcetera. macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): g_conversion1:={}: # RootOf syntax -> rootof syntax g_conversion2:={}: # rootof syntax -> RootOf syntax # evala(Normal( )) g_normal:=proc(aa) global g_conversion1,g_conversion2; local a; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if has(aa,RootOf) then RETURN(evala(Normal(aa),'independent')) fi; a:=subs(g_conversion2,aa); if has(a,RootOf) then subs(g_conversion1,evala(Normal(a),'independent')) else normal(a) fi end: # evala(Expand( )) # Give 3'rd argument x for normalizing the lowest coefficient in x # So then ldegree(result,x) gives the right answer g_expand:=proc(a,ext) global g_conversion1,g_conversion2; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if type(a,polynom) then g_evala(expand(a),ext) elif nargs=3 and type(a,polynom(anything,ext)) then normal_tcoeff(g_evala(expand(a),ext),args[3]) elif nargs=3 then normal_tcoeff(procname(a,ext),args[3]) else subs(g_conversion1,evala(Expand(subs(g_conversion2,a)) ,'independent')) fi end: # Normalize the lowest coefficient normal_tcoeff:=proc(a,x) local c,i,r,d; # if not testeq(subs(g_conversion2,tcoeff(a,x))) then RETURN(a) fi; # Strangely enough this testeq line doesn't seem to speed things up in # the test examples option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; r:=a; do d:=ldegree(r,x); c:=g_normal(coeff(r,x,d)); if coeff(r,x,d)=c then RETURN(r) fi; r:=convert([c*x^d,seq(coeff(r,x,i)*x^i,i=d+1..degree(r,x))],`+`); if c<>0 then RETURN(r) fi od end: # Input: A polynomial in the rootofs in ext # Output: a evala normalization g_evala:=proc(a,ext) local i,e; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if nops(ext)=0 then RETURN(a) elif nops(ext)=1 then e:=ext[1]; expand(convert([seq(coeff(a,e,i)*g_evala_rem(e^i) ,i=0..degree(a,e))],`+`)) else e:=g_evala(a,[op(2..nops(ext),ext)]); g_evala(expand(convert([seq(coeff(e,ext[1],i)*g_evala_rem(ext[1]^i) ,i=0..degree(e,ext[1]))],`+`)),[op(2..nops(ext),ext)]) fi end: g_evala_rem:=proc() global g_conversion1,g_conversion2; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; expand(subs(g_conversion1,evala(Expand(subs(g_conversion2,args)) ,'independent'))) end: # Input : an irreducible polynomial kk in x, not necessarily monic # Output: a root of k # If an algebraic extension is needed it will be placed in ext. g_zero_of:=proc(k,x,ext) global g_conversion1,g_conversion2; local a,v; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if degree(k,x)=1 then ext:=NULL; RETURN(g_normal(-coeff(k,x,0)/coeff(k,x,1))) fi; a:=RootOf(subs(g_conversion2,k),x); if not member(_Z,indets(op(subs(g_conversion1,a)))) then a:=subs(g_conversion1,a); ext:=a; RETURN(a) fi; if g_conversion1={} then g_conversion1:=NULL fi; v:=nops(g_conversion2); g_conversion1:=a=rootof.v,g_conversion1; g_conversion2:={rootof.v=a,op(g_conversion2)}; ext:=rootof.v end: # Factorization over the groundfield # Both RootOf and rootof syntax # ext: describes the groundfield g_factors:=proc(f,x,ext) global g_conversion1,g_conversion2; local v,w,i; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if not has(f,x) then RETURN([]) elif ext<>[] and not has(ext,RootOf) then # rootof syntax input, so rootof syntax output RETURN(subs(g_conversion1, procname(op(subs(g_conversion2,[args]))))) fi; userinfo(5,algcurves,` Factorizing`,f); v:=traperror(factors(numer(f),{op(ext)})); if v=lasterror then # Now try evala(Factor()), it can handle more # cases than factors g_ext([args]); v:=subs(g_conversion2,factors(numer(subs(g_conversion1 ,evala(Factor(numer(f),{op(ext)})))))) fi; v:=v[2]; userinfo(5,algcurves,`Done factorization`); w:=NULL; for i in v do if has(i,x) then w:=w,i fi od; [w] end: # Gives the algebraic extensions appearing in aa. g_ext:=proc(aa) global g_conversion1,g_conversion2; local v,i,result,ii,vv; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; v:=g_ext_r(subs(g_conversion2,aa)); vv:=subs(g_conversion1,v); result:=NULL; for i from 0 to nops(g_conversion2) do if member(rootof.i,vv) then result:=rootof.i,result fi od; for i from nops(v) by -1 to 1 do if not member(subs(g_conversion1,v[i]) ,{seq(rootof.ii,ii=0..nops(g_conversion2))}) then if g_conversion1={} then g_conversion1:=NULL fi; vv:=nops(g_conversion2); g_conversion1:=v[i]=rootof.vv,g_conversion1; g_conversion2:={rootof.vv=v[i],op(g_conversion2)}; result:=subs(g_conversion1,v[i]),result fi od; [result] end: # used by g_ext g_ext_r:=proc(a) local v,vv,i,tail; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; v:=indets(a,RootOf); if nops(v)=0 then RETURN([]) fi; vv:={}; for i in v do vv:=vv union indets(op(i),RootOf) od; tail:=g_ext_r(vv); v:=[op(v minus vv)]; [op(v),op(tail)] end: truncate:=proc(aa,n,y,ext) local dummy,a; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; a:=collect(aa,y); a:=expand(convert([seq(y^dummy*coeff(a,y,dummy) ,dummy=ldegree(a,y)..n-1)],`+`)); g_evala(a,ext) end: #savelib('g_conversion1','g_conversion2','g_normal','g_expand',\ 'normal_tcoeff','g_evala','g_evala_rem','g_zero_of','g_factors',\ 'g_ext','g_ext_r','truncate','`algcurves/g_expand.m`'): #------------------------ filename: algcurves/puiseux --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/puiseux,v $ # $Notify: hoeij@sci.kun.nl $ macro( puiseux=`algcurves/puiseux`, v_ext_m=`algcurves/v_ext_m`, lift_exp=`algcurves/lift_exp`, lift_exp_m1=`algcurves/lift_exp_m1`, truncate_subs=`algcurves/truncate_subs`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, `integral_basis/bound`=`algcurves/ib_bound`, Newtonpolygon=`algcurves/Newtonpolygon` ): # This code uses the procedures in `algcurves/g_expand`. This # file will be read in the puiseux procedure macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): # Gives the zeros of the factors, their multiplicities and algebraic extensions # The input of this procedure is the output of g_factors v_ext_m:=proc(f,x) local ext,nulp,i,res; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; res:=NULL; for i in f do nulp:=g_zero_of(i[1],x,'ext'); ext:=eval(ext); res:=res,[nulp,i[2],[ext],degree(i[1],x)] od; [res] end: # This procedure lifts an expansion until it has multiplicity 1 # Input is a list v # v[1] = The expansion so far # v[2] = The accuracy, i.e. the lowest not yet determined coefficient. This # procedure starts with computing the coefficient of x^v[2] # v[3] = Ramification, if v[2]=x^2, then x stands for sqrt(x). # v[4] = The algebraic extensions in this expansion # v[5] = The multiplicity of this factor # v[6] = The degree of the algebraic extensions above the groundfield # v[7] = The number Int_p, which is the sum of the valuations of the differences # of this expansion with the other expansions. This sum is needed for # computing an integral basis in an algebraic function field, cf JSC 94. # For an approximation r of a puiseux expansion p we have v(p-r)+Int_p=v(f(r)) # v[8] = list of the previous expansions, needed for determining the valuations # of the differences with the other Puiseux expansions # f = minimal polynomial of y over L(x), must be monic # Start computation with [0,0,x,ext,degree(f,y),1,1-degree(f,y),[]],f,x,y lift_exp:=proc(v,f,x,y) local i,ii,r,res,v7,vv7,v3,ext,a,j,n,np,ram,j3; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if v[5]=1 then RETURN({v}) fi; v3:=degree(v[3],x); # the ramification index. res:={}; r:=v[1]+y*x^v[2]; vv7:=v[7]*v3+v[2]-1; # The ldegree of the previous step vv7:=vv7+v[5]; # Upperbound for the ldegree of this step # we try to find an equation for the unknown y ii:=truncate_subs(subs(x=v[3],f),x,y,r,vv7+1,v[4]); if ii=0 then ERROR(`degree estimate was wrong`) fi; v7:=(ldegree(ii,x)-v[2])/v3; # The number Int_p=v(f(r))-v(p-r) r:=v_ext_m(g_factors(tcoeff(ii,x),y,v[4]),y); for i in r do res:=res union lift_exp([v[1]+x^v[2]*i[1],v[2]+1,v[3] ,[op(i[3]),op(v[4])],i[2],v[6]*i[4],v7 ,[op(v[8]),[op(1..4,v)]]],f,x,y) od; if convert([seq(i[5]*i[6]/v[6]*degree(i[3],x)/v3,i=res)],`+`) <>degree(tcoeff(ii,x),y) then ERROR(`found wrong number of expansions`) fi; if v[5]=degree(tcoeff(ii,x),y) then if ldegree(ii,x)<>vv7 then ERROR(`degree estimate was wrong`) fi; RETURN(res) fi; # Make sure that all valuations of the coefficients are correct ii:=collect(ii,y); ii:=convert([seq(normal_tcoeff(coeff(ii,y,i),x)*y^i ,i=0..degree(ii,y))],`+`); # Now compute all slopes s between 0 and 1 np:=Newtonpolygon(ii,x,y); if nops(np)=2 and np[1][3]=0 then ERROR(`found wrong number of expansions`) fi; for j in np do if nops(j)>2 and j[3]>0 and j[3]<1 then r:=g_factors(j[4],x,v[4]); r:=v_ext_m(r,x); for i in r do j3:=j[3]-v[2]; ext:=[op(i[3]),op(v[4])]; n:=mods(1/numer(j3),denom(j3)); ram:=g_normal(i[1]^n)*x^denom(j3); a:=v[2]*denom(j3)-numer(j[3]); res:=res union lift_exp([ g_expand(subs(x=ram,v[1])+x^a* g_normal( i[1]^( (1-n*numer(j3))/denom(j3) ) ) ,ext),a+1, g_expand(subs(x=ram,v[3]),ext),ext,i[2],i[4]*v[6], (j[2]-j[1]*j[3]-a/degree(ram,x))/degree(v[3],x), [op(v[8]),[op(1..4,v)]]],f,x,y) od fi od; res end: # This procedure lifts (= computes more terms of) expansions that have # multiplicity 1. This could also be done with lift_exp, but this # procedure is faster. lift_exp_m1:=proc(v,f,x,y,nnk) local nk,res,powx,r,rr,v3,i; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; v3:=degree(v[3],x); # The ramification index nk:=ceil(nnk*v3)/v3; powx:=v[2]; if nkv[2]+v3*v[7] then ERROR(`wrong input`) fi; r:=collect(r/x^ldegree(r,x),x); while nk>powx/v3 do rr:=collect(coeff(r,x,0),y); if degree(rr,y)<>1 then ERROR(`degree is not 1`) fi; rr:=g_normal(-coeff(rr,y,0)/coeff(rr,y,1)); res:=collect(subs(y=rr+x*y,res),x); powx:=powx+1; if nk<=powx/v3 then break fi; r:=subs(y=rr+x*y,r); r:=g_expand(r,v[4],x); if ldegree(r,x)<>1 then ERROR(`ldegree is not 1`) fi; r:=convert([seq(coeff(r,x,i)*x^(i-1),i=1..v3*nk-powx)],`+`) od; [subs(y=0,res),powx,op(3..7,v)] end: # Input: f in k[x,y] # Output f(y=y_value) mod x^n truncate_subs:=proc(f,x,y,y_value,n,ext) local ym,i,res,F; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if n=0 then RETURN(0) elif degree(f,x)>=n then F:=collect(f,x); F:=convert([seq(coeff(F,x,i)*x^i,i=0..n-1)],`+`) else F:=f fi; res:=0; for i from 0 to degree(F,y) do if i=0 then ym:=1 else ym:=truncate(ym*y_value,n,x,ext) fi; res:=res+coeff(f,y,i)*ym od; res:=truncate(res,n,x,ext); # Make sure that ldegree(res,x) gives a right answer: normal_tcoeff(res,x) end: monic:=proc(f,x,y,ext,q) global g_conversion1,g_conversion2; local j,ff,lc,qq; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; ff:=expand(numer(g_normal(f))); lc:=lcoeff(ff,y); if not has(lc,x) then q:=1; RETURN(g_expand(g_normal(ff/lc),ext)) fi; lc:=subs(g_conversion1,evala(Sqrfree(subs(g_conversion2,lc)))); lc:=lc[2]; lc:=g_expand(product(lc[j][1],j=1..nops(lc)),ext); ff:=monic(subs(y=y/lc,ff),x,y,ext,'qq'); qq:=eval(qq); q:=g_expand(qq*lc,ext); ff end: # See ?puiseux puiseux:=proc(aa,eqn) global g_conversion1,g_conversion2; local x,v,a,f,y,ext,q,verz,i,res,n; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`); if nargs<3 then ERROR(`too few arguments`) fi; if nargs=3 and not type(aa,RootOf) then ERROR(`wrong arguments`) fi; x:=eval(op(1,eqn)); v:=op(2,eqn); if v=infinity then a:=subs(x=1/x,aa) else a:=subs(x=x+v,aa) fi; ext:=g_ext([op(a),op(2..nargs,[args])]); if nargs=3 then f:=subs(g_conversion1,op(a)); f:=subs(_Z=y,f); n:=args[3] else f:=subs(g_conversion1,a); y:=args[3]; n:=args[4] fi; if n=evaln(minimal) then n:=0 fi; if not type(n,numeric) then ERROR(`wrong arguments`) fi; f:=monic(f,x,y,ext,'q'); q:=expand(eval(q)); verz:=`puiseux/technical_answer`(f,x,y,`if`(n=0,0,n+ldegree(q,x)),ext); res:={}; for i in verz do if nargs>4 then res:=res union {[x=subs(x=args[5],i[3]), y=subs(x=args[5],i[1])/q]} else res:=res union {subs(x=(x/lcoeff(i[3],x))^(1/degree(i[3],x)), i[1])/q} fi od; if v=infinity then res:=subs(x=1/x,res) else res:=subs(x=x-v,res) fi; subs(g_conversion2,res) end: # Note: the options remember won't work because of the local y in puiseux # This procedure computes Puiseux expansions. # It's input must be in internal format. # Output: a list, see lift_exp `puiseux/technical_answer`:=proc(f,x,y,n,ext) local i,v; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; userinfo(5,algcurves,`Computing Puiseux expansions...`); # It works for degree >=1, but not degree 0, so: if not has(f,y) then ERROR(`wrong input`) fi; if n=0 then i:=degree(collect(subs(x=0,f),y,g_normal),y); i:=lift_exp([0,0,x,ext,i,1,1-i,[]],f,x,y) else v:=procname(f,x,y,0,ext); if n=`integral basis` then i:={seq(lift_exp_m1(i,f,x,y,`integral_basis/bound`(i,v,x)),i=v)} else i:={seq(lift_exp_m1(i,f,x,y,n),i=v)} fi fi; userinfo(5,algcurves,`Done computing Puiseux expansions.`); i end: # Give an bound upper bound for the accuracy of the Puiseux # expansion p that is required for the integral basis. # Bound=max of (Int_i - Int_p + v(p-i)) taken over all Puiseux expansions i `integral_basis/bound`:=proc(p,v,x) local V,i,m,j,p1,p2,k; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; m:=p[2]/degree(p[3],x); p1:=[op(p[8]),p]; for i in v minus {p} do # Now determine the valuation of the difference p2:=[op(i[8]),i]; j:=1; while p1[j]=p2[j] do j:=j+1 od; V:=min(seq((k[2]-1)/degree(k[3],x),k=[p1[j],p2[j]])); m:=max(m,i[7]-p[7]+V) od; min(m,floor(max(seq(i[7],i=v)))) end: # Input: a monic polynomial f in k((x))[y] # Output: sort of a Newton polygon # Note: only the last element in the output is correct Newtonpolygon:=proc(f,x,y) local i,n,ff,dummy,m,npg,powd,res,slope,val_powd,vals; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; n:=degree(f,y); ff:=expand(f+(y+1)^n*x^(degree(f,x)+1)); if n<2 then ERROR(`nothing to be factored`) fi; vals:=[seq(ldegree(coeff(ff,y,i),x),i=0..n)]; powd:=0; val_powd:=min(op(vals)); npg:=[[powd,val_powd]]; while powd(i-powd)*m do i:=i-1 od; powd:=i;val_powd:=vals[i+1]; npg:=[op(npg),[i,vals[i+1]]] od; res:=NULL; for i from 1 to nops(npg)-1 do # Now we compute the Newton polynomial of slope number i slope:=(npg[i+1][2]-npg[i][2])/(npg[i+1][1]-npg[i][1]); res:=res,[op(npg[i]),slope,convert([seq( coeff(coeff(ff,y,denom(slope)*dummy+npg[i][1]) ,x,dummy*numer(slope)+npg[i][2]) *x^dummy,dummy=0..(npg[i+1][1]-npg[i][1])/denom(slope) )],`+`)] od; [res,npg[nops(npg)]] end: #savelib ('v_ext_m','lift_exp','lift_exp_m1','truncate_subs', \ 'monic','puiseux','`puiseux/technical_answer`', \ '`integral_basis/bound`','Newtonpolygon','`algcurves/puiseux.m`'): #------------------------ filename: algcurves/plot_knot --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/plot_knot,v $ # $Notify: hoeij@sci.kun.nl $ ######################### # Singularity knots # ######################### # This code uses the procedures in `algcurves/g_expand` and # `algcurves/puiseux`. These files will be read in the plot_knot procedure macro( plot_knot=`algcurves/plot_knot`, all_substitutions=`algcurves/pl_allsub`, curve=`algcurves/pl_curve`, `puiseux/technical_answer`=`algcurves/puiseux_te`, g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_expand=`algcurves/g_expand`, g_ext=`algcurves/g_ext` ): # All possible floating point substitutions for algebraic numbers. # v: a list of RootOf's all_substitutions:=proc(v) local w,i,j; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if v=[] then RETURN({v}) fi; w:=evalf(allvalues(v)); # Note: syntax of allvalues has changed from 5.3 to 5.4 {seq([seq(v[i]=j[i],i=1..nops(v))],j=w)} end: # f is a polynomial in x en y plot_knot:=proc(f,x,y) global color,colours,epsilon,g_conversion1,g_conversion2,radius,tubepoints; local ext_f,pui,ext_pui,i,pui2,curven,t,opt,eps,cols; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`): readlib(`algcurves/puiseux`): if has(evala(Content(f,y)),x) then # If f has a component which is independent of y, this component # will be ignored by the puiseux algorithm, but I want to # include it in the computation. This is a fix for this case RETURN(procname(subs(x=x+y,f),args[2..nargs])) fi; opt:={op(4..nargs,[args])}; # default options if not member(numpoints,indets(opt)) then opt:=opt union {numpoints=150} fi; if not member(radius,indets(opt)) then opt:=opt union {radius=0.05} fi; if not member(tubepoints,indets(opt)) then opt:=opt union {tubepoints=5} fi; if not member(scaling,indets(opt)) then opt:=opt union {scaling=CONSTRAINED} fi; if not member(style,indets(opt)) then opt:=opt union {style=PATCHNOGRID} fi; eps:=1; if member(epsilon,indets(opt)) then eps:=subs(opt,epsilon); opt:=opt minus {epsilon=eps} fi; cols:=[]; if member(colours,indets(opt)) then cols:=subs(opt,colours); opt:=opt minus {colours=cols} fi; ext_f:=g_ext(f); pui2:=`puiseux/technical_answer`(g_expand(subs(g_conversion1,f),ext_f) ,x,y,0,ext_f); pui:={}; for i in pui2 do if coeff(i[1],x,0)=0 and ldegree(i[1],x)>=0 then pui:=pui union {i} fi od; ext_f:=subs(g_conversion2,ext_f); ext_f:=all_substitutions(ext_f); ext_f:=ext_f[1]; ext_pui:=subs(ext_f,subs(g_conversion2,g_ext(pui))); ext_pui:=all_substitutions(ext_pui); pui:=subs(ext_f,subs(g_conversion2,pui)); # Now replace the RootOf's by all corresponding floating point numbers pui2:={}; for i in ext_pui do pui2:=pui2 union subs(i,pui) od; curven:={}; for i in pui2 do curven:=curven union {curve(subs(x=x*evalf( lcoeff(i[3],x)^(-1/degree(i[3],x))),i[1]) ,degree(i[3],x),t,eps,x)} od; if printlevel>1 then print(`Number of branches:`,nops(curven)) fi; if cols=[] then plots[tubeplot](curven,t=0..2*Pi,op(opt)) else curven:=[op(curven)]; plots[display]([seq(plots[tubeplot](curven[i],t=0..2*Pi,op(opt), color=cols[1+irem(i,nops(cols))]),i=1..nops(curven))]) fi end: curve:=proc(functie_in_x,d,t,eps,x) local f,fre,fim,factor; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; f:=expand(subs(x=eps*(cos(t)+I*sin(t)),functie_in_x)); fre:=coeff(f,I,0)/eps;fim:=coeff(f,I,1)/eps; factor:=1/(sqrt(1+fre^2+fim^2)-fre); [factor*cos(d*t),factor*sin(d*t),factor*fim]; end: #savelib('`algcurves/pl_allsub`','`algcurves/pl_curve`',\ 'plot_knot','`algcurves/plot_knot.m`'): #------------------------ filename: algcurves/genus --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/genus,v $ # $Notify: hoeij@sci.kun.nl $ macro( genus=`algcurves/genus` ): macro( g_normal=`algcurves/g_normal` ): # f is a polynomial in x and y genus:=proc(f,x,y) local F,d,v,i; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if args[nargs]=`irr check` then F:=evala(AFactor(args[1])); if not type(F,`+`) then i:=2; if type(F,`*`) then i:=0; for d in [op(F)] do if indets(d) intersect {x,y} <> {} then i:=i+1 fi od; fi; if i>1 then RETURN(args[1],`factored as`,F) fi; fi; RETURN(procname(args[1..nargs-1])) fi; if nargs<3 or (not type(x,name)) or (not type(y,name)) or (not type(f,polynom(anything,[x,y]))) then ERROR(`wrong number or type of arguments`) fi; readlib(`algcurves/g_expand`): # to load g_normal # The following is to ensure that degree(..) gives the right answer. F:=collect(numer(f),{x,y},distributed,g_normal); if degree(F,{x,y})=0 then ERROR(`degree is zero`) elif degree(F,x)=1 or degree(F,y)=1 then 0 elif degree(F,x)*degree(F,y)=0 then ERROR(`Input has degree >1 in one variable and degree 0 in the other`) elif degree(F,y)=2 then # We could call `algcurves/singularities` but the following is # likely to be faster: F:=collect(F/lcoeff(F,y),y,g_normal); F:=collect(subs(y=y-coeff(F,y,1)/2,F),y,g_normal); F:=coeff(F,y,0); # Take the factors with odd multiplicity F:=numer(F)*denom(F); floor(convert([-1,seq(`if`(irem(i[2],2)=0,0, degree(collect(i[1],x),x)),i=evala(Sqrfree(F))[2])],`+`)/2) elif degree(F,x)=2 or degree(F,x)+degree(lcoeff(F,x),y) 4 and what=`ratfunction` then i:=collect(numer(f),{x,y},distributed); j:=collect(denom(f),{x,y},distributed); n:=degree(i,{x,y})-degree(j,{x,y}); if n<0 then i:=i*z^(-n) elif n>0 then j:=j*z^n fi; RETURN(procname(i,x,y,z,`polynom`) /procname(j,x,y,z,`polynom`)) fi; n:=degree(f,{x,y}); convert([seq(seq(coeff(coeff(f,x,i),y,j)*x^i*y^j*z^(n-i-j) ,i=0..n-j),j=0..n)],`+`) end: #savelib ('`algcurves/homogeneous`','`algcurves/homogeneous.m`'): #------------------------ filename: algcurves/integral_basis --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/integral_basis,v $ # $Notify: hoeij@sci.kun.nl $ macro( integral_basis=`algcurves/integral_basis`, local_intbasis23=`algcurves/ib23`, double_factors=`algcurves/db_factors`, local_intbasis=`algcurves/local_ib`, ext_to_coeffs=`algcurves/e_to_coeff`, g_gcdex=`algcurves/gcdex` ): # This code uses the procedures in `algcurves/g_expand` and # `algcurves/puiseux` . These files will be read in the # integral_basis procedure # In Maple 5.4 there is the following bug in solve: # solve(y+1-RootOf(-3-_Z+2*_Z^2+_Z^3)-RootOf(-3-_Z+2*_Z^2+_Z^3)^2=0,y); # The following is a fix for this: macro( solve=readlib(`solve/linear`) ): macro( v_ext_m=`algcurves/v_ext_m`, lift_exp=`algcurves/lift_exp`, lift_exp_m1=`algcurves/lift_exp_m1`, truncate_subs=`algcurves/truncate_subs`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, `integral_basis/bound`=`algcurves/ib_bound`, Newtonpolygon=`algcurves/Newtonpolygon` ): macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): # Purpose: efficiently compute the local integral # basis at simple singularities # Input: f in Q(ext)[x,y], not necessarily monic in y. # subs(x=p,lcoeff(f,y))<>0 # p is algebraic over Q(ext), it is the root of P. # Input in rootof syntax # D is the valuation of the discriminant at x=p # D is 2 or 3 # Output: local integral basis, in RootOf syntax local_intbasis23:=proc(f,x,y,ext,p,P,D) global g_conversion1,g_conversion2,g_normal,algcurves; local f_local,v,m,n,i,N; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; n:=degree(f,y); f_local:=subs(g_conversion2,numer(subs(x=p,g_expand(rem(f,P,x),ext)))); userinfo(5,algcurves,`Doing a squarefree factorization of`,f_local); v:=evala(Sqrfree(f_local),'independent')[2]; userinfo(5,algcurves,`Done squarefree factorization`); for i in v do if has(i,y) then if type(m[i[2]],list) or i[2]>4 then ERROR(`Wrong input`) fi; m[i[2]]:=i fi od; if (D=2 and type(m[3],list)) or (D=2 and type(m[2],list) and degree(m[2][1],y)>1) or (D=3 and type(m[4],list)) or (D=3 and type(m[3],list) and type(m[2],list)) or (D=3 and type(m[2],list) and degree(m[2][1],y)>2) then userinfo(4,algcurves, `There are ramification points but no singularities at`, x=subs(g_conversion2,p)); RETURN([seq(y^i,i=0..n-1)]) fi; for i from 1 to 3 do if not type(m[i],list) then m[i]:=[1] fi od; if D=2 then N:=m[2][1] # location of the singularity elif m[3]<>[1] then N:=m[3][1] elif degree(m[2][1],y)=1 then N:=m[2][1] else # there is 1 ramification point and 1 singularity f_local:=subs(g_conversion2,numer(subs(x=p,g_expand(rem( collect(diff(f,x),x),P,x),ext)))); N:=numer(g_normal(m[2][1]/evala(Gcd(f_local,m[2][1])))) fi; if degree(N,y)<>1 then ERROR(`wrong degree`) fi; userinfo(4,algcurves,`Found one singularity at`, x=subs(g_conversion2,p)); if nargs>7 then # This option is used in the algorithm for finding the singularities RETURN(N) fi; if D<>2 and m[3]<>[1] then N:=N^2 fi; if D<>2 and m[3]=[1] and degree(m[2][1],y)<>1 then N:=N*m[2][1] fi; N:=subs(g_conversion1,evala(Expand(N*m[1][1]))); # Numerator of the integral basis N:=collect(N/lcoeff(N,y),y,g_normal); if degree(P,x)>1 then N:=subs(p=x,N) fi; [seq(y^i,i=0..n-2),N/P] end: # Input: a polynomial f in x # Both RootOf and rootof syntax. # ext: gives the ground field # Output: the factors that appear with multiplicity > 1, and their # multiplicities. Plus the multiplicity 1 factor, not factorized. double_factors:=proc(ff,x,ext) global algcurves,g_conversion1,g_conversion2; local f,v,i,df,j; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if ext<>[] and not has(ext,RootOf) then RETURN(subs(g_conversion1, procname(op(subs(g_conversion2,[args]))))) fi; userinfo(5,algcurves,`Doing a squarefree factorization of`,ff); f:=numer(ff); if v=[] then v:=sqrfree(f) else v:=evala(Sqrfree(f),'independent') fi; userinfo(5,algcurves,`Done squarefree factorization`); df:=NULL; for i in v[2] do if i[2]>1 then for j in g_factors(i[1],x,ext) do df:=df,[j[1],i[2]] od else df:=df,i fi od; [df] end: # Syntax: integral_basis(r) or integral_basis(r,x) where r is a RootOf # Or: integral_basis(f,x,y) # Can add an extra option specifying which factors of the discriminant # need to be considered (for computing a local integral basis) integral_basis:=proc() global algcurves,g_conversion1,g_conversion2,g_normal; local f,x,y,B,R,alfa,df,ext,extl,i,k,nulp,q,res,v; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`): readlib(`algcurves/puiseux`): alfa:=eval(args[1]); if nargs=1 then x:=op(indets(alfa,`name`) minus {_Z}) else x:=args[2] fi; if not (type(alfa,algfun(rational)) and type(x,`name`)) then ERROR(`Wrong arguments`) fi; if nargs>=3 and type(args[2],`name`) and type(args[3],`name`) then f:=args[1]; x:=args[2]; y:=args[3]; alfa:=y; extl:=g_ext(f); f:=subs(g_conversion1,f) else extl:=g_ext(op(alfa)); f:=subs(g_conversion1,op(alfa)); f:=subs(_Z=y,f) fi; f:=monic(f,x,y,extl,'q'); q:=eval(q); # f is the minimal polynomial of y over L(x), we made it monic if nargs=4 then # The purpose of the 4'th argument is to compute a local integral # basis. This argument must be a set of factors of the discriminant. # Each element of this set must be either an irreducible polynomial, # or must be of the form [irr poly, multiplicity of that poly in the # discriminant]. If the multiplicity is not specified, we set it # equal to 4: df:={seq(`if`(type(i,list),i,[i,4]),i=args[4])} # so that the procedure local_intbasis23 is not called. That procedure # can only be used if the multiplicity is either 2 or 3. elif nargs=3 and not type(args[3],`name`) then df:={seq(`if`(type(i,list),i,[i,4]),i=args[3])} else df:=discrim(f,y); if testeq(subs(g_conversion2,df)) and g_normal(df)=0 then ERROR(`input should be squarefree`) fi; df:=double_factors(df,x,extl) # df contains those factors k that appear more than once in the # discriminant discrim(f,y) fi; res:=[seq(y^i,i=0..degree(f,y)-1)]; for k in df do if k[2]>1 then nulp:=g_zero_of(k[1],x,'ext'); ext:=eval(ext); ext:=[ext,op(extl)]; userinfo(3,algcurves,`Computing local integral basis for the factor`, subs(g_conversion2,k[1])^k[2]); if k[2]=2 or k[2]=3 then B:=local_intbasis23(f,x,y,extl,nulp,op(k)) else B:=local_intbasis(f,x,y,ext,extl,nulp,k[1]) fi; userinfo(3,algcurves,`Done computing local integral basis.`); R:=NULL; for i from 1 to degree(f,y) do if not has(lcoeff(res[i],y),x) then R:=R,B[i] elif not has(lcoeff(B[i],y),x) then R:=R,res[i] else v:=g_gcdex(1/lcoeff(res[i],y),1/lcoeff(B[i],y),1,x,extl); R:=R,collect(v[2]*B[i]+v[3]*res[i],y,g_normal) fi od; res:=[R] fi od; subs(g_conversion2,subs(y=alfa*q,res)) end: local_intbasis:=proc(f,x,y,ext,extl,nulp,k) global algcurves,g_conversion1,g_conversion2; local j,j2,b,basis,d,equations,f_translated,found_something,i, max_power_k,max_v7,pl_transl,places,power_of_k,value_new_one, values_basis_in_places; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; f_translated:=g_expand(subs(x=x+nulp,f),ext); pl_transl:=[op(`puiseux/technical_answer`( f_translated,x,y,`integral basis`,ext))]; max_v7:=max(seq(i[7],i=pl_transl)); max_power_k:=floor(max_v7); userinfo(2,algcurves,`Maximum number of factors`,k, `in the denominator is`,max_power_k); if max_power_k=0 then RETURN([seq(y^i,i=0..degree(f,y)-1)]) fi; places:=[seq([i[1],degree(i[3],x),i[3],i[4]],i=pl_transl)]; # places[i] contains: (i=1..number of puiseux expansions) # 1) The puiseux expansion (but translated x=x+nulp) # 2) The ramification index # 3) The ramification # In order to interpret a normal x we substitute # subs(x=i[3]+nulp,..) # 4) The algebraic extension belonging to this puiseux expansion values_basis_in_places[1]:=[seq(1,i=pl_transl)]; power_of_k:=1; basis:=[1]; for d from 2 to degree(f,y) do basis:=[op(basis),y*basis[d-1]]; # This basis is our first guess # Now we compute the values of this new basis-element in all # places (that is, we substitute y=puiseux expansions) values_basis_in_places[d]:=[seq(truncate (values_basis_in_places[d-1][j]*places[j][1] ,places[j][2]*max_power_k,x ,places[j][4]),j=1..nops(places))]; # using the values of this new basis-element in the places, # we will see if we can find an integral, with a bigger denominator found_something:=true; while found_something and power_of_k<=max_power_k do for i from 1 to d-1 do b[i]:=evaln(b[i]) od; b[d]:=1; # Now we compute the values of basis[1]*b[1]+..+basis[d]*b[d] # in the places # and we try to put an extra factor k in the denominator value_new_one:=[seq(truncate(convert([seq(b[j2] *values_basis_in_places[j2][j],j2=1..d)],`+`) ,places[j][2]*power_of_k,x ,places[j][4]),j=1..nops(places))]; # All coefficients of powers of x less than power_of_k must # be zero, in order for this new_one to be dividable by # k^power_of_k. So we find equations by taking the # remainder equations:={seq(ext_to_coeffs(j,ext),j=value_new_one)}; equations:={seq(coeffs(j,x),j=equations)}; equations:=subs(g_conversion1,evala({solve( subs(g_conversion2,{seq(expand(numer(i)),i=equations)}) ,{seq(b[i],i=1..d-1)})},'independent')); # Now we know what values b[1] .. b[d] must have if equations={} then found_something:=false else if ext<>extl then equations:=subs(nulp=x,equations) fi; assign(op(equations)); # Now basis[1]*b[1]+..basis[d]*b[d] is dividable by k # In the following for_do_od we compute the values of # basis[1]*b[1]+..+basis[d]*b[d] in all places. values_basis_in_places[d]:=[seq( truncate(convert([seq(subs(x=places[j][3]+nulp ,eval(b[j2]))*values_basis_in_places[j2][j],j2= 1..d)],`+`),places[j][2] *max_power_k,x,places[j][4]) ,j=1..nops(places))]; # Now we will put (basis[1]*b[1]+..+basis[d]*b[d])/k # in the basis basis:=subsop(d=normal( convert([seq(eval(b[j2])*basis[j2],j2=1..d)],`+`) /k),basis); # Now we should divide values_basis_in_places[.. ,d] # by k, but instead we multiply all other # values_basis_in_places[.., <>d] by k, this will # be done in the following nested for_do_od for i from 1 to d-1 do values_basis_in_places[i]:=[seq( truncate(values_basis_in_places[i][ j]*subs(x=places[j][3]+nulp,k),places[j][2]*max_power_k ,x,places[j][4]),j=1..nops(places))] od; # Now we increase power_of_k, so we will try to put # more factors k in the denominator power_of_k:=power_of_k+1 fi od od; basis end: # ext_to_coeffs does basically: coeffs(expression,RootOf( .. )) # We need this procedure because coeffs doesn't always work this way. ext_to_coeffs:=proc(a,ext) global g_conversion2; local i,aa; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; aa:=(indets(a) minus indets(ext)) intersect {seq(rootof.i,i=0..nops(g_conversion2))}; coeffs(a,[op(aa)]) end: g_gcdex:=proc(a,b,c,x,ext) global g_conversion1,g_conversion2; local r,ss,tt,s,t,i,j; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; r:=subs(g_conversion2,[a,b,c]); if ext=[] then r:=gcdex(op(r),x,'ss','tt') else r:=subs(g_conversion1,evala(Gcdex(r[1],r[2],x,'ss','tt'))); if g_normal(r-c)<>0 then ERROR(`found wrong gcd`) fi fi; s:=subs(g_conversion1,ss); t:=subs(g_conversion1,tt); [c,s,t] end: #savelib ('local_intbasis23','double_factors','integral_basis',\ 'local_intbasis','ext_to_coeffs','g_gcdex','`algcurves/integral_basis.m`'): #------------------------ filename: algcurves/singularities --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/singularities,v $ # $Notify: hoeij@sci.kun.nl $ macro( singularities=`algcurves/singularities`, find_points=`algcurves/find_points`, degree_ext=`algcurves/degree_ext` ): # This code uses the procedures in `algcurves/g_expand`, # `algcurves/puiseux` and `algcurves/integral_basis`. These files will # be read in the singularities procedure macro( homogeneous=`algcurves/homogeneous` ): macro( local_intbasis23=`algcurves/ib23`, double_factors=`algcurves/db_factors`, local_intbasis=`algcurves/local_ib`, ext_to_coeffs=`algcurves/e_to_coeff`, g_gcdex=`algcurves/gcdex` ): macro( v_ext_m=`algcurves/v_ext_m`, lift_exp=`algcurves/lift_exp`, lift_exp_m1=`algcurves/lift_exp_m1`, truncate_subs=`algcurves/truncate_subs`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, `integral_basis/bound`=`algcurves/ib_bound`, Newtonpolygon=`algcurves/Newtonpolygon` ): macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): # This is used as a global variable for making the curve homogeneous. The # reason for using a global variable here is that the options remembers in # various places in this package will be more useful this way. macro( z=`algcurves/z` ): # Input: f a polynomial (in expanded form) in x and y. # rootof syntax # Output: the singularities, plus some other info that is used # in odd_singularity_on_C # Furthermore: the multiplicity and the delta invariant of these points find_points:=proc(f,x,y,ext,where) global g_conversion1,z; local ff,n,p,delta,extp,f_translated,i,j,k,mult,pl_transl,res; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; n:=degree(f,{x,y}); if nargs=4 or (nargs=5 and has(args[5],`search points on`)) then j:=discrim(f,y); # First a quick testeq, then if necessary a thorough test g_normal=0 if testeq(subs(g_conversion2,j)) and g_normal(j)=0 then ERROR(`Input should be squarefree`) fi; {seq(op(procname(args[1..4],i,args[5..nargs])),i= [op(double_factors(j,x,ext)),infinity])} elif where=infinity then k:=-1; while k=-1 or degree(subs(z=0,ff),x)0 and not member([0,0],{seq([i[1][2],i[1][3]],i=res)}) then res union {[[1,0,0],1,0]} else res fi elif where[2]<2 then {} else p:=g_zero_of(where[1],x,'extp'); extp:=eval(extp); if member(where[2],{2,3}) and not # If the following is the case, then use the other method because # the other method (although slower) will be efficient enough # (because degree_ext is small) and with some luck it may give # an odd degree place. (has([args],`search points on`) and member(degree_ext(p,ext), args[6][2])) then i:=local_intbasis23(f,x,y,ext,p,op(where),z); if not type(i,list) then # i should be of degree 1 in y subs(g_conversion1, {[[p,g_normal(-coeff(i,y,0)/coeff(i,y,1)),1],2,1]}) else # perhaps this line is useful in odd_singularity_on_C, if # where[2]=3 we have a place of algebraic odd-degree {[[p,0,1],0,0,where[2]]} fi else res:=NULL; extp:=[extp,op(ext)]; f_translated:=g_expand(subs(x=x+p,f),extp); pl_transl:=[op(`puiseux/technical_answer`(f_translated,x,y ,0,extp))]; if irem(degree_ext(p,ext),2)=1 then # add info that could be useful in odd_singularity_on_C for i in pl_transl do if irem(i[6],2)=1 then res:=[[p,0,1],0,0,3]; break fi od fi; ff:={seq(subs(x=0,i[1]),i=pl_transl)}; # the y values for i in ff do # i is y value mult:=0; delta:=0; for j in pl_transl do if subs(x=0,j[1])=i then mult:=mult+degree_ext(j,[i,op(extp)])*min( # Minimum of the valuation of x and y in this place degree(j[3],x), max(ldegree(j[1]-subs(x=0,j[1]),x),1) ); delta:=delta+1/2*degree_ext(j,[i,op(extp)])* (degree(j[3],x)*j[7]-degree(j[3],x)+1) fi od; if not (type(mult,integer) and type(delta,integer)) then ERROR(`found wrong values`) fi; res:=res,[[p,i,1],mult,delta] od; {res} fi fi end: # Gives the algebraic degree of a over b degree_ext:=proc(aa,bb) global g_conversion2; local a,b,v,i,all,d,var; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; a:=subs(g_conversion2,aa); b:=subs(g_conversion2,bb); v:=indets(a,RootOf) minus indets(b,RootOf); all:=[op(indets([a,b],RootOf))]; all:={seq(all[i]=var[i],i=1..nops(all))}; d:=1; for i in v do d:=d*degree(subs(all,op(i)),_Z) od; d end: # Input: f is a polynomial in x and y # Output is a list of the following lists: # [ [x,y,z], multiplicity, delta invariant ] singularities:=proc(f,x,y) global g_conversion1,g_conversion2; local F,v,i,res,ext; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if nargs=3 and has(evala(Content(f,y)),x) then # If f has a component which is independent of y, this component # will be ignored by `algcurves/find_points`, but I want to # include it in the computation. This is a fix for this case v:=procname(subs(x=x+y,f),x,y); RETURN({seq([[i[1][1]+i[1][2],i[1][2],i[1][3]],i[2],i[3]],i=v)}) fi; readlib(`algcurves/g_expand`): readlib(`algcurves/puiseux`): readlib(`algcurves/integral_basis`): readlib(`algcurves/homogeneous`): ext:=g_ext([args]); F:=g_expand(subs(g_conversion1,f),ext); v:=find_points(F,x,y,ext); res:=NULL; for i in v do if i[2]>1 then res:=res,i fi od; subs(g_conversion2,{res}) end: #savelib ('find_points','degree_ext',\ 'singularities','`algcurves/singularities.m`'): #------------------------ filename: algcurves/ratpar --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/ratpar,v $ # $Notify: hoeij@sci.kun.nl $ # #--> 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( solve=readlib(`solve/linear`) ); macro( ratpar=`algcurves/ratpar`, 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( homogeneous=`algcurves/homogeneous` ): # from singularities macro( singularities=`algcurves/singularities`, find_points=`algcurves/find_points`, degree_ext=`algcurves/degree_ext` ): # from integral_basis macro( integral_basis=`algcurves/integral_basis`, local_intbasis23=`algcurves/ib23`, double_factors=`algcurves/db_factors`, local_intbasis=`algcurves/local_ib`, ext_to_coeffs=`algcurves/e_to_coeff`, g_gcdex=`algcurves/gcdex` ): # from puiseux macro( puiseux=`algcurves/puiseux`, v_ext_m=`algcurves/v_ext_m`, lift_exp=`algcurves/lift_exp`, lift_exp_m1=`algcurves/lift_exp_m1`, truncate_subs=`algcurves/truncate_subs`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, `integral_basis/bound`=`algcurves/ib_bound`, Newtonpolygon=`algcurves/Newtonpolygon` ): # From g_expand macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): macro( 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.`; f:=collect(ff,{x,y},distributed,g_normal); 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([g_normal(res[1]/res[2]),g_normal(1/res[2])]) fi; res:=procname(subs(x=x+y,f),args[2..nargs]); RETURN([g_normal(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 g_normal(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([g_normal(res[1]/(res[1]+1)), g_normal(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 convert([seq(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(expand(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))[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 g_gcd(subs(y=i0,f),d1)<>1 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:=convert([seq(co[i]*evala(Rem(collect( convert([seq(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) ,convert([seq(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:=convert([seq(denom(i),i=integral_basis(f,x,y))],`*`); # The following contains the relevant factors of the # discriminant for frac_integral_a rf:=numer(g_normal(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:=convert([seq(seq(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(g_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.`; 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; # This may need some modification of charpol option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; 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]; 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.`; n:=degree(f,y); # Not necessarily true: a[n,2]:=1; ansatz:=convert([seq(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; if g_gcd(op(subs(x=j,[denom(x0v),f])))<>1 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:=expand(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.`; 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() # = bug 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.`; 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:=convert([seq(seq(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.`; 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.`; 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.`; 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:=convert([seq(seq(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.`; P:=rat_point_on_C2(args); if P=failed then 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.`; 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() 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.`; 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; # Take a point at infinity and compute a parametrization option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; s:=RootOf(g_factors(subs(t=0,x=1,homogeneous(f,x,y,t,`polynom`)) ,y,g_ext_r(f))[1][1],y); express_in_p(f,x,y,t,y-x*s,1) 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.`; 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.`; 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:=convert([seq(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 g_gcd(op(R))<>1 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.`; if not member(y,indets(f)) then RETURN(g_solve(f)) 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:=convert([seq(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,expand(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.`; 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:=convert([seq(qc[i]*y^i,i=0..d-1)],`+`); q:=g_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; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; 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:=convert([seq(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.`; 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',\ 'L_inf','g_gcd','`algcurves/ratpar.m`'): #------------------------ filename: algcurves/iss94 --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/iss94,v $ # $Notify: hoeij@sci.kun.nl $ macro( solve=readlib(`solve/linear`) ); macro( iss94=`algcurves/iss94`, function_with_one_pole=`algcurves/f_with_1_p` ): macro( homogeneous=`algcurves/homogeneous` ): macro( singularities=`algcurves/singularities`, find_points=`algcurves/find_points`, degree_ext=`algcurves/degree_ext` ): macro( integral_basis=`algcurves/integral_basis`, local_intbasis23=`algcurves/ib23`, double_factors=`algcurves/db_factors`, local_intbasis=`algcurves/local_ib`, ext_to_coeffs=`algcurves/e_to_coeff`, g_gcdex=`algcurves/gcdex` ): macro( puiseux=`algcurves/puiseux`, v_ext_m=`algcurves/v_ext_m`, lift_exp=`algcurves/lift_exp`, lift_exp_m1=`algcurves/lift_exp_m1`, truncate_subs=`algcurves/truncate_subs`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, `integral_basis/bound`=`algcurves/ib_bound`, Newtonpolygon=`algcurves/Newtonpolygon` ): macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): macro( ratpar=`algcurves/ratpar`, 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` ): macro( z=`algcurves/z` ): # f is a polynomial in x and y # f must be irreducible over Qbar # If f is reducible this procedure does not work properly. iss94:=proc(f,x,y,t,point) global z; local d,fu; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; d:=degree(f,{x,y}); if coeff(f,y,d)=0 then if coeff(f,x,d)<>0 then # switch x and y d:=procname(f,y,x,t,[point[2],point[1],point[3]]); RETURN([d[2],d[1]]) elif coeff(coeff(f,x,0),y,0)<>0 then # switch z and y d:=procname(subs(y=1,z=y,homogeneous(f,x,y,z, `polynom`)),x,y,t,[point[1],point[3],point[2]]); RETURN(map(g_normal,[d[1]/d[2],1/d[2]])) else d:=procname(collect(subs(x=x+y,f),{x,y},distributed,g_normal) ,x,y,t,[point[1]-point[2],point[2],point[3]]); RETURN([g_normal(d[1]+d[2]),d[2]]) fi fi; # coeff(f,y,d)<>0, so [0,1,0] is not a point on the curve anymore. # Now we compute a function that has a pole in point. If the # point is finite, this function has no other poles in finite # otherwise it has no other poles in infinity. # It may still have other poles, they will be removed by # the procedure function_with_one_pole if point[3]=0 then fu:=[(homogeneous(evala(Quo(subs({t=0,x=1} ,homogeneous(f,x,y,t,`polynom`)),y-point[2]/point[1],y)) ,y,t,x,`polynom`)),`infinity`]; else fu:=[evala(Quo(evala(Expand(subs(x=point[1]/point[3],f))) ,y-point[2]/point[3],y))/(x-point[1]/point[3]),`finite`] fi; fu:=function_with_one_pole(f,x,y,fu); userinfo(2,algcurves,`Computed a generator of the function field`); express_in_p(f,x,y,t,numer(fu),denom(fu)) end: # The following procedure computes a parameter. A parameter is a function # that has only 1 pole, with multiplicity 1. Then k(p)=k(x,y), p generates # the function field. The pole will be in point. In fact, as the syntax # is now, point is not a point, but a function with a pole in a certain point. # # point=[function,`finite`] or [function,`infinity`]. This function must # have one pole in finite or infinity. This function is the start for our # search for a parameter. We add an ansatz to this point. Since this function # is already good in a part of the plane (either `finite` or `infinity`) # we will change no poles in that part of the plane, and try to remove the # poles in the other part. # # The point [0,1,0] must not be a point on the curve. The reason is that # we have divided the projective plane in 2 affine parts, a line and # a plane. Doing this, there remains 1 point, and therefore that point must # not be on the curve. # # Note that the way we divide the plane in parts is not really relevant # for the method. The fact that one part is only a line does not give much # asymmetry in the algorithm. The important thing about these parts is that # they are affine, and not projective. We need this to apply integral basis. # # We will compute an integral basis for the finite part of the plane, and a # local integral basis for the line `infinity`. Using these integral basis # we can find linear equations stating that an ansatz, or an ansatz + a # given function is integral (i.e. no poles) in a part of the plane (the # part `finite` and `infinity`. function_with_one_pole:=proc(f,x,y,point) local d,j,f_infty,ib1,ib2,den1,den2,d1,d2,ansatz,a,f1,f2,equations, ext,i,zero,point1,deg_inf; global g_conversion1,g_conversion2,`algcurves/residue`,z; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; d:=degree(f,{x,y}); f_infty:=subs(x=1,homogeneous(f,x,y,z,`polynom`)); ib1:=integral_basis(f,x,y); ib2:=integral_basis(f_infty,z,y,{[z,4]}); den1:=expand(denom(ib1[d])); den2:=expand(denom(ib2[d])); if point[2]=`infinity genus1` then # Syntax for genus1, the function point1 may not be changed # so we only determine the upper bound for the denominators now point1:=point[1]; # den1:=evala(Lcm(den1,denom(point1))) den1:=numer(g_normal(den1*denom(point1) /evala(Gcm(den1,denom(point1))))) elif point[2]=`infinity` then deg_inf:=degree(den1,x)-degree(point[1],{x,y})+1; point1:=x^deg_inf*point[1]/den1; if deg_inf<0 then den1:=evala(Expand(den1*x^(-deg_inf))) fi else point1:=point[1] fi; # den1 is an upperbound for denominators d1:=degree(den1,x); d2:=ldegree(den2,z); d2:=max(d2,degree(expand(numer(point1)),{x,y}) -degree(expand(denom(point1)),{x,y})); # d2 is an upperbound for the denominators in infinity ansatz:=convert([seq(seq(a[i,j]*x^i*y^j,i=0..d1+d2-j) ,j=0..d-1)],`+`)/den1; if point[2]=`finite` then f1:=ansatz; f2:=ansatz+point1 else f1:=ansatz+point1; f2:=ansatz fi; # Now we solve the following set of linear equations: # {f1 has no poles in finite, and f2 has no poles in infinity} ext:=g_ext([f,f1,f2]); f1:=g_normal(f1); den1:=expand(denom(f1)); f1:=subs(g_conversion1,expand(numer(f1))); f2:=g_normal(subs(x=1,homogeneous(f2,x,y,z,`ratfunction`))); den2:=expand(denom(f2)); den2:=z^ldegree(den2,z); f2:=subs(g_conversion1,expand(numer(f2))); ib1:=subs(g_conversion1,[seq(g_expand(g_normal(i*den1),ext),i=ib1)]); ib2:=subs(g_conversion1,[seq(g_expand(g_normal(i*den2),ext),i=ib2)]); zero:=f1; for i from d-1 by -1 to 0 do while degree(coeff(zero,y,i),x) >=degree(lcoeff(ib1[i+1],y),x) and coeff(zero,y,i)<>0 do zero:=g_expand(zero-ib1[i+1]* g_normal(lcoeff(coeff(zero,y,i),x)/lcoeff(lcoeff(ib1[i+1],y),x)) *x^(degree(coeff(zero,y,i),x)-degree(lcoeff(ib1[i+1],y),x)) ,ext); if (indets([f,point],`name`) minus {x,y})<>{} then zero:=expand(collect(zero,{x,y},distributed,g_normal)) fi od od; equations:={a[d1,0],coeffs(zero,[x,y])}; # These equations state that f1 is integral. We've added and 1 extra # linear condition, to make the solution unique. We still have to # add the equations for infinity: zero:=f2; for i from d-1 by -1 to 0 do while degree(coeff(zero,y,i),z) >=degree(lcoeff(ib2[i+1],y),z) and coeff(zero,y,i)<>0 do zero:=g_expand(zero-ib2[i+1]* g_normal(lcoeff(coeff(zero,y,i),z)/lcoeff(lcoeff(ib2[i+1],y),z)) *z^(degree(coeff(zero,y,i),z)-degree(lcoeff(ib2[i+1],y),z)) ,ext); if (indets([f,point],`name`) minus {x,y})<>{} then zero:=expand(collect(zero,{x,y},distributed,g_normal)) fi od od; equations:=subs(g_conversion2,{op(equations),coeffs(zero,[y,z])}); d:=g_normal(subs(solve(equations,indets(equations,`name`) minus indets(f)),ansatz+point1)); if indets(d,`name`) minus indets([args]) <> {} then ERROR(`can not find parameter`,f=evala(AFactor(f))) fi; d end: #savelib ('iss94','function_with_one_pole','`algcurves/iss94.m`'): #------------------------ filename: algcurves/parametrization --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/parametrization,v $ # $Notify: hoeij@sci.kun.nl $ macro( parametrization=`algcurves/parametrization`, find_rat_point=`algcurves/find_rat_point` ): # Input: an irreducible algebraic curve f in x and y with genus 0. # Output: a parametrization [X(s),Y(s)]. Here X(s) and Y(s) are rational # functions in s, giving a bijective morphims from the projective line to # the algebraic curve f. parametrization:=proc(F,x,y,s) local f,t,p; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`): readlib(`algcurves/puiseux`): readlib(`algcurves/integral_basis`): readlib(`algcurves/singularities`): readlib(`algcurves/ratpar`): readlib(`algcurves/homogeneous`): # Make sure that the coefficients are evala normalized: f:=collect(F,{x,y},distributed,`algcurves/g_normal`); if degree(f,{x,y})<4 or degree(f,x)<2 or degree(f,y)<2 then # there is special code for these cases in the file ratpar p:=`algcurves/ratpar`(f,x,y,s,t,`parametrize by line`) else p:=find_rat_point(f,x,y); if p=`genus is not 0` then ERROR(p) fi; # by giving an extra argument one can choose for the second method if p<>failed and nargs=4 then # when given a rational regular point, the method in # ISSAC'94, is usually faster, so: userinfo(1,algcurves,`Using the rational point`,p); p:=readlib(`algcurves/iss94`)(f,x,y,s,p) else # In the cases where this heuristic find_rat_point can not # find a rational regular point ratpar is usually quicker. p:=`algcurves/ratpar`(f,x,y,s,t,`parametrize by line`) fi fi; # verifying correctness... t:=12345; while testeq(subs(s=t,denom(p[1])*denom(p[2]))) do t:=t+1 od; if not testeq(subs(x=p[1],y=p[2],s=t,f)) then ERROR(`found wrong answer`) fi; userinfo(1,algcurves,`verified correctness OK.`); p # See also http://www-math.sci.kun.nl/math/compalg/IntBasis/ # for papers, more documentation and examples. end: # Try to find a regular point on the curve # For use only in the procedure parametrization, for deciding # which method to take. find_rat_point:=proc(F,x,y) local vp,i,ext,f,d; global `algcurves/g_conversion1`,`algcurves/g_conversion2`; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; ext:=`algcurves/g_ext`(F); f:=subs(`algcurves/g_conversion1`,F); vp:=`algcurves/find_points`(f,x,y,ext,[`search points on`,{1}]); d:=degree(f,{x,y}); if (d-1)*(d-2)/2-convert([seq(i[3]*`algcurves/degree_ext`(i,f),i=vp)] ,`+`)<>0 then RETURN(`genus is not 0`) fi; d:={}; for i in vp do if `algcurves/degree_ext`(i,f)=1 and i[2]=1 then d:=d union {i[1]} fi od; if d<>{} then # Take the smallest one, that gives a nicer parametrization vp:=min(seq(length(i),i=d)); for i in d do if length(i)=vp then RETURN(subs(`algcurves/g_conversion2`,i)) fi od fi; failed end: #savelib ('parametrization','find_rat_point',\ '`algcurves/parametrization.m`'): #------------------------ filename: algcurves/genus1 --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/genus1,v $ # $Notify: hoeij@sci.kun.nl $ macro( genus1=`algcurves/genus1`, genus1_compute_x=`algcurves/g1_comp_x` ): macro( solve=readlib(`solve/linear`) ); macro( iss94=`algcurves/iss94`, function_with_one_pole=`algcurves/f_with_1_p` ): macro( homogeneous=`algcurves/homogeneous` ): macro( singularities=`algcurves/singularities`, find_points=`algcurves/find_points`, degree_ext=`algcurves/degree_ext` ): macro( integral_basis=`algcurves/integral_basis`, local_intbasis23=`algcurves/ib23`, double_factors=`algcurves/db_factors`, local_intbasis=`algcurves/local_ib`, ext_to_coeffs=`algcurves/e_to_coeff`, g_gcdex=`algcurves/gcdex` ): macro( puiseux=`algcurves/puiseux`, v_ext_m=`algcurves/v_ext_m`, lift_exp=`algcurves/lift_exp`, lift_exp_m1=`algcurves/lift_exp_m1`, truncate_subs=`algcurves/truncate_subs`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, `integral_basis/bound`=`algcurves/ib_bound`, Newtonpolygon=`algcurves/Newtonpolygon` ): macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): macro( ratpar=`algcurves/ratpar`, 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` ): # Input: a polynomial f in x and y with genus 1, the variables x,y,x0,y0 # and a regular point on the curve. # f must be irreducible over Qbar # If the last argument is `no inverse` then the inverse of the isomorphism # is not computed (i.e. the image of x and y is not computed). # If [0,1,0] is a point on the curve then: # - with the option `no inverse` the output is not necessarily # a polynomial in y. # - without the option `no inverse` it does become a polynomial # by doing a division in the function field # Output: a list containing: # f0 such that Qbar(x)[y]/(f) is isomorphic with Qbar(x0)[y0]/(f0) # The image of x0 under this isomorphism # The image of y0 # The image of x under the inverse isomorphism # The image of y # f0 is of the form y0^2 + polynomial of degree 3 in x0^2 genus1:=proc(ff,x,y,x0,y0,point) global `algcurves/residue`; local a,d,t,v,i,f,f0,x0v,y0v,ansatz,p2,p3,p; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; f:=collect(ff,{x,y},distributed,g_normal); #if nargs<6 or not type(point,list) then # # No point was specified, compute one: # RETURN(procname(args[1..5],find_a_point(f,x,y),args[6..nargs])) if not member(`no inverse`,{args}) then v:=procname( # remove the `j invariant` argument, otherwise the options remember # in this procedure does not help. `if`(args[nargs]=`j invariant`,args[1..nargs-1],args) ,`no inverse`); f0:=v[1]; if member(`j invariant`,{args}) then f0:=subs(y0=0,f0); f0:=f0/lcoeff(f0,x0); a:=coeff(f0,x0,1); d:=coeff(f0,x0,0); RETURN(evala(6912*a^3/(27*d^2+4*a^3))) fi; # Remove possible y's in the denominators: v:=seq(FFDIV(f,x,y,v[i]),i=2..3); # v=x0v,y0v RETURN([f0,v,genus1_compute_x(v,x,y,f,f0,x0,y0), genus1_compute_x(v,y,x,f,f0,x0,y0)]) fi; d:=degree(f,{x,y}); if coeff(f,y,d)=0 then if coeff(f,x,d)=0 then # [0,1,0] is a point on the curve, we want to avoid this v:=[point[1]-point[2],point[2],point[3]],args[7..nargs]; # apply transformation on f to remove the point [0,1,0] v:=procname(subs(x=x+y,f),args[2..5],v); if nops(v)>1 then # Do the reverse transformation i:=subs(x=x-y,[v[2],v[3]]); if nops(v)>3 then v:=[v[1],op(i),v[4]+v[5],v[5]] else v:=[v[1],op(i)] fi fi; RETURN(v) else # [0,1,0] is a point on the curve, we want to avoid this v:=[point[2],point[1],point[3]],args[7..nargs]; v:=procname(subs({x=y,y=x},f),args[2..5],v); if nops(v)>1 then # Do the reverse transformation i:=subs({x=y,y=x},[v[2],v[3]]); if nops(v)>3 then v:=[v[1],op(i),v[5],v[4]] else v:=[v[1],op(i)] fi fi; RETURN(v) fi fi; # Now coeff(f,y,d)<>0, so [0,1,0] is not a point on the curve # Compute a function with a pole in the given point: if point[3]=0 then # point in infinity p:=[(homogeneous(evala(Quo(subs({t=0,x=1} ,homogeneous(f,x,y,t,`polynom`)),y-point[2]/point[1],y)) ,y,t,x,`polynom`))/x^(degree(f,y)-2),`infinity genus1`] else p:=[evala(Quo(evala(Expand(subs(x=point[1]/point[3],f))) ,y-point[2]/point[3],y))/(x-point[1]/point[3]),`finite`] fi; # Now make it a double pole with an indeterminate as residue p2:=[evala(Normal(rem(expand(-p[1]^2+`algcurves/residue` *p[1]),f,y))),p[2]]; # Compute a function with this pole and no other poles: x0v:=function_with_one_pole(f,x,y,p2); p3:=[evala(Normal(rem(expand(-p[1]*x0v+`algcurves/residue` *p[1]),f,y))),p[2]]; # Now compute a function with pole order 3: y0v:=function_with_one_pole(f,x,y,p3); # look for a relation f0 between x0v and y0v using an ansatz: i:=0; ansatz:=y0^2+x0^3+a[1]*y0+a[2]*x0*y0+a[3]+a[4]*x0+a[5]*x0^2; while has(ansatz,{a[1],a[2],a[3],a[4],a[5]}) do i:=i+1; # avoid certain values i for x while g_normal(subs(x=i,denom(x0v)*denom(y0v)))=0 do i:=i+1 od; # find linear equations in a[1] .. a[5] ansatz:=subs(solve({coeffs(expand(rem(expand(numer( subs(x0=x0v,y0=y0v,x=i,ansatz))),expand(subs(x=i,f)),y)),y)} ,{a[1],a[2],a[3],a[4],a[5]}),ansatz) od; # repeat until all a[i] are determined (usually in the 1'st step) f0:=ansatz; # Now normalize f0 a little further: y0v:=g_normal(y0v+subs(x0=x0v,coeff(f0,y0,1))/2); # clear the coefficient of y0^1: f0:=collect(subs(y0=y0-coeff(f0,y0,1)/2,f0), [x0,y0],distributed,g_normal); x0v:=g_normal(x0v+coeff(f0,x0,2)/3); # clear the coefficient of x0^2: f0:=collect(subs(x0=x0-coeff(f0,x0,2)/3,f0), [x0,y0],distributed,g_normal); if has(`Weierstrass`,{args}) then # Weierstrass normal form x0v:=-x0v; y0v:=2*y0v; f0:=y0^2+4*subs(x0=-x0,y0=0,f0) fi; [f0,x0v,y0v] end: # Express x as an expression in x0 and y0. # f0 = algebraic relation between x0 and y0 # f = algebraic relation between x and y # x0v = image of x0 in Qbar(x,y) # y0v = image of x0 in Qbar(x,y) genus1_compute_x:=proc(x0v,y0v,x,y,f,f0,x0,y0) local cp,result,v,i; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; cp:=genus1_charpol(x0v,x,y,f,x0); result:=-coeff(cp,x,1)/2; cp:=g_normal((x^2-subs(x=x+result,cp))/coeff(f0,y0,0)); if cp<>0 then # Now add y0*sqrt(cp) to result: v:=[seq(evala(Sqrfree(i))[2],i=[numer(cp),denom(cp)])]; v:=g_normal(convert([seq(i[1]^(i[2]/2),i=v[1])],`*`)/ convert([seq(i[1]^(i[2]/2),i=v[2])],`*`)); cp:=factors(numer(x^2+g_normal(cp/v^2)) ,indets([args],RootOf))[2]; if has(cp,x0) then bug() fi; for i in cp do if has(i,x) and degree(i[1],x)=1 then v:=v*subs(solve({i[1]},{x}),x)*y0; result:=result+v,result-v; break fi od; fi; i:=0; result:={result}; while nops(result)>1 do for v in result do i:=i+1; # avoid certain values i for x while g_normal(subs(x=i,denom(x0v)*denom(y0v)))=0 do i:=i+1 od; if g_normal(rem(expand(numer(subs(x0=x0v,y0=y0v,x=i,numer(v)) -i*subs(x0=x0v,y0=y0v,x=i,denom(v)))),subs(x=i,f),y))<>0 then result:=result minus {v}; break fi od od; op(result) end: #savelib ('genus1','genus1_compute_x','`algcurves/genus1.m`'): #------------------------ filename: algcurves/Weierstrassform --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/Weierstrassform,v $ # $Notify: hoeij@sci.kun.nl $ macro( Weierstrassform=`algcurves/Weierstrassform`, find_a_point=`algcurves/find_a_point` ): macro( singularities=`algcurves/singularities`, find_points=`algcurves/find_points`, degree_ext=`algcurves/degree_ext` ): macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, normal_tcoeff=`algcurves/normal_tcoeff`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, rootof=`algcurves/rootof`, g_ext=`algcurves/g_ext`, g_ext_r=`algcurves/g_ext_r`, truncate=`algcurves/truncate` ): Weierstrassform:=proc(F,x,y,x0,y0) global `algcurves/x0`,`algcurves/y0`; local f,p; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`): readlib(`algcurves/puiseux`): readlib(`algcurves/integral_basis`): readlib(`algcurves/singularities`): readlib(`algcurves/homogeneous`): readlib(`algcurves/iss94`): readlib(`algcurves/genus1`): readlib(`algcurves/ratpar`): # Make sure that the coefficients are evala normalized: f:=collect(F,{x,y},distributed,g_normal); if nargs<6 or not type(args[6],list) then p:=find_a_point(f,x,y); if p=`genus is not 1` then ERROR(p) fi else p:=NULL fi; # The reason for using this global `algcurves/x0` is that the # options remember in `algcurves/genus1` will be more effective # this way. subs(`algcurves/x0`=x0,`algcurves/y0`=y0,`algcurves/genus1`( f,x,y,`algcurves/x0`,`algcurves/y0`,p,args[6..nargs])) end: # Find a regular point on the curve # For use only in the procedure Weierstrassform, for taking # a regular point and checking g=1. find_a_point:=proc(F,x,y) local vp,i,ext,f,d,m; global g_conversion1,g_conversion2; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; ext:=g_ext(F); f:=subs(g_conversion1,F); vp:=find_points(f,x,y,ext,[`search points on`,{1,2,3}]); d:=degree(f,{x,y}); if (d-1)*(d-2)/2-convert([seq(i[3]*degree_ext(i,f),i=vp)] ,`+`)<>1 then RETURN(`genus is not 1`) fi; # Add a regular point, just in case we found no regular points yet. i:=0; while testeq(discrim(subs(x=i,F),y)) do i:=i+1 od; vp:=vp union {seq([[i,RootOf(d[1],y),1],1,0], d=g_factors(subs(x=i,f),y,g_ext_r(F)))}; d:={}; m:=min(seq(`if`(i[2]=1,degree_ext(i,f),infinity),i=vp)); for i in vp do if degree_ext(i,f)=m and i[2]=1 then d:=d union {i[1]} fi od; # Take the smallest one. m:=min(seq(length(i),i=d)); for i in d do if length(i)=m then RETURN(subs(g_conversion2,i)) fi od; # The following can't happen, d can not be empty because # there should at least be one regular point in there (we # placed it there in the beginning of this procedure. ERROR() end: #savelib ('Weierstrassform','find_a_point','`algcurves/Weierstrassform.m`'): #------------------------ filename: algcurves/j_invariant --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/j_invariant,v $ # $Notify: hoeij@sci.kun.nl $ macro( j_invariant=`algcurves/j_invariant` ): macro( g_normal=`algcurves/g_normal` ): # Compute the j invariant for curves with genus 1. j_invariant:=proc(F,x,y) local f,i,c; global `algcurves/x0`,`algcurves/y0`; # global, because then the options remember in the other # procedures work options `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`): # to read g_normal f:=collect(F,{x,y},distributed,g_normal); if degree(f,y)<2 then ERROR(`genus is not 1`) elif degree(f,y)=2 then # We could call `algcurves/Weierstrassform` but the following # is likely to be faster: f:=collect(f/lcoeff(f,y),y,g_normal); f:=collect(subs(y=y-coeff(f,y,1)/2,f),y,g_normal); f:=coeff(f,y,0); # Take the factors with odd multiplicity f:=numer(f)*denom(f); f:=collect(convert([seq(`if`(irem(i[2],2)=0,1,i[1]), i=evala(Sqrfree(f))[2])],`*`),x,g_normal); if not member(degree(f,x),{3,4}) then ERROR(`genus is not 1`) fi; for i from 0 to 4 do c[i]:=coeff(f,x,i) od; g_normal( 256*(c[2]^2+12*c[4]*c[0]-3*c[3]*c[1])^3/(-4*c[2]^3*c[0]*c[3]^2-4*c[2]^3*c[4]*c[ 1]^2+16*c[2]^4*c[4]*c[0]-80*c[3]*c[2]^2*c[1]*c[4]*c[0]+144*c[4]^2*c[1]^2*c[0]*c [2]-192*c[4]^2*c[0]^2*c[3]*c[1]+18*c[0]*c[3]^3*c[2]*c[1]-6*c[0]*c[3]^2*c[4]*c[1 ]^2+144*c[0]^2*c[3]^2*c[4]*c[2]+c[3]^2*c[2]^2*c[1]^2-128*c[4]^2*c[0]^2*c[2]^2- 27*c[0]^2*c[3]^4-27*c[4]^2*c[1]^4+256*c[4]^3*c[0]^3-4*c[3]^3*c[1]^3+18*c[3]*c[2 ]*c[1]^3*c[4]) ) elif degree(f,x)=2 then procname(f,y,x) else readlib(`algcurves/Weierstrassform`)(f,x,y, `algcurves/x0`,`algcurves/y0`,`j invariant`) fi end: #savelib ('j_invariant','`algcurves/j_invariant.m`'): #------------------------ filename: algcurves/genus_fx_plus_gy --------- # $Source: /u/maple/research/lib/algcurves/src/RCS/genus_fx_plus_gy,v $ # $Notify: hoeij@sci.kun.nl $ # Purpose of this file: compute the genus in an efficient way for curves # of the form (polynomial in x) + (polynomial in y). macro( genus_fx_plus_gy=`algcurves/genus_fx_plus_gy`, pq_formula=`algcurves/pq_formula` ): macro( g_factors=`algcurves/g_factors`, g_ext_r=`algcurves/g_ext_r` ): # Input: a polynomial f(x) and a polynomial g(y) # Output: consider the curve F(x,y)=f(x)+g(y), and the morphism # (x,y,z) -> (x,z) to P^1. Then compute S = Sum (e_P-1) where e_P are # the ramification indices of this morphism. genus_fx_plus_gy:=proc(f,g,x,y) local S,v,w,i,j,m; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; readlib(`algcurves/g_expand`): # to load g_factors # F:=f+g: # The sum of the e_P's at the line at infinity S:=pq_formula(degree(g,y),degree(g,y)-degree(f,x)); # Now for the finite points, determine the possible y coordinates # of ramification points: v:=g_factors(diff(g,y),y,g_ext_r([args])); for i in v do # # Determine the x coordinates: # w:=evala(Sqrfree(numer(subs(y=RootOf(i[1],y),F)),x))[2]; # S:=S+convert([seq(pq_formula(i[2]+1,j[2]) # *degree(j[1],x)*degree(i[1],y),j=w)],`+`) # the following is quicker: m:=sqrfree(numer(resultant(x-g,i[1],y)),x)[2][1]; w:=sqrfree(numer(resultant(f+y,subs(x=y,m[1]),y)),x)[2]; S:=S+m[2]*convert([seq(pq_formula(i[2]+1,j[2]) *degree(j[1],x),j=w)],`+`) od; # Now apply the Hurwitz formula and return the genus: S/2+1-degree(g,y) end: # If a curve looks locally at 0 like y^p - x^q, then the Puiseux expansions are # of the form y=x^{q/p} * (the p different p'th roots of unity) + higher terms. # Then the ramification indices are denom(q/p), and the number of places at # (x,y)=(0,0) is p/denom(q/p). So the sum of (e_P-1), where e_P is the # ramification index, this sum taken over the places P at (0,0) is the # following: pq_formula:=proc(p,q) option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; (denom(q/p)-1)*p/denom(q/p) end: #savelib('genus_fx_plus_gy','pq_formula','`algcurves/genus_fx_plus_gy.m`'): #------------------------ filename: algcurves/help_5.3_ervoor --------- lprint(`Help page is available for the following procedures:`): lprint(`genus,integral_basis,plot_knot,puiseux,singularities,parametrization`); lprint(`Weierstrassform,homogeneous and j_invariant.`): macro( genus=genus, integral_basis=integral_basis, plot_knot=plot_knot, puiseux=puiseux, singularities=singularities, parametrization=parametrization, Weierstrassform=Weierstrassform, j_invariant=j_invariant, homogeneous=homogeneous ): genus:=proc() `algcurves/genus`(args) end: integral_basis:=proc() `algcurves/integral_basis`(args) end: plot_knot:=proc() `algcurves/plot_knot`(args) end: puiseux:=proc() `algcurves/puiseux`(args) end: singularities:=proc() `algcurves/singularities`(args) end: parametrization:=proc() `algcurves/parametrization`(args) end: Weierstrassform:=proc() `algcurves/Weierstrassform`(args) end: j_invariant:=proc() `algcurves/j_invariant`(args) end: homogeneous:=proc() `algcurves/homogeneous`(args) end: #------------------------ filename: algcurves/help_5.3 --------- `help/text/genus` := TEXT(`Function: genus - Compute the genus of an algebraic curve`,` ` ,`Calling Sequence:`,` genus(f,x,y,opt)`,` `,`Parameters:`, ` f - a squarefree polynomial specifying an algebraic curve`, ` x,y - variables`,` opt - (optional) a sequence of options`,` `, `Description:`,` `, `- The genus of an irreducible algebraic curve is a nonnegative integer.`, ` However, if f is reducible this procedure could give a negative outcome. By` ,` default, this procedure does not check if f is irreducible. This default c\ an`, ` be altered by giving the optional argument ``irr check``. The polynomial f`, ` must be squarefree, otherwise an error message follows. `,` `, `Examples: `,`> `,`> f:=x^2+x*y+y^2;`, ` 2 2`, ` f := x + x y + y`,` `, `> factor(f); # It is irreducible in Q[x,y]`, ` 2 2`, ` x + x y + y`,` `, `> genus(f,x,y); # The result is impossible for an irreducible curve`, ` -1`,` `, `> evala(AFactor(f)); # So f must be reducible over the algebraic numbers`, ` 2 2`, `(y + (1 + 1/4 RootOf(16 + 4 _Z + _Z )) x) (y - 1/4 RootOf(16 + 4 _Z + _Z ) x)` ,` `,`> f:=761328152*x^6*z^4-5431439286*x^2*y^8+2494*x^2*z^8+`, `> 228715574724*x^6*y^4+9127158539954*x^10-15052058268*x^6*y^2*z^2+`, `> 3212722859346*x^8*y^2-134266087241*x^8*z^2-202172841*y^8*z^2`, `> -34263110700*x^4*y^6-6697080*y^6*z^4-2042158*x^4*z^6-201803238*y^10+`, `> 12024807786*x^4*y^4*z^2-128361096*x^4*y^2*z^4+506101284*x^2*z^2*y^6+`, `> 47970216*x^2*z^4*y^4+660492*x^2*z^6*y^2-z^10-474*z^8*y^2-84366*z^6*y^4:`, `# This f is a polynomial of degree 10 having a maximal number of cusps`, `# according to the Plucker formulas. It was found by Rob Koelman. It has 26`, `# cusps and no other singularities, hence the genus is 10. `, `> f:=subs(z=1,f);`, ` 6 2 8 2 6 4`, `f := 761328152 x - 5431439286 x y + 2494 x + 228715574724 x y`,` `, ` 10 6 2 8 2`, ` + 9127158539954 x - 15052058268 x y + 3212722859346 x y`,` `, ` 8 8 4 6 6`, ` - 134266087241 x - 202172841 y - 34263110700 x y - 6697080 y`,` `, ` 4 10 4 4 4 2`, ` - 2042158 x - 201803238 y + 12024807786 x y - 128361096 x y`,` `, ` 2 6 2 4 2 2 2 4`, ` + 506101284 x y + 47970216 x y + 660492 x y - 1 - 474 y - 84366 y`, ` `,`> genus(f,x,y);`,` 10`,` `, `See Also: singularities,parametrization `): `help/text/integral_basis` := TEXT( `Function: integral_basis - the integral basis of an algebraic function field`, ` `,`Calling Sequence:`,` integral_basis(a)`,` integral_basis(a,x`, ` integral_basis(f,x,y)`,` `,`Parameters:`, ` a - an algebraic function in RootOf form`,` x - a variable`, ` y - a variable`,` f - a squarefree polynomial in x and y`,` `, `Description:`,` `, `- This procedure computes an integral basis for an algebraic function field.`, ` The algebraic function field is specified either by a RootOf a giving an`, ` algebraic extension L(x) subset L(x)[a], or else by a minimal polynomial f`, ` of a in the variable y. This polynomial must be squarefree. `,` `, `- The output is a list of elements of the function field which generate (as a` ,` L[x] module) the integral closure of L[x] in the function field L(x)[a]. `, ` `,`- If a contains no other parameters than x (so then L is an algebraic`, ` extension of Q) then the second argument x is not necessary.`,` `, `- By setting infolevel[integral_basis] to a positive value some information`, ` will be printed during the computation. `,` `,`- The method used to compu\ te an integral basis is based on Puiseux expansions.`,` See the paper: "An al\ gorithm for computing an integral basis in an algebraic`, ` function field", J. of Symbolic Computation, 1994, 18, p. 353-363.`,` `, `Examples: `,`> `,`> alpha:=RootOf(x^3+7,x);`, ` 3`, ` alpha := RootOf(_Z + 7)`,` `, `> f:=y^8+3*x*y^5*alpha+5*x^4+x^6*alpha;`, ` 8 5 3 4 6 3`, ` f := y + 3 x y RootOf(_Z + 7) + 5 x + x RootOf(_Z + 7)`,` `, `> integral_basis(f,x,y);`, ` 3 4 2 3 4`, ` 2 y y y %1 y %1 y %1`, ` [1, y, y , ----, ----, -----, -----, -----]`, ` x x 2 2 3`, ` x x x`, ` `, ` 3 3`, ` %1 := 3 x RootOf(_Z + 7) + y`,` `, `# integral_basis(RootOf(f,y)) does the same, except that the output will be`, `# polynomials in RootOf(f,y) instead of in y.`,`See Also: puiseux,maxorder `): `help/text/plot_knot` := TEXT( `Function: plot_knot - make a tubeplot for a singularity knot`,` `, `Calling Sequence:`,` plot_knot(f,x,y,opt)`,` `,`Parameters:`, ` f - an algebraic curve with a singularity at the point 0`, ` x,y - variables`,` opt - (optional) a sequence of options`,` `, `Description:`,` `, `- Let f be a polynomial in x and y giving an algebraic curve in the plane C^2` ,` with a singularity at the point (x,y)=(0,0). The output of this procedure \ is`, ` called the singularity knot of this singularity. This knot is defined as`, ` follows: By identifying C^2 with R^4 the curve can be viewed as a 2`, ` dimensional surface over the real numbers. This procedure computes the`, ` intersection of this surface with a sphere in R^4 with radius epsilon and`, ` center 0. The intersection consists of a number of closed curves over the`, ` real numbers. After applying a projection from the sphere (is 3 dimensional` ,` over R) to R^3 these curves can be plotted by the tubeplot command in the`, ` plots package. Such a plot gives information about the singularity of f in`, ` the point 0. See also: E. Brieskorn, H. Kn"orrer: Ebene Algebraische Kurven\ ,`,` Birkh"auser 1981. `,` `, `- The curve given by f need not be irreducible, but f must be squarefree`, ` otherwise this procedure does not work. `,` `,`- If printlevel > 1 the nu\ mber of branches will be printed to the screen. Each`, ` branch (i.e. place above the point 0) corresponds to one component in the`, ` knot. `,` `,`- OPTIONS: `, ` epsilon: the radius of the sphere. The default is 1. If there are other`, ` singularities near 0 then a smaller number must be chosen, because the`, ` picture is not correct if there is more than 1 singularity inside the`, ` sphere. `, ` colours: Specifying a list of colours results in a plot where each branch`, ` gets its own colour.`, ` The options for tubeplot can be used as well. In plot_knot these options`, ` have the following default values: numpoints=150, radius=0.05, tubepoints=5\ ,`,` scaling=CONSTRAINED and style=PATCHNOGRID. `,` `,`Examples: `, `> `,`> printlevel:=2:`,`> plot_knot(y^2-x^3,x,y);`, `> f:=(y^3-x^7)*(y^2-2*x^5);`, `> plot_knot(f,x,y,epsilon=0.8,radius=0.03,colours=[blue,red]);`, `> plot_knot(f+y^3,x,y,epsilon=0.8);`, `> plot_knot(f+y^3,y,x,epsilon=0.8); # Same knot, but it looks`, `> # different because the projection point is different now`, `> # that x and y are switched.`,`> f:=(y^3-x^7)*(y^2-2*x^5)*(y^2+2*x^5);`, `> plot_knot(f,y,x,epsilon=0.8,radius=0.03,colours=[blue,red,pink]);`, `> g:=(y^3-x^7)*(y^3-x^7+100*x^13)*(y^3-x^7-100*x^13);`, `> plot_knot(g,x,y,epsilon=0.8,radius=0.03,numpoints=250,`, `> colours=[blue,red,green]);`,`See Also: tubeplot `): `help/text/puiseux` := TEXT( `Function: puiseux - determines the Puiseux expansions of an algebraic`, `function `,` `,`Calling Sequence:`,` puiseux(a,x=p,n)`, ` puiseux(f,x=p,y,n)`,` `,`Parameters:`, ` a - an algebraic function in RootOf form`, ` x=p - gives the point around which the expansions are computed`, ` n - a number specifying the desired accuracy of the expansions`, ` y - a variable`,` f - a squarefree polynomial in x and y`,` `, `Description:`,` `, `- A polynomial f of degree N in the variable y with coefficients in a field`, ` L(x) has N roots in the algebraic closure of the field L((x-p)). These root\ s`, ` are called the Puiseux expansions of f at x=p. Each Puiseux expansion is of` ,` the form Sum(a_i * (x-p)^(i/r),i=-N..infinity) for some integer r (r is`, ` called the ramification index of the Puiseux expansion), and some integer N` ,` and elements a_i in the algebraic closure of L. `,` `, `- The polynomial f must be squarefree, otherwise the puiseux procedure does`, ` not work. `,` `,`- If f is irreducible then f gives an algebraic extensio\ n L(x)[y]/(f) of L(x).`, ` Instead of giving f, this algebraic extension can also be specified with a`, ` RootOf a of f. This a can be viewed as a multivalued function of x. The`, ` Puiseux expansions give the local expansions of this multivalued function. ` ,` `, `- The procedure puiseux determines the field L from a or f. The groundfield L` ,` of the computation is the smallest field such that f is in L(x)[y]. `,` ` ,`- The Puiseux expansions are only computed up to conjugation over L((x-p)). \ So`, ` if a number of expansions are algebraically conjugated over L((x-p)) then`, ` only one of these expansions is given. `,` `, `- The Puiseux expansions are computed modulo x^n. So if for instance n=10,`, ` then the term x^(49/5) would be computed, but not the term x^10.`, ` If n=0 then the expansions are not computed modulo x^0, but in this case`, ` the number of terms that is computed is precisely the number that is needed` ,` to distinguish the expansions from the other expansions. `,` `, `- NOTE: The Maple "alias" function does not recognize an alias in terms of`, ` another alias. Therefore, you must not use nested aliasses for algebraic`, ` numbers because then the puiseux algorithm is not able to construct the`, ` field L over which the algebraic function is defined. `,` `, `Examples: `,`> `,`> alpha:=RootOf(x^3+7);`, ` 3`, ` alpha := RootOf(_Z + 7)`,` `, `> f:=y^8+x*y^5+x^4-x^6*alpha;`, ` 8 5 4 6 3`, ` f := y + x y + x - x RootOf(_Z + 7)`,` `, `> puiseux(f,x=0,y,0); # is the same as puiseux(RootOf(f,y),x=0,0);`, ` 1/3 3/5`, ` {(-x) , -x }`,` `,`> puiseux(f,x=0,y,5);` ,` 170 13/3 11/3 3 5/3 1/3`, `{- --- %2 + 1/3 %1 %2 + 2/3 x - 1/3 %2 + %2 ,`,` 81 `,` `, `/ 25389 2\\ 23/5 42 21/5 483 19/5 17/5`, `|- ----- + 2/25 %1 | x + --- %1 x - --- x + 4/25 %1 x`, `\\ 15625 / 125 625`,` `, ` 3 13/5 11/5 7/5 3/5`, ` - 2/5 x + 1/5 %1 x - 6/25 x - 1/5 x - x }`,` `, ` 3`,`%1 := RootOf(_Z + 7)`,` `,`%2 := -x `,` `, `# To get all 8 Puiseux expansions substitute all conjugates over L((x)) of`, `# x^(1/5) and (-x)^(1/3) in this set.`,`See Also: series `): `help/text/singularities` := TEXT( `Function: singularities - Compute the singularities of an algebraic curve`, ` `,`Calling Sequence:`,` singularities(f,x,y)`,` `,`Parameters:`, ` f - a polynomial specifying an algebraic curve`,` x,y - variables`, ` `,`Description:`,` `,`- Let f be squarefree polynomial in x and y. Then \ f defines an algebraic curve`, ` in the plane C^2, and also in the projective plane P^2 by making f`,` homo\ geneous. This procedure computes the singular points of the curve in the`, ` projective plane. The points are given by homogeneous co-ordinates [X,Y,Z].` ,` For each singularity this procedure also computes the multiplicity and the` ,` delta invariant. `,` `, `- The output of this procedure is a set consisting of lists of the following`, ` form [ point , multiplicity , delta invariant]. `,` `, `- This procedure computes all singularities up to conjugation. So if a`, ` singularity [RootOf(_Z^2-2),1,1] is given in the output, and if`,` RootOf(\ _Z^2-2) does not appear in the input, then [-RootOf(_Z^2-2),1,1] is a`, ` singular point as well but will not be given in the output. `,` `, `- The genus of a curve is the number (d-1)*(d-2)/2 - Sum(delta invariants)`, ` where d is the degree of the curve. Note that if we apply this formula to`, ` compute the genus, then for each singularity we must multiply the delta`, ` invariant by the degree of the algebraic extension over which the`, ` singularity is defined, because only 1 singularity of each conjugacy class`, ` is given in the output. `,` `,`Examples: `,`> `, `> f:=-8*y^5-207*y*x^4+180*x^5-35*y^4-128*y^3*x+621*y*x^3-450*x^4+`, `> 82*y^3-521*y*x^2+369*x^3-19*y^2+135*y*x-100*x^2-28*y-7*x+8;`, ` 5 4 5 4 3 3 4 3` ,`f := -8 y - 207 y x + 180 x - 35 y - 128 y x + 621 y x - 450 x + 82 y` ,` `,` 2 3 2 2`, ` - 521 y x + 369 x - 19 y + 135 y x - 100 x - 28 y - 7 x + 8`,` `, `> degree(f,{x,y});`,` 5`,` `, `> singularities(f,x,y);`, `{[[0, 1, 1], 2, 1], [[1, 0, 1], 2, 1], [[1, -1, 1], 2, 1],`,` `, ` 2 2 `, `[[RootOf(_Z + 1), RootOf(_Z + 1), 1], 2, 1]}`,` `, `# Note that the conjugate (replace RootOf(_Z^2+1) by -RootOf(_Z^2+1)) is also` ,`# a singularity. So the genus is (5-1)*(5-2)/2-1-1-1-2*1=1`, `See Also: puiseux,genus,singular `): `help/text/parametrization` := TEXT( `Function: parametrization - Compute a parametrization for a curve with genus`, `0 `,` `,`Calling Sequence:`,` parametrization(f,x,y,t)`,` `, `Parameters:`,` f - a squarefree polynomial in x and y, with genus 0`, ` x,y,t - variables`,` `,`Description:`,` `, `- This procedure computes, if it exists, a parametrization of an algebraic`, ` curve f. A parametrization is a birational equivalence from a projective`, ` line to the given curve f. Such a parametrization exists if and only if the` ,` genus is 0 and the curve is irreducible (which can be checked by AFactor).` ,` `, `- The output of the procedure is a list [X(t),Y(t)] of rational functions in`, ` t, such that [X(t),Y(t)] is a point on the curve f for every value of t`, ` (although one can not use substitution for those values of t for which the`, ` denominator of X(t) or Y(t) vanishes). `,` `,`Examples: `, `> `,`> f:=y^5+2*x*y^2+2*x*y^3+x^2*y-4*x^3*y+2*x^5:`, `> v:=parametrization(f,x,y,t);`, ` 13500 + 8775 t + 2070 %1 + 212 %2 + 8 %3`, `v := [-----------------------------------------------------,`,` 5 `, ` 4 t + 160 %3 + 2480 %2 + 18920 %1 + 71590 t + 107897`,` `, ` 960 + 848 t + 276 %1 + 39 %2 + 2 %3`, `-2 -----------------------------------------------------]`,` 5 `, ` 4 t + 160 %3 + 2480 %2 + 18920 %1 + 71590 t + 107897`,` `,` 2 `, `%1 := t `,` `,` 3 `,`%2 := t `,` `,` 4 `,`%3 := t `, ` `, `# Now subs(t=any number,v) should be a point on the curve. Test the result`, `# (this should be 0): `,`> normal(subs(x=v[1],y=v[2],f));`, ` 0`, ` `, `> parametrization(x^4+y^4+a*x^2*y^2+b*y^3,x,y,t);`, ` 3 3 4`, ` t b t b`, ` [- -----------------, - -----------------]`, ` 8 2 4 4 8 2 4 4`, ` b + t b a + t b + t b a + t`,` `, `See Also: genus,AFactor `): `help/text/Weierstrassform` := TEXT( `FUNCTION: Weierstrassform - compute an isomorphism for function fields with genus 1`, ` `,`CALLING SEQUENCE: Weierstrassform(f,x,y,x0,y0,opt)`,` `,`PARAMETERS:`, ` f - a polynomial in x and y describing the curve of genus 1`, ` x,y,x0,y0 - variables`,` opt - (optional) a sequence of options`, ` `,`SYNOPSIS: `, ` - An algebraic function field Qbar(x)[y]/(f) where Qbar is the algebraic`, ` closure of Q is isomorphic to the field Qbar(x0)[y0]/(f0) where f0 is of`, ` the form y0^2 + square free polynomial in x0 of degree 3`, ` if and only if the genus is equal to 1.`, ` This procedure computes such an isomorphism.`,` `, ` - The output of genus is a list of 5 items containing:`,` The curve f0`, ` The image of x0 under this isomorphism`, ` The image of y0 under this isomorphism`, ` The image of x under the inverse isomorphism`, ` The image of y under the inverse isomorphism`,` `, ` - A regular point [x,y,z] on the curve can be specified as a 6'th argument.`, ` Then x0 and y0 will have their poles in this point on the curve. If the`, ` option ``no inverse`` is used then the inverse isomorphism is not`, ` computed. The option ``Weierstrass`` results in a Weierstrass normal form`, ` i.e. y0^2 - 4*x0^3 - a*x0 - b. If the option ``j invariant`` is specified`, ` then only the j-invariant is given.`,` `,`EXAMPLES: `, `f:=x^4+y^4-2*x^3+x^2*y-2*y^3+x^2-x*y+y^2;`,`v:=Weierstrassform(f,x,y,x0,y0);`, ` 3 2 - 3 y + 2 x`, `v := [x0 + 2/3 x0 + 1/108 + y0 , 1/3 -----------,`, ` x`,` `, ` 2 2 3 3 2`, ` 2 x y - 2 y - x + x + 2 y - 2 x y`, `- 1/2 --------------------------------------,`,` 2`, ` x (x - 1)`,` `, ` 2 3`, ` - 36 y0 + 54 y0 x0 + 45 x0 - 108 x0 - 52 + 54 x0`, `- 3/2 --------------------------------------------------,`, ` 2 3 4`, ` 97 - 96 x0 + 216 x0 - 216 x0 + 81 x0`,` `,` \ 2 2 4 3`,` 72 y0 - 21\ 6 y0 x0 + 162 y0 x0 - 246 x0 + 351 x0 + 162 x0 + 104 - 432 x0`,`1/2 -------\ --------------------------------------------------------------------`, ` 2 3 4`, ` 97 - 96 x0 + 216 x0 - 216 x0 + 81 x0`,` `,`] `, ` `,`# check if the image of x and y still satisfy the relation f in`, `# the field Qbar(x0)[y0]/(f0)`,`subs({x=v[4],y=v[5]},f):`, `evala(subs(y0=RootOf(v[1],y0),"));`,` 0` ,` `,`# check if the image of x0 and y0 still satisfy the relation f0 in`, `# the field Qbar(x)[y]/(f)`,`subs({x0=v[2],y0=v[3]},v[1]):`, `evala(subs(y=RootOf(f,y),"));`,` 0`, `SEE ALSO: genus, j_invariant`): `help/text/j_invariant` := TEXT( `FUNCTION: j_invariant - compute the j invariant of an elliptic curve.`,` `, `CALLING SEQUENCE: j_invariant(f,x,y)`,` `,`PARAMETERS:`, ` f - a polynomial in x and y representing a curve of genus 1`, ` x,y - variables`,` `,`SYNOPSIS: `, ` - For algebraic curves with genus 1 one can compute a number called the`, ` j invariant. The main property of this j invariant is the following: two`, ` elliptic (this means: genus = 1) curves are birationally`, ` equivalent (i.e. can be transformed to each other with rational`, ` transformations over an algebraically closed field) if and only if`, ` their j invariants are the same.`,` `, ` - The curve must be irreducible and have genus 1, otherwise the j invariant`, ` is not defined and this procedure will fail.`,` `,`EXAMPLE: `, `f:=y^5+4/3-23/3*y^2+11*y^3-17/3*y^4-16/3*x^2+16/3*x^3-4/3*x^4:`, `# Check that the genus is 1, because only then is the j invariant defined.`, `genus(f,x,y);`,` `,` 1`,` `, `j_invariant(f,x,y);`,` `, ` 1404928`, ` - -------`, ` 171`,`SEE ALSO: genus, Weierstrassform` ): `help/text/homogeneous`:=TEXT( `Function: homogeneous - Make a polynomial in two variables homogeneous in`, `three variables.`,` `,`Calling Sequence:`,` homogeneous(f,x,y,z)`,` `, `Parameters:`,` f - a polynomial in x and y`,` x,y,z - variables`,` ` ,`Description:`,` `, `- Algebraic curves are often represented using a polynomial in two variables,` ,` or as a homogeneous polynomial in three variables. This procedure converts` ,` to the homogeneous representation.`,` `, `- If f is a polynomial in x and y of degree d, and f can be written as Sum`, ` c[i,j] * x^i * y^j for some coefficients c[i,j], then the output of this`, ` procedure is Sum c[i,j] * x^i * y^j * z^(d-i-j).`,` `,`Examples: `, `> `,`> f:=y^2-x^3;`, ` 2 3`, ` f := y - x`,` `,`> homogeneous(f,x,y,z);`, ` 3 2`, ` -x + y z`,` `, `> # Now one can convert back as follows:`,`> subs(z=1,");`, ` 2 3`, ` y - x`,` `,`See Also: algcurves `): `algcurves/help_5.3`:=proc() end: # To avoid: Error, ineffective readlib of help_5.3 macro( readlib=readlib ); macro( solve=solve );