--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 category NTSCAT NormalizedTriangularSetCategory
++ Author: Marc Moreno Maza
++ Date Created: 10/07/1998
++ Date Last Updated: 12/12/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: essai Graphisme
++ AMS Classifications:
++ Keywords: polynomial, multivariate, ordered variables set
++ Description:
++ The category of normalized triangular sets. A triangular
++ set \spad{ts} is said normalized if for every algebraic
++ variable \spad{v} of \spad{ts} the polynomial \spad{select(ts,v)}
++ is normalized w.r.t. every polynomial in \spad{collectUnder(ts,v)}.
++ A polynomial \spad{p} is said normalized w.r.t. a non-constant 
++ polynomial \spad{q} if \spad{p} is constant or \spad{degree(p,mdeg(q)) = 0}
++ and \spad{init(p)} is normalized w.r.t. \spad{q}. One of the important
++ features of normalized triangular sets is that they are regular sets.\newline
++ References :
++  [1] D. LAZARD "A new method for solving algebraic systems of 
++      positive dimension" Discr. App. Math. 33:147-160,1991
++  [2] P. AUBRY, D. LAZARD and M. MORENO MAZA "On the Theories
++      of Triangular Sets" Journal of Symbol. Comp. (to appear)
++  [3] M. MORENO MAZA and R. RIOBOO "Computations of gcd over
++      algebraic towers of simple extensions" In proceedings of AAECC11
++      Paris, 1995.
++  [4] M. MORENO MAZA "Calculs de pgcd au-dessus des tours
++      d'extensions simples et resolution des systemes d'equations
++      algebriques" These, Universite P.etM. Curie, Paris, 1997.


NormalizedTriangularSetCategory(R:GcdDomain,E:OrderedAbelianMonoidSup,_
 V:OrderedSet,P:RecursivePolynomialCategory(R,E,V)):
         Category ==  RegularTriangularSetCategory(R,E,V,P) 

)abbrev package NORMPK NormalizationPackage
++ Author: Marc Moreno Maza
++ Date Created: 09/23/1998
++ Date Last Updated: 12/16/1998
++ Basic Functions:
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords:
++ Description: 
++ A package for computing normalized assocites of univariate polynomials
++ with coefficients in a tower of simple extensions of a field.\newline
++ References :
++  [1] D. LAZARD "A new method for solving algebraic systems of 
++      positive dimension" Discr. App. Math. 33:147-160,1991
++  [2] M. MORENO MAZA and R. RIOBOO "Computations of gcd over
++      algebraic towers of simple extensions" In proceedings of AAECC11
++      Paris, 1995.
++  [3] M. MORENO MAZA "Calculs de pgcd au-dessus des tours
++      d'extensions simples et resolution des systemes d'equations
++      algebriques" These, Universite P.etM. Curie, Paris, 1997.
++ Version: 1. 

