pfo.spad

)abbrev package FORDER FindOrderFinite
++ Finds the order of a divisor over a finite field
++ Author: Manuel Bronstein
++ Date Created: 1988
FindOrderFinite(F, UP, UPUP, R) : Exports == Implementation where
  F   : Join(Finite, Field)
  UP  : UnivariatePolynomialCategory F
  UPUP : UnivariatePolynomialCategory Fraction UP
  R   : FunctionFieldCategory(F, UP, UPUP)

  Exports ==> with
    order : FiniteDivisor(F, UP, UPUP, R) -> NonNegativeInteger
        ++ order(x) \undocumented
  Implementation ==> add
    order d ==
      dd := d := reduce d
      for i in 1.. repeat
        principal? dd => return(i::NonNegativeInteger)
        dd := reduce(d + dd)

)abbrev package RDIV ReducedDivisor
++ Finds the order of a divisor over a finite field
++ Author: Manuel Bronstein
++ Date Created: 1988
ReducedDivisor(F1, UP, UPUP, R, F2) : Exports == Implementation where
  F1    : Field
  UP    : UnivariatePolynomialCategory F1
  UPUP  : UnivariatePolynomialCategory Fraction UP
  R     : FunctionFieldCategory(F1, UP, UPUP)
  F2    : Join(Finite, Field)

  N     ==> NonNegativeInteger
  FD    ==> FiniteDivisor(F1, UP, UPUP, R)
  UP2   ==> SparseUnivariatePolynomial F2
  UPUP2 ==> SparseUnivariatePolynomial Fraction UP2

  Exports ==> with
    order : (FD, UPUP, F1 -> F2) -> N
        ++ order(f, u, g) \undocumented

  Implementation ==> add
    algOrder : (FD, UPUP, F1 -> F2)  -> N
    rootOrder : (FD, UP, N, F1 -> F2) -> N

