# $Source: /u/maple/research/lib/algcurves/src/RCS/integral_basis,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ macro( integral_basis1=`algcurves/integral_basis`, integral_basis2=`algcurves/integral_basis2`, 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( RAD={radical,nonreal}, solve=`solve/linear`, monic=`algcurves/monic`, `puiseux/technical_answer`=`algcurves/puiseux_te`, g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, g_normal=`algcurves/g_normal`, g_expand=`algcurves/g_expand`, g_zero_of=`algcurves/g_zero_of`, g_factors=`algcurves/g_factors`, g_ext=`algcurves/g_ext`, truncate=`algcurves/truncate` ): ########################################################################### # # Checks types and calls the appropriate routine. MJR 03/97 # macro( RADFIELD = radfield, RADEXT = `algcurves/integral_basis/radext` ); integral_basis1 := proc(alg) local i,j,f,r,x,y,z,rofs,rads,rel,back,back1,back2,forw,tm, ground,B,rof,d,minpol,F,s; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; userinfo(1,{'algcurves','integral_basis'}, `computing an integral basis: start time`,time()); if nargs=1 then x := select(type,indets(alg),name); if nops(x)>1 then error "choose a variable: %1", x else x := x[1] fi elif nargs>=2 then if args[2]::name then x := args[2] else error "second argument should be a name" fi; fi; if nargs>=4 then return integral_basis2(args) fi; if alg::{radext,algext,set({radext,algext})} then if alg::{algfun(rational),set(algfun(rational))} then rofs := indets(alg,RootOf); elif alg::{'radalgfun',set('radalgfun')} then rofs := indets(alg,{radical,nonreal,function}); rads := select(type,rofs,radnum); if rads <> {} then s := RADFIELD(rads); forw := s[1]; back := s[2]; rofs := eval(subs(forw,rofs)); fi; rads := indets(rofs,{radical,nonreal}); if rads <> {} then forw := [seq(tm=convert(tm,RootOf),tm=rads)]; back1 := [seq(rhs(tm)=lhs(tm),tm=forw)]; F := proc(t) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if op(1,t)::(RootOf^integer) then op([1,1],t)=op([2,1],t)^(1/denom(op([2,2],t))) elif type(op(1,t),`*`) and nops(op(1,t))=2 then if op([1,1],t) = x or op([1,2],t) = x then procname(t/x) elif type(op([1,1],t),`^`) and op([1,1,1],t)=x and type(op([1,1,2],t),integer) then procname(t/op([1,1],t)); elif type(op([1,2],t),`^`) and op([1,2,1],t)=x and type(op([1,2,2],t),integer) then procname(t/op([1,2],t)); fi; else t fi; end: back1 := map(F,{op(back1)}); rofs := indets(eval(subs(forw,rofs)),function) fi; else error "invalid first argument" fi; ground := remove(has,rofs,x); rofs := rofs minus ground; if not evala(Indep(ground,'rel')) then error "the RootOfs should be independent. Relations are %1",rel fi; if rofs={} then error "first argument is not a function of %1" ,x elif nops(rofs)=1 then if frontend(numboccur,[op(1,rofs[1]),_Z])=1 then r := RADEXT(rofs[1],x) # Use Van Hoeij's code. else r := integral_basis2(rofs[1],x); fi; if assigned(back1) then r := subs(back1,r) fi; if assigned(back) then r := subs(back,r) fi; userinfo(1,{'algcurves','integral_basis'}, `computing an integral basis: end time`,time()); return r fi; if not evala(Indep(rofs,'rel')) then error "the RootOfs should be independent. Relations are %1",rel fi; else if nargs>=3 then if args[3]::name then y := args[3]; if x=y then error "the second argument should be " "different than the third one" fi; else error "third argument should be a name" fi else error "missing third argument or invalid first argument" fi; if not alg::polynom('radalgfun',{x,y}) then error "invalid arguments" fi; ground := indets(alg,{radical,nonreal,function}); f := alg; if ground::{radnum,set(radnum)} then s := RADFIELD(ground); forw := s[1]; back := s[2]; f := eval(subs(forw,f)); elif ground::{'radalgfun',set('radalgfun')} then rads := select(type,ground,radnum); if rads <> {} then s := RADFIELD(rads); forw := s[1]; back := s[2]; f := eval(subs(forw,f)); fi; rads := indets(f,{radical,nonreal}); if rads <> {} then # There's not much we can do for now. forw := [seq(tm=convert(tm,RootOf),tm=rads)]; back1 := [seq(rhs(tm)=lhs(tm),tm=forw)]; f := eval(subs(forw,f)) fi; if not evala(Indep(rads,'rel')) then error "the RootOfs should be independent. Relations are %1",rel fi; else error "invalid first argument" fi; if not has(f,y) then error "first argument should depend on %1",y fi; if numboccur(f,y)=1 then r := RADEXT(f,x,y) else r := integral_basis2(f,x,y,args[4..nargs]); fi; if assigned(back1) then r := subs(back1,r) fi; if assigned(back) then r := subs(back,r) fi; return r fi; forw := `evala/mkmonic`(rofs,'back2'); rofs := indets(subs(forw,rofs),'RootOf'); ground := remove(has,rofs,x); rofs := rofs minus ground; rofs := sort(rofs,(a,b) -> evalb(`evala/algdeg`(a) < `evala/algdeg`(b))); # build a basis to start with # B := [1]; for rof in rofs do minpol := frontend(subs,['_Z'=z,op(1,rof)],[{`+`,`*`},{'_Z'=z}]); minpol := collect(minpol,z); d := degree(minpol,z); B := [seq(seq(B[i]*rof^j,i=1..nops(B)),j=0..d-1)]; od; userinfo(2,{'algcurves','integral_basis'}, `degree of the field is`, nops(B)); r := `algcurves/integral_basis/intff`(B,x,rofs,ground); r := convert(evalm(transpose(r) &* B), 'list'); r := map(normal,subs(back2,r)); if assigned(back1) then r := eval(subs(back1,r)) fi; if assigned(back) then r := eval(subs(back,r)) fi; userinfo(1,{'algcurves','integral_basis'}, `computing an integral basis: end time`,time()); r end: # # End of add-on. # ############################################################################ # 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. Author: M. van Hoeij`; n:=degree(f,y); f_local:=subs(g_conversion2,numer(subs(x=p,g_expand(rem(f,P,x),ext)))); userinfo(5,{'algcurves','integral_basis'}, `Doing a squarefree factorization of`,f_local); v:=evala(Sqrfree(f_local,y),{'expanded','independent'})[2]; userinfo(5,{'algcurves','integral_basis'}, `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','integral_basis'}, `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] # singularity with a ramification elif degree(m[2][1],y)=1 then N:=m[2][1] # a cusp 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:=evala(Gcd(f_local,m[2][1])) fi; if degree(N,y)<>1 then error "wrong degree" fi; userinfo(4,{'algcurves','integral_basis'},`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 elif D<>2 and m[3]=[1] and degree(m[2][1],y)<>1 then N:=numer(g_normal(m[2][1]/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,y. # Both RootOf and rootof syntax. # ext: gives the ground field # Output: those factors (and their multiplicities) of discrim(f,y) # of which the roots could be x-coordinates of singularities. # Note that the root should have at least multiplicity 2. double_factors:=proc(f,x,y,ext) global algcurves,g_conversion1,g_conversion2; local df,v,R,i,j; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved. Author: M. van Hoeij`; if ext<>[] and not has(ext,RootOf) then subs(g_conversion1,procname(op(subs(g_conversion2,[args])))) else df:=discrim(f,y); if testeq(df) and evala(df)=0 then error "Curve should be squarefree" fi; userinfo(5,algcurves, `Doing a squarefree factorization of`,df); v:=sqrfree(df,x); userinfo(5,algcurves, `Done squarefree factorization`); df:=NULL; # Compute gcd's with the following polynomial, because if # alpha is not a root of R, then there's no singularity at # x=alpha. R:=evala(Primpart(resultant(diff(f,x),diff(f,y),y),x)); for i in v[2] do if i[2]>1 then df:=df,seq([collect(j[1]/lcoeff(j[1],x),x,evala),i[2]], j=g_factors(evala(Gcd(R,i[1])),x,ext)) fi od; [df] fi 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_basis2 :=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. Author: M. van Hoeij`; if indets([args],RAD)<>{} then f:=[args]; v:=radfield(indets(f,RAD)); return subs(v[2],procname(op(subs(v[1],f)))) fi; 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: extl:=g_ext([args]); df:=subs(g_conversion1,{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:=double_factors(f,x,y,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)]; Normalizer:=`if`(extl=[],normal,`evala/Normal`); 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','integral_basis'}, `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','integral_basis'}, `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(subs(g_conversion2,v[2]*B[i]+v[3]*res[i]) ,y,Normalizer) fi od; res:=[R] fi od; res:=subs(g_conversion2,subs(y=alfa*q,res)); if has(alfa,x) then # entries need to be normalized: map(Normalizer,res) else res fi 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. Author: M. van Hoeij`; 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','integral_basis'}, `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(add(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(add(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( add(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: #savelib ('local_intbasis23','double_factors','integral_basis1',\ 'local_intbasis','integral_basis2'):