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

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

)abbrev package FORDER FindOrderFinite
++ Finds the order of a divisor over a finite field
++ Author: Manuel Bronstein
++ Date Created: 1988
++ Date Last Updated: 11 Jul 1990
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
++ Date Last Updated: 8 November 1994
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
++ Date Last Updated: 11 Jul 1990
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 IntegerPrimesPackage(Z)
    import 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 ==
      lcm [commonDenominator(retract(c)@UP) for c in coefficients r] * 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
++ Date Last Updated: 3 August 1993
++ 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
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 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
      one? n => n
      m   := rat(modulus, d, getGoodPrime(p * bad))
      n = m => n
      0

    rat(pp, d, p) ==
      gf := InnerPrimeField p
      order(d, pp,
        numer(#1)::gf / denom(#1)::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 one?(#(numer ideal d)) => 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
++ Date Last Updated: 11 Jul 1990
++ 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 new().
++ Keywords: function, space, redcution.
FunctionSpaceReduce(R, F): Exports == Implementation where
  R: Join(IntegralDomain, RetractableTo Integer)
  F: FunctionSpace R

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

  Exports ==> with
    bringDown: F -> Q
	++ bringDown(f) \undocumented
    bringDown: (F, K) -> UP
	++ bringDown(f,k) \undocumented
    newReduc : () -> Void
	++ newReduc() \undocumented

  Implementation ==> add
    macro ALGOP  == '%alg
    import SparseUnivariatePolynomialFunctions2(F, Q)
    import PolynomialCategoryQuotientFunctions(IndexedExponents K,
                         K, R, SparseMultivariatePolynomial(R, K), F)

    K2Z : K -> F

    redmap := table()$AssociationList(K, Z)

    newReduc() ==
      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()$Z)::F
      u::Z::F

)abbrev package PFO PointsOfFiniteOrder
++ Finds the order of a divisor on a curve
++ Author: Manuel Bronstein
++ Date Created: 1988
++ Date Last Updated: 22 July 1998
++ 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(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
  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:Z, dfpoly:UPQ)
  RC0 ==> Record(ncurve:UPUPQ, disc:Z)
  ID  ==> FractionalIdeal(UP, QF, UPUP, R)
  SMP ==> SparseMultivariatePolynomial(R0,K)

  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
    macro ALGOP == '%alg
    import IntegerPrimesPackage(Z)
    import PointsOfFiniteOrderTools(UPQ, UPUPQ)
    import UnivariatePolynomialCommonDenominator(Z, Q, UPQ)

    cmult: List SMP -> SMP
    raise         : (UPQ, K) -> F
    raise2        : (UP2, K) -> UP
    qmod          : F     -> Q
    fmod          : UPF   -> UPQ
    rmod          : UP    -> UPQ
    pmod          : UPUP  -> UPUPQ
    kqmod         : (F,    K) -> UPQ
    krmod         : (UP,   K) -> UP2
    kpmod         : (UPUP, K) -> UP3
    selectIntegers: K   -> REC
    selIntegers:    ()  -> RC0
    possibleOrder : FD -> N
    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
    toQ1          : (UP2, UPQ) -> UP
    toQ2          : (UP3, UPQ) -> R
    Q2F           : Q -> F
    Q2UPUP        : UPUPQ -> UPUP

    q := FunctionSpaceReduce(R0, F)

    torsion? d == order(d) case N
    Q2F x      == numer(x)::F / denom(x)::F
    qmod x     == bringDown(x)$q
    kqmod(x,k) == bringDown(x, k)$q
    fmod p     == map(qmod, p)$SparseUnivariatePolynomialFunctions2(F, Q)
    pmod p     == map(qmod, p)$MultipleMap(F, UP, UPUP, Q, UPQ, UPUPQ)
    Q2UPUP p   == map(Q2F, p)$MultipleMap(Q, UPQ, UPUPQ, F, UP, UPUP)
    klist d    == "setUnion"/[kernels c for c in coefficients d]
    notIrr? d  == #(factors 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(UPQ2F(#1, k), p)$UnivariatePolynomialCategoryFunctions2(UPQ,UP2,F,UP)

    UP32UPUP(p:UP3, k:K):UPUP ==
      map(UP22UP(#1,k)::QF,
          p)$UnivariatePolynomialCategoryFunctions2(UP2, UP3, QF, UPUP)

    cmult(l:List SMP):SMP ==
      R0 has GcdDomain => lcm l
      */l

    doubleDisc(f:UP3):Z ==
      d := discriminant f
      g := gcd(d, differentiate d)
      d := (d exquo g)::UP2
      zero?(e := discriminant d) => 0
      gcd [retract(c)@Z for c in coefficients e]

    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)

    alglist d ==
      n := numer(i := ideal d)
      select!(has?(operator #1, ALGOP),
               setUnion(klist denom i,
                 "setUnion"/[aklist qelt(n,i) for i in minIndex n..maxIndex n]))

    krmod(p,k) ==
       map(kqmod(#1, k),
           p)$UnivariatePolynomialCategoryFunctions2(F, UP, UPQ, UP2)

    rmod p ==
       map(qmod, p)$UnivariatePolynomialCategoryFunctions2(F, UP, Q, UPQ)

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

    raise2(p, k) ==
      map(raise(#1, k),
          p)$UnivariatePolynomialCategoryFunctions2(UPQ, UP2, F, UP)

    algcurve(d, rc, k) ==
      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), k)
      s  := kqmod(resultant(primitivePart separate(raise2(krmod(
                    retract(norm hh)@UP, k), k), b2).primePart, b2), k)
      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(retract(#1)@Z :: gf,
                         mm)$SparseUnivariatePolynomialFunctions2(Q, gf)
      one? degree m =>
        alpha := - coefficient(m, 0) / leadingCoefficient m
        order(d, pp,
           (map(numer(#1)::gf / denom(#1)::gf,
            kqmod(#1,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,
           reduce(map(numer(#1)::gf / denom(#1)::gf,
            kqmod(#1,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 [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
      one? n => 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, (qq := qmod #1;numer(qq)::gf / denom(qq)::gf)
                                    )$ReducedDivisor(F, UP, UPUP, R, gf)

-- returns the potential order of d, 0 if d is of infinite order
    possibleOrder d ==
      zero?(genus()) or one?(#(numer ideal d)) => 1
      empty?(la := alglist d) => ratcurve(d, selIntegers())
      not(empty? rest la) =>
           error "PFO::possibleOrder: more than 1 algebraic constant"
      algcurve(d, selectIntegers first la, first la)

    selIntegers():RC0 ==
      f := definingPolynomial()$R
      repeat
        r := polyred pmod f
        d := doubleDisc r
        if zero? d then
          newReduc()$q
        else
          return [r,d]

    selectIntegers(k:K):REC ==
      g := polyred(f := definingPolynomial()$R)
      p := minPoly k
      repeat
        r := kpmod(g, k)
        d := doubleDisc r
        if zero? d or notIrr? fmod p then
          newReduc()$q
        else
          return [r, d, splitDenominator(fmod p).num]

    toQ1(p, d) ==
      map(Q2F(retract(#1 rem d)@Q),
          p)$UnivariatePolynomialCategoryFunctions2(UPQ, UP2, F, UP)

    toQ2(p, d) ==
      reduce map(toQ1(#1, d)::QF,
        p)$UnivariatePolynomialCategoryFunctions2(UP2, UP3, QF, UPUP)

    kpmod(p, k) ==
      map(krmod(retract(#1)@UP, k),
        p)$UnivariatePolynomialCategoryFunctions2(QF, UPUP, UP2, UP3)

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

    kgetGoodPrime(rec, res, h, b, d) ==
      p:PI := 3
      u : Union(UPQ, "failed")
      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) ==
      zero?(rec.disc rem p) => "failed"
      gf := InnerPrimeField p
      l  := [f.factor for f in factors factor(map(retract(#1)@Z :: gf,
               rec.dfpoly)$SparseUnivariatePolynomialFunctions2(Q,
                 gf))$DistinctDegreeFactorize(gf,
                   SparseUnivariatePolynomial gf) | one?(f.exponent)]
      empty? l => "failed"
      mdg := first l
      for ff in rest l repeat
        if degree(ff) < degree(mdg) then mdg := ff
      md := map(convert(#1)@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))