-- pp is not necessarily monic
    order(d, pp, f) ==
      (r := retractIfCan(reductum pp)@Union(Fraction UP, "failed"))
        case "failed" => algOrder(d, pp, f)
      rootOrder(d, - retract(r@Fraction(UP) / leadingCoefficient pp)@UP,
                degree pp, f)

    algOrder(d, modulus, reduce) ==
      redmod := map(reduce, modulus)$MultipleMap(F1, UP, UPUP, F2, UP2, UPUP2)
      curve  := AlgebraicFunctionField(F2, UP2, UPUP2, redmod)
      order(map(reduce,
              d)$FiniteDivisorFunctions2(F1, UP, UPUP, R, F2, UP2, UPUP2, curve)
                                 )$FindOrderFinite(F2, UP2, UPUP2, curve)

    rootOrder(d, radicand, n, reduce) ==
      redrad := map(reduce,
           radicand)$UnivariatePolynomialCategoryFunctions2(F1, UP, F2, UP2)
      curve  := RadicalFunctionField(F2, UP2, UPUP2, redrad::Fraction UP2, n)
      order(map(reduce,
              d)$FiniteDivisorFunctions2(F1, UP, UPUP, R, F2, UP2, UPUP2, curve)
                                 )$FindOrderFinite(F2, UP2, UPUP2, curve)

)abbrev package PFOTOOLS PointsOfFiniteOrderTools
++ Utilities for PFOQ and PFO
++ Author: Manuel Bronstein
++ Date Created: 25 Aug 1988
PointsOfFiniteOrderTools(UP, UPUP) : Exports == Implementation where
  UP   : UnivariatePolynomialCategory Fraction Integer
  UPUP : UnivariatePolynomialCategory Fraction UP

  PI  ==> PositiveInteger
  N   ==> NonNegativeInteger
  Z   ==> Integer
  Q   ==> Fraction Integer

  Exports ==> with
    getGoodPrime : Z -> PI
      ++ getGoodPrime n returns the smallest prime not dividing n
    badNum       : UP   -> Record(den : Z, gcdnum : Z)
        ++ badNum(p) \undocumented
    badNum       : UPUP -> Z
        ++ badNum(u) \undocumented
    mix          : List Record(den : Z, gcdnum : Z) -> Z
        ++ mix(l) \undocumented
    doubleDisc   : UPUP -> Z
        ++ doubleDisc(u) \undocumented
    polyred      : UPUP -> UPUP
        ++ polyred(u) \undocumented

  Implementation ==> add
    import from IntegerPrimesPackage(Z)
    import from UnivariatePolynomialCommonDenominator(Z, Q, UP)

    mix l          == lcm(lcm [p.den for p in l], gcd [p.gcdnum for p in l])
    badNum(p : UPUP) == mix [badNum(retract(c)@UP) for c in coefficients p]

    polyred r ==
        ld := [commonDenominator(retract(c)@UP) for c in coefficients r]
        lcm(ld)*r

    badNum(p : UP) ==
      cd := splitDenominator p
      [cd.den, gcd [retract(c)@Z for c in coefficients(cd.num)]]

    getGoodPrime n ==
      p : PI := 3
      while zero?(n rem p) repeat
        p := nextPrime(p::Z)::PI
      p

    doubleDisc r ==
      d := retract(discriminant r)@UP
      retract(discriminant((d exquo gcd(d, differentiate d))::UP))@Z

)abbrev package PFOQ PointsOfFiniteOrderRational
++ Finds the order of a divisor on a rational curve
++ Author: Manuel Bronstein
++ Date Created: 25 Aug 1988
++ Description:
++ This package provides function for testing whether a divisor on a
++ curve is a torsion divisor.
++ Keywords: divisor, algebraic, curve.
++ Examples: )r PFOQ INPUT
-- Currently unused.
PointsOfFiniteOrderRational(UP, UPUP, R) : Exports == Implementation where
  UP   : UnivariatePolynomialCategory Fraction Integer
  UPUP : UnivariatePolynomialCategory Fraction UP
  R    : FunctionFieldCategory(Fraction Integer, UP, UPUP)

  PI  ==> PositiveInteger
  N   ==> NonNegativeInteger
  Z   ==> Integer
  Q   ==> Fraction Integer
  FD  ==> FiniteDivisor(Q, UP, UPUP, R)

  Exports ==> with
    order       : FD -> Union(N, "failed")
        ++ order(f) \undocumented
    torsion?    : FD -> Boolean
        ++ torsion?(f) \undocumented
    torsionIfCan : FD -> Union(Record(order:N, function:R), "failed")
        ++ torsionIfCan(f) \undocumented

  Implementation ==> add
    import from PointsOfFiniteOrderTools(UP, UPUP)

    possibleOrder : FD -> N
    ratcurve     : (FD, UPUP, Z) -> N
    rat          : (UPUP, FD, PI) -> N

    torsion? d  == order(d) case N

-- returns the potential order of d, 0 if d is of infinite order
    ratcurve(d, modulus, disc) ==
      mn  := minIndex(nm := numer(i := ideal d))
      h   := lift(hh := nm(mn + 1))
      s   := separate(retract(norm hh)@UP,
               b := retract(retract(nm.mn)@Fraction(UP))@UP).primePart
      bd  := badNum denom i
      r   := resultant(s, b)
      bad := lcm [disc, numer r, denom r, bd.den * bd.gcdnum, badNum h]$List(Z)
      n   := rat(modulus, d, p := getGoodPrime bad)
-- if n > 1 then it is cheaper to compute the order modulo a second prime,
-- since computing n * d could be very expensive
      (n = 1) => n
      m   := rat(modulus, d, getGoodPrime(p * bad))
      n = m => n
      0

    rat(pp, d, p) ==
      gf := InnerPrimeField p
      order(d, pp,
        (z1 : Q) : gf +-> numer(z1)::gf / denom(z1)::gf)$ReducedDivisor(Q, UP, UPUP, R, gf)