NormalizationPackage(R,E,V,P,TS): Exports == Implementation where

  R : GcdDomain
  E : OrderedAbelianMonoidSup
  V : OrderedSet
  P : RecursivePolynomialCategory(R,E,V)
  TS : RegularTriangularSetCategory(R,E,V,P)
  N ==> NonNegativeInteger
  Z ==> Integer
  B ==> Boolean
  S ==> String
  K ==> Fraction R
  LP ==> List P
  PWT ==> Record(val : P, tower : TS)

  BWT ==> Record(val : Boolean, tower : TS)
  LpWT ==> Record(val : (List P), tower : TS)
  Split ==> List TS
  --KeyGcd ==> Record(arg1: P, arg2: P, arg3: TS, arg4: B)
  --EntryGcd ==> List PWT
  --HGcd ==> TabulatedComputationPackage(KeyGcd, EntryGcd)
  --KeyInvSet ==> Record(arg1: P, arg3: TS)
  --EntryInvSet ==> List TS
  --HInvSet ==> TabulatedComputationPackage(KeyInvSet, EntryInvSet)
  polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P)
  regsetgcdpack ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,TS)

  Exports ==  with

     recip: (P, TS) -> Record(num:P,den:P)
       ++ \axiom{recip(p,ts)} returns the inverse of \axiom{p} w.r.t \spad{ts}
       ++ assuming that \axiom{p} is invertible w.r.t \spad{ts}.
     normalizedAssociate: (P, TS) -> P
       ++ \axiom{normalizedAssociate(p,ts)} returns a normalized polynomial 
       ++ \axiom{n} w.r.t. \spad{ts} such that \axiom{n} and \axiom{p} are
       ++ associates w.r.t \spad{ts} and assuming that \axiom{p} is invertible 
       ++ w.r.t \spad{ts}.
     normalize: (P, TS) -> List PWT
       ++ \axiom{normalize(p,ts)} normalizes \axiom{p} w.r.t \spad{ts}.
     outputArgs: (S, S, P, TS) -> Void
       ++ \axiom{outputArgs(s1,s2,p,ts)} 
       ++ is an internal subroutine, exported only for developement.
     normInvertible?: (P, TS) -> List BWT
       ++ \axiom{normInvertible?(p,ts)} 
       ++ is an internal subroutine, exported only for developement.

  Implementation == add

     if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P)
     then

       normInvertible?(p:P, ts:TS): List BWT ==
         stoseInvertible?_sqfreg(p,ts)$regsetgcdpack

     else

       normInvertible?(p:P, ts:TS): List BWT ==
         stoseInvertible?_reg(p,ts)$regsetgcdpack

     if (R has RetractableTo(Integer)) and (V has ConvertibleTo(Symbol))
     then 

       outputArgs(s1:S, s2: S, p:P,ts:TS): Void ==
         if not empty? s1 then output(s1, p::OutputForm)$OutputPackage
         if not empty? s1 then output(s1,(convert(p)@String)::OutputForm)$OutputPackage
         output(" ")$OutputPackage
         if not empty? s2 then output(s2, ts::OutputForm)$OutputPackage       
         empty? s2 => void()
         output(s2,("[")::OutputForm)$OutputPackage
         lp: List P := members(ts)
         for q in lp repeat
            output((convert(q)@String)::OutputForm)$OutputPackage
         output("]")$OutputPackage
         output(" ")$OutputPackage

     else

       outputArgs(s1:S, s2: S, p:P,ts:TS): Void ==
         if not empty? s1 then output(s1, p::OutputForm)$OutputPackage
         output(" ")$OutputPackage
         if not empty? s2 then output(s2, ts::OutputForm)$OutputPackage       
         output(" ")$OutputPackage

     recip(p:P,ts:TS): Record(num:P, den:P) ==
     -- ASSUME p is invertible w.r.t. ts
     -- ASSUME mvar(p) is algebraic w.r.t. ts
       v := mvar(p)
       ts_v := select(ts,v)::P
       d : P
       n : P
       if mdeg(p) < mdeg(ts_v)
         then
           hesrg2 := halfExtendedSubResultantGcd2(ts_v,p)$P
           d := hesrg2.gcd
           n := hesrg2.coef2
         else
           hesrg1 := halfExtendedSubResultantGcd1(p,ts_v)$P
           d :=  hesrg1.gcd
           n := hesrg1.coef1
       g := gcd(n,d)
       (n, d) := ((n exquo g)::P, (d exquo g)::P)
       remn, remd: Record(rnum:R,polnum:P,den:R)
       remn := remainder(n,ts); remd := remainder(d,ts)
       cn := remn.rnum; pn := remn.polnum; dn := remn.den
       cd := remd.rnum; pd := remd.polnum; dp := remd.den
       k: K := (cn / cd) * (dp / dn)
       pn := removeZero(pn,ts)
       pd := removeZero(pd,ts)
       [numer(k) * pn, denom(k) * pd]$Record(num:P, den:P)

     normalizedAssociate(p:P,ts:TS): P ==
     -- ASSUME p is invertible or zero w.r.t. ts
       empty? ts => p
       zero?(p) => p
       ground?(p) => 1
       zero? initiallyReduce(init(p),ts) =>
         error "in normalizedAssociate$NORMPK: bad #1"
       vp := mvar(p)
       ip: P := p
       mp: P := 1
       tp: P := 0
       while not ground?(ip) repeat
         v := mvar(ip)
         if algebraic?(v,ts)
           then
             if v = vp
               then
                 ts_v := select(ts,v)::P
                 ip := lastSubResultant(ip,ts_v)$P
                 ip := remainder(ip,ts).polnum
                 -- ip := primitivePart stronglyReduce(ip,ts)
                 ip := primitivePart initiallyReduce(ip,ts)
               else
                 qr := recip(ip,ts)
                 ip := qr.den
                 tp := qr.num * tp
                 zero? ip =>
                     outputArgs("p = ", " ts = ",p,ts)
                     error "in normalizedAssociate$NORMPK: should never happen !"
           else
             tp := tail(ip) * mp + tp
             mp := mainMonomial(ip) * mp
             ip := init(ip)
       r := ip * mp + tp
       r := remainder(r,ts).polnum
       -- primitivePart stronglyReduce(r,ts)
       primitivePart initiallyReduce(r,ts)

     normalize(p: P, ts: TS): List PWT ==
       zero? p => [[p,ts]$PWT]
       ground? p => [[1,ts]$PWT]
       zero? initiallyReduce(init(p),ts) =>
         error "in normalize$NORMPK: init(#1) reduces to 0 w.r.t. #2"
       --output("Entering  normalize")$OutputPackage
       --outputArgs("p = ", " ts = ",p,ts)
       --output("Calling  normInvertible?")$OutputPackage
       lbwt: List BWT := normInvertible?(p,ts)
       --output("Result is: ")$OutputPackage
       --output(lbwt::OutputForm)$OutputPackage
       lpwt: List PWT := []
       for bwt in lbwt repeat
         us := bwt.tower
         q := remainder(p,us).polnum
         q := removeZero(q,us)
         bwt.val =>
           --output("Calling  normalizedAssociate")$OutputPackage
           --outputArgs("q = ", " us = ",q,us)
           lpwt := cons([normalizedAssociate(q,us)@P,us]$PWT, lpwt)
           --output("Leaving  normalizedAssociate")$OutputPackage
         zero? q => lpwt := cons([0$P,us]$PWT, lpwt)
         lpwt := concat(normalize(q,us)@(List PWT),lpwt)
       lpwt

