Initial population.

1 files changed, 571 insertions, 0 deletions

diff --git a/src/algebra/padiclib.spad.pamphlet b/src/algebra/padiclib.spad.pamphlet
new file mode 100644
index 00000000..9b3e79cc
--- /dev/null
+++ b/src/algebra/padiclib.spad.pamphlet
@@ -0,0 +1,571 @@
+\documentclass{article}
+\usepackage{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) =>
+      not (ee.generator = 1) =>
+        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}