-- returns the potential order of d, 0 if d is of infinite order
    possibleOrder d ==
      zero?(genus()) or (#(numer ideal d) = 1) => 1
      r := polyred definingPolynomial()$R
      ratcurve(d, r, doubleDisc r)

    order d ==
      zero?(n := possibleOrder(d := reduce d)) => "failed"
      principal? reduce(n::Z * d) => n
      "failed"

    torsionIfCan d ==
      zero?(n := possibleOrder(d := reduce d)) => "failed"
      (g := generator reduce(n::Z * d)) case "failed" => "failed"
      [n, g@R]

)abbrev package FSRED FunctionSpaceReduce
++ Reduction from a function space to the rational numbers
++ Author: Manuel Bronstein
++ Date Created: 1988
++ Description:
++ This package provides function which replaces transcendental kernels
++ in a function space by random integers. The correspondence between
++ the kernels and the integers is fixed between calls to newReduc().
++ Keywords: function, space, redcution.
FunctionSpaceReduce(R, F) : Exports == Implementation where
  R : Join(Comparable, IntegralDomain, RetractableTo Integer)
  F : FunctionSpace R

  Z   ==> Integer
  Q   ==> Fraction Integer
  UP  ==> SparseUnivariatePolynomial Q
  K   ==> Kernel F

  Exports ==> with
    bringDown : F -> Q
        ++ bringDown(f) replaces transcendental kernels in f
        ++ by random integers.  Error if f contains algebraic
        ++ kernel.
    bringDown : (F, K) -> UP
        ++ bringDown(f, k) is like bringDown(f) but k is
        ++ algebraic kernel and the result is reduced modulo
        ++ minimal polynomial of k.
    newReduc : Integer -> Void
        ++ newReduc(n) empties reduction state.  n is size
        ++ of random integers to use during next reduction.

  Implementation ==> add
    import from SparseUnivariatePolynomialFunctions2(F, Q)
    import from PolynomialCategoryQuotientFunctions(IndexedExponents K,
                         K, R, SparseMultivariatePolynomial(R, K), F)

    ALGOP  := '%alg

    K2Z : K -> F

    redmap := table()$AssociationList(K, Z)
    nn : Integer := 10^6

    newReduc(n) ==
      nn := n
      for k in keys redmap repeat remove!(k, redmap)

    bringDown(f, k) ==
      ff := univariate(f, k)
      (bc := extendedEuclidean(map(bringDown, denom ff),
                m := map(bringDown, minPoly k), 1)) case "failed" =>
                     error "denominator is 0"
      (map(bringDown, numer ff) * bc.coef1) rem m

    bringDown f ==
      retract(eval(f, lk := kernels f, [K2Z k for k in lk]))@Q

    K2Z k ==
      has?(operator k, ALGOP) => error "Cannot reduce constant field"
      (u := search(k, redmap)) case "failed" =>
                                      setelt!(redmap, k, random(nn)$Z)::F
      u@Z::F

)abbrev package PFO PointsOfFiniteOrder
++ Finds the order of a divisor on a curve
++ Author: Manuel Bronstein
++ Date Created: 1988
++ Description:
++ This package provides function for testing whether a divisor on a
++ curve is a torsion divisor.
++ Keywords: divisor, algebraic, curve.
++ Examples: )r PFO INPUT
PointsOfFiniteOrder(R0, F, UP, UPUP, R) : Exports == Implementation where
  R0   : Join(Comparable, IntegralDomain, RetractableTo Integer)
  F    : FunctionSpace R0
  UP   : UnivariatePolynomialCategory F
  UPUP : UnivariatePolynomialCategory Fraction UP
  R    : FunctionFieldCategory(F, UP, UPUP)

  PI  ==> PositiveInteger
  N   ==> NonNegativeInteger
  Z   ==> Integer
  Q   ==> Fraction Integer
  UPF ==> SparseUnivariatePolynomial F
  UPQ ==> SparseUnivariatePolynomial Q
  UPZ ==> SparseUnivariatePolynomial(Z)
  QF  ==> Fraction UP
  UPUPQ ==> SparseUnivariatePolynomial Fraction UPQ
  UP2 ==> SparseUnivariatePolynomial UPQ
  UP3 ==> SparseUnivariatePolynomial UP2
  FD  ==> FiniteDivisor(F, UP, UPUP, R)
  K   ==> Kernel F
  REC ==> Record(ncurve : UP3, disc : UPZ, dfpoly : UPQ)
  RC0 ==> Record(ncurve : UPUPQ, disc : Z)
  ID  ==> FractionalIdeal(UP, QF, UPUP, R)
  SMP ==> SparseMultivariatePolynomial(R0, K)
  IREC ==> Record(ncurve : UPUP, n_div_numer : Vector(UPUP),
                  n_div_denom : UPUP, need_change : Boolean)

  Exports ==> with
    order        : FD -> Union(N, "failed")
        ++ order(fd) returns order of divisor \spad{fd} or \spad{"failed"} if
        ++ \spad{fd} is not of finite order.
    torsion?     : FD -> Boolean
        ++ torsion?(fd) returns \spad{true} if \spad{fd} is of finite order,
        ++ \spad{false} otherwise.
    torsionIfCan : FD -> Union(Record(order:N, function:R), "failed")
        ++ torsionIfCan(fd) returns \spad{[n, f]} such that
        ++ \spad{n*fd = div(f)} when \spad{fd} is of finite order.
        ++ torsionIfCan(fd) returns \spad{"failed"} if
        ++ \spad{fd} is not of finite order.
    cmult        : List SMP -> SMP
        ++ cmult(x) should be local but conditional
    possibleOrder : FD -> N
        ++ possibleOrder(d) returns n such that d is of order n,
        ++ or of infinite order.
    handle_imaginary : FD -> IREC
        ++ handle_imaginary(d) should be local but conditional.
    if R0 has CharacteristicZero and F has AlgebraicallyClosedField then
        simplifyCoeffs : (IREC, List K) -> N
          ++ simplifyCoeffs(d, la) should be local but conditional.

  Implementation ==> add
    import from IntegerPrimesPackage(Z)
    import from PointsOfFiniteOrderTools(UPQ, UPUPQ)
    import from UnivariatePolynomialCommonDenominator(Z, Q, UPQ)

    ALGOP  := '%alg

    selectIntegers : K   -> REC
    selIntegers :    ()  -> RC0
    ratcurve      : (FD, RC0)    -> N
    algcurve      : (FD, REC, K) -> N
    kbad3Num      : (UP3, UPQ) -> Z
    kbadBadNum    : (UP2, UPQ) -> Z
    kgetGoodPrime : (REC, UPQ, UP3, UP2, UP2) -> Record(prime : PI, poly : UPQ)
    goodRed       : (REC, UPQ, UP3, UP2, UP2, PI) -> Union(UPQ, "failed")
    good?         : (UPQ, UP3, UP2, UP2, PI, UPQ) -> Boolean
    klist         : UP -> List K
    aklist        : R  -> List K
    alglist       : FD -> List K
    notIrr?       : UPQ -> Boolean
    rat           : (UPUP, FD, PI) -> N

    q := FunctionSpaceReduce(R0, F)

    torsion? d == order(d) case N

    Q2F(x : Q) : F     == numer(x)::F / denom(x)::F

    Q2UPUP(p : UPUPQ) : UPUP ==
        map(Q2F, p)$MultipleMap(Q, UPQ, UPUPQ, F, UP, UPUP)

    raise(p : UPQ, k : K) : F ==
        (map(Q2F, p)$SparseUnivariatePolynomialFunctions2(Q, F)) (k::F)

    raise2(p : UP2, k : K) : UP ==
        map(z1 +-> raise(z1, k),
            p)$UnivariatePolynomialCategoryFunctions2(UPQ, UP2, F, UP)

    -- reduction functions
    qmod(x : F) : Q == bringDown(x)$q
    fmod(p : UPF) : UPQ ==
        map(qmod, p)$UnivariatePolynomialCategoryFunctions2(F, UPF, Q, UPQ)
    pmod(p : UPUP) : UPUPQ ==
        map(qmod, p)$MultipleMap(F, UP, UPUP, Q, UPQ, UPUPQ)
    rmod(p : UP) : UPQ ==
        map(qmod, p)$UnivariatePolynomialCategoryFunctions2(F, UP, Q, UPQ)

    kqmod(x : F, k : K) : UPQ == bringDown(x, k)$q
    krmod(p : UP, k : K) : UP2 ==
        map(z1 +-> kqmod(z1, k),
            p)$UnivariatePolynomialCategoryFunctions2(F, UP, UPQ, UP2)
    kpmod(p : UPUP, k : K) : UP3 ==
        map((p1 : QF) : UP2 +-> krmod(retract(p1)@UP, k),
            p)$UnivariatePolynomialCategoryFunctions2(QF, UPUP, UP2, UP3)


    klist d    == algtower(coefficients(d))
    notIrr? d  == #(factorList factor(d)$RationalFactorize(UPQ)) > 1
    kbadBadNum(d, m) == mix [badNum(c rem m) for c in coefficients d]
    kbad3Num(h, m)   == lcm [kbadBadNum(c, m) for c in coefficients h]

    torsionIfCan d ==
      zero?(n := possibleOrder(d := reduce d)) => "failed"
      (g := generator reduce(n::Z * d)) case "failed" => "failed"
      [n, g@R]

    UPQ2F(p : UPQ, k : K) : F ==
      map(Q2F, p)$UnivariatePolynomialCategoryFunctions2(Q, UPQ, F, UP) (k::F)

    UP22UP(p : UP2, k : K) : UP ==
      map((p1 : UPQ) : F +-> UPQ2F(p1, k), p)$UnivariatePolynomialCategoryFunctions2(UPQ, UP2, F, UP)

    UP32UPUP(p : UP3, k : K) : UPUP ==
      map((p1 : UP2) : QF +-> UP22UP(p1, k)::QF,
          p)$UnivariatePolynomialCategoryFunctions2(UP2, UP3, QF, UPUP)

    if R0 has GcdDomain then
       cmult(l : List SMP) : SMP == lcm l
    else
       cmult(l : List SMP) : SMP == */l

    commonDen(p : UP) : SMP ==
      l1 : List F := coefficients p
      l2 : List SMP := [denom c for c in l1]
      cmult l2

    polyred(f : UPUP) : UPUP ==
      cmult([commonDen(retract(c)@UP) for c in coefficients f])::F::UP::QF * f

    aklist f ==
      (r := retractIfCan(f)@Union(QF, "failed")) case "failed" =>
        "setUnion"/[klist(retract(c)@UP) for c in coefficients lift f]
      klist(retract(r@QF)@UP)

    alglist0(d : FD) : List K ==
      n := numer(i := ideal d)
      select!((k1 : K) : Boolean +-> has?(operator k1, ALGOP),
               setUnion(klist denom i,
                 "setUnion"/[aklist qelt(n,i) for i in minIndex n..maxIndex n]))

    alglist(d : FD) : List K ==
        f := definingPolynomial()$R
        l1 : List K := "setUnion"/[klist(retract(c)@UP) for c in
                                   coefficients f]
        setUnion(l1, alglist0(d))

    alglist2(f : UPUP) : List K ==
        "setUnion"/[klist(retract(c)@UP) for c in coefficients f]

    alglist1(nc : UPUP, nva : Vector(UPUP), nd : UPUP) : List K ==
        l1 : List(UPUP) := entries(nva)
        l2 := cons(nc, cons(nd, l1))
        "setUnion"/[alglist2(f) for f in l2]

    if R0 is Complex(Integer) and F is Expression(R0) then

        subst_mp_ni(p : SMP, ni : F) : F ==
            map((c : R0) : R0 +-> real(c)::R0, p)::F +
              ni*map((c : R0) : R0 +-> imag(c)::R0, p)::F

        subst_upup_ni(upup : UPUP, ni : F) : UPUP ==
            map((a : QF) : QF +->
              map((b : UP) : UP +->
                map((c : F) : F +-> subst_mp_ni(numer(c), ni)/
                                    subst_mp_ni(denom(c), ni), b), a), upup)

        handle_imaginary(fd) ==
            -- p0 : F := (new()$Symbol::F)^2 + 1
            p0 := monomial(1, 2)$UPF + 1$UPF
            ni := rootOf(p0)
            f := definingPolynomial()$R
            nf := subst_upup_ni(f, ni)
            id := ideal(fd)$FD
            n := numer(id)
            nva := vector([subst_upup_ni(lift(vi), ni)
                           for vi in entries(n)])$Vector(UPUP)
            nd := subst_upup_ni(denom(id)::QF::UPUP, ni)
            [nf, nva, nd, f ~= nf]

    else

        handle_imaginary(fd) ==
            id := ideal(fd)$FD
            n := numer(id)
            d := denom(id)
            nva := vector([lift(vi) for vi in entries(n)])$Vector(UPUP)
            [definingPolynomial()$R, nva, d::QF::UPUP, false]

    UPQ_to_UPF(p : UPQ) : UPF ==
        map(Q2F, p)$SparseUnivariatePolynomialFunctions2(Q, F)

    algcurve(d, rc, k) ==
      k1 : K :=
          mp := minPoly k
          every?((x : F) : Boolean +->
                   retractIfCan(x)@Union(Q, "failed") case Q,
                 coefficients(mp)) => k
          F has AlgebraicallyClosedField =>
              first(kernels(rootOf(UPQ_to_UPF(rc.dfpoly))$F))
          error "cannot reduce root"
      mn := minIndex(n := numer(i := minimize ideal d))
      h  := kpmod(lift(hh := n(mn + 1)), k)
      b2 := primitivePart
                 raise2(b := krmod(retract(retract(n.mn)@QF)@UP, k), k1)
      s  := kqmod(resultant(primitivePart separate(raise2(krmod(
                    retract(norm hh)@UP, k), k1), b2).primePart, b2), k1)
      pr := kgetGoodPrime(rc, s, h, b, dd := krmod(denom i, k))
      p  := pr.prime
      pp := UP32UPUP(rc.ncurve, k)
      mm := pr.poly
      gf := InnerPrimeField p
      m  := map((z1 : Q) : gf +-> retract(z1)@Z :: gf,
                         mm)$SparseUnivariatePolynomialFunctions2(Q, gf)
      (degree m = 1) =>
        alpha := - coefficient(m, 0) / leadingCoefficient m
        order(d, pp,
           (z1 : F) : gf +-> (map((q1 : Q) : gf +-> numer(q1)::gf / denom(q1)::gf,
            kqmod(z1, k))$SparseUnivariatePolynomialFunctions2(Q, gf))(alpha)
                                   )$ReducedDivisor(F, UP, UPUP, R, gf)
        -- d1 := toQ1(dd, mm)
        -- rat(pp, divisor ideal([(toQ1(b, mm) / d1)::QF::R,
                                       -- inv(d1::QF) * toQ2(h, mm)])$ID, p)
      sae:= SimpleAlgebraicExtension(gf, SparseUnivariatePolynomial gf, m)
      order(d, pp,
           (z1 : F) : sae +-> reduce(map((q1 : Q) : gf +-> numer(q1)::gf / denom(q1)::gf,
            kqmod(z1, k))$SparseUnivariatePolynomialFunctions2(Q, gf))$sae
                                   )$ReducedDivisor(F, UP, UPUP, R, sae)

-- returns the potential order of d, 0 if d is of infinite order
    ratcurve(d, rc) ==
      mn  := minIndex(nm := numer(i := minimize ideal d))
      h   := pmod lift(hh := nm(mn + 1))
      b   := rmod(retract(retract(nm.mn)@QF)@UP)
      s   := separate(rmod(retract(norm hh)@UP), b).primePart
      bd  := badNum rmod denom i
      r   := resultant(s, b)
      bad := lcm [factorial(rank()$R), rc.disc, numer r, denom r, bd.den*bd.gcdnum, badNum h]$List(Z)
      pp  := Q2UPUP(rc.ncurve)
      n   := rat(pp, d, p := getGoodPrime bad)
-- if n > 1 then it is cheaper to compute the order modulo a second prime,
-- since computing n * d could be very expensive
      (n = 1) => n
      m   := rat(pp, d, getGoodPrime(p * bad))
      n = m => n
      0

-- returns the order of d mod p
    rat(pp, d, p) ==
      gf := InnerPrimeField p
      order(d, pp, (q1 : F) : gf +-> (qq := qmod q1;numer(qq)::gf / denom(qq)::gf)
                                    )$ReducedDivisor(F, UP, UPUP, R, gf)

    evalup(upup : UPUP, la : List K, lrhs : List F) : UPUP ==
        map((a : QF) : QF +->
            map((b : UP) : UP +->
                map((c : F) : F +-> eval(c, la, lrhs), b), a), upup)

    if R0 has CharacteristicZero and F has AlgebraicallyClosedField then
        simplifyCoeffs(irec : IREC, la : List K) : N ==
            rec := primitiveElement([k::F for k in la]
                               )$FunctionSpacePrimitiveElement(R0, F)
            y := rootOf(rec.prim)
            lrhs := [pol(y) for pol in rec.poly]$List(F)
            nf := evalup(irec.ncurve, la, lrhs)
            nR := AlgebraicFunctionField(F, UP, UPUP, nf)
            nFD := FiniteDivisor(F, UP, UPUP, nR)
            nFR := FractionalIdeal(UP, Fraction(UP), UPUP, nR)
            nPFO := PointsOfFiniteOrder(R0, F, UP, UPUP, nR)
            nd := reduce(evalup(irec.n_div_denom, la, lrhs))$nR
            va := [reduce(evalup(ni, la, lrhs))$nR/nd
                    for ni in entries(irec.n_div_numer)]$Vector(nR)
            ndiv := divisor(ideal(va)$nFR)$nFD
            possibleOrder(ndiv)$nPFO

