\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} <>= )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?(infOrder(#1) < 0, 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 (n := infOrder f) > 0 => 0 zero? n => (leadingCoefficient numer f) / (leadingCoefficient denom f) error "f not locally integral at infinity" rfmonom n == n < 0 => 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} <>= )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} <>= )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} <>= )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} <>= )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; discPoly()) radcand() == (INIT; 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 void 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 == 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] discPoly() := primitivePart(numer dsc) / denom(dsc) void char0StartUp() == rp := rootPoly(radicnd, n) rp.exponent ~= n => error "RadicalFunctionField: curve is not irreducible" newrad() := rp.radicand ib := iBasis(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) void 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} <>= )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; infBr?()) discriminant() == (INIT; 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) void 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) void startUp b == 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" infBr?() := degree(numer dsc0) < degree(denom dsc0) dsc := dsc * determinant(ibasis) ** 2 discPoly() := primitivePart(numer dsc) / denom(dsc) void 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} <>= --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. @ <<*>>= <> -- 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 <> <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}