########################## # nfintbas ### # Convert from 5.2 to Release 3 # Authors: Michael Beckker (beckker@sci.kun.nl) # Mark van Hoeij (hoeij@sci.kun.nl) # Description: Computation of an integral basis in an algebraic number field # Methods: * Round two: implemented by Beckker, cf. "Het berekenen van een # basis van gehelen in een getallenlichaam" (written in Dutch), # graduation work, University of Nijmegen. # * A method described in "An algorithm for computing an integral # basis in an algebraic function field", to be published in # JSC. implemented by Mark van Hoeij # # For the case of a wild ramification the method used by nfintbas # is roundtwo, for the other cases the default is the second # method. This default can be overruled by an argument roundtwo # # Bug reports: Please send them to hoeij@sci.kun.nl # Note: The code is poorly documented at the moment, but if someone # is interested in reading the code this could be improved. # date: Aug 22 1994 `help/text/nfintbas`:=TEXT( ``, `FUNCTION: nfintbas - number field integral basis`, ``, `CALLING SEQUENCE: nfintbas(f,opt)`, ``, `PARAMETERS:`, ` f - a polynomial in one variable`, ` opt - (optional) a sequence of options`, ``, `SYNOPSIS:`, `- This procedure computes an integral basis for the algebraic number field`, ` Q[x]/(f), if x is the variable used in f.`, ``, `- A local integral basis in a set S of prime numbers can be obtained by`, ` including the arguments ``local``,S. The default computation method can be`, ` changed to the round two method by including an argument ``roundtwo``.`, ` If infolevel[nfintbas] > 0 some information will be printed during the`, ` computation.`, ``, `EXAMPLES:`, ` f:=x^4+20*x^2+20;`, ` nfintbas(f);`, ` 2 3`, ` [1, x, 3/4 + 1/8 x , 3/4 x + 1/8 x ]`, ``, ` g:=x^26-x^3+125;`, ` nfintbas(g,``local``,{5});`, ` 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17`, `[1, x, x , x , x , x , x , x , x , x , x , x , x , x , x , x , x , x ,`, ``, ` 18 19 20 21 22 23 24 24 2 25`, ` x , x , x , x , x , x , 1/5 x + 4/5 x, ---- x + 1/25 x ]`, ` 25`, ``, ` # The following 2 may take some time:`, ` f:=y^20+5*y^18+6*y^16+6*y^15+3*y^14+3*y^13+6*y^11+6*y^10+3*y^9+3*y^6+3*y^5`, ` +6*y^4+5*y^2+6*y+1;`, ` nfintbas(f);`, ` nfintbas(y^20+8*y^19+9*y^16+6*y^12+6*y^8+9*y^4+8*y^3+1);`, ``, `SEE ALSO: maxorder` ): # (: # Dutch to English translation of Beckker's procedure names: macro( basis_wortelideaal=`nfintbas/basis_radical`, controleer=`nfintbas/check`, corrigeer_basis=`nfintbas/correct_basis`, druk_uit_in=`nfintbas/write_as`, goede_rooster=`nfintbas/lattice`, idealizer=`nfintbas/idealizer`, kwadraat=`nfintbas/square`, locale_ring_van_gehelen=`nfintbas/local_intbasis_r2`, maal=`nfintbas/multiply`, macht_modulof=`nfintbas/exponent_mod_f`, meng_bases=`nfintbas/mix_bases`, monisch=`nfintbas/monic`, pvaluatie=`nfintbas/p_valuation`, reduceer=`nfintbas/reduce`, ring_van_gehelen=`nfintbas/roundtwo`, verz_priemfac_waarvan_kwadr_disc_deelt=`nfintbas/square_factors`, zoek_minimale_ai=`nfintbas/minimal_a_i` ): ################################################################# # Round two implementation macro( matrix=linalg[matrix], transpose=linalg[transpose], inverse=linalg[inverse], alfa=`gehelen/alfa` ): ring_van_gehelen:=proc(f,x) local disc,n,P,B,R,i,j,g; if not (((monisch(f,x) and irreduc(f)) and type(f,polynom(integer,x))) and degree(f,x)>1) then ERROR(`f geen monisch irreducibel polynoom van graad groter dan 1`, `f moet gehele coefficienten bevatten`) fi; disc:=discrim(f,x); n:=degree(f,x); readlib(ifactors): #R:=[seq(alfa^i,i=0..n-1)]; R:=[seq(x^i,i=0..n-1)]; if nargs=2 then P:=verz_priemfac_waarvan_kwadr_disc_deelt(ifactors(disc)[2]) elif nargs=4 and args[3]=locaal then P:={[args[4],pvaluatie(disc,args[4])]}; if P[1][2]<=1 then P:={} fi else ERROR(`wrong arguments`) fi; if P={} then RETURN(R) #(subs(alfa=x,R)) fi; B:={}; # g:=subs(x=alfa,f); g:=f; for i to nops(P) do B:=B union {locale_ring_van_gehelen(P[i][1],trunc(P[i][2]/2),R,g,n,x)} od; for j to nops(B) do R:=meng_bases(R,B[j],x) od; # subs(alfa=x,reduceer(R)) reduceer(R,x) end: monisch:=proc(f,x) lcoeff(f,x)=1 end: verz_priemfac_waarvan_kwadr_disc_deelt:=proc(V) local i,verz; verz:={}; for i to nops(V) do if V[i][2] >= 2 then verz:=verz union {V[i]} fi od; verz end: pvaluatie:=proc(a,p) if a=0 then infinity elif irem(a,p)=0 then 1+pvaluatie(a/p,p) else 0 fi end: locale_ring_van_gehelen:=proc(p,t,R,f,n,x) local V,W,IdW,i,k,d,oude_d; oude_d:=1; V:=R; k:=1; while p^k < n do k:=k+1 od; do W:=basis_wortelideaal(p,t,V,f,k,n,x); IdW:=idealizer(W,f,p,t,n,x); d:=1; for i to nops(IdW) do d:=d*denom(IdW[i]) od; if goede_rooster(d,p,oude_d,t) then RETURN(IdW) else V:=eval(IdW); oude_d:=d fi od end: meng_bases:=proc(basis1,basis2,x) local i,result,s,t; result:=NULL; for i to nops(basis1) do # if igcdex(denom(lcoeff(basis1[i],alfa)),denom(lcoeff(basis2[i],alfa)),'s','t')<>1 if igcdex(denom(lcoeff(basis1[i],x)),denom(lcoeff(basis2[i],x)),'s','t')<>1 then ERROR(`noemers hebben niet ggd 1`) fi; # Deze igcdex heeft nu s en t bepaald. result:=result,s*basis2[i]+t*basis1[i] od; [result] end: reduceer:=proc(V,x) local W,w,d,i,j,k,l,s; W:=V; for i to nops(W) do if denom(W[i]) <> 1 then w:=(denom(W[i])*W[i] mod denom(W[i]))/denom(W[i]); W:=[W[1..i-1],w,W[i+1..nops(W)]] fi od; for j from nops(V) by -1 to 2 do l:=W[j]; d:=denom(W[j]); for k from j-1 by -1 to 1 do # s:=iquo(d*coeff(l,alfa,k-1),d*lcoeff(W[k])); s:=iquo(d*coeff(l,x,k-1),d*lcoeff(W[k])); l:=expand(l-s*W[k]); if denom(l) <> 1 then l:=(denom(l)*l mod denom(l))/denom(l) fi od; W:=[W[1..j-1],l,W[j+1..nops(W)]] od; W end: basis_wortelideaal:=proc(p,t,R,f,k,n,x) local a,b,i,j,l,m,M,B,D,J,pt; a:=0; pt:=p^t; for i to nops(R) do #a:=collect(plus(a,z.i*macht_modulof(pt*R[i],p^k,f,p,t)),alfa) a:=collect(a+z.i*macht_modulof(pt*R[i],p^k,f,p,t,x),x) od; a:=a/pt; b:=druk_uit_in(a,R,x); M:=transpose(matrix(n,n,[seq([seq(coeff(b[j],z.i),j=1..n)],i=1..n)])); B:=Nullspace(M) mod p; D:=[seq(0,i=1..nops(B))]; J:=[seq(p*R[i],i=1..nops(R))]; for l to nops(B) do for m to n do D:=[D[1..l-1],D[l]+B[l][m]*R[m],D[l+1..nops(B)]] od; od; D:=controleer(D,p,x); corrigeer_basis(J,D,x) end: idealizer:=proc(V,f,p,t,n,x) local W,Z,M,d,g,h,i,j,k,l,q,pivot,Y,u,z,pt; pt:=p^t; #g:=sum('a.z*alfa^z','z'=0..n-1); g:=sum('a.z*x^z','z'=0..n-1); W:={}; for i to n do #h:=collect(maal(g,pt*V[i],f,p,t),alfa); h:=collect(maal(g,pt*V[i],f,p,t,x),x); for j to n do #c.(n+1-j):=coeff(h,alfa,n-j)/lcoeff(V[n+1-j],alfa); c.(n+1-j):=coeff(h,x,n-j)/lcoeff(V[n+1-j],x); W:=W union {c.(n+1-j) mod pt}; #h:=collect(h-c.(n+1-j)*V[n+1-j],alfa) h:=collect(h-c.(n+1-j)*V[n+1-j],x) od od; Z:=[]; for k to n do q:=zoek_minimale_ai(W,k-1,p,t); if q[2] <> 0 then W:=W minus {q[1]};Y:={}; pivot:=((1/(q[2]/q[3])mod pt)*q[1]) mod pt; for l in W do u:=(l-(coeff(l,a.(k-1))/q[3])*pivot) mod pt; Y:=Y union {u} od; W:=eval(Y); Z:=[op(Z),pivot] else Z:=[op(Z),q[1]] fi; od; M:=matrix(n,n,[seq([seq(coeff(Z[i],a.j),j=0..n-1)],i=1..n)]); M:=transpose(inverse(M)); #d:=convert([seq(alfa^i,i=0..n-1)],list); d:=convert([seq(x^i,i=0..n-1)],list); RETURN([seq((i*pt mod pt*p)/pt,i=convert(evalm(M&*d),list))]) end: goede_rooster:=proc(d,p,oude_d,t) d=p^t or d=oude_d end: macht_modulof:=proc(getal,exp,pol,p,t,x) local pt; pt:=p^(t+1); if exp = 2 then kwadraat(getal,pol,p,t,x) mod pt elif exp mod 2 = 0 then kwadraat(macht_modulof(getal,exp/2,pol,p,t,x),pol,p,t,x) mod pt else maal(macht_modulof(getal,exp-1,pol,p,t,x),getal,pol,p,t,x) mod pt fi end: druk_uit_in:=proc(f,V,x) local d,g,i,l,m; g:=f; l:=[seq(0,i=1..nops(V))]; while g <> 0 do #d:=degree(g,alfa)+1; d:=degree(g,x)+1; #m:=lcoeff(g,alfa)/lcoeff(V[d],alfa); m:=lcoeff(g,x)/lcoeff(V[d],x); l:=[l[1..d-1],m,l[d+1..nops(V)]]; #g:=collect(g-m*V[d],alfa) g:=collect(g-m*V[d],x) od; l end: controleer:=proc(V,p,x) local W,d1,d2,i,j,k; W:=V; for i to nops(V)-1 do #d1:=degree(V[i],alfa); d1:=degree(V[i],x); for j from i+1 to nops(V) do #d2:=degree(V[j],alfa); d2:=degree(V[j],x); if d1 = d2 #then k:=lcoeff(V[j],alfa)/lcoeff(V[i],alfa) mod p; then k:=lcoeff(V[j],x)/lcoeff(V[i],x) mod p; W:=[W[1..j-1],V[j]-k*V[i] mod p,W[j+1..nops(V)]] fi od od; RETURN(convert(W,list)) end: corrigeer_basis:=proc(V,W,x) local d,h,i,R; R:=V; for i to nops(W) do h:=W[i]; #d:=degree(h,alfa)+1; d:=degree(h,x)+1; R:=[R[1..d-1],h,R[d+1..nops(V)]] od; R end: zoek_minimale_ai:=proc(W,k,p,t) local i,j,c; for i to t do for j to nops(W) do c:=coeff(W[j],a.k); if c mod p^i <> 0 then RETURN(W[j],c,p^(i-1)) fi od od; RETURN(p^t*a.k,0,0) end: maal:=proc(a,b,f,p,t,x) #rem(expand(a*b/p^t),f,alfa) rem(expand(a*b/p^t),f,x) end: kwadraat:=proc(getal,f,p,t,x) maal(getal,getal,f,p,t,x) end: ########################################################## # Factorization in finite fields macro( factorise_modp=`nfintbas/factorise_modp`, degree_ext=`nfintbas/degree_ext`, gcd_modp=`nfintbas/gcd_modp`, rem_modp=`nfintbas/rem_modp` ): # Input: f is a polynomial in x # ext=[RootOf(..),RootOf(..),...] gives an extension on Fp # ext[nops(ext)] is the last extension # p the characteristic # Output: the factors + multiplicities of f # Comments: rather inefficient code factorise_modp:=proc(f,x,ext,p) local ff,i,ii,result,ra,g,dummy; if degree(f,x)=1 then RETURN([[evala(f/lcoeff(f,x)) mod p,1]]) fi; if nops(ext)<2 then ff:=Factors(f,op(ext)) mod p; RETURN(ff[2]) fi; ff:=factorise_modp(evala(Norm(f,ext,[ext[1..nops(ext)-1]])) mod p, x,[ext[1..nops(ext)-1]],p); result:=NULL; for i in ff do ii:=gcd_modp(f,i[1],x,p); if degree(i[1],x)=degree(ii,x)*degree_ext(ext,nops(ext)) then result:=result,[ii,i[2]]; else ra:=evala(randpoly(ext[nops(ext)],degree=2)) mod p; ii:=gcd_modp(f,evala(Expand(i[1]^i[2])) mod p,x,p); g:=factorise_modp(evala(Expand(subs(x=x+ra,ii))) mod p,x,ext,p); result:=result,seq([evala(Expand(subs(x=x-ra, dummy[1]))) mod p,dummy[2]],dummy=g) fi od; [result] end: # degree_ext gives the degree of the extension degree_ext:=proc(ext) local a,j,i,b; options remember; if nargs=2 then b:=args[2] else b:=1 fi; convert([seq(degree(subs(seq(ext[i-j]=a[i-j],j=1..i-1), op(ext[i])),_Z),i=b..nops(ext))],`*`) end: # input must be monic gcd_modp:=proc(a,b,x,p) if b=0 then RETURN(a) fi; gcd_modp(b,rem_modp(a,b,x,p),x,p) end: # input must be monic # output will also be monic (so this is no ordinary remainder) rem_modp:=proc(a,b,x,p) local aa,db; if b=1 then RETURN(0) fi; aa:=a mod p; db:=degree(b,x); while degree(aa,x)>=db do aa:=evala(Expand(aa-lcoeff(aa,x)*b*x^(degree(aa,x)-db))) mod p od; if aa<>0 then aa:=evala(aa/lcoeff(aa,x)) mod p fi; aa end: ########################################################## # algebraic numbers macro( g_evala=`nfintbas/g_evala`, g_evala_rem=`nfintbas/g_evala_rem`, g_expand=`nfintbas/g_expand`, truncate_subs=`nfintbas/truncate_subs`, v_ext_mult=`nfintbas/v_ext_mul`, g_solve=`nfintbas/g_solve`, truncate=`nfintbas/truncate`, g_conversion1=`nfintbas/g_conversion1`, g_conversion2=`nfintbas/g_conversion2`, g_normal=`nfintbas/g_normal` ): # evala replacement, ext contains a description of the algebraic extension # See IntBasis for a description of g_evala, g_expand etc. # Note the order of ext=[RootOf,..] is reversed g_evala:=proc(a,ext) local dummy,e,en; if nops(ext)=0 then RETURN(a) elif not type(a,polynom(anything,ext)) then ERROR(`bug alert, please mail hoeij@sci.kun.nl`) elif nops(ext)=1 then e:=ext[1]; expand(convert([seq(coeff(a,e,dummy)*g_evala_rem(e^dummy) ,dummy=0..degree(a,e))],`+`)) else e:=g_evala(a,[ext[1..nops(ext)-1]]); en:=ext[nops(ext)]; g_evala(expand(convert([seq(coeff(e,en,dummy)*g_evala_rem( en^dummy),dummy=ldegree(e,en)..degree(e,en))],`+`)), [ext[1..nops(ext)-1]]) fi end: # used by g_evala g_evala_rem:=proc() options remember; expand(subs(g_conversion1,evala(Expand(subs(g_conversion2,args))))) end: g_expand:=proc(a,ext) g_evala(expand(a),ext) end: # f is a polynomial in y. y_val is being substituted for y # This is being computed modulo p^n # Here x stands for ap^(1/ram) where ap has p-valuation 1. # Note ap^(1/ram) may not be simplified by expand, otherwise bug. truncate_subs:=proc(f,x,y,y_val,p,n,ext,ram,ap) local ym,i,result,y_value,pn; pn:=p^n; y_value:=subs(x=ap^(1/ram),y_val) mod pn; result:=0; for i from 0 to degree(f,y) do if i=0 then ym:=1 else ym:=g_expand(ym*y_value,ext) mod pn fi; result:=result+coeff(f,y,i)*ym od; subs({seq(ap^(i/ram)=x^i,i=1..ram-1)},result) mod pn end: g_conversion1:={}: # RootOf syntax -> my own syntax g_conversion2:={}: # my syntax -> RootOf syntax # Gives the zeros, multiplicities and extensions needed, for a # monic irreducible polynomial v_ext_mult:=proc(f,x,ext,p) global g_conversion1, g_conversion2; local r,result,i,ro,ii,vv; r:=subs(g_conversion1,factorise_modp(op(subs( g_conversion2,[f,x,ext,p])))); result:=NULL; for i in r do if degree(i[1],x)=1 then result:=result,[x-i[1],ext,i[2]]; else ro:=RootOf(i[1],x); # fix global variables for g_evala if not member(subs(g_conversion1,ro), {seq(`nfintbas/rootof`.ii,ii=0..nops(g_conversion2))}) then if g_conversion1={} then g_conversion1:=NULL fi; vv:=nops(g_conversion2); g_conversion1:=ro=`nfintbas/rootof`.vv,g_conversion1; g_conversion2:={`nfintbas/rootof`.vv=ro,op(g_conversion2)}; fi; ro:=subs(g_conversion1,ro); result:=result,[ro,[op(ext),ro],i[2]] fi od; [result] end: g_solve:=proc(eqn,p) msolve(eqn,p) end: # Expands modulo m truncate:=proc(aa,ext,ram,x,m,p,ap) g_expand(rem(expand(aa),x^ram-ap,x),ext) mod m end: g_normal:=proc(aa) local a; if indets(aa,RootOf)<>{} then RETURN(evala(Normal(aa))) fi; a:=subs(g_conversion2,aa); if indets(a,RootOf)={} then normal(a) else subs(g_conversion1,evala(Normal(a))) fi end: ########################################################## # p-adic numbers macro( padic=`nfintbas/padic`, valuation=`nfintbas/valuation`, cont_exp=`nfintbas/cont_exp`, el_of=`nfintbas/el_of` ): # [valuation,lowest coefficient] # valuation(p)=1 # valuation(xx)=1/ram valuation:=proc(a,xx,ram,p) local x,aa,result; if ram>1 then x:=xx fi; aa:=a; result:=0; while aa mod p=0 do aa:=aa/p;result:=result+1 od; aa:=aa mod p; [result+ldegree(aa,x)/ram,tcoeff(aa,x)] end: # Input f: monic polynomial over the integers in y. # p: prime number # ac: desired accuracy, ac is a number or the string `intbasis` # Output: the zeros of f over the algebraic closure of the p-adic numbers, # computed accurate enough for integral basis computation if ac=`intbasis` padic:=proc(f,y,p,ac) local x,v,i,result,dummy,m; # x denotes p^(1/ramification), ramification is now 1 if ac<>`intbasis` then RETURN(cont_exp([0,0,[],degree(f,y),1,x,0],f,y,p,ac)) fi; userinfo(1,nfintbas,`computing p-adic factorization in the prime`,p); v:=cont_exp([0,0,[],degree(f,y),1,x,0],f,y,p,0); result:=NULL; m:=max(seq(dummy[7]+(dummy[2]-1)/dummy[5],dummy=v)); for i in v do result:=result,op(cont_exp(i,f,y,p,m-i[7]+1/i[5])) od; userinfo(1,nfintbas,`upper bound for the multiplicity of`,p,`is `, trunc(max(seq(dummy[7],dummy=v)))); [result] end: # This procedure continues an expansion until it has multiplicity 1 # Input is a list v # v[1] = The expansion so far # v[2] = This expansion is determined modulo x^v[2] # where x=p^(1/ramification) # v[3] = The algebraic extensions in this expansion # v[4] = The multiplicity of this factor # v[5] = the ramification # v[6] = the x used for ramification (x=p^(1/ramification)) # This will also return (not needed in the input) # v[7] = sum of the valuations ( v(x^(1/d))=1/d ) of the differences # of this expansion with the other expansions + the accuracy # f = minimal polynomial of y over Q, must be monic with integer coefficients # ac = the desired accuracy (0 means splitting) # Comments: this code is derived from the puiseux algorithm in IntBasis cont_exp:=proc(v,f,y,p,ac) local x,t,n,a,i,ii,r,som,result,vv7,rp,tabel,ld; if v[5]>degree(f,y) then ERROR(`wild ramification: case not implemented`) fi; if (ac=0 and v[4]=1) or (ac>0 and v[2]/v[5]>=ac) then RETURN({v}) fi; result:={}; x:=v[6]; r:=rem(v[1]+a*x^v[2],x^v[5]-p,x); # I will find an equation for the unknown a i:=1; ii:=0; while ii=0 do i:=i+2; ii:=truncate_subs(f,x,y,r,p,ceil(v[2]/v[5]+v[7]+i),v[3],v[5],p) od; vv7:=valuation(ii,x,v[5],p); r:=vv7[2]; vv7:=vv7[1]-v[2]/v[5]; rp:=r mod p; ld:=ldegree(rp,a); if degree(rp,a)>0 then if v[4]=0 then rp:=expand(rp/a^ld) fi; r:=v_ext_mult(rp,a,v[3],p); for i in r do result:=result union cont_exp([v[1]+x^v[2]*i[1],v[2]+1,i[2], i[3],v[5],v[6],vv7],f,y,p,ac) od fi; som:=0; for ii in result do som:=som+ii[4]*degree_ext(subs(g_conversion2,ii[3])) od; som:=som/degree_ext(subs(g_conversion2,v[3])); n:=2; if som{} do for t from 1 to n-1 do if igcd(t,n)=1 and el_of(t/n,tabel) then rp:=cont_exp([subs(x=x^n,v[1]),(v[2]-1)*n+t,v[3],0 ,v[5]*n,v[6],vv7],f,y,p,ac); result:=result union rp[2]; som:=0; for i in rp[2] do som:=som+i[4]*degree_ext(subs(g_conversion2,i[3]))/ degree_ext(subs(g_conversion2,v[3])) od; tabel:=el_of(t/n,tabel,rp[1],som) fi od; n:=n+1; od; if v[4]=0 then result:=[ld,result] fi; result end: # used for deciding which t/n to use in cont_exp el_of:=proc(q,tabel,rp,som) local i,v,r; if nargs=2 then for i in tabel do if i[1]q then RETURN(true) fi od; RETURN(false) fi; for i in tabel do if i[1]q then v:=i fi od; r:=tabel minus {v}; if v[3]>rp+som then r:=r union {[v[1],q,v[3],rp+som]} fi; if rp>v[4] then r:=r union {[q,v[2],rp,v[4]]} fi; r end: ###################################################################### # nfintbas macro( local_intbas=`nfintbas/local_intbas`, monic=`nfintbas/monic` ): # Input, Output: see ?nfintbas # This procedure computes a set df, either obtained from the input # if an argument `local` is given, or df is the set of multiplicity > 1 # factors of the discrimint. For each prime in df a local integral # basis is computed by local_intbas. These bases are combined by # the procedure meng_bases. nfintbas:=proc(ff) local y,f,q,disc,tmp,df,i,opties,basis; options remember; y:=indets(ff); if nops(y)<>1 then ERROR(`wrong arguments`) fi; y:=op(y); f:=monic(ff,y,'q'); q:=eval(q); disc:=discrim(f,y); if not member(`local`,{args}) then userinfo(1,nfintbas,`factorizing the discriminant`,disc); tmp:=readlib(ifactors)(disc); df:=NULL; for i in tmp[2] do if i[2]>1 then df:=df,i[1] fi od; df:={df}; # df contains those factors k that appear more than once in the # discriminant discrim(f,y) else df:={}; for i in {args} do if type(i,integer) then df:=df union {i} elif type(i,set) then df:=df union i fi od fi; userinfo(1,nfintbas,`set of primes`,df); if member(roundtwo,{args}) then # Use the round two algorithm to compute local integral bases opties:=roundtwo fi; basis:=[seq(y^i,i=0..degree(f,y)-1)]; for i in df do # Compute local integral basis and mix it with basis. basis:=meng_bases(basis,local_intbas(f,y,i,disc,opties),y) od; subs(y=y*q,reduceer(basis,y)) end: # Input: a monic polynomial f in y, a prime number p, the discriminant # and a variable opties describing the options (which method to choose). # Output: a local integral basis # Comments: poorly documented, but the code and method resembles the # IntBasis code which is (a bit) better documented. See also IntBasis.tex local_intbas:=proc(f,y,p,disc,opties) local n,defect,i,basis,d,power,max_power,b,value_new_one,dummy ,values_basis_in_places,places,dummy2,equations,found_something; userinfo(1,nfintbas,`computing local integral basis in the prime`,p); n:=degree(f,y); if opties<>roundtwo then places:=traperror(padic(f,y,p,`intbasis`)) fi; # A wild ramification causes places=lasterror if places=lasterror or opties=roundtwo then # Call Beckker's round two userinfo(1,nfintbas,`calling roundtwo method`); RETURN(locale_ring_van_gehelen(p, floor(pvaluatie(disc,p)/2),[seq(y^i,i=0..n-1)],f,n,y)) fi; power:=1; # We try to put p^power in the denominator, if this succeeds # then increase power. max_power:=floor(max(seq(i[7],i=places))); defect:=p^max_power; # maximal possible denominator if defect=1 then RETURN([seq(y^i,i=0..n-1)]) fi; places:=[seq([i[1],i[3],i[5],i[6],defect,p],i=places)]; basis:=[1]; values_basis_in_places[1]:=[seq(1,i=places)]; for d from 2 to degree(f,y) do basis:=[op(basis),y*basis[d-1]]; # This basis is our first guess values_basis_in_places[d]:=[seq(truncate( values_basis_in_places[d-1][dummy]*places[dummy][1] ,places[dummy][2..6],p),dummy=1..nops(places))]; found_something:=true; while found_something and power<=max_power 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[dummy2]* values_basis_in_places[dummy2][dummy],dummy2=1..d)],`+`) ,places[dummy][2..6],p),dummy=1..nops(places))]; equations:=g_solve({seq(coeffs(dummy mod p,indets(dummy) minus {seq(b[dummy2],dummy2=1..d-1)}),dummy=value_new_one)},p); # Now we know what values b[1] .. b[d] must have if equations=NULL then found_something:=false else values_basis_in_places[d]:=[seq(dummy/p,dummy=subs( equations,value_new_one))]; basis:=[basis[1..d-1],subs(equations,convert([seq( b[dummy]*basis[dummy],dummy=1..d)],`+`))/p]; power:=power+1 fi od od; expand(basis) end: # Makes f monic, returns also q, the factor needed to multiply y by, to # get it integral. monic:=proc(f,y,q) local ff,qq; ff:=numer(f); ff:=ff/icontent(ff); qq:=lcoeff(ff,y); q:=qq; ff:=subs(y=y/qq,ff); ff/lcoeff(ff,y) end: