--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.
 
import Field
import UnivariatePolynomialCategory
import CharacteristicZero
import MonogenicAlgebra
import Factored
)abbrev package IALGFACT InnerAlgFactor
++ Factorisation in a simple algebraic extension
++ Author: Patrizia Gianni
++ Date Created: ???
++ Date Last Updated: 20 Jul 1988
++ Description:
++ Factorization of univariate polynomials with coefficients in an
++ algebraic extension of a field over which we can factor UP's;
++ Keywords: factorization, algebraic extension, univariate polynomial

InnerAlgFactor(F, UP, AlExt, AlPol): Exports == Implementation where
  F    : Field
  UP   : UnivariatePolynomialCategory F
  AlPol: UnivariatePolynomialCategory AlExt
  AlExt : Join(Field, CharacteristicZero, MonogenicAlgebra(F,UP))
  NUP   ==> SparseUnivariatePolynomial UP
  N     ==> NonNegativeInteger
  Z     ==> Integer
  FR    ==> Factored UP
  UPCF2 ==> UnivariatePolynomialCategoryFunctions2


  Exports ==> with
    factor: (AlPol, UP -> FR)  ->  Factored AlPol
      ++ factor(p, f) returns a prime factorisation of p;
      ++ f is a factorisation map for elements of UP;
 
  Implementation ==> add
    pnorm        : AlPol -> UP
    convrt       : AlPol -> NUP
    change       : UP    -> AlPol
    perturbfactor: (AlPol, Z, UP -> FR) -> List AlPol
    irrfactor    : (AlPol, Z, UP -> FR) -> List AlPol
 
 
    perturbfactor(f, k, fact) ==
      pol   := monomial(1$AlExt,1)-
               monomial(reduce monomial(k::F,1)$UP ,0)
      newf  := elt(f, pol)
      lsols := irrfactor(newf, k, fact)
      pol   := monomial(1, 1) +
               monomial(reduce monomial(k::F,1)$UP,0)
      [elt(pp, pol) for pp in lsols]
 
    ---  factorize the square-free parts of f  ---
    irrfactor(f, k, fact) ==
      degree(f) =$N 1 => [f]
      newf := f
      nn   := pnorm f
      --newval:RN:=1
      --pert:=false
      --if ^ SqFr? nn then
      --  pert:=true
      --  newterm:=perturb(f)
      --  newf:=newterm.ppol
      --  newval:=newterm.pval
      --  nn:=newterm.nnorm
      listfact := factors fact nn
      #listfact =$N 1 =>
        first(listfact).exponent =$Z 1 => [f]
        perturbfactor(f, k + 1, fact)
      listerm:List(AlPol):= []
      for pelt in listfact repeat
        g    := gcd(change(pelt.factor), newf)
        newf := (newf exquo g)::AlPol
        listerm :=
          pelt.exponent =$Z 1 => cons(g, listerm)
          append(perturbfactor(g, k + 1, fact), listerm)
      listerm
 
    factor(f, fact) ==
      sqf := squareFree f
      unit(sqf) * _*/[_*/[primeFactor(pol, sqterm.exponent)
                          for pol in irrfactor(sqterm.factor, 0, fact)]
                                            for sqterm in factors sqf]
 
    p := definingPolynomial()$AlExt
    newp := map(#1::UP, p)$UPCF2(F, UP, UP, NUP)
 
    pnorm  q == resultant(convrt q, newp)
    change q == map(coerce, q)$UPCF2(F,UP,AlExt,AlPol)
 
    convrt q ==
      swap(map(lift, q)$UPCF2(AlExt, AlPol,
           UP, NUP))$CommuteUnivariatePolynomialCategory(F, UP, NUP)

import UnivariatePolynomialCategory
import CharacteristicZero
import Field
import MonogenicAlgebra
import Fraction
import Integer
import Factored
)abbrev package SAEFACT SimpleAlgebraicExtensionAlgFactor
++ Factorisation in a simple algebraic extension;
++ Author: Patrizia Gianni
++ Date Created: ???
++ Date Last Updated: ???
++ Description:
++ Factorization of univariate polynomials with coefficients in an
++ algebraic extension of the rational numbers (\spadtype{Fraction Integer}).
++ Keywords: factorization, algebraic extension, univariate polynomial
 
