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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra padiclib.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package IBPTOOLS IntegralBasisPolynomialTools}
<<package IBPTOOLS IntegralBasisPolynomialTools>>=
)abbrev package IBPTOOLS IntegralBasisPolynomialTools
++ Author: Clifton Williamson
++ Date Created: 13 August 1993
++ Date Last Updated: 17 August 1993
++ Basic Operations: mapUnivariate, mapBivariate
++ Related Domains: PAdicWildFunctionFieldIntegralBasis(K,R,UP,F)
++ Also See: WildFunctionFieldIntegralBasis, FunctionFieldIntegralBasis
++ AMS Classifications:
++ Keywords: function field, finite field, integral basis
++ Examples:
++ References:
++ Description:  IntegralBasisPolynomialTools provides functions for
++   mapping functions on the coefficients of univariate and bivariate
++   polynomials.

IntegralBasisPolynomialTools(K,R,UP,L): Exports == Implementation where
  K  : Ring
  R  : UnivariatePolynomialCategory K
  UP : UnivariatePolynomialCategory R
  L  : Ring

  MAT ==> Matrix
  SUP ==> SparseUnivariatePolynomial

  Exports ==> with
    mapUnivariate: (L -> K,SUP L) -> R
      ++ mapUnivariate(f,p(x)) applies the function \spad{f} to the
      ++ coefficients of \spad{p(x)}.

    mapUnivariate: (K -> L,R) -> SUP L
      ++ mapUnivariate(f,p(x)) applies the function \spad{f} to the
      ++ coefficients of \spad{p(x)}.

    mapUnivariateIfCan: (L -> Union(K,"failed"),SUP L) -> Union(R,"failed")
      ++ mapUnivariateIfCan(f,p(x)) applies the function \spad{f} to the
      ++ coefficients of \spad{p(x)}, if possible, and returns
      ++ \spad{"failed"} otherwise.

    mapMatrixIfCan: (L -> Union(K,"failed"),MAT SUP L) -> Union(MAT R,"failed")
      ++ mapMatrixIfCan(f,mat) applies the function \spad{f} to the
      ++ coefficients of the entries of \spad{mat} if possible, and returns
      ++ \spad{"failed"} otherwise.

    mapBivariate: (K -> L,UP) -> SUP SUP L
      ++ mapBivariate(f,p(x,y)) applies the function \spad{f} to the
      ++ coefficients of \spad{p(x,y)}.

  Implementation ==> add

    mapUnivariate(f:L -> K,poly:SUP L) ==
      ans : R := 0
      while not zero? poly repeat
        ans := ans + monomial(f leadingCoefficient poly,degree poly)
        poly := reductum poly
      ans

    mapUnivariate(f:K -> L,poly:R) ==
      ans : SUP L := 0
      while not zero? poly repeat
        ans := ans + monomial(f leadingCoefficient poly,degree poly)
        poly := reductum poly
      ans

    mapUnivariateIfCan(f,poly) ==
      ans : R := 0
      while not zero? poly repeat
        (lc := f leadingCoefficient poly) case "failed" => return "failed"
        ans := ans + monomial(lc :: K,degree poly)
        poly := reductum poly
      ans

    mapMatrixIfCan(f,mat) ==
      m := nrows mat; n := ncols mat
      matOut : MAT R := new(m,n,0)
      for i in 1..m repeat for j in 1..n repeat
        (poly := mapUnivariateIfCan(f,qelt(mat,i,j))) case "failed" =>
          return "failed"
        qsetelt!(matOut,i,j,poly :: R)
      matOut

    mapBivariate(f,poly) ==
      ans : SUP SUP L := 0
      while not zero? poly repeat
        ans :=
          ans + monomial(mapUnivariate(f,leadingCoefficient poly),degree poly)
        poly := reductum poly
      ans

@
\section{package IBACHIN ChineseRemainderToolsForIntegralBases}
<<package IBACHIN ChineseRemainderToolsForIntegralBases>>=
)abbrev package IBACHIN ChineseRemainderToolsForIntegralBases
++ Author: Clifton Williamson
++ Date Created: 9 August 1993
++ Date Last Updated: 3 December 1993
++ Basic Operations: chineseRemainder, factorList
++ Related Domains: PAdicWildFunctionFieldIntegralBasis(K,R,UP,F)
++ Also See: WildFunctionFieldIntegralBasis, FunctionFieldIntegralBasis
++ AMS Classifications:
++ Keywords: function field, finite field, integral basis
++ Examples:
++ References:
++ Description:

