--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.
 
--% Polynomial composition and decomposition functions
--  If f = g o h then  g = leftFactor(f, h)  &  h = rightFactor(f, g)
--  SMW Dec 86

)abbrev package PCOMP PolynomialComposition
++ Description:
++ This package \undocumented
PolynomialComposition(UP: UnivariatePolynomialCategory(R), R: Ring): with
        compose: (UP, UP) -> UP
		++ compose(p,q) \undocumented
    == add
        compose(g, h) ==
            r: UP := 0
            while g ~= 0 repeat
                r := leadingCoefficient(g)*h**degree(g) + r
                g := reductum g
            r

)abbrev package PDECOMP PolynomialDecomposition
++ Description:
++ This package \undocumented
--  Ref: Kozen and Landau, Cornell University  TR 86-773
PolynomialDecomposition(UP, F): PDcat == PDdef where
    F:Field
    UP:UnivariatePolynomialCategory F
    NNI ==> NonNegativeInteger
    LR  ==> Record(left: UP, right: UP)
 
    PDcat == with
        decompose: UP -> List UP
		++ decompose(up) \undocumented
        decompose: (UP, NNI, NNI) -> Union(LR, "failed")
		++ decompose(up,m,n) \undocumented
        leftFactor: (UP, UP) -> Union(UP, "failed") 
		++ leftFactor(p,q) \undocumented
        rightFactorCandidate:  (UP, NNI) -> UP
		++ rightFactorCandidate(p,n) \undocumented
    PDdef == add
        leftFactor(f, h) ==
             g: UP := 0
             for i in 0.. while f ~= 0 repeat
                 fr := divide(f, h)
                 f := fr.quotient; r := fr.remainder
                 positive? degree r => return "failed"
                 g := g + r * monomial(1, i)
             g
 
        decompose(f, dg, dh) ==
            df := degree f
            dg*dh ~= df => "failed"
            h := rightFactorCandidate(f, dh)
            g := leftFactor(f, h)
            g case "failed" => "failed"
            [g::UP, h]
 
        decompose f ==
            df := degree f
            for dh in 2..df-1 | df rem dh = 0 repeat
                h := rightFactorCandidate(f, dh)
                g := leftFactor(f, h)
                g case UP => return
                    append(decompose(g::UP), decompose h)
            [f]
        rightFactorCandidate(f, dh) ==
            f  := f/leadingCoefficient f
            df := degree f
            dg := df quo dh
            h  := monomial(1, dh)
            for k in 1..dh repeat
                hdg:= h**dg
                c  := (coefficient(f,(df-k)::NNI)-coefficient(hdg,(df-k)::NNI))/(dg::F)
                h  := h + monomial(c, (dh-k)::NNI)
            h - monomial(coefficient(h, 0), 0) -- drop constant term

