open-axiom repository from github
\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra curve.spad}
\author{Manuel Bronstein}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category FFCAT FunctionFieldCategory}
<<category FFCAT FunctionFieldCategory>>=
)abbrev category FFCAT FunctionFieldCategory
++ Function field of a curve
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 19 Mai 1993
++ Description: This category is a model for the function field of a
++ plane algebraic curve.
++ Keywords: algebraic, curve, function, field.
FunctionFieldCategory(F, UP, UPUP): Category == Definition where
F : UniqueFactorizationDomain
UP : UnivariatePolynomialCategory F
UPUP: UnivariatePolynomialCategory Fraction UP
Z ==> Integer
Q ==> Fraction F
P ==> Polynomial F
RF ==> Fraction UP
QF ==> Fraction UPUP
SY ==> Symbol
REC ==> Record(num:$, den:UP, derivden:UP, gd:UP)
Definition ==> MonogenicAlgebra(RF, UPUP) with
numberOfComponents : () -> NonNegativeInteger
++ numberOfComponents() returns the number of absolutely irreducible
++ components.
genus : () -> NonNegativeInteger
++ genus() returns the genus of one absolutely irreducible component
absolutelyIrreducible? : () -> Boolean
++ absolutelyIrreducible?() tests if the curve absolutely irreducible?
rationalPoint? : (F, F) -> Boolean
++ rationalPoint?(a, b) tests if \spad{(x=a,y=b)} is on the curve.
branchPointAtInfinity? : () -> Boolean
++ branchPointAtInfinity?() tests if there is a branch point at infinity.
branchPoint? : F -> Boolean
++ branchPoint?(a) tests whether \spad{x = a} is a branch point.
branchPoint? : UP -> Boolean
++ branchPoint?(p) tests whether \spad{p(x) = 0} is a branch point.
singularAtInfinity? : () -> Boolean
++ singularAtInfinity?() tests if there is a singularity at infinity.
singular? : F -> Boolean
++ singular?(a) tests whether \spad{x = a} is singular.
singular? : UP -> Boolean
++ singular?(p) tests whether \spad{p(x) = 0} is singular.
ramifiedAtInfinity? : () -> Boolean
++ ramifiedAtInfinity?() tests if infinity is ramified.
ramified? : F -> Boolean
++ ramified?(a) tests whether \spad{x = a} is ramified.
ramified? : UP -> Boolean
++ ramified?(p) tests whether \spad{p(x) = 0} is ramified.
integralBasis : () -> Vector $
++ integralBasis() returns the integral basis for the curve.
integralBasisAtInfinity: () -> Vector $
++ integralBasisAtInfinity() returns the local integral basis at infinity.
integralAtInfinity? : $ -> Boolean
++ integralAtInfinity?() tests if f is locally integral at infinity.
integral? : $ -> Boolean
++ integral?() tests if f is integral over \spad{k[x]}.
complementaryBasis : Vector $ -> Vector $
++ complementaryBasis(b1,...,bn) returns the complementary basis
++ \spad{(b1',...,bn')} of \spad{(b1,...,bn)}.
normalizeAtInfinity : Vector $ -> Vector $
++ normalizeAtInfinity(v) makes v normal at infinity.
reduceBasisAtInfinity : Vector $ -> Vector $
++ reduceBasisAtInfinity(b1,...,bn) returns \spad{(x**i * bj)}
++ for all i,j such that \spad{x**i*bj} is locally integral at infinity.
integralMatrix : () -> Matrix RF
++ integralMatrix() returns M such that
++ \spad{(w1,...,wn) = M (1, y, ..., y**(n-1))},
++ where \spad{(w1,...,wn)} is the integral basis of
++ \spadfunFrom{integralBasis}{FunctionFieldCategory}.
inverseIntegralMatrix : () -> Matrix RF
++ inverseIntegralMatrix() returns M such that
++ \spad{M (w1,...,wn) = (1, y, ..., y**(n-1))}
++ where \spad{(w1,...,wn)} is the integral basis of
++ \spadfunFrom{integralBasis}{FunctionFieldCategory}.
integralMatrixAtInfinity : () -> Matrix RF
++ integralMatrixAtInfinity() returns M such that
++ \spad{(v1,...,vn) = M (1, y, ..., y**(n-1))}
++ where \spad{(v1,...,vn)} is the local integral basis at infinity
++ returned by \spad{infIntBasis()}.
inverseIntegralMatrixAtInfinity: () -> Matrix RF
++ inverseIntegralMatrixAtInfinity() returns M such
++ that \spad{M (v1,...,vn) = (1, y, ..., y**(n-1))}
++ where \spad{(v1,...,vn)} is the local integral basis at infinity
++ returned by \spad{infIntBasis()}.
yCoordinates : $ -> Record(num:Vector(UP), den:UP)
++ yCoordinates(f) returns \spad{[[A1,...,An], D]} such that
++ \spad{f = (A1 + A2 y +...+ An y**(n-1)) / D}.
represents : (Vector UP, UP) -> $
++ represents([A0,...,A(n-1)],D) returns
++ \spad{(A0 + A1 y +...+ A(n-1)*y**(n-1))/D}.
integralCoordinates : $ -> Record(num:Vector(UP), den:UP)
++ integralCoordinates(f) returns \spad{[[A1,...,An], D]} such that
++ \spad{f = (A1 w1 +...+ An wn) / D} where \spad{(w1,...,wn)} is the
++ integral basis returned by \spad{integralBasis()}.
integralRepresents : (Vector UP, UP) -> $
++ integralRepresents([A1,...,An], D) returns
++ \spad{(A1 w1+...+An wn)/D}
++ where \spad{(w1,...,wn)} is the integral
++ basis of \spad{integralBasis()}.
integralDerivationMatrix:(UP -> UP) -> Record(num:Matrix(UP),den:UP)
++ integralDerivationMatrix(d) extends the derivation d from UP to $
++ and returns (M, Q) such that the i^th row of M divided by Q form
++ the coordinates of \spad{d(wi)} with respect to \spad{(w1,...,wn)}
++ where \spad{(w1,...,wn)} is the integral basis returned
++ by integralBasis().
integral? : ($, F) -> Boolean
++ integral?(f, a) tests whether f is locally integral at \spad{x = a}.
integral? : ($, UP) -> Boolean
++ integral?(f, p) tests whether f is locally integral at \spad{p(x) = 0}.
differentiate : ($, UP -> UP) -> $
++ differentiate(x, d) extends the derivation d from UP to $ and
++ applies it to x.
primitivePart : $ -> $
++ primitivePart(f) removes the content of the denominator and
++ the common content of the numerator of f.
elt : ($, F, F) -> F
++ elt(f,a,b) or f(a, b) returns the value of f at the point \spad{(x = a, y = b)}
++ if it is not singular.
elliptic : () -> Union(UP, "failed")
++ elliptic() returns \spad{p(x)} if the curve is the elliptic
++ defined by \spad{y**2 = p(x)}, "failed" otherwise.
hyperelliptic : () -> Union(UP, "failed")
++ hyperelliptic() returns \spad{p(x)} if the curve is the hyperelliptic
++ defined by \spad{y**2 = p(x)}, "failed" otherwise.
algSplitSimple : ($, UP -> UP) -> REC
++ algSplitSimple(f, D) returns \spad{[h,d,d',g]} such that \spad{f=h/d},
++ \spad{h} is integral at all the normal places w.r.t. \spad{D},
++ \spad{d' = Dd}, \spad{g = gcd(d, discriminant())} and \spad{D}
++ is the derivation to use. \spad{f} must have at most simple finite
++ poles.
if F has Field then
nonSingularModel: SY -> List Polynomial F
++ nonSingularModel(u) returns the equations in u1,...,un of
++ an affine non-singular model for the curve.
if F has Finite then
rationalPoints: () -> List List F
++ rationalPoints() returns the list of all the affine rational points.
add
import InnerCommonDenominator(UP, RF, Vector UP, Vector RF)
import UnivariatePolynomialCommonDenominator(UP, RF, UPUP)
repOrder: (Matrix RF, Z) -> Z
Q2RF : Q -> RF
infOrder: RF -> Z
infValue: RF -> Fraction F
intvalue: (Vector UP, F, F) -> F
rfmonom : Z -> RF
kmin : (Matrix RF,Vector Q) -> Union(Record(pos:Z,km:Z),"failed")
Q2RF q == numer(q)::UP / denom(q)::UP
infOrder f == (degree denom f)::Z - (degree numer f)::Z
integral? f == ground?(integralCoordinates(f).den)
integral?(f:$, a:F) == (integralCoordinates(f).den)(a) ~= 0
absolutelyIrreducible? == one? numberOfComponents()
yCoordinates f == splitDenominator coordinates f
hyperelliptic() ==
degree(f := definingPolynomial()) ~= 2 => "failed"
(u:=retractIfCan(reductum f)@Union(RF,"failed")) case "failed" => "failed"
(v := retractIfCan(-(u::RF) / leadingCoefficient f)@Union(UP, "failed"))
case "failed" => "failed"
odd? degree(p := v::UP) => p
"failed"
algSplitSimple(f, derivation) ==
cd := splitDenominator lift f
dd := (cd.den exquo (g := gcd(cd.den, derivation(cd.den))))::UP
[reduce(inv(g::RF) * cd.num), dd, derivation dd,
gcd(dd, retract(discriminant())@UP)]
elliptic() ==
(u := hyperelliptic()) case "failed" => "failed"
degree(p := u::UP) = 3 => p
"failed"
rationalPoint?(x, y) ==
zero?((definingPolynomial() (y::UP::RF)) (x::UP::RF))
if F has Field then
import PolyGroebner(F)
import MatrixCommonDenominator(UP, RF)
UP2P : (UP, P) -> P
UPUP2P: (UPUP, P, P) -> P
UP2P(p, x) ==
(map(#1::P, p)$UnivariatePolynomialCategoryFunctions2(F, UP,
P, SparseUnivariatePolynomial P)) x
UPUP2P(p, x, y) ==
(map(UP2P(retract(#1)@UP, x),
p)$UnivariatePolynomialCategoryFunctions2(RF, UPUP,
P, SparseUnivariatePolynomial P)) y
nonSingularModel u ==
d := commonDenominator(coordinates(w := integralBasis()))::RF
n := #w
vars := [concat(string u, string i)::SY for i in 1..n]
x := '%%dummy1
y := '%%dummy2
select!(zero?(degree(#1, x)) and zero?(degree(#1, y)),
lexGroebner([v::P - UPUP2P(lift(d * w.i), x::P, y::P)
for v in vars for i in 1..n], concat([x, y], vars)))
if F has Finite then
ispoint: (UPUP, F, F) -> List F
-- must use the 'elt function explicitely or the compiler takes 45 mins
-- on that function MB 5/90
-- still takes ages : I split the expression up. JHD 6/Aug/90
ispoint(p, x, y) ==
jhd:RF:=p(y::UP::RF)
zero?(jhd (x::UP::RF)) => [x, y]
empty()
rationalPoints() ==
p := definingPolynomial()
concat [[pt for y in 1..size()$F | not empty?(pt :=
ispoint(p, index(x::PositiveInteger)$F,
index(y::PositiveInteger)$F))]$List(List F)
for x in 1..size()$F]$List(List(List F))
intvalue(v, x, y) ==
singular? x => error "Point is singular"
mini := minIndex(w := integralBasis())
rec := yCoordinates(+/[qelt(v, i)::RF * qelt(w, i)
for i in mini .. maxIndex w])
n := +/[(qelt(rec.num, i) x) *
(y ** ((i - mini)::NonNegativeInteger))
for i in mini .. maxIndex w]
zero?(d := (rec.den) x) =>
zero? n => error "0/0 -- cannot compute value yet"
error "Shouldn't happen"
(n exquo d)::F
elt(f, x, y) ==
rec := integralCoordinates f
n := intvalue(rec.num, x, y)
zero?(d := (rec.den) x) =>
zero? n => error "0/0 -- cannot compute value yet"
error "Function has a pole at the given point"
(n exquo d)::F
primitivePart f ==
cd := yCoordinates f
d := gcd([content qelt(cd.num, i)
for i in minIndex(cd.num) .. maxIndex(cd.num)]$List(F))
* primitivePart(cd.den)
represents [qelt(cd.num, i) / d
for i in minIndex(cd.num) .. maxIndex(cd.num)]$Vector(RF)
reduceBasisAtInfinity b ==
x := monomial(1, 1)$UP ::RF
concat([[f for j in 0.. while
integralAtInfinity?(f := x**j * qelt(b, i))]$Vector($)
for i in minIndex b .. maxIndex b]$List(Vector $))
complementaryBasis b ==
m := inverse(traceMatrix b)::Matrix(RF)
[represents row(m, i) for i in minRowIndex m .. maxRowIndex m]
integralAtInfinity? f ==
not any?(negative? infOrder(#1),
coordinates(f) * inverseIntegralMatrixAtInfinity())$Vector(RF)
numberOfComponents() ==
count(integralAtInfinity?, integralBasis())$Vector($)
represents(v:Vector UP, d:UP) ==
represents
[qelt(v, i) / d for i in minIndex v .. maxIndex v]$Vector(RF)
genus() ==
ds := discriminant()
d := degree(retract(ds)@UP) + infOrder(ds * determinant(
integralMatrixAtInfinity() * inverseIntegralMatrix()) ** 2)
dd := (((d exquo 2)::Z - rank()) exquo numberOfComponents())::Z
(dd + 1)::NonNegativeInteger
repOrder(m, i) ==
nostart:Boolean := true
ans:Z := 0
r := row(m, i)
for j in minIndex r .. maxIndex r | qelt(r, j) ~= 0 repeat
ans :=
nostart => (nostart := false; infOrder qelt(r, j))
min(ans, infOrder qelt(r,j))
nostart => error "Null row"
ans
infValue f ==
zero? f => 0
positive?(n := infOrder f) => 0
zero? n =>
(leadingCoefficient numer f) / (leadingCoefficient denom f)
error "f not locally integral at infinity"
rfmonom n ==
negative? n => inv(monomial(1, (-n)::NonNegativeInteger)$UP :: RF)
monomial(1, n::NonNegativeInteger)$UP :: RF
kmin(m, v) ==
nostart:Boolean := true
k:Z := 0
ii := minRowIndex m - (i0 := minIndex v)
for i in minIndex v .. maxIndex v | qelt(v, i) ~= 0 repeat
nk := repOrder(m, i + ii)
if nostart then (nostart := false; k := nk; i0 := i)
else
if nk < k then (k := nk; i0 := i)
nostart => "failed"
[i0, k]
normalizeAtInfinity w ==
ans := copy w
infm := inverseIntegralMatrixAtInfinity()
mhat := zero(rank(), rank())$Matrix(RF)
ii := minIndex w - minRowIndex mhat
repeat
m := coordinates(ans) * infm
r := [rfmonom repOrder(m, i)
for i in minRowIndex m .. maxRowIndex m]$Vector(RF)
for i in minRowIndex m .. maxRowIndex m repeat
for j in minColIndex m .. maxColIndex m repeat
qsetelt!(mhat, i, j, qelt(r, i + ii) * qelt(m, i, j))
sol := first nullSpace transpose map(infValue,
mhat)$MatrixCategoryFunctions2(RF, Vector RF, Vector RF,
Matrix RF, Q, Vector Q, Vector Q, Matrix Q)
(pr := kmin(m, sol)) case "failed" => return ans
qsetelt!(ans, pr.pos,
+/[Q2RF(qelt(sol, i)) * rfmonom(repOrder(m, i - ii) - pr.km)
* qelt(ans, i) for i in minIndex sol .. maxIndex sol])
integral?(f:$, p:UP) ==
(r:=retractIfCan(p)@Union(F,"failed")) case F => integral?(f,r::F)
(integralCoordinates(f).den exquo p) case "failed"
differentiate(f:$, d:UP -> UP) ==
differentiate(f, differentiate(#1, d)$RF)
@
\section{package MMAP MultipleMap}
<<package MMAP MultipleMap>>=
)abbrev package MMAP MultipleMap
++ Lifting a map through 2 levels of polynomials
++ Author: Manuel Bronstein
++ Date Created: May 1988
++ Date Last Updated: 11 Jul 1990
++ Description: Lifting of a map through 2 levels of polynomials;
MultipleMap(R1,UP1,UPUP1,R2,UP2,UPUP2): Exports == Implementation where
R1 : IntegralDomain
UP1 : UnivariatePolynomialCategory R1
UPUP1: UnivariatePolynomialCategory Fraction UP1
R2 : IntegralDomain
UP2 : UnivariatePolynomialCategory R2
UPUP2: UnivariatePolynomialCategory Fraction UP2
Q1 ==> Fraction UP1
Q2 ==> Fraction UP2
Exports ==> with
map: (R1 -> R2, UPUP1) -> UPUP2
++ map(f, p) lifts f to the domain of p then applies it to p.
Implementation ==> add
import UnivariatePolynomialCategoryFunctions2(R1, UP1, R2, UP2)
rfmap: (R1 -> R2, Q1) -> Q2
rfmap(f, q) == map(f, numer q) / map(f, denom q)
map(f, p) ==
map(rfmap(f, #1),
p)$UnivariatePolynomialCategoryFunctions2(Q1, UPUP1, Q2, UPUP2)
@
\section{package FFCAT2 FunctionFieldCategoryFunctions2}
<<package FFCAT2 FunctionFieldCategoryFunctions2>>=
)abbrev package FFCAT2 FunctionFieldCategoryFunctions2
++ Lifts a map from rings to function fields over them
++ Author: Manuel Bronstein
++ Date Created: May 1988
++ Date Last Updated: 26 Jul 1988
++ Description: Lifts a map from rings to function fields over them.
FunctionFieldCategoryFunctions2(R1, UP1, UPUP1, F1, R2, UP2, UPUP2, F2):
Exports == Implementation where
R1 : UniqueFactorizationDomain
UP1 : UnivariatePolynomialCategory R1
UPUP1: UnivariatePolynomialCategory Fraction UP1
F1 : FunctionFieldCategory(R1, UP1, UPUP1)
R2 : UniqueFactorizationDomain
UP2 : UnivariatePolynomialCategory R2
UPUP2: UnivariatePolynomialCategory Fraction UP2
F2 : FunctionFieldCategory(R2, UP2, UPUP2)
Exports ==> with
map: (R1 -> R2, F1) -> F2
++ map(f, p) lifts f to F1 and applies it to p.
Implementation ==> add
map(f, f1) ==
reduce(map(f, lift f1)$MultipleMap(R1, UP1, UPUP1, R2, UP2, UPUP2))
@
\section{package CHVAR ChangeOfVariable}
<<package CHVAR ChangeOfVariable>>=
)abbrev package CHVAR ChangeOfVariable
++ Sends a point to infinity
++ Author: Manuel Bronstein
++ Date Created: 1988
++ Date Last Updated: 22 Feb 1990
++ Description:
++ Tools to send a point to infinity on an algebraic curve.
ChangeOfVariable(F, UP, UPUP): Exports == Implementation where
F : UniqueFactorizationDomain
UP : UnivariatePolynomialCategory F
UPUP: UnivariatePolynomialCategory Fraction UP
N ==> NonNegativeInteger
Z ==> Integer
Q ==> Fraction Z
RF ==> Fraction UP
Exports ==> with
mkIntegral: UPUP -> Record(coef:RF, poly:UPUP)
++ mkIntegral(p(x,y)) returns \spad{[c(x), q(x,z)]} such that
++ \spad{z = c * y} is integral.
++ The algebraic relation between x and y is \spad{p(x, y) = 0}.
++ The algebraic relation between x and z is \spad{q(x, z) = 0}.
radPoly : UPUP -> Union(Record(radicand:RF, deg:N), "failed")
++ radPoly(p(x, y)) returns \spad{[c(x), n]} if p is of the form
++ \spad{y**n - c(x)}, "failed" otherwise.
rootPoly : (RF, N) -> Record(exponent: N, coef:RF, radicand:UP)
++ rootPoly(g, n) returns \spad{[m, c, P]} such that
++ \spad{c * g ** (1/n) = P ** (1/m)}
++ thus if \spad{y**n = g}, then \spad{z**m = P}
++ where \spad{z = c * y}.
goodPoint : (UPUP,UPUP) -> F
++ goodPoint(p, q) returns an integer a such that a is neither
++ a pole of \spad{p(x,y)} nor a branch point of \spad{q(x,y) = 0}.
eval : (UPUP, RF, RF) -> UPUP
++ eval(p(x,y), f(x), g(x)) returns \spad{p(f(x), y * g(x))}.
chvar : (UPUP,UPUP) -> Record(func:UPUP,poly:UPUP,c1:RF,c2:RF,deg:N)
++ chvar(f(x,y), p(x,y)) returns
++ \spad{[g(z,t), q(z,t), c1(z), c2(z), n]}
++ such that under the change of variable
++ \spad{x = c1(z)}, \spad{y = t * c2(z)},
++ one gets \spad{f(x,y) = g(z,t)}.
++ The algebraic relation between x and y is \spad{p(x, y) = 0}.
++ The algebraic relation between z and t is \spad{q(z, t) = 0}.
Implementation ==> add
import UnivariatePolynomialCommonDenominator(UP, RF, UPUP)
algPoly : UPUP -> Record(coef:RF, poly:UPUP)
RPrim : (UP, UP, UPUP) -> Record(coef:RF, poly:UPUP)
good? : (F, UP, UP) -> Boolean
infIntegral?: (UPUP, UPUP) -> Boolean
eval(p, x, y) == map(#1 x, p) monomial(y, 1)
good?(a, p, q) == p(a) ~= 0 and q(a) ~= 0
algPoly p ==
ground?(a:= retract(leadingCoefficient(q:=clearDenominator p))@UP)
=> RPrim(1, a, q)
c := d := squareFreePart a
q := clearDenominator q monomial(inv(d::RF), 1)
while not ground?(a := retract(leadingCoefficient q)@UP) repeat
c := c * (d := gcd(a, d))
q := clearDenominator q monomial(inv(d::RF), 1)
RPrim(c, a, q)
RPrim(c, a, q) ==
one? a => [c::RF, q]
[(a * c)::RF, clearDenominator q monomial(inv(a::RF), 1)]
-- always makes the algebraic integral, but does not send a point to infinity
-- if the integrand does not have a pole there (in the case of an nth-root)
chvar(f, modulus) ==
r1 := mkIntegral modulus
f1 := f monomial(r1inv := inv(r1.coef), 1)
infIntegral?(f1, r1.poly) =>
[f1, r1.poly, monomial(1,1)$UP :: RF,r1inv,degree(retract(r1.coef)@UP)]
x := (a:= goodPoint(f1,r1.poly))::UP::RF + inv(monomial(1,1)::RF)
r2c:= retract((r2 := mkIntegral map(#1 x, r1.poly)).coef)@UP
t := inv((monomial(1, 1)$UP - a::UP)::RF)
[- inv(monomial(1, 2)$UP :: RF) * eval(f1, x, inv(r2.coef)),
r2.poly, t, r1.coef * r2c t, degree r2c]
-- returns true if y is an n-th root, and it can be guaranteed that p(x,y)dx
-- is integral at infinity
-- expects y to be integral.
infIntegral?(p, modulus) ==
(r := radPoly modulus) case "failed" => false
ninv := inv(r.deg::Q)
degy:Q := degree(retract(r.radicand)@UP) * ninv
degp:Q := 0
while p ~= 0 repeat
c := leadingCoefficient p
degp := max(degp,
(2 + degree(numer c)::Z - degree(denom c)::Z)::Q + degree(p) * degy)
p := reductum p
degp <= ninv
mkIntegral p ==
(r := radPoly p) case "failed" => algPoly p
rp := rootPoly(r.radicand, r.deg)
[rp.coef, monomial(1, rp.exponent)$UPUP - rp.radicand::RF::UPUP]
goodPoint(p, modulus) ==
q :=
(r := radPoly modulus) case "failed" =>
retract(resultant(modulus, differentiate modulus))@UP
retract(r.radicand)@UP
d := commonDenominator p
for i in 0.. repeat
good?(a := i::F, q, d) => return a
good?(-a, q, d) => return -a
radPoly p ==
(r := retractIfCan(reductum p)@Union(RF, "failed")) case "failed"
=> "failed"
[- (r::RF), degree p]
-- we have y**m = g(x) = n(x)/d(x), so if we can write
-- (n(x) * d(x)**(m-1)) ** (1/m) = c(x) * P(x) ** (1/n)
-- then z**q = P(x) where z = (d(x) / c(x)) * y
rootPoly(g, m) ==
zero? g => error "Should not happen"
pr := nthRoot(squareFree((numer g) * (d := denom g) ** (m-1)::N),
m)$FactoredFunctions(UP)
[pr.exponent, d / pr.coef, */(pr.radicand)]
@
\section{domain RADFF RadicalFunctionField}
<<domain RADFF RadicalFunctionField>>=
)abbrev domain RADFF RadicalFunctionField
++ Function field defined by y**n = f(x)
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 27 July 1993
++ Keywords: algebraic, curve, radical, function, field.
++ Description: Function field defined by y**n = f(x);
++ Examples: )r RADFF INPUT
RadicalFunctionField(F, UP, UPUP, radicnd, n): Exports == Impl where
F : UniqueFactorizationDomain
UP : UnivariatePolynomialCategory F
UPUP : UnivariatePolynomialCategory Fraction UP
radicnd : Fraction UP
n : NonNegativeInteger
N ==> NonNegativeInteger
Z ==> Integer
RF ==> Fraction UP
QF ==> Fraction UPUP
UP2 ==> SparseUnivariatePolynomial UP
REC ==> Record(factor:UP, exponent:Z)
MOD ==> monomial(1, n)$UPUP - radicnd::UPUP
INIT ==> if (deref brandNew?) then startUp false
Exports ==> FunctionFieldCategory(F, UP, UPUP)
Impl ==> SimpleAlgebraicExtension(RF, UPUP, MOD) add
import ChangeOfVariable(F, UP, UPUP)
import InnerCommonDenominator(UP, RF, Vector UP, Vector RF)
import UnivariatePolynomialCategoryFunctions2(RF, UPUP, UP, UP2)
diag : Vector RF -> Vector $
startUp : Boolean -> Void
fullVector : (Factored UP, N) -> PrimitiveArray UP
iBasis : (UP, N) -> Vector UP
inftyBasis : (RF, N) -> Vector RF
basisvec : () -> Vector RF
char0StartUp: () -> Void
charPStartUp: () -> Void
getInfBasis : () -> Void
radcand : () -> UP
charPintbas : (UPUP, RF, Vector RF, Vector RF) -> Void
brandNew?:Reference(Boolean) := ref true
discPoly:Reference(RF) := ref(0$RF)
newrad:Reference(UP) := ref(0$UP)
n1 := (n - 1)::N
modulus := MOD
ibasis:Vector(RF) := new(n, 0)
invibasis:Vector(RF) := new(n, 0)
infbasis:Vector(RF) := new(n, 0)
invinfbasis:Vector(RF):= new(n, 0)
mini := minIndex ibasis
discriminant() == (INIT; deref discPoly)
radcand() == (INIT; deref newrad)
integralBasis() == (INIT; diag ibasis)
integralBasisAtInfinity() == (INIT; diag infbasis)
basisvec() == (INIT; ibasis)
integralMatrix() == diagonalMatrix basisvec()
integralMatrixAtInfinity() == (INIT; diagonalMatrix infbasis)
inverseIntegralMatrix() == (INIT; diagonalMatrix invibasis)
inverseIntegralMatrixAtInfinity()==(INIT;diagonalMatrix invinfbasis)
definingPolynomial() == modulus
ramified?(point:F) == zero?(radcand() point)
branchPointAtInfinity?() == (degree(radcand()) exquo n) case "failed"
elliptic() == (n = 2 and degree(radcand()) = 3 => radcand(); "failed")
hyperelliptic() == (n=2 and odd? degree(radcand()) => radcand(); "failed")
diag v == [reduce monomial(qelt(v,i+mini), i) for i in 0..n1]
integralRepresents(v, d) ==
ib := basisvec()
represents
[qelt(ib, i) * (qelt(v, i) /$RF d) for i in mini .. maxIndex ib]
integralCoordinates f ==
v := coordinates f
ib := basisvec()
splitDenominator
[qelt(v,i) / qelt(ib,i) for i in mini .. maxIndex ib]$Vector(RF)
integralDerivationMatrix d ==
dlogp := differentiate(radicnd, d) / (n * radicnd)
v := basisvec()
cd := splitDenominator(
[(i - mini) * dlogp + differentiate(qelt(v, i), d) / qelt(v, i)
for i in mini..maxIndex v]$Vector(RF))
[diagonalMatrix(cd.num), cd.den]
-- return (d0,...,d(n-1)) s.t. (1/d0, y/d1,...,y**(n-1)/d(n-1))
-- is an integral basis for the curve y**d = p
-- requires that p has no factor of multiplicity >= d
iBasis(p, d) ==
pl := fullVector(squareFree p, d)
d1 := (d - 1)::N
[*/[pl.j ** ((i * j) quo d) for j in 0..d1] for i in 0..d1]
-- returns a vector [a0,a1,...,a_{m-1}] of length m such that
-- p = a0^0 a1^1 ... a_{m-1}^{m-1}
fullVector(p, m) ==
ans:PrimitiveArray(UP) := new(m, 0)
ans.0 := unit p
l := factors p
for i in 1..maxIndex ans repeat
ans.i :=
(u := find(#1.exponent = i, l)) case "failed" => 1
(u::REC).factor
ans
-- return (f0,...,f(n-1)) s.t. (f0, y f1,..., y**(n-1) f(n-1))
-- is a local integral basis at infinity for the curve y**d = p
inftyBasis(p, m) ==
rt := rootPoly(p(x := inv(monomial(1, 1)$UP :: RF)), m)
m ~= rt.exponent =>
error "Curve not irreducible after change of variable 0 -> infinity"
a := (rt.coef) x
b:RF := 1
v := iBasis(rt.radicand, m)
w:Vector(RF) := new(m, 0)
for i in mini..maxIndex v repeat
qsetelt!(w, i, b / ((qelt(v, i)::RF) x))
b := b * a
w
charPintbas(p, c, v, w) ==
degree(p) ~= n => error "charPintbas: should not happen"
q:UP2 := map(retract(#1)@UP, p)
ib := integralBasis()$FunctionFieldIntegralBasis(UP, UP2,
SimpleAlgebraicExtension(UP, UP2, q))
not diagonal?(ib.basis)=> error "charPintbas: integral basis not diagonal"
a:RF := 1
for i in minRowIndex(ib.basis) .. maxRowIndex(ib.basis)
for j in minColIndex(ib.basis) .. maxColIndex(ib.basis)
for k in mini .. maxIndex v repeat
qsetelt!(v, k, (qelt(ib.basis, i, j) / ib.basisDen) * a)
qsetelt!(w, k, qelt(ib.basisInv, i, j) * inv a)
a := a * c
charPStartUp() ==
r := mkIntegral modulus
charPintbas(r.poly, r.coef, ibasis, invibasis)
x := inv(monomial(1, 1)$UP :: RF)
invmod := monomial(1, n)$UPUP - (radicnd x)::UPUP
r := mkIntegral invmod
charPintbas(r.poly, (r.coef) x, infbasis, invinfbasis)
startUp b ==
setref(brandNew?,b)
if zero?(p := characteristic$F) or p > n then char0StartUp()
else charPStartUp()
dsc:RF := ((-1)$Z ** ((n *$N n1) quo 2::N) * (n::Z)**n)$Z *
radicnd ** n1 *
*/[qelt(ibasis, i) ** 2 for i in mini..maxIndex ibasis]
setref(discPoly,primitivePart(numer dsc) / denom(dsc))
char0StartUp() ==
rp := rootPoly(radicnd, n)
rp.exponent ~= n => error "RadicalFunctionField: curve is not irreducible"
setref(newrad,rp.radicand)
ib := iBasis(deref newrad, n)
infb := inftyBasis(radicnd, n)
invden:RF := 1
for i in mini..maxIndex ib repeat
qsetelt!(invibasis, i, a := qelt(ib, i) * invden)
qsetelt!(ibasis, i, inv a)
invden := invden / rp.coef -- always equals 1/rp.coef**(i-mini)
qsetelt!(infbasis, i, a := qelt(infb, i))
qsetelt!(invinfbasis, i, inv a)
ramified?(p:UP) ==
(r := retractIfCan(p)@Union(F, "failed")) case F =>
singular?(r::F)
(radcand() exquo p) case UP
singular?(p:UP) ==
(r := retractIfCan(p)@Union(F, "failed")) case F =>
singular?(r::F)
(radcand() exquo(p**2)) case UP
branchPoint?(p:UP) ==
(r := retractIfCan(p)@Union(F, "failed")) case F =>
branchPoint?(r::F)
((q := (radcand() exquo p)) case UP) and
((q::UP exquo p) case "failed")
singular?(point:F) ==
zero?(radcand() point) and
zero?(((radcand() exquo (monomial(1,1)$UP-point::UP))::UP) point)
branchPoint?(point:F) ==
zero?(radcand() point) and not
zero?(((radcand() exquo (monomial(1,1)$UP-point::UP))::UP) point)
@
\section{domain ALGFF AlgebraicFunctionField}
<<domain ALGFF AlgebraicFunctionField>>=
)abbrev domain ALGFF AlgebraicFunctionField
++ Function field defined by f(x, y) = 0
++ Author: Manuel Bronstein
++ Date Created: 3 May 1988
++ Date Last Updated: 24 Jul 1990
++ Keywords: algebraic, curve, function, field.
++ Description: Function field defined by f(x, y) = 0.
++ Examples: )r ALGFF INPUT
AlgebraicFunctionField(F, UP, UPUP, modulus): Exports == Impl where
F : Field
UP : UnivariatePolynomialCategory F
UPUP : UnivariatePolynomialCategory Fraction UP
modulus: UPUP
N ==> NonNegativeInteger
Z ==> Integer
RF ==> Fraction UP
QF ==> Fraction UPUP
UP2 ==> SparseUnivariatePolynomial UP
SAE ==> SimpleAlgebraicExtension(RF, UPUP, modulus)
INIT ==> if (deref brandNew?) then startUp false
Exports ==> FunctionFieldCategory(F, UP, UPUP) with
knownInfBasis: N -> Void
++ knownInfBasis(n) \undocumented{}
Impl ==> SAE add
import ChangeOfVariable(F, UP, UPUP)
import InnerCommonDenominator(UP, RF, Vector UP, Vector RF)
import MatrixCommonDenominator(UP, RF)
import UnivariatePolynomialCategoryFunctions2(RF, UPUP, UP, UP2)
startUp : Boolean -> Void
vect : Matrix RF -> Vector $
getInfBasis: () -> Void
brandNew?:Reference(Boolean) := ref true
infBr?:Reference(Boolean) := ref true
discPoly:Reference(RF) := ref 0
n := degree modulus
n1 := (n - 1)::N
ibasis:Matrix(RF) := zero(n, n)
invibasis:Matrix(RF) := copy ibasis
infbasis:Matrix(RF) := copy ibasis
invinfbasis:Matrix(RF):= copy ibasis
branchPointAtInfinity?() == (INIT; deref infBr?)
discriminant() == (INIT; deref discPoly)
integralBasis() == (INIT; vect ibasis)
integralBasisAtInfinity() == (INIT; vect infbasis)
integralMatrix() == (INIT; ibasis)
inverseIntegralMatrix() == (INIT; invibasis)
integralMatrixAtInfinity() == (INIT; infbasis)
branchPoint?(a:F) == zero?((retract(discriminant())@UP) a)
definingPolynomial() == modulus
inverseIntegralMatrixAtInfinity() == (INIT; invinfbasis)
vect m ==
[represents row(m, i) for i in minRowIndex m .. maxRowIndex m]
integralCoordinates f ==
splitDenominator(coordinates(f) * inverseIntegralMatrix())
knownInfBasis d ==
if deref brandNew? then
alpha := [monomial(1, d * i)$UP :: RF for i in 0..n1]$Vector(RF)
ib := diagonalMatrix
[inv qelt(alpha, i) for i in minIndex alpha .. maxIndex alpha]
invib := diagonalMatrix alpha
for i in minRowIndex ib .. maxRowIndex ib repeat
for j in minColIndex ib .. maxColIndex ib repeat
infbasis(i, j) := qelt(ib, i, j)
invinfbasis(i, j) := invib(i, j)
getInfBasis() ==
x := inv(monomial(1, 1)$UP :: RF)
invmod := map(#1 x, modulus)
r := mkIntegral invmod
degree(r.poly) ~= n => error "Should not happen"
ninvmod:UP2 := map(retract(#1)@UP, r.poly)
alpha := [(r.coef ** i) x for i in 0..n1]$Vector(RF)
invalpha := [inv qelt(alpha, i)
for i in minIndex alpha .. maxIndex alpha]$Vector(RF)
invib := integralBasis()$FunctionFieldIntegralBasis(UP, UP2,
SimpleAlgebraicExtension(UP, UP2, ninvmod))
for i in minRowIndex ibasis .. maxRowIndex ibasis repeat
for j in minColIndex ibasis .. maxColIndex ibasis repeat
infbasis(i, j) := ((invib.basis)(i,j) / invib.basisDen) x
invinfbasis(i, j) := ((invib.basisInv) (i, j)) x
ib2 := infbasis * diagonalMatrix alpha
invib2 := diagonalMatrix(invalpha) * invinfbasis
for i in minRowIndex ib2 .. maxRowIndex ib2 repeat
for j in minColIndex ibasis .. maxColIndex ibasis repeat
infbasis(i, j) := qelt(ib2, i, j)
invinfbasis(i, j) := invib2(i, j)
startUp b ==
setref(brandNew?,b)
nmod:UP2 := map(retract, modulus)
ib := integralBasis()$FunctionFieldIntegralBasis(UP, UP2,
SimpleAlgebraicExtension(UP, UP2, nmod))
for i in minRowIndex ibasis .. maxRowIndex ibasis repeat
for j in minColIndex ibasis .. maxColIndex ibasis repeat
qsetelt!(ibasis, i, j, (ib.basis)(i, j) / ib.basisDen)
invibasis(i, j) := ((ib.basisInv) (i, j))::RF
if zero?(infbasis(minRowIndex infbasis, minColIndex infbasis))
then getInfBasis()
ib2 := coordinates normalizeAtInfinity vect ibasis
invib2 := inverse(ib2)::Matrix(RF)
for i in minRowIndex ib2 .. maxRowIndex ib2 repeat
for j in minColIndex ib2 .. maxColIndex ib2 repeat
ibasis(i, j) := qelt(ib2, i, j)
invibasis(i, j) := invib2(i, j)
dsc := resultant(modulus, differentiate modulus)
dsc0 := dsc * determinant(infbasis) ** 2
degree(numer dsc0) > degree(denom dsc0) =>error "Shouldn't happen"
setref(infBr?,degree(numer dsc0) < degree(denom dsc0))
dsc := dsc * determinant(ibasis) ** 2
setref(discPoly,primitivePart(numer dsc) / denom(dsc))
integralDerivationMatrix d ==
w := integralBasis()
splitDenominator(coordinates([differentiate(w.i, d)
for i in minIndex w .. maxIndex w]$Vector($))
* inverseIntegralMatrix())
integralRepresents(v, d) ==
represents(coordinates(represents(v, d)) * integralMatrix())
branchPoint?(p:UP) ==
INIT
(r:=retractIfCan(p)@Union(F,"failed")) case F =>branchPoint?(r::F)
not ground? gcd(retract(discriminant())@UP, p)
@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--Copyright (C) 2007-2010, Gabriel Dos Reis.
--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.
@
<<*>>=
<<license>>
-- SPAD files for the integration world should be compiled in the
-- following order:
--
-- intaux rderf intrf CURVE curvepkg divisor pfo
-- intalg intaf efstruc rdeef intef irexpand integrat
<<category FFCAT FunctionFieldCategory>>
<<package MMAP MultipleMap>>
<<package FFCAT2 FunctionFieldCategoryFunctions2>>
<<package CHVAR ChangeOfVariable>>
<<domain RADFF RadicalFunctionField>>
<<domain ALGFF AlgebraicFunctionField>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}