My goal is to implement the algebraic closure of Q along the lines of
\cite{Steel:AlgebraicallyClosedFields:2010} (bibtex item below).
Such an implementation seemingly needs the algebraic closure of a finite
field. This is what I present in this mail.
The main idea is that not a true algebraic closure is implemented, but a
domain that is close to it in the following sense.
At each point in time the domain represents a finite extension of the
ground field (that's currently GF(p) in my attached files).
So the factorization works exactly as in a finite field extension.
However, if one uses rootOf(somepoly), then the polynomial is internally
factored and if there is no linear factor, then one of the irreducible
factors is used to extend the current field.
Technically, this means that a new variable is created and the
irreducible polynomial (considered in this new variable) is added to an
internal triangular set. The new variable is then returned as a new root.
So in fact, this domain DynamicAlgebraicClosurePrimeField behaves
like FiniteFieldExtensionByPolynomial, just that instead of one
polynomial, I use a triangular set.
As can be seen from the code, I actually only build on the functions of
PrimeField and some packages, in particular, the factorization in
ddfact.spad, but not on any FiniteFieldExtension* domains.
Luckily, most things I needed were already implemented as category defaults.
However, there is one important change that would have to be done if
DynamicAlgebraicClosurePrimeField should make it into FriCAS. It
concerns ModMonic. Seemingly for efficiency reasons, ModMonic checks in
"setpoly" for certain conditions and then avoids some recomputations.
It took me some time to realize that ModMonic behaves badly for my
dynamic domain. ModMonic(F, SUP F) silently assumes that for finite
fields size()$F is a constant. Well, a reasonable assumption, but bad
for me, since DynamicAlgebraicClosurePrimeField grows in size with each
new root that is not yet in the field. So I had to patch modmon.spad in
a few places. Since these changes are not too costly, I think they can
also be used for the other cases where size()$F is a true constant.
Attached is the code and a small session to demonstrate the things that
I want to work. Maybe there are still some things that do not work
properly, so just take the code as a "work-in-progress".
Comments and bug reports are welcome.
Ralf
@article{Steel:AlgebraicallyClosedFields:2010,
title = {Computing with algebraically closed fields},
journal = JSC,
volume = 45,
number = 3,
pages = {342 - 372},
year = 2010,
issn = {0747-7171},
doi = {10.1016/j.jsc.2009.09.005},
url =
{http://www.sciencedirect.com/science/article/pii/S0747717109001497},
author = {Allan K. Steel},
keywords = {Algebraic closure, Algebraic number field, Algebraic
function field, Field extension, Inseparability,
Non-perfect field, Polynomial factorization, Root
finding},
}
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/fricas-devel/173d86e0-acb5-3b26-d073-476c9831724c%40hemmecke.org.
-------------------------------------------------------------------
---
--- DynamicAlgebraicClosurePrimeField
--- Copyright (C) 2019, Ralf Hemmecke <[email protected]>
---
-------------------------------------------------------------------
-- This program is free software: you can redistribute it and/or modify
-- it under the terms of the GNU General Public License as published by
-- the Free Software Foundation, either version 3 of the License, or
-- (at your option) any later version.
--
-- This program is distributed in the hope that it will be useful,
-- but WITHOUT ANY WARRANTY; without even the implied warranty of
-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-- GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with this program. If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------
)if LiterateDoc
\documentclass{article}
\usepackage{url}
\DeclareUrlCommand\code{\urlstyle{tt}}
\begin{document}
\title{Dynamic Algebraic Closure Prime Field}
\author{Ralf Hemmecke}
\date{22-Oct-2019}
\maketitle
\begin{abstract}
\begin{verbatim}
@article{Steel:AlgebraicallyClosedFields:2010,
title = {Computing with algebraically closed fields},
journal = JSC,
volume = 45,
number = 3,
pages = {342 - 372},
year = 2010,
issn = {0747-7171},
doi = {10.1016/j.jsc.2009.09.005},
url =
{http://www.sciencedirect.com/science/article/pii/S0747717109001497},
author = {Allan K. Steel},
keywords = {Algebraic closure, Algebraic number field, Algebraic
function field, Field extension, Inseparability,
Non-perfect field, Polynomial factorization, Root
finding},
}
\end{verbatim}
\end{abstract}
\tableofcontents
\section{Overview}
For the implementation of an ACF in the sense of
\cite{Steel:AlgebraicallyClosedFields:2010} we need finite fields and
algebraic extensions of them. Since in FriCAS we the type of the
finite field must be known at compile time, we create a finite field
that behaves like the algebraic closure of $GF(p)$.
Maybe, eventually, we must create a domain that behaves like the union
of such fields for all primes $p$.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
)endif
OF==>OutputForm
dbgPrint(x,y) ==> print(([":> "::Symbol::OF, x::Symbol::OF,
y::OF]$List(OF)::OF))
rep x ==> (x@%) pretend Rep
per x ==> (x@Rep) pretend %
P ==> PositiveInteger
N ==> NonNegativeInteger
Z ==> Integer
Q ==> Fraction Z
SI ==> SingleInteger
SY ==> Symbol
asN x ==> x pretend N
asP x ==> x pretend P
asZ x ==> x pretend Z
)if LiterateDoc
%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{IndexedVariable}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
)endif
)abbrev domain IVAR IndexedVariable
++ Author: Ralf Hemmecke
++ Description:
++ IndexedVariable(x) represents the variables x1, x2, x3. etc.
++ and orders them like the natural order of their indices.
IndexedVariable(v: String): Exports == Implementation where
Exports ==> Join(OrderedSet, ConvertibleTo Symbol) with
new: () -> %
++ new() returns a new (not yet used) variable.
rank: () -> Z
++ rank() returns the highest index of the currently known
++ variables. It returns 0 if no variables are known.
lookup: % -> P
++ lookup(x) returns the index of the variable.
index: P -> %
++ index(p) returns the variable v such that lookup(v) = p. It
++ is an error to call index(p) for p > rank().
Implementation ==> P add
next: Reference Z := ref(1$Z)
syms: FlexibleArray Symbol := empty()
Rep ==> Z -- we only use positive integers (x>0)
rank(): Z ==
-- free next
deref(next) - 1
lookup(x: %): P == asP rep x
index(p: P): % ==
z: Z := asZ p
z < deref next => per z
error "No variable with that index is known."
new(): % ==
-- free next, syms
m: Z := deref next
s: String := concat(v, string m)
sym: Symbol := s::Symbol
concat!(syms, sym) -- destructively changes syms
setelt!(next, m+1)$Reference(Z)
per m
convert(x: %): Symbol == syms(rep x)
coerce(x: %): OutputForm ==
import from Symbol
(convert(x)@Symbol)::OutputForm
)if LiterateDoc
%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Dynamic Algebraic Closure Prime Field}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The constructor \code{DynamicAlgebraicClosurePrimeField} takes as
input a prime $p$ and returns the algebraic closure of
$GF(p)$ in the sense of \cite{Steel:AlgebraicallyClosedFields:2010}.
)endif
)abbrev domain DACPF DynamicAlgebraicClosurePrimeField
++ Author: Ralf Hemmecke
++ Description:
++ Implement an ACF according to \cite{Steel:AlgebraicallyClosedFields:2010}
++ but here for a prime field of characteristic p.
DynamicAlgebraicClosurePrimeField(p: P): Exports == Implementation where
V ==> IndexedVariable("f")
K ==> PrimeField p
EX ==> IndexedExponents V
R ==> NewSparseMultivariatePolynomial(K, V)
TS ==> RegularTriangularSet(K, EX, V, R)
NUP ==> NewSparseUnivariatePolynomial
SUP ==> SparseUnivariatePolynomial
Exports ==> Join(AlgebraicallyClosedField, FiniteFieldCategory) with
rank: () -> Z
show: () -> Void()
extendBy!: SUP % -> Void()
getGenerator: P -> %
Implementation ==> R add
-- We assume that all elements are always reduced wrt. the
-- triangular set ts.
Rep ==> R
import from Rep
xx ==> rep x
yy ==> rep y
RZ ==> SUP R
G1 ==> Record(gcd: NUP %, coef1: NUP %) -- for "inv" computation
FUNION ==> Union("nil", "sqfr", "irred", "prime")
FREC ==> Record(flag: FUNION, factor: SUP %, exponent: N)
FL ==> List FREC
univar(r, v) ==> (univariate(r, v)$R) pretend NUP % -- (%, V) -> NUP %
multivar(s, v) ==> per(multivariate(s pretend RZ, v)$R) -- (NUP %, V) -> %
-- internally used vars, sorted decreasingly f4,f3,f2,f1
mvars: Reference List V := ref(empty()$List(V))
-- assign(sy: Symbol, x: %): Void() == (assignSymbol(sy, x, %)$Lisp; x)
newVariable(): V ==
import from List V
v: V := new()
setelt!(mvars, cons(v, deref mvars))$Reference(List V)
v
-- In contrast to \cite{Steel:AlgebraicallyClosedFields:2010} we
-- only allow irreducible polynomial to extend. In other words,
-- if ts belongs to the current representation of the extension field
-- and f is the polynomial with the highest variable, then f is
-- irreducible over the field that belongs to ts without f.
-- We assume that we can construct all roots of f in the current
-- field belonging to ts. Not that the highest variable is such a root.
rank(): Z ==
import from List V
empty? deref mvars => 0
asZ(lookup(first deref mvars)$V)
-- Domain variable ts represents the ideal I (\subseteq J)
-- Note that the first entry in ts it the one with the highest variable.
ts: Reference TS := ref(empty()$TS)
getRelation(n: P): R ==
--assert(not empty? deref gamma)
import from TS
triset: List R := deref(ts)::List(R)
m: Z := n pretend Z
triset(rank() - m + 1)
getGenerator(n: P): % ==
import from N
v: V := index(n)$V
multivar(monomial(1, 1)$RZ, v)
extendBy!(f: SUP %): Void() ==
v: V := newVariable()
minpoly: R := rep multivar(f, v)
triset: TS := deref ts
triset := extend(triset, minpoly)
setelt!(ts, triset)
--assert(# deref ts = # deref gamma)
--assert(# deref ts = # deref mvars)
--assert((vs := deref mvars; if not empty? vs then #vs = lookup first vs))
show(): Void() ==
import from V
triset: TS := deref ts
print(triset::OutputForm)$OF
vars: List V := deref mvars
l: List(OF) := [v::OutputForm for v in vars]
print(bracket(l)$OF)$OF
reduce(r: Rep): % ==
z: Record(rnum: K, polnum: R, den: K) := remainder(r, deref ts)$TS
per((z.rnum/z.den) * z.polnum)
_-(x: %): % == reduce(- xx)
((x: %) + (y: %)): % == reduce(xx + yy)
((x: %) - (y: %)): % == reduce(xx - yy)
((x: %) * (y: %)): % == reduce(xx * yy)
((x: %) ^ (n: P)): % == expt(x, n)$RepeatedSquaring(%)
((x: %) ^ (n: N)): % ==
zero? n => 1
x^asP(n)
((x: %) ^ (n: Z)): % ==
n < 0 => inv(x) ^ asP(-n)
x^asN(n)
zero?(x: %): Boolean == ground? xx and zero?(ground xx)$K
((x: %) = (y: %)): Boolean == zero?(x-y)
((x: %) / (y: %)): % == x * inv y
inv(x: %): % ==
zero? x => error "cannot invert zero element"
--assert(not zero? x)
ground? xx => per((inv(ground xx)$K)::R)
v: V := mvar xx
f: NUP % := univar(xx, v)
n: P := lookup v
gn: R := getRelation n
g: NUP % := univar(gn, v)
cu: G1 := halfExtendedSubResultantGcd1(f, g)
c: NUP % := cu.gcd
import from N
not zero? degree c => error "ffalgclos: inv: gcd is not a constant"
lcc: % := leadingCoefficient c
inv(lcc) * multivar(cu.coef1, v)
rootOf(poly: SUP %, sy: Symbol): % ==
-- We currently ignore the second parameter.
fl: FL := factorList(factor poly)$Factored(SUP %)
-- Note that fl cannot be empty.
for frec in fl repeat
f: SUP % := frec.factor
if one? degree f then return (- coefficient(f, 0))
-- Arriving here means that there was no linear factor, i.e.,
-- we must extend our field by an irreducible polynomial.
-- For simplicity, we simply take the first element in fl.
extendBy!(first(fl).factor)
-- After extension, there must be a root, namely the value
-- corresponding to the varlable v.
v: V := first(deref mvars)$List(V)
multivar(monomial(1, 1)$RZ, v)
((x: %) ^ (q: Q)): % ==
n: Z := numer q
d: N := asN denom q
f: SUP % := monomial(1, d)$SUP(%) - (x^n)::SUP(%)
rootOf f
sqrt(x: %): % ==
import from N
f: SUP % := monomial(1, 2)$SUP(%) - x::SUP(%)
rootOf f
representationType(): Union("prime","polynomial","normal","cyclic") ==
-- representationType() returns the type of the representation,
-- one of: prime, polynomial, normal, or cyclic.
error "NOT YET: representationType"
primitiveElement(): % == error "NOT YET: primitiveElement"
tableForDiscreteLogarithm(n: Z): Table(P, N) ==
error "NOT YET: tableForDiscreteLogarithm"
factorsOfCyclicGroupSize(): List Record(factor: Z, exponent: N) ==
factors(factor(asZ(size())-1)$Z)$Factored(Z)
levelSizes(triset: List R): List Z == -- local function
empty? triset => [asZ(size()$K)]$List(Z)
d: N := mdeg first triset
triset := rest triset
ls: List Z := levelSizes triset
sz: Z := first ls
cons(sz^d, ls)
allLevelSizes(): List Z == -- local function
import from TS
import from R
triset: List R := deref(ts)::List(R)
levelSizes triset
size(): N == asN(first(allLevelSizes())$List(Z))
indexRec(n: Z, ls: List Z): % ==
--assert(not empty? ls)
--assert(n < first ls)
--assert(n >= 0)
zero? n => 0
empty? ls => per((index(asP n)$K)::R) -- Base case; element of K
b: Z := first ls
n < b => indexRec(n, rest ls)
k: P := asP(#ls)
v: V := index(k)$V
-- Now find the "digits" of n in base b
s: NUP % := 0
j: N := 0
while n > 0 repeat
rec: Record(quotient: Z, remainder: Z) := divide(n, b)
r: Z := rec.remainder
if not zero? r then
c: % := indexRec(r, ls)
s := s + monomial(c, j)$NUP(%)
n := rec.quotient
j := j + 1
multivar(s, v)
index(n: P): % ==
-- Important: Even after extending the field, i.e. growing ts,
-- index(n) should give the same element as before.
-- The basic ideas is that we consider the constant of the
-- polynomial first.
-- The following algorithm works recursively.
ls: List Z := allLevelSizes()
nn: Z := first ls
-- Note that index expects n in the range 1..size(), but
-- indexRec in the range 0..size()-1.
indexRec(positiveRemainder(asZ n - 1, nn), rest ls)
lookupRec(r: R, triset: List R, ls: List Z): Z == -- local function
zero? r => 0
ground? r => asZ(lookup(ground r)$K)
-- assert(#triset + 1 = #ls
-- lookupRec returns values in the range 0..size()-1.
v: V := mvar r
while v < mvar first triset repeat
triset := rest triset
ls := rest ls
triset := rest triset
ls := rest ls
sz: Z := first ls
result: Z := lookupRec(leadingCoefficient(r, v), triset, ls)
d: Z := asZ mdeg r
while not zero?(r := reductum r) repeat
e: N := asN(d - (d := asZ mdeg r))
c: Z := lookupRec(leadingCoefficient(r, v), triset, ls)
result := result * sz^e + c
result * sz^asN(d)
lookup(x: %): P ==
--assert(x = reduce x)
-- lookup returns values in the range 1..size().
import from TS
import from R
triset: List R := deref(ts)::List(R)
ls: List Z := levelSizes triset
asP(1 + lookupRec(xx, triset, ls))
--: *: (%, Fraction Integer) -> %
--: *: (Fraction Integer, %) -> %
--: *: (Integer, %) -> %
--: *: (NonNegativeInteger, %) -> %
--: *: (PositiveInteger, %) -> %
--: ^: (%, Integer) -> %
--: ^: (%, NonNegativeInteger) -> %
--: ^: (%, PositiveInteger) -> %
--: ~=: (%, %) -> Boolean
--: annihilate?: (%, %) -> Boolean
--: antiCommutator: (%, %) -> %
--: associates?: (%, %) -> Boolean
--: associator: (%, %, %) -> %
--: characteristic: () -> NonNegativeInteger
--: charthRoot: % -> %
--: -- charthRoot(a) takes the characteristic’th root of a. Note: such
--: -- a root is alway defined in finite fields.
--: charthRoot: % -> Union(%, failed)
--: coerce: % -> %
--: coerce: Fraction Integer -> %
--: coerce: Integer -> %
--: commutator: (%, %) -> %
--: conditionP: Matrix % -> Union(Vector %, failed)
--: convert: % -> InputForm
--: createPrimitiveElement: () -> %
--: -- createPrimitiveElement() computes a generator of the (cyclic)
--: -- multiplicative group of the field.
--: D: % -> %
--: D: (%, NonNegativeInteger) -> %
--: differentiate: % -> %
--: differentiate: (%, NonNegativeInteger) -> %
--: discreteLog: % -> NonNegativeInteger
--: -- discreteLog(a) computes the discrete logarithm of a with
--: -- respect to primitiveElement() of the field.
--: discreteLog: (%, %) -> Union(NonNegativeInteger, failed)
--: divide: (%, %) -> Record(quotient: %, remainder: %)
--: enumerate: () -> List %
--: euclideanSize: % -> NonNegativeInteger
--: expressIdealMember: (List %, %) -> Union(List %, failed)
--: exquo: (%, %) -> Union(%, failed)
--: extendedEuclidean: (%, %) -> Record(coef1: %, coef2: %, generator: %)
--: extendedEuclidean: (%, %, %) -> Union(Record(coef1: %, coef2: %),
failed)
--: factor: % -> Factored %
--: factorPolynomial: SparseUnivariatePolynomial % -> Factored
SparseUnivariatePolynomial %
--: factorsOfCyclicGroupSize: () -> List Record(factor: Integer, exponent:
NonNegativeInteger)
--: -- factorsOfCyclicGroupSize() returns the factorization of
--: -- size()-1
--: factorSquareFreePolynomial: SparseUnivariatePolynomial % -> Factored
SparseUnivariatePolynomial %
--: gcd: (%, %) -> %
--: gcd: List % -> %
--: gcdPolynomial: (SparseUnivariatePolynomial %,
SparseUnivariatePolynomial %) -> SparseUnivariatePolynomial %
--: hash: % -> SingleInteger
--: hashUpdate!: (HashState, %) -> HashState
--: index: PositiveInteger -> %
--: init: %
--: from StepThrough
--: inv: % -> %
--: latex: % -> String
--: lcm: (%, %) -> %
--: lcm: List % -> %
--: lcmCoef: (%, %) -> Record(llcm_res: %, coeff1: %, coeff2: %)
--: leftPower: (%, NonNegativeInteger) -> %
--: leftPower: (%, PositiveInteger) -> %
--: leftRecip: % -> Union(%, failed)
--: lookup: % -> PositiveInteger
--: multiEuclidean: (List %, %) -> Union(List %, failed)
--: nextItem: % -> Union(%, failed)
--: one?: % -> Boolean
--: opposite?: (%, %) -> Boolean
--: order: % -> OnePointCompletion PositiveInteger
--: order: % -> PositiveInteger
--: -- order(b) computes the order of an element b in the
--: -- multiplicative group of the field. Error: if b equals 0.
--: prime?: % -> Boolean
--: primeFrobenius: % -> %
--: primeFrobenius: (%, NonNegativeInteger) -> %
--: primitive?: % -> Boolean
--: -- primitive?(b) tests whether the element b is a generator of the
--: -- (cyclic) multiplicative group of the field, i.e. is a primitive
--: -- element. Implementation Note: see ch.IX.1.3, th.2 in D. Lipson.
--: primitiveElement: () -> %
--: -- primitiveElement() returns a primitive element stored in a
--: -- global variable in the domain. At first call, the primitive
--: -- element is computed by calling createPrimitiveElement.
--: principalIdeal: List % -> Record(coef: List %, generator: %)
--: quo: (%, %) -> %
--: random: () -> %
--: recip: % -> Union(%, failed)
--: rem: (%, %) -> %
--: representationType: () -> Union(prime, polynomial, normal, cyclic)
--: -- representationType() returns the type of the representation,
--: -- one of: prime, polynomial, normal, or cyclic.
--: rightPower: (%, NonNegativeInteger) -> %
--: rightPower: (%, PositiveInteger) -> %
--: rightRecip: % -> Union(%, failed)
--: sample: %
--: size: () -> NonNegativeInteger
--: sizeLess?: (%, %) -> Boolean
--: smaller?: (%, %) -> Boolean
--: solveLinearPolynomialEquation: (List SparseUnivariatePolynomial %,
SparseUnivariatePolynomial %) -> Union(List SparseUnivariatePolynomial %,
failed)
--: squareFree: % -> Factored %
--: squareFreePart: % -> %
--: squareFreePolynomial: SparseUnivariatePolynomial % -> Factored
SparseUnivariatePolynomial %
--: subtractIfCan: (%, %) -> Union(%, failed)
--: tableForDiscreteLogarithm: Integer -> Table(PositiveInteger,
NonNegativeInteger)
--: -- tableForDiscreteLogarithm(a, n) returns a table of the discrete
--: -- logarithms of a^0 up to a^(n-1) which, called with key
--: -- lookup(a^i) returns i for i in 0..n-1. Error: if not called for
--: -- prime divisors of order of multiplicative group.
--: unit?: % -> Boolean
--: unitCanonical: % -> %
--: unitNormal: % -> Record(unit: %, canonical: %, associate: %)
)if LiterateDoc
\bibliography{qeta}
\end{document}
)endif
)abbrev domain MODMON ModMonic
++ Description:
++ This package \undocumented
-- following line prevents caching ModMonic
)bo PUSH('ModMonic, $mutableDomains)
ModMonic(R, Rep) : C == T
where
R : Ring
Rep : UnivariatePolynomialCategory(R)
C == UnivariatePolynomialCategory(R) with
--operations
setPoly : Rep -> Rep
++ setPoly(x) \undocumented
modulus : -> Rep
++ modulus() \undocumented
reduce : Rep -> %
++ reduce(x) \undocumented
lift : % -> Rep --reduce lift = identity
++ lift(x) \undocumented
coerce : Rep -> %
++ coerce(x) \undocumented
Vectorise : % -> Vector(R)
++ Vectorise(x) \undocumented
UnVectorise : Vector(R) -> %
++ UnVectorise(v) \undocumented
An : % -> Vector(R)
++ An(x) \undocumented
pow : -> PrimitiveArray(%)
++ pow() \undocumented
computePowers : -> PrimitiveArray(%)
++ computePowers() \undocumented
if R has FiniteFieldCategory then
frobenius : % -> %
++ frobenius(x) computes \spad{x^q} where \spad{q} is
++ the size of \spad{R}
--LinearTransf: (%, Vector(R)) -> SquareMatrix<deg> R
--assertions
if R has Finite then Finite
T == add
--constants
m : Rep := monomial(1, 1)$Rep --| degree(m) > 0 and LeadingCoef(m) = R$1
d := degree(m)$Rep
d1 := (d-1)::NonNegativeInteger
twod := 2*d1+1
frobenius? : Boolean := R has FiniteFieldCategory
sz: NonNegativeInteger := size()$R -- to figure out when size() changed
--VectorRep := DirectProduct(d: NonNegativeInteger, R)
--declarations
x, y : %
p : Rep
d, n : Integer
e, k1, k2 : NonNegativeInteger
c : R
--vect: Vector(R)
power : PrimitiveArray(%)
frobeniusPower : PrimitiveArray(%)
computeFrobeniusPowers : () -> PrimitiveArray(%)
--representations
--mutable m --take this out??
--define
power := new(0, 0)
frobeniusPower := new(0, 0)
setPoly (mon : Rep) ==
oldsz := sz
sz := size()$R
sz=oldsz and mon =$Rep m => mon
oldm := m
leadingCoefficient mon ~= 1 => error "polynomial must be monic"
-- following copy code needed since FFPOLY can modify mon
copymon : Rep := 0
while not zero? mon repeat
copymon := monomial(leadingCoefficient mon, degree mon)$Rep + copymon
mon := reductum mon
m := copymon
d := degree(m)$Rep
d1 := (d-1)::NonNegativeInteger
twod := 2*d1+1
power := computePowers()
if frobenius? then
if sz = oldsz then
degree(oldm)>1 and not((oldm exquo$Rep m) case "failed") =>
for i in 1..d1 repeat
frobeniusPower(i) := reduce lift frobeniusPower(i)
frobeniusPower := computeFrobeniusPowers()
m
modulus == m
if R has Finite then
size == d * size$R
random == UnVectorise([random()$R for i in 0..d1])
0 == 0$Rep
1 == 1$Rep
c * x == c *$Rep x
n * x == (n::R) *$Rep x
coerce(c : R) : % == monomial(c, 0)$Rep
coerce(x : %) : OutputForm == coerce(x)$Rep
coefficient(x, e) : R == coefficient(x, e)$Rep
reductum(x) == reductum(x)$Rep
leadingCoefficient x == (leadingCoefficient x)$Rep
degree x == (degree x)$Rep
lift(x) == x pretend Rep
reduce(p) == (monicDivide(p, m)$Rep).remainder
coerce(p) == reduce(p)
x = y == x =$Rep y
x + y == x +$Rep y
- x == -$Rep x
x * y ==
p := x *$Rep y
ans := 0$Rep
while (n := degree p)>d1 repeat
ans := ans + leadingCoefficient(p)*power.(n-d)
p := reductum p
ans+p
Vectorise(x) == [coefficient(lift(x), i) for i in 0..d1]
UnVectorise(vect) ==
reduce(+/[monomial(vect.(i+1), i) for i in 0..d1])
computePowers ==
mat : PrimitiveArray(%) := new(d, 0)
mat.0 := reductum(-m)$Rep
w : % := monomial$Rep (1, 1)
for i in 1..d1 repeat
mat.i := w *$Rep mat.(i-1)
if degree mat.i = d then
mat.i := reductum mat.i + leadingCoefficient mat.i * mat.0
mat
if frobenius? then
computeFrobeniusPowers() ==
mat : PrimitiveArray(%) := new(d, 1)
mat.1 := mult := monomial(1, size$R)$%
for i in 2..d1 repeat
mat.i := mult * mat.(i-1)
mat
frobenius(a : %) : % ==
aq : % := 0
while a ~= 0 repeat
aq := aq + leadingCoefficient(a)*frobeniusPower(degree a)
a := reductum a
aq
pow == power
monomial(c, e)==
if e<d then monomial(c, e)$Rep
else
if e<=twod then
c * power.(e-d)
else
k1 := e quo twod
k2 := (e-k1*twod)::NonNegativeInteger
reduce((power.d1 ^k1)*monomial(c, k2))
if R has Field then
(x:% exquo y:%):Union(%, "failed") ==
uv := extendedEuclidean(y, modulus(), x)$Rep
uv case "failed" => "failed"
return reduce(uv.coef1)
recip(y:%):Union(%, "failed") == 1 exquo y
divide(x : %, y : %) ==
(q := (x exquo y)) case "failed" => error "not divisible"
[q, 0]
-- An(MM) == Vectorise(-(reduce(reductum(m))::MM))
-- LinearTransf(vect, MM) ==
-- ans := 0::SquareMatrix<d>(R)
-- for i in 1..d do setelt!(ans, i, 1, vect.i)
-- for j in 2..d do
-- setelt!(ans, 1, j, elt(ans, d, j-1) * An(MM).1)
-- for i in 2..d do
-- setelt!(ans, i, j, elt(ans, i-1, j-1) + elt(ans, d, j-1) *
An(MM).i)
-- ans
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the
-- names of its contributors may be used to endorse or promote products
-- derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
)clear complete
)co src/algebra/modmon.spad
)co src/algebra/ffalgclos.spad
P ==> PositiveInteger
Z ==> Integer
V ==> IndexedVariable("f")
E ==> IndexedExponents V
R ==> NewSparseMultivariatePolynomial(K, V)
TS ==> RegularTriangularSet(Q, E, V, R)
SI ==> SingleInteger
maxi2 ==> shift(max()$SI, -1)::Z
p := qcoerce(prevPrime(maxi2)$IntegerPrimesPackage(Z))@P
--p := 53
K ==> PrimeField p
A ==> DynamicAlgebraicClosurePrimeField(p)
F ==> FiniteFieldExtensionByPolynomial(K, kx)
KZ ==> SUP K
AZ ==> SUP(A)
FZ ==> SUP(F)
az d ==> monomial(1, d)$AZ
kz d ==> monomial(1, d)$KZ
fz d ==> monomial(1, d)$FZ
kx := kz 2 + 18 * kz 1 + 6
ax := az 2 + 18 * az 1 + 6
fx := fz 2 + 18 * fz 1 + 6
rootOf(ax, 'v)
show()$A
-- factor fx -- use only when kx is irreducible over K
factor ax
ax2 := monomial(1,5)$AZ-3
factor ax2
irreducible? ax2
rts := rootsOf(ax2)
[x^5 for x in rts]
[ax2 x for x in rts]
irreducible? ax2
factor(ax2)
ax3 := monomial(1,3)$AZ-7
factor ax3
fl := factorList factor(ax3)
reduce(_*, [x.factor for x in fl])
a7: A := -7
show()$A
f1 := getGenerator(1)$A
ax3 := 9*f1 + 35
inv ax3
r := f1^(1/7)
r^7
r := f1^(1/9)
r^9
show()$A
(1) -> P ==> PositiveInteger
Type: Void
(2) -> Z ==> Integer
Type: Void
(3) -> V ==> IndexedVariable("f")
Type: Void
(4) -> E ==> IndexedExponents V
Type: Void
(5) -> R ==> NewSparseMultivariatePolynomial(K, V)
Type: Void
(6) -> TS ==> RegularTriangularSet(Q, E, V, R)
Type: Void
(7) -> SI ==> SingleInteger
Type: Void
(8) -> maxi2 ==> shift(max()$SI, -1)::Z
Type: Void
(9) -> p := qcoerce(prevPrime(maxi2)$IntegerPrimesPackage(Z))@P
(9) 2305843009213693921
Type: PositiveInteger
(10) -> --p := 53
(10) -> K ==> PrimeField p
Type: Void
(11) -> A ==> DynamicAlgebraicClosurePrimeField(p)
Type: Void
(12) -> F ==> FiniteFieldExtensionByPolynomial(K, kx)
Type: Void
(13) -> KZ ==> SUP K
Type: Void
(14) -> AZ ==> SUP(A)
Type: Void
(15) -> FZ ==> SUP(F)
Type: Void
(16) -> az d ==> monomial(1, d)$AZ
Type: Void
(17) -> kz d ==> monomial(1, d)$KZ
Type: Void
(18) -> fz d ==> monomial(1, d)$FZ
Type: Void
(19) -> kx := kz 2 + 18 * kz 1 + 6
2
(19) ? + 18 ? + 6
Type: SparseUnivariatePolynomial(PrimeField(2305843009213693921))
(20) -> ax := az 2 + 18 * az 1 + 6
2
(20) ? + 18 ? + 6
Type:
SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(21) -> fx := fz 2 + 18 * fz 1 + 6
2
(21) ? + 18 ? + 6
Type:
SparseUnivariatePolynomial(FiniteFieldExtensionByPolynomial(PrimeField(2305843009213693921),?^2+18*?+6))
(22) -> rootOf(ax, 'v)
(22) 2175765076733634914
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(23) -> show()$A
{}
[]
Type: Void
(24) -> -- factor fx -- use only when kx is irreducible over K
(24) -> factor ax
(24) (? + 130077932480059007)(? + 2175765076733634932)
Type:
Factored(SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921)))
(25) -> ax2 := monomial(1,5)$AZ-3
5
(25) ? + 2305843009213693918
Type:
SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(26) -> factor ax2
5
(26) ? + 2305843009213693918
Type:
Factored(SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921)))
(27) -> irreducible? ax2
(27) true
Type: Boolean
(28) -> rts := rootsOf(ax2)
(28)
[f1, 1916729304821096802 f1, 1074224509982630437 f1, 378638512013051338 f1,
1242093691610609264 f1]
Type: List(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(29) -> [x^5 for x in rts]
(29) [3, 3, 3, 3, 3]
Type: List(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(30) -> [ax2 x for x in rts]
(30) [0, 0, 0, 0, 0]
Type: List(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(31) -> irreducible? ax2
(31) false
Type: Boolean
(32) -> factor(ax2)
(32)
(? + 389113704392597119 f1)(? + 1063749317603084657 f1)
*
(? + 1231618499231063484 f1)(? + 1927204497200642583 f1)
*
(? + 2305843009213693920 f1)
Type:
Factored(SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921)))
(33) -> ax3 := monomial(1,3)$AZ-7
3
(33) ? + 2305843009213693914
Type:
SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(34) -> factor ax3
(34)
(? + 8387903196959538)(? + 431301322456251594)(? + 1866153783560482789)
Type:
Factored(SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921)))
(35) -> fl := factorList factor(ax3)
(35)
[[flag = "prime", factor = ? + 8387903196959538, exponent = 1],
[flag = "prime", factor = ? + 431301322456251594, exponent = 1],
[flag = "prime", factor = ? + 1866153783560482789, exponent = 1]]
Type: List(Record(flag: Union("nil","sqfr","irred","prime"),factor:
SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921)),exponent:
NonNegativeInteger))
(36) -> reduce(_*, [x.factor for x in fl])
3
(36) ? + 2305843009213693914
Type:
SparseUnivariatePolynomial(DynamicAlgebraicClosurePrimeField(2305843009213693921))
(37) -> a7: A := -7
(37) 2305843009213693914
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(38) -> show()$A
5
{f1 + 2305843009213693918}
[f1]
Type: Void
(39) -> f1 := getGenerator(1)$A
(39) f1
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(40) -> ax3 := 9*f1 + 35
(40) 9 f1 + 35
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(41) -> inv ax3
(41)
4 3 2
1906771754488978363 f1 + 14715095564764556 f1 + 967593743565335136 f1
+
848821460117751202 f1 + 1823123231128065150
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(42) -> r := f1^(1/7)
3
(42) 1208701501887878640 f1
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(43) -> r^7
(43) f1
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(44) -> r := f1^(1/9)
(44) f2
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(45) -> r^9
(45) f1
Type: DynamicAlgebraicClosurePrimeField(2305843009213693921)
(46) -> show()$A
5 9
{f1 + 2305843009213693918, f2 + 2305843009213693920 f1}
[f2, f1]
Type: Void
(47) ->
