--* From postmaster%watson.vnet.ibm.com@yktvmv.watson.ibm.com  Thu Jul 14 07:25:05 1994
--* Received: from yktvmv-ob.watson.ibm.com by asharp.watson.ibm.com (AIX 3.2/UCB 5.64/930311)
--*           id AA30200; Thu, 14 Jul 1994 07:25:05 -0400
--* Received: from watson.vnet.ibm.com by yktvmv.watson.ibm.com (IBM VM SMTP V2R3)
--*    with BSMTP id 0669; Thu, 14 Jul 94 07:25:08 EDT
--* Received: from YKTVMV by watson.vnet.ibm.com with "VAGENT.V1.0"
--*           id <A.BMT.NOTE.YKTVMV.3469.Jul.14.07:25:07.-0400>
--*           for asbugs@watson; Thu, 14 Jul 94 07:25:08 -0400
--* Received: from ICNUCEVX.CNUCE.CNR.IT by watson.ibm.com (IBM VM SMTP V2R3)
--*    with TCP; Thu, 14 Jul 94 07:25:01 EDT
--* Received: from cardano (cardano.dm.unipi.it) by icnucevx.cnuce.cnr.it (PMDF
--*  V4.2-11 #4369) id <01HEP517UTUOB7ASQY@icnucevx.cnuce.cnr.it>; Thu,
--*  14 Jul 1994 13:24:51 MET
--* Received: by cardano (4.1/SMI-4.1) id AA07754; Thu, 14 Jul 94 13:26:10 +0200
--* Date: Thu, 14 Jul 1994 13:26:10 +0200
--* From: bmt@posso.dm.unipi.it (Barry Trager)
--* Subject: [3] complains of syntax errors but cannot locate them
--*  [pdecomp.as][36.0]
--* To: asbugs@watson.ibm.com
--* Message-Id: <9407141126.AA07754@cardano>
--* Content-Transfer-Encoding: 7BIT

--@ Fixed  by:  SSD   Fri Jul 15 12:01:01 EDT 1994 
--@ Tested by:  none 
--@ Summary:    Not a bug. Error messages correctly points to an unbalanced close parenthesis. 

#pile
#include "axiom.as"
-- (c) Copyright International Business Machines Corporation 1986.
-- All rights reserved.
--
-- Address comments and questions to the
--   Computer Algebra Group, Mathematical Sciences Department
--   IBM Thomas J. Watson Research Center, Box 218
--   Yorktown Heights, New York  10598  USA

--% Polynomial composition and decomposition functions
--  If f = g o h then  g = leftFactor(f, h)  &  h = rightFactor(f, g)
--  SMW Dec 86

--% PolynomialComposition
--)abbrev package PCOMP PolynomialComposition
--)abbrev package PDECOMP PolynomialDecomposition

PolynomialComposition(UP: UnivariatePolynomialCategory(R), R: Ring): with
        compose: (UP, UP) -> UP

    == add
        compose(g:UP, h:UP):UP ==
            r: UP := 0
            while g ~= 0 repeat
                r := leadingCoefficient(g)*h**degree(g) + r
                g := reductum g
            r


--  Ref: Kozen and Landau, Cornell University  TR 86-773

--% PolynomialDecomposition

PolynomialDecomposition(UP, F): PDcat == PDdef where
    default F:Field
    default UP:UnivariatePolynomialCategory F
    NNI ==> NonNegativeInteger
    LR  ==> Record(left: UP, right: UP)

    PDcat ==> with
        decompose: UP -> List UP
        decompose: (UP, NNI, NNI) -> Union(LR, "failed")
        leftFactor: (UP, UP) -> Union(UP, "failed")
        rightFactorCandidate:  (UP, NNI) -> UP
    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
                 degree r > 0 => 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) ==
            df := degree f
            dg := df quo dh
            h  := monomial(1, dh:NNI)
            for k in 1..dh repeat
                hdg:= h**(dg:NNI)
                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