\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra updecomp.spad}
\author{Frederic Lehobey}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package UPDECOMP UnivariatePolynomialDecompositionPackage}
<<package UPDECOMP UnivariatePolynomialDecompositionPackage>>=
)abbrev package UPDECOMP UnivariatePolynomialDecompositionPackage
++ Author: Frederic Lehobey
++ Date Created: 17 June 1996
++ Date Last Updated: 4 June 1997
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keyword:
++ Exemples:
++ References:
++ [1] Peter Henrici, Automatic Computations with Power Series,
++ Journal of the Association for Computing Machinery, Volume 3, No. 1,
++ January 1956, 10-15
++ [2] Dexter Kozen and Susan Landau, Polynomial Decomposition
++ Algorithms, Journal of Symbolic Computation (1989) 7, 445-456
-- Decomposition would be speeded up (O(n log n) instead of O(n^2)) by
-- implementing the algorithm described in [3] based on [4] and [5]. 
++ [3] Joachim von zur Gathen, Functional Decomposition Polynomials:
++ the Tame Case, Journal of Symbolic Computation (1990) 9, 281-299
++ [4] R. P. Brent and H. T. Kung, Fast Algorithms for Manipulating
++ Formal Power Series, Journal of the Association for Computing
++ Machinery, Vol. 25, No. 4, October 1978, 581-595
++ [5] R. P. Brent, Multiple-Precision Zero-Finding Methods and the
++ Complexity of Elementary Function Evaluation, Analytic
++ Computational Complexity, J. F. Traub, Ed., Academic Press,
++ New York 1975, 151-176 
++ Description: UnivariatePolynomialDecompositionPackage implements
++ functional decomposition of univariate polynomial with coefficients
++ in an \spad{IntegralDomain} of \spad{CharacteristicZero}.
UnivariatePolynomialDecompositionPackage(R,UP): Exports == Implementation where
  R : Join(IntegralDomain,CharacteristicZero)
  UP : UnivariatePolynomialCategory(R)
  N ==> NonNegativeInteger
  LR ==> Record(left: UP, right: UP)
  QR ==> Record(quotient: UP, remainder: UP)


  Exports ==> with

    monicRightFactorIfCan: (UP,N) -> Union(UP,"failed")
      ++ monicRightFactorIfCan(f,d) returns a candidate to be the
      ++ monic right factor (h in f = g o h) of degree d of a
      ++ functional decomposition of the polynomial f or
      ++ \spad{"failed"} if no such candidate.
    rightFactorIfCan: (UP,N,R) -> Union(UP,"failed")
      ++ rightFactorIfCan(f,d,c) returns a candidate to be the
      ++ right factor (h in f = g o h) of degree d with leading
      ++ coefficient c of a functional decomposition of the
      ++ polynomial f or \spad{"failed"} if no such candidate. 
    leftFactorIfCan: (UP,UP) -> Union(UP,"failed")
      ++ leftFactorIfCan(f,h) returns the left factor (g in f = g o h)
      ++ of the functional decomposition of the polynomial f with
      ++ given h or \spad{"failed"} if g does not exist. 
    monicDecomposeIfCan: UP -> Union(LR,"failed")
      ++ monicDecomposeIfCan(f) returns a functional decomposition
      ++ of the monic polynomial f of "failed" if it has not found any.
    monicCompleteDecompose: UP -> List UP
      ++ monicCompleteDecompose(f) returns a list of factors of f for
      ++ the functional decomposition ([ f1, ..., fn ] means 
      ++ f = f1 o ... o fn).

  Implementation ==> add

    rightFactorIfCan(p,dq,lcq) ==
      dp := degree p
      zero? lcq =>
       error "rightFactorIfCan: leading coefficient may not be zero"
      (zero? dp) or (zero? dq) => "failed"
      nc := dp exquo dq
      nc case "failed" => "failed"
      n := nc::N
      s := subtractIfCan(dq,1)::N
      lcp := leadingCoefficient p
      q: UP := monomial(lcq,dq)
      k: N 
      for k in 1..s repeat
        c: R := 0
        i: N
        for i in 0..subtractIfCan(k,1)::N repeat
         c := c+(k::R-(n::R+1)*(i::R))*
          coefficient(q,subtractIfCan(dq,i)::N)*
           coefficient(p,subtractIfCan(dp+i,k)::N)
        cquo := c exquo ((k*n)::R*lcp)
        cquo case "failed" => return "failed"
        q := q+monomial(cquo::R,subtractIfCan(dq,k)::N)
      q

    monicRightFactorIfCan(p,dq) == rightFactorIfCan(p,dq,1$R)

    import UnivariatePolynomialDivisionPackage(R,UP)

    leftFactorIfCan(f,h) ==
      g: UP := 0
      zero? degree h => "failed"
      for i in 0.. while not zero? f repeat
        qrf := divideIfCan(f,h)
        qrf case "failed" => return "failed"
        qr := qrf :: QR
        r := qr.remainder
        not ground? r => return "failed"
        g := g+monomial(ground(r),i)
        f := qr.quotient
      g

    monicDecomposeIfCan f ==
      df := degree f
      zero? df => "failed"  
      for dh in 2..subtractIfCan(df,1)::N | zero?(df rem dh) repeat
        h := monicRightFactorIfCan(f,dh)
        h case UP =>
         g := leftFactorIfCan(f,h::UP)
         g case UP => return [g::UP,h::UP]
      "failed"

    monicCompleteDecompose f ==
      cf := monicDecomposeIfCan f
      cf case "failed" => [ f ]
      lr := cf :: LR
      append(monicCompleteDecompose lr.left,[lr.right])

@
\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 UPDECOMP UnivariatePolynomialDecompositionPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}