ChineseRemainderToolsForIntegralBases(K,R,UP): Exports == Implementation where
  K  : FiniteFieldCategory
  R  : UnivariatePolynomialCategory K
  UP : UnivariatePolynomialCategory R

  DDFACT ==> DistinctDegreeFactorize
  I      ==> Integer
  L      ==> List
  L2     ==> ListFunctions2
  Mat    ==> Matrix R
  NNI    ==> NonNegativeInteger
  PI     ==> PositiveInteger
  Q      ==> Fraction R
  SAE    ==> SimpleAlgebraicExtension
  SUP    ==> SparseUnivariatePolynomial
  SUP2   ==> SparseUnivariatePolynomialFunctions2
  Result ==> Record(basis: Mat, basisDen: R, basisInv: Mat)

  Exports ==> with
    factorList: (K,NNI,NNI,NNI) -> L SUP K
	++ factorList(k,n,m,j) \undocumented

    listConjugateBases: (Result,NNI,NNI) -> List Result
      ++ listConjugateBases(bas,q,n) returns the list
      ++ \spad{[bas,bas^Frob,bas^(Frob^2),...bas^(Frob^(n-1))]}, where
      ++ \spad{Frob} raises the coefficients of all polynomials
      ++ appearing in the basis \spad{bas} to the \spad{q}th power.

    chineseRemainder: (List UP, List Result, NNI) -> Result
	++ chineseRemainder(lu,lr,n) \undocumented

  Implementation ==> add
    import ModularHermitianRowReduction(R)
    import TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)

    applyFrobToMatrix: (Matrix R,NNI) -> Matrix R
    applyFrobToMatrix(mat,q) ==
      -- raises the coefficients of the polynomial entries of 'mat'
      -- to the qth power
      m := nrows mat; n := ncols mat
      ans : Matrix R := new(m,n,0)
      for i in 1..m repeat for j in 1..n repeat
        qsetelt!(ans,i,j,map(#1 ** q,qelt(mat,i,j)))
      ans

    listConjugateBases(bas,q,n) ==
      outList : List Result := list bas
      b := bas.basis; bInv := bas.basisInv; bDen := bas.basisDen
      for i in 1..(n-1) repeat
        b := applyFrobToMatrix(b,q)
        bInv := applyFrobToMatrix(bInv,q)
        bDen := map(#1 ** q,bDen)
        newBasis : Result := [b,bDen,bInv]
        outList := concat(newBasis,outList)
      reverse! outList

    factorList(a,q,n,k) ==
      coef : SUP K := monomial(a,0); xx : SUP K := monomial(1,1)
      outList : L SUP K := list((xx - coef)**k)
      for i in 1..(n-1) repeat
        coef := coef ** q
        outList := concat((xx - coef)**k,outList)
      reverse! outList

    basisInfoToPolys: (Mat,R,R) -> L UP
    basisInfoToPolys(mat,lcm,den) ==
      n := nrows(mat) :: I; n1 := n - 1
      outList : L UP := empty()
      for i in 1..n repeat
        pp : UP := 0
        for j in 0..n1 repeat
          pp := pp + monomial((lcm quo den) * qelt(mat,i,j+1),j)
        outList := concat(pp,outList)
      reverse! outList

    basesToPolyLists: (L Result,R) -> L L UP
    basesToPolyLists(basisList,lcm) ==
      [basisInfoToPolys(b.basis,lcm,b.basisDen) for b in basisList]

    OUT ==> OutputForm

    approximateExtendedEuclidean: (UP,UP,R,NNI) -> Record(coef1:UP,coef2:UP)
    approximateExtendedEuclidean(f,g,p,n) ==
      -- f and g are monic and relatively prime (mod p)
      -- function returns [coef1,coef2] such that
      -- coef1 * f + coef2 * g = 1 (mod p^n)
      sae := SAE(K,R,p)
      fSUP : SUP R := makeSUP f; gSUP : SUP R := makeSUP g
      fBar : SUP sae := map(convert(#1)@sae,fSUP)$SUP2(R,sae)
      gBar : SUP sae := map(convert(#1)@sae,gSUP)$SUP2(R,sae)
      ee := extendedEuclidean(fBar,gBar)
      not one?(ee.generator) =>
        error "polynomials aren't relatively prime"
      ss1 := ee.coef1; tt1 := ee.coef2
      s1 : SUP R := map(convert(#1)@R,ss1)$SUP2(sae,R); s := s1
      t1 : SUP R := map(convert(#1)@R,tt1)$SUP2(sae,R); t := t1
      pPower := p
      for i in 2..n repeat
        num := 1 - s * fSUP - t * gSUP
        rhs := (num exquo pPower) :: SUP R
        sigma := map(#1 rem p,s1 * rhs); tau := map(#1 rem p,t1 * rhs)
        s := s + pPower * sigma; t := t + pPower * tau
        quorem := monicDivide(s,gSUP)
        pPower := pPower * p
        s := map(#1 rem pPower,quorem.remainder)
        t := map(#1 rem pPower,t + fSUP * (quorem.quotient))
      [unmakeSUP s,unmakeSUP t]

    --mapChineseToList: (L SUP Q,L SUP Q,I) -> L SUP Q
    --mapChineseToList(list,polyList,i) ==
    mapChineseToList: (L UP,L UP,I,R) -> L UP
    mapChineseToList(list,polyList,i,den) ==
      -- 'polyList' consists of MONIC polynomials
      -- computes a polynomial p such that p = pp (modulo polyList[i])
      -- and p = 0 (modulo polyList[j]) for j ~= i for each 'pp' in 'list'
      -- create polynomials
      q : UP := 1
      for j in 1..(i-1) repeat
        q := q * first polyList
        polyList := rest polyList
      p := first polyList
      polyList := rest polyList
      for j in (i+1).. while not empty? polyList repeat
        q := q * first polyList
        polyList := rest polyList
      --p := map((numer(#1) rem den)/1, p)
      --q := map((numer(#1) rem den)/1, q)
      -- 'den' is a power of an irreducible polynomial
      --!! make this computation more efficient!!
      factoredDen := factor(den)$DDFACT(K,R)
      prime := nthFactor(factoredDen,1)
      n := nthExponent(factoredDen,1) :: NNI
      invPoly := approximateExtendedEuclidean(q,p,prime,n).coef1
      -- monicDivide may be inefficient?
      [monicDivide(pp * invPoly * q,p * q).remainder for pp in list]

    polyListToMatrix: (L UP,NNI) -> Mat
    polyListToMatrix(polyList,n) ==
      mat : Mat := new(n,n,0)
      for i in 1..n for poly in polyList repeat
        while not zero? poly repeat
          mat(i,degree(poly) + 1) := leadingCoefficient poly
          poly := reductum poly
      mat

    chineseRemainder(factors,factorBases,n) ==
      denLCM : R := reduce("lcm",[base.basisDen for base in factorBases])
      denLCM = 1 => [scalarMatrix(n,1),1,scalarMatrix(n,1)]
      -- compute local basis polynomials with denominators cleared
      factorBasisPolyLists := basesToPolyLists(factorBases,denLCM)
      -- use Chinese remainder to compute basis polynomials w/o denominators
      basisPolyLists : L L UP := empty()
      for i in 1.. for pList in factorBasisPolyLists repeat
        polyList := mapChineseToList(pList,factors,i,denLCM)
        basisPolyLists := concat(polyList,basisPolyLists)
      basisPolys := concat reverse! basisPolyLists
      mat := squareTop rowEchelon(polyListToMatrix(basisPolys,n),denLCM)
      matInv := UpTriBddDenomInv(mat,denLCM)
      [mat,denLCM,matInv]

@
\section{package PWFFINTB PAdicWildFunctionFieldIntegralBasis}
<<package PWFFINTB PAdicWildFunctionFieldIntegralBasis>>=
)abbrev package PWFFINTB PAdicWildFunctionFieldIntegralBasis
++ Author: Clifton Williamson
++ Date Created: 5 July 1993
++ Date Last Updated: 17 August 1993
++ Basic Operations: integralBasis, localIntegralBasis
++ Related Domains: WildFunctionFieldIntegralBasis(K,R,UP,F)
++ Also See: FunctionFieldIntegralBasis
++ AMS Classifications:
++ Keywords: function field, finite field, integral basis
++ Examples:
++ References:
++ Description:
++ In this package K is a finite field, R is a ring of univariate
++ polynomials over K, and F is a monogenic algebra over R.
++ We require that F is monogenic, i.e. that \spad{F = K[x,y]/(f(x,y))},
++ because the integral basis algorithm used will factor the polynomial
++ \spad{f(x,y)}.  The package provides a function to compute the integral
++ closure of R in the quotient field of F as well as a function to compute
++ a "local integral basis" at a specific prime.

PAdicWildFunctionFieldIntegralBasis(K,R,UP,F): Exports == Implementation where
  K  : FiniteFieldCategory
  R  : UnivariatePolynomialCategory K
  UP : UnivariatePolynomialCategory R
  F  : MonogenicAlgebra(R,UP)

  I        ==> Integer
  L        ==> List
  L2       ==> ListFunctions2
  Mat      ==> Matrix R
  NNI      ==> NonNegativeInteger
  PI       ==> PositiveInteger
  Q        ==> Fraction R
  SAE      ==> SimpleAlgebraicExtension
  SUP      ==> SparseUnivariatePolynomial
  CDEN     ==> CommonDenominator
  DDFACT   ==> DistinctDegreeFactorize
  WFFINTBS ==> WildFunctionFieldIntegralBasis
  Result   ==> Record(basis: Mat, basisDen: R, basisInv:Mat)
  IResult  ==> Record(basis: Mat, basisDen: R, basisInv:Mat,discr: R)
  IBPTOOLS ==> IntegralBasisPolynomialTools
  IBACHIN  ==> ChineseRemainderToolsForIntegralBases
  IRREDFFX ==> IrredPolyOverFiniteField
  GHEN     ==> GeneralHenselPackage

  Exports ==> with
    integralBasis : () -> Result
      ++ \spad{integralBasis()} returns a record
      ++ \spad{[basis,basisDen,basisInv] } containing information regarding
      ++ the integral closure of R in the quotient field of the framed
      ++ algebra F.  F is a framed algebra with R-module basis
      ++ \spad{w1,w2,...,wn}.
      ++ If 'basis' is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of 'basis' contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix 'basisInv' contains the coordinates of \spad{wi} with respect
      ++ to the basis \spad{v1,...,vn}: if 'basisInv' is the matrix
      ++ \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
    localIntegralBasis : R -> Result
      ++ \spad{integralBasis(p)} returns a record
      ++ \spad{[basis,basisDen,basisInv] } containing information regarding
      ++ the local integral closure of R at the prime \spad{p} in the quotient
      ++ field of the framed algebra F.  F is a framed algebra with R-module
      ++ basis \spad{w1,w2,...,wn}.
      ++ If 'basis' is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the local integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of 'basis' contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix 'basisInv' contains the coordinates of \spad{wi} with respect
      ++ to the basis \spad{v1,...,vn}: if 'basisInv' is the matrix
      ++ \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
    reducedDiscriminant: UP -> R
	++ reducedDiscriminant(up) \undocumented

  Implementation ==> add
    import IntegralBasisTools(R, UP, F)
    import GeneralHenselPackage(R,UP)
    import ModularHermitianRowReduction(R)
    import TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)

    reducedDiscriminant f ==
      ff : SUP Q := mapUnivariate(#1 :: Q,f)$IBPTOOLS(R,UP,SUP UP,Q)
      ee := extendedEuclidean(ff,differentiate ff)
      cc := concat(coefficients(ee.coef1),coefficients(ee.coef2))
      cden := splitDenominator(cc)$CDEN(R,Q,L Q)
      denom := cden.den
      gg := gcd map(numer,cden.num)$L2(Q,R)
      (ans := denom exquo gg) case "failed" =>
        error "PWFFINTB: error in reduced discriminant computation"
      ans :: R

    compLocalBasis: (UP,R) -> Result
    compLocalBasis(poly,prime) ==
      -- compute a local integral basis at 'prime' for k[x,y]/(poly(x,y)).
      sae := SAE(R,UP,poly)
      localIntegralBasis(prime)$WFFINTBS(K,R,UP,sae)

    compLocalBasisOverExt: (UP,R,UP,NNI) -> Result
    compLocalBasisOverExt(poly0,prime0,irrPoly0,k) ==
      -- poly0 = irrPoly0**k (mod prime0)
      n := degree poly0; disc0 := discriminant poly0
      (disc0 exquo prime0) case "failed" =>
        [scalarMatrix(n,1), 1, scalarMatrix(n,1)]
      r := degree irrPoly0
      -- extend scalars:
      -- construct irreducible polynomial of degree r over K
      irrPoly := generateIrredPoly(r :: PI)$IRREDFFX(K)
      -- construct extension of degree r over K
      E := SAE(K,SUP K,irrPoly)
      -- lift coefficients to elements of E
      poly := mapBivariate(#1 :: E,poly0)$IBPTOOLS(K,R,UP,E)
      redDisc0 := reducedDiscriminant poly0
      redDisc := mapUnivariate(#1 :: E,redDisc0)$IBPTOOLS(K,R,UP,E)
      prime := mapUnivariate(#1 :: E,prime0)$IBPTOOLS(K,R,UP,E)
      sae := SAE(E,SUP E,prime)
      -- reduction (mod prime) of polynomial of which poly is the kth power
      redIrrPoly :=
        pp := mapBivariate(#1 :: E,irrPoly0)$IBPTOOLS(K,R,UP,E)
        mapUnivariate(reduce,pp)$IBPTOOLS(SUP E,SUP SUP E,SUP SUP SUP E,sae)
      -- factor the reduction
      factorListSAE := factors factor(redIrrPoly)$DDFACT(sae,SUP sae)
      -- list the 'primary factors' of the reduction of poly
      redFactors : List SUP sae := [(f.factor)**k for f in factorListSAE]
      -- lift these factors to elements of SUP SUP E
      primaries : List SUP SUP E :=
        [mapUnivariate(lift,ff)$IBPTOOLS(SUP E,SUP SUP E,SUP SUP SUP E,sae) _
             for ff in redFactors]
      -- lift the factors to factors modulo a suitable power of 'prime'
      deg := (1 + order(redDisc,prime) * degree(prime)) :: PI
      henselInfo := HenselLift(poly,primaries,prime,deg)$GHEN(SUP E,SUP SUP E)
      henselFactors := henselInfo.plist
      psi1 := first henselFactors
      FF := SAE(SUP E,SUP SUP E,psi1)
      factorIb := localIntegralBasis(prime)$WFFINTBS(E,SUP E,SUP SUP E,FF)
      bs := listConjugateBases(factorIb,size()$K,r)$IBACHIN(E,SUP E,SUP SUP E)
      ib := chineseRemainder(henselFactors,bs,n)$IBACHIN(E,SUP E,SUP SUP E)
      b : Matrix R :=
        bas := mapMatrixIfCan(retractIfCan,ib.basis)$IBPTOOLS(K,R,UP,E)
        bas case "failed" => error "retraction of basis failed"
        bas :: Matrix R
      bInv : Matrix R :=
        --bas := mapMatrixIfCan(ric,ib.basisInv)$IBPTOOLS(K,R,UP,E)
        bas := mapMatrixIfCan(retractIfCan,ib.basisInv)$IBPTOOLS(K,R,UP,E)
        bas case "failed" => error "retraction of basis inverse failed"
        bas :: Matrix R
      bDen : R :=
        p := mapUnivariateIfCan(retractIfCan,ib.basisDen)$IBPTOOLS(K,R,UP,E)
        p case "failed" => error "retraction of basis denominator failed"
        p :: R
      [b,bDen,bInv]

    padicLocalIntegralBasis: (UP,R,R,R) -> IResult
    padicLocalIntegralBasis(p,disc,redDisc,prime) ==
      -- polynomials in x modulo 'prime'
      sae := SAE(K,R,prime)
      -- find the factorization of 'p' modulo 'prime' and lift the
      -- prime powers to elements of UP:
      -- reduce 'p' modulo 'prime'
      reducedP := mapUnivariate(reduce,p)$IBPTOOLS(R,UP,SUP UP,sae)
      -- factor the reduced polynomial
      factorListSAE := factors factor(reducedP)$DDFACT(sae,SUP sae)
      -- if only one prime factor, perform usual integral basis computation
      (# factorListSAE) = 1 =>
        ib := localIntegralBasis(prime)$WFFINTBS(K,R,UP,F)
        index := diagonalProduct(ib.basisInv)
        [ib.basis,ib.basisDen,ib.basisInv,disc quo (index * index)]
      -- list the 'prime factors' of the reduced polynomial
      redPrimes : List SUP sae :=
        [f.factor for f in factorListSAE]
      -- lift these factors to elements of UP
      primes : List UP :=
        [mapUnivariate(lift,ff)$IBPTOOLS(R,UP,SUP UP,sae) for ff in redPrimes]
      -- list the exponents
      expons : List NNI := [((f.exponent) :: NNI) for f in factorListSAE]
      -- list the 'primary factors' of the reduced polynomial
      redPrimaries : List SUP sae :=
        [(f.factor) **((f.exponent) :: NNI) for f in factorListSAE]
      -- lift these factors to elements of UP
      primaries : List UP :=
        [mapUnivariate(lift,ff)$IBPTOOLS(R,UP,SUP UP,sae) for ff in redPrimaries]
      -- lift the factors to factors modulo a suitable power of 'prime'
      deg := (1 + order(redDisc,prime) * degree(prime)) :: PI
      henselInfo := HenselLift(p,primaries,prime,deg)
      henselFactors := henselInfo.plist
      -- compute integral bases for the factors
      factorBases : List Result := empty(); degPrime := degree prime
      for pp in primes for k in expons for qq in henselFactors repeat
        base :=
          degPp := degree pp
          degPp > 1 and gcd(degPp,degPrime) = 1 =>
            compLocalBasisOverExt(qq,prime,pp,k)
          compLocalBasis(qq,prime)
        factorBases := concat(base,factorBases)
      factorBases := reverse! factorBases
      ib := chineseRemainder(henselFactors,factorBases,rank()$F)$IBACHIN(K,R,UP)
      index := diagonalProduct(ib.basisInv)
      [ib.basis,ib.basisDen,ib.basisInv,disc quo (index * index)]

    localIntegralBasis prime ==
      p := definingPolynomial()$F; disc := discriminant p
      --disc := determinant traceMatrix()$F
      redDisc := reducedDiscriminant p
      ib := padicLocalIntegralBasis(p,disc,redDisc,prime)
      [ib.basis,ib.basisDen,ib.basisInv]

    listSquaredFactors: R -> List R
    listSquaredFactors px ==
      -- returns a list of the factors of px which occur with
      -- exponent > 1
      ans : List R := empty()
      factored := factor(px)$DistinctDegreeFactorize(K,R)
      for f in factors(factored) repeat
        if f.exponent > 1 then ans := concat(f.factor,ans)
      ans

    integralBasis() ==
      p := definingPolynomial()$F; disc := discriminant p; n := rank()$F
      --traceMat := traceMatrix()$F; n := rank()$F
      --disc := determinant traceMat        -- discriminant of current order
      singList := listSquaredFactors disc -- singularities of relative Spec
      redDisc := reducedDiscriminant p
      runningRb := runningRbinv := scalarMatrix(n,1)$Mat
      -- runningRb    = basis matrix of current order
      -- runningRbinv = inverse basis matrix of current order
      -- these are wrt the original basis for F
      runningRbden : R := 1
      -- runningRbden = denominator for current basis matrix
      empty? singList => [runningRb, runningRbden, runningRbinv]
      for prime in singList repeat
        lb := padicLocalIntegralBasis(p,disc,redDisc,prime)
        rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
        disc := lb.discr
        mat := vertConcat(rbden * runningRb,runningRbden * rb)
        runningRbden := runningRbden * rbden
        runningRb := squareTop rowEchelon(mat,runningRbden)
        --runningRb := squareTop rowEch mat
        runningRbinv := UpTriBddDenomInv(runningRb,runningRbden)
      [runningRb, runningRbden, runningRbinv]

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<package IBPTOOLS IntegralBasisPolynomialTools>>
<<package IBACHIN ChineseRemainderToolsForIntegralBases>>
<<package PWFFINTB PAdicWildFunctionFieldIntegralBasis>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}