-- returns the potential order of d, 0 if d is of infinite order
    possibleOrder d ==
      zero?(genus()) or (#(numer ideal d) = 1) => 1
      d := divisor(decompose(d).id)
      irec := handle_imaginary(d)
      la := alglist1(irec.ncurve, irec.n_div_numer, irec.n_div_denom)
      empty?(la) => ratcurve(d, selIntegers())
      not(empty? rest la) or irec.need_change =>
          if R0 has CharacteristicZero and F has AlgebraicallyClosedField then
               return simplifyCoeffs(irec, la)
          error "PFO::possibleOrder: more than 1 algebraic constant"
      nden := reduce(irec.n_div_denom)$R
      va := [reduce(ni)$R/nden for ni in entries(irec.n_div_numer)]$Vector(R)
      d := divisor(ideal(va)$FractionalIdeal(UP, Fraction(UP), UPUP, R))$FD
      algcurve(d, selectIntegers first la, first la)

    double_disc(f : UPUP) : F ==
        d := retract(discriminant(f))@UP
        g := gcd(d, differentiate(d))
        discriminant((d exquo g)::UP)

    selIntegers() : RC0 ==
      f := definingPolynomial()$R
      g := polyred(f)
      d0 := double_disc(g)
      n : Integer := 10^6
      repeat
          newReduc(n)$q
          d := retract(qmod(d0))@Z
          if not(zero?(d)) then return [pmod(g), d]

    selectIntegers(k : K) : REC ==
      g := polyred(f := definingPolynomial()$R)
      d0 := double_disc(g)
      p := minPoly k
      n : Integer := 10^6
      repeat
          newReduc(n)$q
          d1 := kqmod(d0, k)
          d := map(retract, d1)$SparseUnivariatePolynomialFunctions2(
                                Q, Z)
          zero?(d) => iterate
          r := kpmod(g, k)
          return [r, d, splitDenominator(fmod p).num]

    order d ==
      zero?(n := possibleOrder(d := reduce d)) => "failed"
      principal? reduce(n::Z * d) => n
      "failed"

    kgetGoodPrime(rec, res, h, b, d) ==
      p : PI := nextPrime(rank()$R)::PI
      while (u := goodRed(rec, res, h, b, d, p)) case "failed" repeat
        p := nextPrime(p::Z)::PI
      [p, u::UPQ]

    goodRed(rec, res, h, b, d, p) ==
      every?((c : Z) : Boolean +-> zero?(c rem p),
             coefficients(rec.disc)) => "failed"
      gf := InnerPrimeField p
      l  := [f.factor for f in factorList factor(map((z1 : Q) : gf +-> retract(z1)@Z :: gf,
               rec.dfpoly)$SparseUnivariatePolynomialFunctions2(Q,
                 gf))$UnivariateFiniteFieldFactorize(gf,
                   SparseUnivariatePolynomial gf) | (f.exponent = 1)]
      empty? l => "failed"
      discp := map(c +-> c::gf, rec.disc
                  )$SparseUnivariatePolynomialFunctions2(Z, gf)
      mdg := first l
      is_good := false
      for ff in l repeat
          discp rem ff = 0 => iterate
          is_good and degree(ff) < degree(mdg) => mdg := ff
          is_good := true
      not(is_good) => "failed"
      md := map((z1 : gf) : Q +-> convert(z1)@Z :: Q,
                 mdg)$SparseUnivariatePolynomialFunctions2(gf, Q)
      good?(res, h, b, d, p, md) => md
      "failed"

    good?(res, h, b, d, p, m) ==
      bd := badNum(res rem m)
      not (zero?(bd.den rem p) or zero?(bd.gcdnum rem p) or
        zero?(kbadBadNum(b, m) rem p) or zero?(kbadBadNum(d, m) rem p)
              or zero?(kbad3Num(h, m) rem p))

--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.