SimpleAlgebraicExtensionAlgFactor(UP,SAE,UPA):Exports==Implementation where
  UP : UnivariatePolynomialCategory Fraction Integer
  SAE : Join(Field, CharacteristicZero,
                         MonogenicAlgebra(Fraction Integer, UP))
  UPA: UnivariatePolynomialCategory SAE
 
  Exports ==> with
    factor: UPA -> Factored UPA
      ++ factor(p) returns a prime factorisation of p.
 
  Implementation ==> add
    factor q ==
      factor(q, factor$RationalFactorize(UP)
                       )$InnerAlgFactor(Fraction Integer, UP, SAE, UPA)

import UnivariatePolynomialCategory
import Factored
import Polynomial
import Integer
)abbrev package RFFACT RationalFunctionFactor
++ Factorisation in UP FRAC POLY INT
++ Author: Patrizia Gianni
++ Date Created: ???
++ Date Last Updated: ???
++ Description:
++ Factorization of univariate polynomials with coefficients which
++ are rational functions with integer coefficients.

RationalFunctionFactor(UP): Exports == Implementation where
  UP: UnivariatePolynomialCategory Fraction Polynomial Integer
 
  SE ==> Symbol
  P  ==> Polynomial Integer
  RF ==> Fraction P
  UPCF2 ==> UnivariatePolynomialCategoryFunctions2
 
  Exports ==> with
    factor: UP -> Factored UP
      ++ factor(p) returns a prime factorisation of p.
 
  Implementation ==> add
    likuniv: (P, SE, P) -> UP
 
    dummy := new()$SE
 
    likuniv(p, x, d) ==
      map(#1 / d, univariate(p, x))$UPCF2(P,SparseUnivariatePolynomial P,
                                          RF, UP)
 
    factor p ==
      d  := denom(q := elt(p,dummy::P :: RF))
      map(likuniv(#1,dummy,d),
          factor(numer q)$MultivariateFactorize(SE,
               IndexedExponents SE,Integer,P))$FactoredFunctions2(P, UP)

import UnivariatePolynomialCategory
import Field
import CharacteristicZero
import Polynomial
import Fraction
import Integer
)abbrev package SAERFFC SAERationalFunctionAlgFactor
++ Factorisation in UP SAE FRAC POLY INT
++ Author: Patrizia Gianni
++ Date Created: ???
++ Date Last Updated: ???
++ Description:
++ Factorization of univariate polynomials with coefficients in an
++ algebraic extension of \spadtype{Fraction Polynomial Integer}.
++ Keywords: factorization, algebraic extension, univariate polynomial
 
SAERationalFunctionAlgFactor(UP, SAE, UPA): Exports == Implementation where
  UP : UnivariatePolynomialCategory Fraction Polynomial Integer
  SAE : Join(Field, CharacteristicZero,
                      MonogenicAlgebra(Fraction Polynomial Integer, UP))
  UPA: UnivariatePolynomialCategory SAE
 
  Exports ==> with
    factor: UPA -> Factored UPA
      ++ factor(p) returns a prime factorisation of p.
 
  Implementation ==> add
    factor q ==
      factor(q, factor$RationalFunctionFactor(UP)
              )$InnerAlgFactor(Fraction Polynomial Integer, UP, SAE, UPA)

import UnivariatePolynomialCategory
import AlgebraicNumber
import Boolean
)abbrev package ALGFACT AlgFactor
++ Factorization of UP AN;
++ Author: Manuel Bronstein
++ Date Created: ???
++ Date Last Updated: ???
++ Description:
++ Factorization of univariate polynomials with coefficients in 
++ \spadtype{AlgebraicNumber}.
 
AlgFactor(UP): Exports == Implementation where
  UP: UnivariatePolynomialCategory AlgebraicNumber
 
  N   ==> NonNegativeInteger
  Z   ==> Integer
  Q   ==> Fraction Integer
  AN  ==> AlgebraicNumber
  K   ==> Kernel AN
  UPQ ==> SparseUnivariatePolynomial Q
  SUP ==> SparseUnivariatePolynomial AN
  FR  ==> Factored UP
 
  Exports ==> with
    factor: (UP, List AN) -> FR
      ++ factor(p, [a1,...,an]) returns a prime factorisation of p
      ++ over the field generated by its coefficients and a1,...,an.
    factor: UP            -> FR
      ++ factor(p) returns a prime factorisation of p
      ++ over the field generated by its coefficients.
    split : UP            -> FR
      ++ split(p) returns a prime factorisation of p
      ++ over its splitting field.
    doublyTransitive?: UP -> Boolean
      ++ doublyTransitive?(p) is true if p is irreducible over
      ++ over the field K generated by its coefficients, and
      ++ if \spad{p(X) / (X - a)} is irreducible over 
      ++ \spad{K(a)} where \spad{p(a) = 0}.
 
  Implementation ==> add
    import PolynomialCategoryQuotientFunctions(IndexedExponents K,
                           K, Z, SparseMultivariatePolynomial(Z, K), AN)

    UPCF2 ==> UnivariatePolynomialCategoryFunctions2

    fact    : (UP,  List K) -> FR
    ifactor : (SUP, List K) -> Factored SUP
    extend  : (UP, Z) -> FR
    allk    : List AN -> List K
    downpoly: UP  -> UPQ
    liftpoly: UPQ -> UP
    irred?  : UP  -> Boolean
 
    allk l       == removeDuplicates concat [kernels x for x in l]
    liftpoly p   == map(#1::AN,  p)$UPCF2(Q, UPQ, AN, UP)
    downpoly p   == map(retract(#1)@Q, p)$UPCF2(AN, UP ,Q, UPQ)
    ifactor(p,l) == (fact(p pretend UP, l)) pretend Factored(SUP)
    factor p     == fact(p, allk coefficients p)
 
    factor(p, l) ==
      fact(p, allk removeDuplicates concat(l, coefficients p))
 
    split p ==
      fp := factor p
      unit(fp) *
            _*/[extend(fc.factor, fc.exponent) for fc in factors fp]
 
    extend(p, n) ==
      one? degree p => primeFactor(p, n)
      q := monomial(1, 1)$UP - zeroOf(p pretend SUP)::UP
      primeFactor(q, n) * split((p exquo q)::UP) ** (n::N)
 
    doublyTransitive? p ==
      irred? p and irred?((p exquo
        (monomial(1, 1)$UP - zeroOf(p pretend SUP)::UP))::UP)
 
    irred? p ==
      fp := factor p
      one? numberOfFactors fp and one? nthExponent(fp, 1)
 
    fact(p, l) ==
      one? degree p => primeFactor(p, 1)
      empty? l =>
        dr := factor(downpoly p)$RationalFactorize(UPQ)
        (liftpoly unit dr) *
          _*/[primeFactor(liftpoly dc.factor,dc.exponent)
            for dc in factors dr]
      q   := minPoly(alpha := "max"/l)$AN
      newl  := remove(alpha = #1, l)
      sae := SimpleAlgebraicExtension(AN, SUP, q)
      ups := SparseUnivariatePolynomial sae
      fr  := factor(map(reduce univariate(#1, alpha, q),
                     p)$UPCF2(AN, UP, sae, ups),
                      ifactor(#1, newl))$InnerAlgFactor(AN, SUP, sae, ups)
      newalpha := alpha::AN
      map((lift(#1)$sae) newalpha, unit fr)$UPCF2(sae, ups, AN, UP) *
            _*/[primeFactor(map((lift(#1)$sae) newalpha,
                      fc.factor)$UPCF2(sae, ups, AN, UP),
                                 fc.exponent) for fc in factors fr]

