# $Source: /u/maple/research/lib/numtheory/integral_basis/src/RCS/nfintbas,v $ # $Notify: mvanhoei@daisy.uwaterloo.ca $ # 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", JSC 1994 # 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 # Calling sequence: # nfintbas(f,x) where f \in Q[x] # gives the integral basis # nfintbas(f,x,`local`,{set of primes}) # this tells nfintbas to look only at that given set of primes, instead # of factoring the discriminant. # An extra option "roundtwo" causes nfintbas to use the first method. This # is useful for testing the code, compute it first with, then without this # option, and compare map(denom, output) (these should be the same). # A test example: # # f:=x^12-5476*x^9-151959*x^8+11244966*x^6-3328509936*x^5+7697179227*x^4- # 10262905636*x^3-2847956313990*x^2-42149753447052*x-126449260341156: # # v1:=nfintbas(f,x); v2:=nfintbas(f,x,roundtwo); # if map(denom,v1)<>map(denom,v2) then not OK. ################################################################# # Round two implementation macro( transpose=LinearAlgebra:-LA_Main:-Transpose, alfa=`numtheory/integral_basis/nfintbas/alfa`, pvaluatie=`numtheory/integral_basis/nfintbas/p_valuation`, meng_bases=`numtheory/integral_basis/nfintbas/mix_bases`, reduceer=`numtheory/integral_basis/nfintbas/reduce`, basis_radical=`numtheory/integral_basis/nfintbas/basis_radical`, idealizer=`numtheory/integral_basis/nfintbas/idealizer`, powmodf=`numtheory/integral_basis/nfintbas/powmodf`, write_in=`numtheory/integral_basis/nfintbas/write_in`, check=`numtheory/integral_basis/nfintbas/check`, a=`numtheory/integral_basis/nfintbas/a`, correct_basis=`numtheory/integral_basis/nfintbas/correct_basis`, minimal_a_i=`numtheory/integral_basis/nfintbas/minimal_a_i`, times=`numtheory/integral_basis/nfintbas/times`, square=`numtheory/integral_basis/nfintbas/square`, local_intbas_R2=`numtheory/integral_basis/nfintbas/local_intbas_R2`, corrigeer_basis=`numtheory/integral_basis/nfintbas/correct_basis`, local_intbas_R4=`numtheory/integral_basis/rd4` ): pvaluatie:=proc(a,p) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if a=0 then infinity elif irem(a,p)=0 then 1+pvaluatie(a/p,p) else 0 fi end: local_intbas_R2:=proc(p,t,R,f,n,x) local V,W,IdW,i,k,d,oude_d; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; oude_d:=1; V:=R; k:=1; while p^k < n do k:=k+1 od; do W:=basis_radical(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 d=p^t or d=oude_d then return IdW else V:=eval(IdW); oude_d:=d fi od end: meng_bases:=proc(basis1,basis2,x) local i,result,s,t; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; result:=NULL; for i to nops(basis1) do if igcdex(denom(lcoeff(basis1[i],x)),denom(lcoeff(basis2[i],x)),'s','t')<>1 then error "denominators do not have gcd 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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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:=[op(W[1..i-1]),w,op(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:=[op(W[1..j-1]),l,op(W[j+1..nops(W)])] od; W end: basis_radical:=proc(p,t,R,f,k,n,x) global a; local b,i,j,l,m,M,B,D,J,pt,z; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; a:=0; pt:=p^t; for i to nops(R) do #a:=collect(plus(a,z[i]*powmodf(pt*R[i],p^k,f,p,t)),alfa) a:=collect(a+z[i]*powmodf(pt*R[i],p^k,f,p,t,x),x) od; a:=a/pt; b:=write_in(a,R,x); M:=Matrix(n,n,[seq([seq(coeff(b[j],z[i]),j=1..n)],i=1..n)],scan=[columns]); 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:=[op(D[1..l-1]),D[l]+B[l][m]*R[m],op(D[l+1..nops(B)])] od; od; D:=check(D,p,x); corrigeer_basis(J,D,x) end: idealizer:=proc(V,f,p,t,n,x) # global a; local a,c,W,Z,M,d,g,h,i,j,k,l,q,pivot,Y,u,z,pt; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; pt:=p^t; g:=sum('a[z]*x^z','z'=0..n-1); W:={}; for i to n do h:=collect(times(g,pt*V[i],f,p,t,x),x); for j to n do 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],x) od od; Z:=[]; for k to n do q:=minimal_a_i(W,k-1,p,t,a); 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 := LinearAlgebra:-LA_Main:-MatrixInverse(M,method=none,conjugate=false, outputoptions=[]); transpose(M,inplace=true,outputoptions=[]); d:=Vector([seq(x^i,i=0..n-1)]); [seq((i*pt mod pt*p)/pt,i=convert(LinearAlgebra:-LA_Main:-MatrixVectorMultiply( M,d,outputoptions=[]),list))] end: powmodf:=proc(getal,exp,pol,p,t,x) local pt; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; pt:=p^(t+1); if exp = 2 then square(getal,pol,p,t,x) mod pt elif exp mod 2 = 0 then square(powmodf(getal,exp/2,pol,p,t,x),pol,p,t,x) mod pt else times(powmodf(getal,exp-1,pol,p,t,x),getal,pol,p,t,x) mod pt fi end: write_in:=proc(f,V,x) local d,g,i,l,m; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; g:=f; l:=[seq(0,i=1..nops(V))]; while g <> 0 do d:=degree(g,x)+1; m:=lcoeff(g,x)/lcoeff(V[d],x); l:=[op(l[1..d-1]),m,op(l[d+1..nops(V)])]; g:=collect(g-m*V[d],x) od; l end: check:=proc(V,p,x) local W,d1,d2,i,j,k; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; W:=V; for i to nops(V)-1 do d1:=degree(V[i],x); for j from i+1 to nops(V) do d2:=degree(V[j],x); if d1 = d2 then k:=lcoeff(V[j],x)/lcoeff(V[i],x) mod p; W:=[op(W[1..j-1]),V[j]-k*V[i] mod p,op(W[j+1..nops(V)])] fi od od; convert(W,list) end: corrigeer_basis:=proc(V,W,x) local d,h,i,R; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; R:=V; for i to nops(W) do h:=W[i]; d:=degree(h,x)+1; R:=[op(R[1..d-1]),h,op(R[d+1..nops(V)])] od; R end: minimal_a_i:=proc(W,k,p,t,a) # global a; local i,j,c; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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; p^t*a[k],0,0 end: times:=proc(a,b,f,p,t,x) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; rem(expand(a*b/p^t),f,x) end: square:=proc(getal,f,p,t,x) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; times(getal,getal,f,p,t,x) end: ########################################################## # Factorization in finite fields macro( factorise_modp=`numtheory/integral_basis/nfintbas/factorise_modp`, degree_ext=`numtheory/integral_basis/nfintbas/degree_ext`, gcd_modp=`numtheory/integral_basis/nfintbas/gcd_modp`, rem_modp=`numtheory/integral_basis/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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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:=procname(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, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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=`numtheory/integral_basis/nfintbas/g_evala`, g_evala_rem=`numtheory/integral_basis/nfintbas/g_evala_rem`, g_expand=`numtheory/integral_basis/nfintbas/g_expand`, truncate_subs=`numtheory/integral_basis/nfintbas/truncate_subs`, v_ext_mult=`numtheory/integral_basis/nfintbas/v_ext_mul`, g_solve=`numtheory/integral_basis/nfintbas/g_solve`, truncate=`numtheory/integral_basis/nfintbas/truncate`, g_conversion1=`numtheory/integral_basis/nfintbas/g_conversion1`, g_conversion2=`numtheory/integral_basis/nfintbas/g_conversion2`, g_normal=`numtheory/integral_basis/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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if nops(ext)=0 then return a elif not type(a,polynom(anything,ext)) then error "assertion failed" 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, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; expand(subs(g_conversion1,evala(Expand(subs(g_conversion2,args))))) end: g_expand:=proc(a,ext) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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(subs(g_conversion2,i[1]),x); # fix global variables for g_evala if not member(subs(g_conversion1,ro), {seq(`numtheory/integral_basis/nfintbas/rootof` || ii,ii=0..nops(g_conversion2))}) then if g_conversion1={} then g_conversion1:=NULL fi; vv:=nops(g_conversion2); g_conversion1:=ro=`numtheory/integral_basis/nfintbas/rootof` || vv,g_conversion1; g_conversion2:={`numtheory/integral_basis/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) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; msolve(eqn,p) end: # Expands modulo m truncate:=proc(aa,ext,ram,x,m,p,ap) option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; g_expand(rem(expand(aa),x^ram-ap,x),ext) mod m end: g_normal:=proc(aa) local a; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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=`numtheory/integral_basis/nfintbas/padic`, valuation=`numtheory/integral_basis/nfintbas/valuation`, cont_exp=`numtheory/integral_basis/nfintbas/cont_exp`, el_of=`numtheory/integral_basis/nfintbas/el_of` ): # [valuation,lowest coefficient] # valuation(p)=1 # valuation(xx)=1/ram valuation:=proc(a,xx,ram,p) local x,aa,result; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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 option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if ac<>` _intbasis` then return cont_exp([0,0,[],degree(f,y),1,x,0],f,y,p,ac) fi; userinfo(1,integral_basis,`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,integral_basis,`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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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=`numtheory/integral_basis/nfintbas/local_intbas`, monic=`numtheory/integral_basis/nfintbas/monic` ): # 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. `numtheory/integral_basis/nfintbas`:=proc(ff) local y,f,q,disc,tmp,df,i,opties,basis,opt; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; 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); ##### Modification to merge D. Ford's code start here. MJR # #if not member(`local`,{args}) then # userinfo(1,integral_basis,`factorizing the discriminant`,disc); # tmp:=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; # # # Process options # for opt in [args[2..nargs]] do if type(opt,identical('primes')=anything) then if assigned(df) then error "redundant option %1",opt fi; if not type(op(2,opt), {posint,set(posint),list(posint)}) then error "%1 are invalid primes", op(2,opt) fi; df := convert(op(2,opt),set); for i in df do if not isprime(i) then error "%1 is not a prime",i fi od; df := select(t->evalb(irem(disc,t^2)=0),df); elif type(opt,identical('method')=anything) then if assigned(opties) then error "%1 is a redundant option", opt fi; if op(2,opt)='roundtwo' or op(2,opt)='round2' then opties := 'roundtwo' elif op(2,opt)='roundfour' or op(2,opt)='round4' then opties := 'roundfour' elif op(2,opt)='p_adic' then opties := 'p_adic' else error "%1 is an unknown method", op(2,opt) fi else error "%1 is an invalid option", opt fi od; # # Default values # if not assigned(df) then userinfo(1,integral_basis,`factorizing the discriminant`,disc); tmp:=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) fi; if not assigned(opties) then opties := 'roundfour' fi; # # # userinfo(1,integral_basis,`set of primes`,df); ### #if member(roundtwo,{args}) then # # Use the round two algorithm to compute local integral bases # opties:=roundtwo #fi; # #### Modifications to merge D. Ford's code end here. MJR 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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; userinfo(1,integral_basis,`computing local integral basis in the prime`,p); n:=degree(f,y); #### Modifications to merge D. Ford's code start here. MJR # # if opties='roundfour' then # Call D. Ford's Round 4 userinfo(1,integral_basis,'`calling Round 4 method`'); return local_intbas_R4(p,f,pvaluatie(disc,p),y) elif opties='roundtwo' then userinfo(1,integral_basis,'`calling Round 2 method`'); return local_intbas_R2(p, floor(pvaluatie(disc,p)/2),[seq(y^i,i=0..n-1)],f,n,y) fi; places:=traperror(padic(f,y,p,` _intbasis`)); # A wild ramification causes places=lasterror if places=lasterror then WARNING(cat("wild ramification detected at ", convert(p,string), ". Switching to Round 4 method.")); userinfo(1,integral_basis,'`calling Round 4 method`'); return local_intbas_R4(p,f,pvaluatie(disc,p),y) fi; userinfo(1,integral_basis,'`calling p_adic method`'); # # #### Modifications to merge D. Ford's code end here. MJR 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] ,op(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)],`+`) ,op(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:=[op(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; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; ff:=numer(f); ff:=ff/icontent(ff); qq:=lcoeff(ff,y); q:=qq; ff:=subs(y=y/qq,ff); ff/lcoeff(ff,y) end: #savelib('`numtheory/integral_basis/nfintbas`', \ 'g_conversion1','g_conversion2', \ 'pvaluatie','local_intbas_R2','meng_bases','reduceer','basis_radical','idealizer', \ 'powmodf','write_in','check','corrigeer_basis','minimal_a_i','times','square', \ 'factorise_modp','degree_ext','gcd_modp','rem_modp','g_evala','g_evala_rem', \ 'g_expand','truncate_subs','v_ext_mult','g_solve','truncate','g_normal', \ 'valuation','padic','cont_exp','el_of','local_intbas,monic', \ '`numtheory/integral_basis/nfintbas.m`'):