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

)abbrev package GENEEZ GenExEuclid
++ Author : P.Gianni.
++ January 1990
++ The equation \spad{Af+Bg=h} and its generalization to n polynomials
++ is solved for solutions over the R, euclidean domain.
++ A table containing the solutions of \spad{Af+Bg=x**k} is used.
++ The operations are performed modulus a prime which are in principle big enough,
++ but the solutions are tested and, in case of failure, a hensel
++ lifting process is used to get to the right solutions.
++ It will be used in the factorization of multivariate polynomials
++ over finite field, with \spad{R=F[x]}.

GenExEuclid(R,BP) : C == T
 where
  R   :   EuclideanDomain
  PI  ==> PositiveInteger
  NNI ==> NonNegativeInteger
  BP  :   UnivariatePolynomialCategory R
  L   ==> List

  C == with
       reduction: (BP,R) -> BP
         ++ reduction(p,prime) reduces the polynomial p modulo prime of R.
         ++ Note: this function is exported only because it's conditional.
       compBound: (BP,L BP) -> NNI
         ++ compBound(p,lp)
         ++ computes a bound for the coefficients of the solution
         ++ polynomials.
         ++ Given a polynomial right hand side p, and a list lp of left hand side polynomials.
         ++ Exported because it depends on the valuation.
       tablePow :    (NNI,R,L BP)     -> Union(Vector(L BP),"failed")
         ++ tablePow(maxdeg,prime,lpol) constructs the table with the
         ++ coefficients of the Extended Euclidean Algorithm for lpol.
         ++ Here the right side is \spad{x**k}, for k less or equal to maxdeg.
         ++ The operation returns "failed" when the elements are not coprime modulo prime.
       solveid  : (BP,R,Vector L BP) -> Union(L BP,"failed")
         ++ solveid(h,table) computes the coefficients of the
         ++ extended euclidean algorithm for a list of polynomials
         ++ whose tablePow is table and with right side h.

       testModulus : (R, L BP)    -> Boolean
         ++ testModulus(p,lp) returns true if the the prime p
         ++ is valid for the list of polynomials lp, i.e. preserves
         ++ the degree and they remain relatively prime.

  T == add
    if R has multiplicativeValuation then
      compBound(m:BP,listpolys:L BP) : NNI ==
        ldeg:=[degree f for f in listpolys]
        n:NNI:= (+/[df for df in ldeg])
        normlist:=[ +/[euclideanSize(u)**2 for u in coefficients f]
                         for f in listpolys]
        nm:= +/[euclideanSize(u)**2 for u in coefficients m]
	normprod := */[g**((n-df)::NNI) for g in normlist for df in ldeg]
        2*(approxSqrt(normprod * nm)$IntegerRoots(Integer))::NNI
    else if R has additiveValuation then
      -- a fairly crude Hadamard-style bound for the solution
      -- based on regarding the problem as a system of linear equations.
      compBound(m:BP,listpolys:L BP) : NNI ==
        "max"/[euclideanSize u for u in coefficients m] +
          +/["max"/[euclideanSize u for u in coefficients p]
             for p in listpolys]
    else
      compBound(m:BP,listpolys:L BP) : NNI ==
        error "attempt to use compBound without a well-understood valuation"
    if R has IntegerNumberSystem then
      reduction(u:BP,p:R):BP ==
        p = 0 => u
        map(symmetricRemainder(#1,p),u)
    else reduction(u:BP,p:R):BP ==
        p = 0 => u
        map(#1 rem p,u)

    merge(p:R,q:R):Union(R,"failed") ==
         p = q => p
         p = 0 => q
         q = 0 => p
         "failed"

    modInverse(c:R,p:R):R ==
        (extendedEuclidean(c,p,1)::Record(coef1:R,coef2:R)).coef1

    exactquo(u:BP,v:BP,p:R):Union(BP,"failed") ==
        invlcv:=modInverse(leadingCoefficient v,p)
        r:=monicDivide(u,reduction(invlcv*v,p))
        not zero? reduction(r.remainder,p) => "failed"
        reduction(invlcv*r.quotient,p)

    FP:=EuclideanModularRing(R,BP,R,reduction,merge,exactquo)

    --make table global variable!
    table:Vector L BP
    import GeneralHenselPackage(R,BP)

                       --local  functions
    makeProducts :    L BP   -> L BP
    liftSol: (L BP,BP,R,R,Vector L BP,BP,NNI) -> Union(L BP,"failed")

    reduceList(lp:L BP,lmod:R): L FP ==[reduce(ff,lmod) for ff in lp]

    coerceLFP(lf:L FP):L BP == [fm::BP for fm in lf]

    liftSol(oldsol:L BP,err:BP,lmod:R,lmodk:R,
           table:Vector L BP,m:BP,bound:NNI):Union(L BP,"failed") ==
      euclideanSize(lmodk) > bound => "failed"
      d:=degree err
      ftab:Vector L FP :=
        map(reduceList(#1,lmod),table)$VectorFunctions2(List BP,List FP)
      sln:L FP:=[0$FP for xx in ftab.1 ]
      for i in 0 .. d |not zero?(cc:=coefficient(err,i)) repeat
        sln:=[slp+reduce(cc::BP,lmod)*pp
              for pp in ftab.(i+1) for slp in sln]
      nsol:=[f-lmodk*reduction(g::BP,lmod) for f in oldsol for g in sln]
      lmodk1:=lmod*lmodk
      nsol:=[reduction(slp,lmodk1) for slp in nsol]
      lpolys:L BP:=table.(#table)
      (fs:=+/[f*g for f in lpolys for g in nsol]) = m => nsol
      a:BP:=((fs-m) exquo lmodk1)::BP
      liftSol(nsol,a,lmod,lmodk1,table,m,bound)

    makeProducts(listPol:L BP):L BP ==
      #listPol < 2 => listPol
      #listPol = 2 => reverse listPol
      f:= first listPol
      ll := rest listPol
      [*/ll,:[f*g for g in makeProducts ll]]

    testModulus(pmod, listPol) ==
         redListPol := reduceList(listPol, pmod)
         for pol in listPol for rpol in redListPol repeat
              degree(pol) ~= degree(rpol::BP) => return false
         while not empty? redListPol repeat
             rpol := first redListPol
             redListPol := rest redListPol
             for rpol2 in redListPol repeat
                not one? gcd(rpol, rpol2) => return false
         true

    if R has Field then
      tablePow(mdeg:NNI,pmod:R,listPol:L BP) ==
        multiE:=multiEuclidean(listPol,1$BP)
        multiE case "failed" => "failed"
        ptable:Vector L BP :=new(mdeg+1,[])
        ptable.1:=multiE
        x:BP:=monomial(1,1)
        for i in 2..mdeg repeat ptable.i:=
            [tpol*x rem fpol for tpol in ptable.(i-1) for fpol in listPol]
        ptable.(mdeg+1):=makeProducts listPol
        ptable

      solveid(m:BP,pmod:R,table:Vector L BP) : Union(L BP,"failed") ==
            -- Actually, there's no possibility of failure
        d:=degree m
        sln:L BP:=[0$BP for xx in table.1]
        for i in 0 .. d | not zero? coefficient(m,i) repeat
          sln:=[slp+coefficient(m,i)*pp
                for pp in table.(i+1) for slp in sln]
        sln

    else

      tablePow(mdeg:NNI,pmod:R,listPol:L BP) ==
        listP:L FP:= [reduce(pol,pmod) for pol in listPol]
        multiE:=multiEuclidean(listP,1$FP)
        multiE case "failed" => "failed"
        ftable:Vector L FP :=new(mdeg+1,[])
        fl:L FP:= [ff::FP for ff in multiE]
        ftable.1:=[fpol for fpol in fl]
        x:FP:=reduce(monomial(1,1),pmod)
        for i in 2..mdeg repeat ftable.i:=
            [tpol*x rem fpol for tpol in ftable.(i-1) for fpol in listP]
        ptable:= map(coerceLFP,ftable)$VectorFunctions2(List FP,List BP)
        ptable.(mdeg+1):=makeProducts listPol
        ptable

      solveid(m:BP,pmod:R,table:Vector L BP) : Union(L BP,"failed") ==
        d:=degree m
        ftab:Vector L FP:=
          map(reduceList(#1,pmod),table)$VectorFunctions2(List BP,List FP)
        lpolys:L BP:=table.(#table)
        sln:L FP:=[0$FP for xx in ftab.1]
        for i in 0 .. d | not zero? coefficient(m,i) repeat
          sln:=[slp+reduce(coefficient(m,i)::BP,pmod)*pp
                for pp in ftab.(i+1) for slp in sln]
        soln:=[slp::BP for slp in sln]
        (fs:=+/[f*g for f in lpolys for g in soln]) = m=> soln
        -- Compute bound
        bound:=compBound(m,lpolys)
        a:BP:=((fs-m) exquo pmod)::BP
        liftSol(soln,a,pmod,pmod,table,m,bound)

