path: root/src/algebra/curve.spad.pamphlet

\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}