Initial population.

1 files changed, 946 insertions, 0 deletions

diff --git a/src/algebra/curve.spad.pamphlet b/src/algebra/curve.spad.pamphlet
new file mode 100644
index 00000000..5ca9495b
--- /dev/null
+++ b/src/algebra/curve.spad.pamphlet
@@ -0,0 +1,946 @@
+\documentclass{article}
+\usepackage{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.
+    represents             : (Vector UP, UP) -> $
+      ++ represents([A0,...,A(n-1)],D) returns
+      ++ \spad{(A0 + A1 y +...+ A(n-1)*y**(n-1))/D}.
+    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()
+    absolutelyIrreducible? == numberOfComponents() = 1
+    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
+        vars := [concat(string u, string i)::SY for i in 1..(n := #w)]
+        x    := "%%dummy1"::SY
+        y    := "%%dummy2"::SY
+        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}
+<<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 = 1) => [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; 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}
+<<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}
+<<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.
+@
+<<*>>=
+<<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}