axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / FFCAT.spad

)abbrev category FFCAT FunctionFieldCategory
++ Author: Manuel Bronstein
++ Date Created: 1987
++ Date Last Updated: 19 Mai 1993
++ Description:
++ Function field of a curve
++ This category is a model for the function field of a
++ plane algebraic curve.

FunctionFieldCategory(F, UP, UPUP) : Category == SIG 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)

  SIG ==> MonogenicAlgebra(RF, UPUP) with

    numberOfComponents : () -> NonNegativeInteger
      ++numberOfComponents() returns the number of absolutely irreducible
      ++ components.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X numberOfComponents()$R

    genus : () -> NonNegativeInteger
      ++genus() returns the genus of one absolutely irreducible component
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X genus()$R

    absolutelyIrreducible? : () -> Boolean
      ++absolutelyIrreducible?() tests if the curve absolutely irreducible?
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R2 := RadicalFunctionField(INT, P0, P1, 2 * x**2, 4)
      ++X absolutelyIrreducible?()$R2

    rationalPoint? : (F, F) -> Boolean
      ++rationalPoint?(a, b) tests if \spad{(x=a,y=b)} is on the curve.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X rationalPoint?(0,0)$R
      ++X R2 := RadicalFunctionField(INT, P0, P1, 2 * x**2, 4)
      ++X rationalPoint?(0,0)$R2

    branchPointAtInfinity? : () -> Boolean
      ++branchPointAtInfinity?() tests if there is a branch point 
      ++ at infinity.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X branchPointAtInfinity?()$R
      ++X R2 := RadicalFunctionField(INT, P0, P1, 2 * x**2, 4)
      ++X branchPointAtInfinity?()$R

    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.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X integralBasis()$R

    integralBasisAtInfinity: () -> Vector $
      ++integralBasisAtInfinity() returns the local integral basis 
      ++ at infinity
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X integralBasisAtInfinity()$R

    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}.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X integralMatrix()$R

    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}.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X inverseIntegralMatrix()$R

    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()}.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X integralMatrixAtInfinity()$R

    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()}.
      ++
      ++X P0 := UnivariatePolynomial(x, Integer)
      ++X P1 := UnivariatePolynomial(y, Fraction P0)
      ++X R := RadicalFunctionField(INT, P0, P1, 1 - x**20, 20)
      ++X inverseIntegralMatrixAtInfinity()$R

    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? == 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((s:F):P +-> s::P, p)_
          $UnivariatePolynomialCategoryFunctions2(F, UP,
                                     P, SparseUnivariatePolynomial P)) x

      UPUP2P(p, x, y) ==
        (map((s:RF):P +-> UP2P(retract(s)@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_!(s+->zero?(degree(s, x)) and zero?(degree(s, 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?(s +-> infOrder(s) < 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, x +-> differentiate(x, d)$RF)