# AlgDyn: MAPLE TOOLS FOR ALGEBRAIC AND DISCRETE DYNAMICS # # F. Vivaldi. Version 2c, May 2006 # # The procedures' names are available both in mixed case # (Mathematica style) and in lower case (Maple style), # e.g. IrrPoly, irrpoly # This library makes use of selected Maple libraries with(linalg): with(numtheory): # ===================================================== CAPTIONS # f mapping [polynomial or procedure] # x indeterminate # L,S,LS a list, a set, a list or set # p characteristic # n period, or degree of algebraic extension # ===================================================== CONTENTS # ----- GENERAL MAPLE (General) # Flat(LS) flatten a list or a set # LShift(L) left shift of a list # RShift(L) right shift of a list # Freq(L) frequency of occurrences of the elements of a list # isMinPeriod(L) TRUE if L -as a periodic list- has minimal period. # Int2Dig(m,b) list of digits of integer m to the base b # Dig2Int(L,b) integer having L as list of digits to the base b # Int2Fibo(n) Fibonacci expansion of a positive integer # Coeff2Poly([i],x) convert list of coefficients into a polynomial in x # ----- NUMBER THEORY (NumThy) # pvalue(r,p) p-adic value of a rational r (or a list/set of rationals) # punit(r,p) unit part of a rational r # Artin(f,p) Artin map: cycle decomposition of f(x) modulo p # the output is a list of inertia degrees # a -ve degree signals ramification # MiniPoly(a,f1,...,fn,x) # minimal polynomial in x # of the polynomial a(x_1,..,x_n), where the # polynomials f_k are such that f_k(x_k)=0 # Eplus(P,Q,E) P+Q on the elliptic curve E=[a,b] => y^2=x^3+a*x+b # ReduForms(d) reduced forms of -ve discriminant d # ReduForm(a,b,c) reduced form equivalent to a +ve def. form (a,b,c) # ----- ALGEBRAIC & DISCRETE DYNAMICS (AlgDyn) # PhiPoly(f,n,x) n-th Phi polynomial of f (roots are cycles of essential period n) # MulPoly(f,n,x) n-th multiplier polynomial (roots are multipliers of n-cycles) # MatMult(f,g,x) matrix of multiplication by a polynomial # LocMatMult(f,g,p,x) same as MatMult, but in characteristic p # Wop(f,g,p,x) local wedge operator of two polynomials f and g # GWop(f,g,x) global wedge operator (simple roots) # GWopBar(f,g,x) global wedge operator (multiple roots) # Orbit(f,ic) orbit of map f, with initial condition ic # SOrbit(f,ic) ditto + handle points where f is undefined. # IOrbit(f,ic) same as Orbit, with f invertible # Orbits(f,L) orbits of f: L -> L # IOrbits(f,L) same as Orbits, with f invertible # Graph(f,L) graph of map f (input for networks[graph]) # Components(f,L) components of map f over L # ----- DATABASE: IRREDUCIBLE POLYNOMIALS OVER FINITE FIELDS (IrrPoly) # - - - - - Database # IrrPoly[p] irreducible polynomials over Fp (all) # IrrPoly[p][n] irreducible polynomials over Fp of degree n # - - - - - Procedures for database # LoadIrrPoly() load database # NumIrrPoly(p,n) number of irreduc. polynomials over Fp of degree n # NumIrrPolyFpn(p,n) number of irreduc. polys over Fp of degree dividing n # PolyRanges(p,n) index ranges for polys over Fp of degree dividing n # Int2Poly(i,p,n,x) convert lexicographical index into a polynomial # Poly2Int(f,p,x) convert a polynomial into a lexicographical index # - - - - - Content of database # Irreducible polynomials of degree n over GF(p) # For each p, the maximum degree corresponds to the # largest n for which there are less than 1000 polys. # p n # ----------- # 2 1-13 # 3 1-8 # 5 1-5 # 7 1-4 # 11 1-3 # 13 1-3 # 17 1-2 # 19 1-2 # 23 1-2 # 29 1-2 # 31 1-2 # 37 1-2 # 41 1-2 # 43 1-2 # LoadIrrPoly_____________________________________________________IrrPoly # load database of irreducible polynomials # > LoadIrrPoly(); LoadIrrPoly:=proc() global IrrPoly: if not assigned(IrrPoly) then read(`/home/fv/public_html/database/irrpoly`) fi end: # IrrPoly_________________________________________________________IrrPoly # extract irreducible polynomials over Fp from database # >IrrPoly[p]; # INPUT: p:prime - the characteristic # OUTPUT: a list of lists of lists of expr sequences: # [[polys of deg 1],[polys of deg 2],...] # where [polys of deg k] = [[poly1],[poly2],...], # where [poly n] = [high coeff, ... ,low coeff] # FUNCTIONS: none # SYNOPSIS: # This is the interface to the IrrPoly database. # List of coefficients of irreducible polynomials over Fp # (ordered lexicographically, as in Lidl and Niederreiter) # EXAMPLES: nops(IrrPoly[2]); -> 13 # Polys over GF(2) have degree up to 13. # NumIrrPoly,NumIrrPolyFpn________________________________________IrrPoly # number of irreduc. polynomials over Fp of degree n (of degree dividing n) # >NumIrrPoly(p, n); # >NumIrrPolyFpn(p, n); # INPUT: p:prime - the characteristic # n:integer - the degree of the extension # OUTPUT: an integer # FUNCTIONS: NumIrrPoly # SYNOPSIS: # Number of irr. polys over Fp of degree n or dividing n. alias(numirrpoly=NumIrrPoly): NumIrrPoly := proc(p::prime, n::posint) local d: 0: for d in divisors(n) do %+mobius(n/d)*p^d od: %/n end: alias(numirrpolyfpn=NumIrrPolyFpn): NumIrrPolyFpn := proc(p::prime, n::posint) global NumIrrPoly: local d: 0: for d in divisors(n) do % + NumIrrPoly(p,d) od end: # PolyRanges______________________________________________________IrrPoly # ranges of lexicographical indices for polynomials over Fp of degree dividing n # >PolyRanges(p, n); # INPUT: p:prime - the characteristic # n:integer - the degree of the extension # OUTPUT: A table containing a list of ranges d = min..max, # for each divisor d of n. # FUNCTIONS: none # SYNOPSIS: # Partition of numbers of irreducible polynomials in Fp^n, # by degree. alias(polyranges=PolyRanges): PolyRanges := proc(p::prime, n::posint) local divn, i, ndivn, npolys, part, ranges, sumpart: #--------------------------------------- initialize divn := [op(divisors(n))]: ndivn := nops(divn): npolys := [seq(NumIrrPoly(p,i),i=1..n)]: #--------------------------------------- partition of the order of Vn part := [seq(npolys[i], i=sort(divn))]: #--------------------------------------- sum of elements of partition sumpart := [part[1]]: for i from 2 to ndivn do sumpart := [op(sumpart), sumpart[i-1] + part[i]] od: #-------------------------------------- ranges ranges := table([1=1..sumpart[1], seq(divn[i]=sumpart[i-1]+1..sumpart[i], i=2..ndivn)]): eval(ranges) end: # Coeff2Poly______________________________________________________General # convert list of coefficients into a polynomial # >Coeff2Poly(polyst, x); # INPUT: polyst:list - a list of coefficients # x:symbol - an indeterminate # OUTPUT: a polynomial # FUNCTIONS: none # SYNOPSIS: # Convert list of coefficients of a polynomial # into a polynomials in the indeterminate x alias(coeff2poly=Coeff2Poly): Coeff2Poly := proc(polylist::list, x::symbol) local i, result, degplus1: result := 0: degplus1 := nops(polylist): for i to degplus1 do result := result + polylist[i] * x^(degplus1 - i) od: result end: # Flat____________________________________________________________General # flatten a list or a set (Written by F J Wright). # >Flat(a); # INPUT: anything (typically a list or a set) # OUTPUT: an expression sequence # FUNCTIONS: none # SYNOPSIS: # Flatten the argument sequence i.e. [[a,[[b]]],c] -> a,b,c alias(flat=Flat): Flat := proc(a) if(nargs=1) then if(type(a,{list,set})) then Flat(op(a)) else a fi else op(map(Flat,[args])) fi end: # Int2Dig,Dig2Int,Int2FIbo________________________________________General # convert from b-adic digits to integer, and vice-versa # >Int2Dig(i,b); # >Int2Dig(i,b,nd); # >Dig2Int(L,b); # >Int2Fibo(n); alias(int2dig=Int2Dig): int2dig:=proc(x::nonnegint,b::posint,nd::nonnegint) local n, i: if x=0 then if nargs = 2 then return [0] else return([0$nd]) fi: fi: if nargs = 2 then n:=max(0,floor(evalf(log[b](abs(x)))) + 1) else n:=nd fi: [seq(irem(iquo(x,b^(i-1)),b),i=1..n)] end: alias(dig2int=Dig2Int): dig2int:=proc(L::list(nonnegint),b::posint) local i: add(op(i,L)*b^(i-1),i=1..nops(L)): end: Int2Fibo:=proc(n::posint) local i, idig, F, FD, x, nFibo: F:=[1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597,2584,4181,6765,10946, 17711,28657,46368,75025,121393,196418,317811,514229,832040,1346269, 2178309,3524578,5702887,9227465,14930352,24157817,39088169,63245986,102334155,165580141]: if n>F[-2]+F[-1] then ERROR("not enough Fibonacci numbers") fi: x:=n: for i while F[i]<=x do od: nFibo:=i-1:FD:=array(1..nFibo,[0$nFibo]): idig:=nFibo:#FD[idig]:=1: while x>0 do FD[idig]:=1: x:=x-F[idig]: for i while F[i]Freq(L); # INPUT: L:list # OUTPUT: a list of lists # FUNCTIONS: none # SYNOPSIS: # Compute the frequency of occurrence of the elements of a list. # The number of operands is equal to the number of distinct # elements of the input, which are sorted. # Each operand consists of the given element, with its frequency. alias(freq=Freq): Freq:=proc(L::list) local sortOpsL, x: sortOpsL:=sort([op({op(L)})]): [seq([x, nops(select((l,x)->evalb(l=x), L, x))], x=sortOpsL)] end: # Int2Poly,Poly2Int_______________________________________________General # convert lexicographical index into a polynomial, and vice-versa # >Int2Poly(i, p, n, x); # >Poly2Int(f, p, x); # INPUT: i:integer - index identifying the polynomial # p:prime - the characteristic # n:integer - degree of the extension. # x:symbol - an indeterminate # f:polynom - a polynomial # OUTPUT: a polynomial (an integer) # FUNCTIONS: Coeff2Poly, PolyRanges # SYNOPSIS: # Convert the integer i into the corresponding irreducible # polynomial of degree n over Fp, in the indeterminate x, # determined according to the lexicographical order of the # irreducible polynomials of degree dividing n. alias(int2poly=Int2Poly): Int2Poly := proc(i::integer, p::prime, n::posint, x::symbol) global Coeff2Poly, PolyRanges: local deg, deglist, coefflist, j, rgtab: # ---------------------------- find range table and list of degrees rgtab := PolyRanges(p,n): deglist := map(op,[indices(rgtab)]): #--------------------------------------------- check input: i if(i < 1 or i > op(2,rgtab[deglist[nops(deglist)]])) then ERROR(`i: out of bounds`) fi: #--------------------------------------------- find divisor for deg in deglist while(not member(i,{$rgtab[deg]})) do od: #--------------------------------------------- load list of coefficients coefflist:=IrrPoly[p][deg]: #--------------------------------------------- find index in list j := i - op(1,rgtab[deg]) + 1: Coeff2Poly(coefflist[j], x) end: alias(poly2int=Poly2Int): Poly2Int := proc(f::polynom, p::prime, x::symbol) local coefflist, floc, i, n, irrpolylocal: floc := collect(f, x) mod p: n := degree(floc, x): coefflist := [1]: for i from n-1 to 0 by -1 do coefflist := [op(coefflist), coeff(floc,x,i)] od: irrpolylocal := IrrPoly[p][n]: for i while(coefflist <> irrpolylocal[i]) do od: i + op(1, PolyRanges(p,n)[n]) - 1 end: # LShift,Rshift___________________________________________________General # left and right shift of a list # >LShift(L); # >RShift(L); # INPUT: lst:list - a list # OUTPUT: a shifted list # FUNCTIONS: none # SYNOPSIS: # Shift the elements of a list to the left/right. alias(lshift=LShift): LShift := proc(L::list) local i, n: n := nops(L): [seq(L[i],i=2..n), L[1]] end: alias(rshift=RShift): RShift := proc(L::list) local i, n: n := nops(L): [L[n], seq(L[i],i=1..n-1)] end: # isMinPeriod_____________________________________________________General # check whether of not a periodic sequence has minimal period # >isMinPeriod(L); alias(isminperiod=isMinPeriod): isMinPeriod:=proc(L::list) local n, t, Lnow: n:=nops(L): Lnow:=[op(L[2..n]),L[1]]: for t while L <> Lnow do Lnow:=[op(Lnow[2..n]),Lnow[1]] od: evalb(t=n) end: # # pvalue__________________________________________________________NumThy # p-adic value pvalue:=proc(r::{integer,fraction,list},p::integer) local n, d, nd, k, pk, segno: if type(r,list) then return min(op(map(pvalue,r,p))) fi: if r=0 then return(infinity) fi: n:=numer(r): d:=denom(r): if igcd(n,p)=1 then nd:=d: segno:=-1 else nd:=n: segno:=1 fi: pk:=p: for k while igcd(nd,pk) = pk do pk:=pk*p od: return (k-1)*segno end: # punit___________________________________________________________NumThy # p-adic unit part of a rational punit:=proc(r::{integer,fraction,list,set},p::integer) global pvalue: if type(r,list) then return [punit(r[1],p),punit(r[2],p)] elif type(r,set) then return {punit(r[1],p),punit(r[2],p)} elif r=0 then return(0) else return(r/p^pvalue(r,p)) fi: end: # Artin___________________________________________________________NumThy # Artin - Artin map: cycle decomposition of f modulo p into irreducbles # >Artin(f,p); # INPUT: f:polynomial - a polynomial over Q # p:prime - the characteristic # OUTPUT: A list of lists containing sorted inertia degrees of the # irreducible factors of f(x) modulo p. # FUNCTIONS: none alias(artin=Artin): Artin:=proc(f::polynom(integer),p::prime) local x, d, fmodp, g, L: x:=indets(f): if not nops(xnow)=1 then ERROR(f,`is not a univariate polynomial`) fi: d:=degree(f,x): fmodp:=Factor(f) mod p: if type(fmodp,`+`) or type(fmodp,symbol) then [d] elif type(fmodp,`*`) then L:=NULL: for g in op(fmodp) do if type(g,`+`) or type(g,symbol) then L:=L,degree(g,x) else L:=L,-degree(op(1,g)) fi od: sort([L]): elif type(fmodp,`^`) then sort([-degree(op(1,fmodp),x)]) fi end: # Eplus___________________________________________________________NumThy # addition on an elliptic curve over Q # >Eplus(P,Q,E); # INPUT: # P,Q: list of two rationals (including infinity) # E: list(integers) the elliptic curve: E=[a,b] => y^2=x^3+a*x+b # OUTPUT: a list of two rationals, representing P+Q Eplus:=proc(P::{list(rational),infinity}, Q::{list(rational),infinity}, E::list(integer)) local x1,x2,x3,y1,y2,y3,r: if P=infinity then Q elif Q=infinity then P else x1:=P[1]: y1:=P[2]: x2:=Q[1]: y2:=Q[2]: if x1<>x2 then r:=(y2-y1)/(x2-x1): x3:=r^2-x1-x2: y3:=-r*x3-(x2*y1-x1*y2)/(x2-x1): [x3,y3] elif y1=y2 then r:=(3*x1^2+E[1])/(2*y1): x3:=r^2-2*x1: y3:=r*(x1-x3)-y1: [x3,y3] else infinity fi fi end: # MiniPoly________________________________________________________NumThy # minimal polynomial of a(x_1,...,x_n), where f_k(x_k)=0 # >MiniPoly(a,f1,...,fn,x); # SYNOPSIS: # Minimal polynomial of an element of an algebraic extension # EXAMPLES: > MiniPoly(a+b,a^2-2,b^2-3,x); # 1-10x^2+x^4 # this is the minimal polynomial of sqrt(2)+sqrt(3). alias(minipoly=MiniPoly): MiniPoly:=proc() local a, jcol, coe, dnow, fnow, inow, qnow, xnow, B, D, E, F, Q, X, d, e, i, j, k, m, n, M: #-------------------------- process arguments if nargs < 2 then ERROR(`USAGE: MiniPoly(a,f1,...,fn,x)`) fi: #------------------------------ 1st argument a:=args[1]: if not type(a,polynom(integer)) then ERROR(`first argument must be polynomial over Z`) fi: #----------------------------- other arguments n:=nargs-2: F:=NULL: D:=NULL: X:=NULL: for k to n do fnow:=args[k+1]: xnow:=indets(fnow): if not nops(xnow)=1 then ERROR(fnow,`is not a univariate polynomial`) elif not type(fnow,polynom(integer)) then ERROR(fnow,`is not a polynomial over Z`) fi: xnow:=op(xnow): if k > 1 and member(xnow,{X}) then ERROR(`two identical indeterminates`) fi: F:=F,fnow: X:=X,xnow: D:=D,degree(fnow,xnow) od: D:=[D]: F:=[F]: X:=[X]: Q:=[1,seq(mul(d,d=D[1..k]),k=1..n-1)]: m:=mul(d,d=D): M:=array(1..m,1..m): #--------------------------------- construct basis B B:=NULL: for i to m do inow:=i-1: E:=NULL: for k from n to 1 by -1 do qnow:=Q[k]: e:=iquo(inow,qnow): inow:=inow-e*qnow: E:=e,E: od: E:=[E]: B:=B,mul(X[k]^E[k],k=1..n) od: B:=[B]: #----- construct matrix M of multiplication by a wrt basis B for j to m do jcol:=a*B[j]: for k to n do jcol:=rem(jcol,F[k],X[k]) od: for i to m do coe:=jcol: for k to n do dnow:=degree(B[i],X[k]): coe:=coeff(coe,X[k],dnow) od: M[i,j]:= coe od od: linalg[minpoly](M,args[nargs]) end: # ReduForm,ReduForms______________________________________________NumThy # the reduced quadratic form equivalent to a given +ve definite form # reduced quadratic forms of -ve discriminant d # >ReduForm(a,b,c); # >ReduForms(d); ReduForm:=proc(A,B,C) local D, zerouno, a, b, c, cnew, k; a:=abs(A): b:=sign(A)*B: c:=sign(A)*C: D:=b^2-4*a*c: zerouno:=D mod 4: if zerouno > 1 or D > 0 then ERROR(`the argument must be a negative discriminant`): fi: while not ((c>a and -a < b and b <= a) or (c=a and 0<=b and b<=a)) do (c-b)/(2*c): k:=floor(%): if k = %% then k:=k-1 fi: cnew:=a+b*k+c*k^2: a:=c: b:=-b-2*k*c: c:=cnew: od: [a,b,c] end: ReduForms:=proc(d::negint) local A, B, C, AC, L, maxb2, Arange, Brange, k, zerouno: zerouno:=d mod 4: if zerouno > 1 then ERROR(`the argument must be a negative discriminant`): fi: #---------------------------- range of B-values maxb2:=abs(d)/3: zerouno: for k while (2*k+zerouno)^2 <= maxb2 do %,2*k+zerouno od: Brange:={%}: L:=NULL: for B in Brange do AC:=(-d+B^2)/4: Arange:=sort([op(numtheory[divisors](AC))]): for A in Arange while A <= AC/A do C:=AC/A: #------------------------ 0 <= B is already satisfied if B <= A then L:=L,[A,B,C]: if 0 < B and B < A and A < C then L:=L,[A,-B,C] fi fi od od: [L] end: # MulPoly_________________________________________________________AlgDyn # multiplier polynomial # >MulPoly(f,n,x); # INPUT: f:polynom - polynomial # n:integer - the period # x:symbol - an indeterminate # OUTPUT: a polynomial. # FUNCTIONS: MiniPoly, PhiPoly # SYNOPSIS: # Polynomial whose roots are the multipliers of the n-cycles of f(x) alias(mulpoly=MulPoly); MulPoly:=proc(f::polynom,n::posint,x::symbol) local df, multiplier, xnow: global MiniPoly, PhiPoly: df:=diff(f,x): multiplier:=df: xnow:=x: to n-1 do xnow:=subs(x=xnow,f): multiplier:=multiplier*subs(x=xnow,df) od: MiniPoly(multiplier,PhiPoly(f,n,x),x) end: # Orbit___________________________________________________________AlgDyn # orbit of map f with initial condition ic # >Orbit(f,ic); # OUTPUT: a list of two lists [[transient],[limit_cycle]] Orbit:=proc(f::procedure,ic) local L, x, t: L:=NULL: x:=ic: while not member(x, [L]) do L:=L,x: x:=f(x) od: L:=[L]: for t while x <> L[t] do od: [L[1..t-1],L[t..nops(L)]] end: # SOrbit__________________________________________________________AlgDyn # orbit of map f with initial condition ic # >SOrbit(f,ic); # OUTPUT: a list of two lists [[transient],[cycle]] # [transient] is empty of the orbit is periodic # [cyle] is empty if the last point is singular # SOrbit:=proc(f::procedure,ic) local L, x, y, t: L:=NULL: x:=ic: while not member(x, [L]) do L:=L,x: y:=f(x): if y <> NULL then x:=y fi od: L:=[L]: if y=NULL then [L,[]] else for t while x <> L[t] do od: [L[1..t-1],L[t..nops(L)]] fi end: # IOrbit_____________________________________________________________AlgDyn # orbit of an invertible map f with initial condition ic # >IOrbit(f,ic); # OUTPUT: a list [cycle] IOrbit:=proc(f::procedure,ic) local T, X, t: t:=1: T:=table([t=ic]): X:=f(ic): while X <> ic do t:=t+1: T[t]:=X: X:=f(X) od: return([seq(T[i],i=1..t)]): end: # Orbits__________________________________________________________AlgDyn # global orbit structure of map f over phase space L # >Orbits(f,L); # OUTPUT: a list of lists [...,[{basin_k},[cycle_k]],...] # basin_k is the basin of attraction of cycle_k Orbits:=proc(f::procedure,L::list) local basin, cycle, Cycles, nCycles, Lnew, Lold, G, compG, Os, i: global IOrbits, Graph: #----------------------- limit Cycles Lold:={op(L)}: Lnew:=map(f,Lold): while Lnew <> Lold do Lold:=Lnew: Lnew:=map(f,Lold): od: Cycles:=IOrbits(f,[op(Lnew)]): #----------------------- components G:=networks[graph](Graph(f,L)): compG:=components(G): #----------------------- basins nCycles:=nops(Cycles): Os:=array(1..nCycles): for i to nCycles do cycle:=Cycles[i]: select((x,c)->member(c[1],x),compG,cycle): basin:=op(%) minus {op(cycle)}: Os[i]:=[basin,cycle] od: [seq(Os[i],i=1..nCycles)] end: # IOrbits_________________________________________________________AlgDyn # global orbit structure of an invertible map f over L # >IOrbits(f,L); # OUTPUT: a list of lists [[cycle_1],...,[cycle_n]] IOrbits:=proc(f::procedure,L::list) local ic, norbits, L_now, orbit_now, T: global IOrbit: T:=table(): norbits:=0: L_now:={op(L)}: while nops(L_now) > 0 do ic:=op(1,L_now): orbit_now:=IOrbit(f,ic): L_now:=L_now minus {op(orbit_now)}: norbits:=norbits+1: T[norbits]:=orbit_now: od: map(x->op(x),[entries(T)]): end: # Graph____________________________________________________________Integ # graph of phase space # >Graph(f,L); # OUTPUT: vset,estet, where vset and esets are sets representing the # vertices and oriented edges of the graph, repectively. Graph:=proc(f::procedure,L::list) local vset, eset: {seq([z,f(z)],z=L)}: eset:=remove(x->evalb(nops(x)=1),%): vset:=map(op,eset): vset,eset end: # PhiPoly_________________________________________________________AlgDyn # Phi-polynomial # >PhiPoly(f, n, x); # INPUT: f:polynom - polynomial # n:integer - the period # x:symbol - an indeterminate # OUTPUT: a polynomial. # FUNCTIONS: none # SYNOPSIS: # Polynomial whose roots are the periodic points of f(x) # of essential period n. alias(phipoly=PhiPoly): PhiPoly := proc(f::polynom, n::posint, x::symbol) local d, fd, result, start: fd := x: result := 1: start:= 1: for d in sort([op(divisors(n))]) do from start to d do fd := subs(x=f, fd) od: result := result * (fd - x)^mobius(n/d): start := d + 1 od: normal(result) end: # MatMult_________________________________________________________AlgDyn # MatMult - matrix of multiplication of g(x) by f(x) in # Z[x]/(g(x)), w.r.t basis 1,x,...,x^(n-1). # LocMatMult - matrix of multiplication by a polynomial (char. p) # > MatMult(f, g, x); # > LocMatMult(f, g, p, x); # INPUT: f,g:polynom - polynomials # x:symbol - indeterminate # p:prime - the characteristic # OUTPUT: a matrix # FUNCTIONS: none # SYNOPSIS: # Matrix of multiplication of g(x) by f(x) in # Z[x]/(g(x)) (or Fp[x]/(g(x)), w.r.t the basis 1,x,x^2,....) alias(matmult=MatMult): MatMult := proc(f::polynom, g::polynom, x::symbol) local i, j, remnow, dg, A: dg := degree(g, x): #----------------------------------------- error control if dg = 0 then ERROR(`the second arguments cannot be a constant polynomial`) fi: A := array(1..dg, 1..dg): #----------- multiply f(x) by basis, then take remainder modulo g(x) for i to dg do remnow := rem(f * x^(i-1), g, x); for j to dg do A[j,i] := coeff(remnow, x, j-1) od od: evalm(A) end: alias(locmatmult=LocMatMult): LocMatMult := proc(f::polynom, g::polynom, p::prime, x::symbol) local i, j, remnow, dg, A: dg := degree(g mod p, x): #----------------------------------------- error control if (dg = 0) then ERROR(`the second arguments cannot be a constant polynomial`) fi: A := array(1..dg, 1..dg): #----------- multiply f(x) by basis, then take remainder modulo g(x) for i to dg do remnow := rem(f * x^(i-1), g, x) mod p; for j to dg do A[j,i] := coeff(remnow, x, j-1): od: od: evalm(A) end: # Wop_____________________________________________________________AlgDyn # local wedge operator of two polynomials # Wop - local wedge operator of two polynomials # >Wop(f, g, p, x); # >Wop2(i, j, p, n); # INPUT: f,g:polynom - sequence of polynomials # i,j:integer - indices of polynomials in lexic. order # x:symbol - an indeterminate # p:prime - the characteristic # n:integer - the degree # OUTPUT: a polynomial. # FUNCTIONS: LocMatMult # SYNOPSIS: # Compute local wedge operator: g(x) > f(x) (simple roots) alias(wop=Wop): Wop := proc(f::polynom, g::polynom, p::prime, x::symbol) global LocMatMult: local poly, diffpoly: poly := linalg[charpoly](LocMatMult(args),x): poly := Factor(poly) mod p: diffpoly := diff(poly, x): normal(poly/gcd(poly,diffpoly)) mod p end: alias(wop2=Wop2): Wop2 := proc(f::posint, g::posint, p::prime, n::posint) local x, fpoly, gpoly; fpoly:=Int2Poly(f,p,n,x): gpoly:=Int2Poly(g,p,n,x): Poly2Int(Wop(fpoly,gpoly,p,x),p,x) end: # GWop,GWopBar____________________________________________________AlgDyn # wedge operator for simple roots, and multiple roots # >GWop(f, g, x); # >GWopBar(f, g, x); # INPUT: f,g:polynom - polynomials # x:symbol - an indeterminate # OUTPUT: a polynomial in the indeterminate x # FUNCTIONS: MatMult # SYNOPSIS: # Compute value of wedge operator: g(x) |> f(x) alias(gwop=GWop): GWop := proc(f::polynom, g::polynom, x::symbol) global MatMult: local poly: poly := linalg[charpoly](MatMult(f,g,x),x): normal(poly/gcd(poly,diff(poly,x))) end: alias(gwopbar=GWopBar): GWopBar := proc(f::polynom, g::polynom, x::symbol) global MatMult: linalg[charpoly](MatMult(f,g,x),x) end: