path: root/src/algebra/moddfact.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra moddfact.spad}
\author{Barry Trager, James Davenport}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package MDDFACT ModularDistinctDegreeFactorizer}
<<package MDDFACT ModularDistinctDegreeFactorizer>>=
)abbrev package MDDFACT ModularDistinctDegreeFactorizer
++ Author: Barry Trager
++ Date Created:
++ Date Last Updated: 20.9.95 (JHD)
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package supports factorization and gcds
++ of univariate polynomials over the integers modulo different
++ primes. The inputs are given as polynomials over the integers
++ with the prime passed explicitly as an extra argument.

ModularDistinctDegreeFactorizer(U):C == T where
  U : UnivariatePolynomialCategory(Integer)
  I ==> Integer
  NNI ==> NonNegativeInteger
  PI ==> PositiveInteger
  V ==> Vector
  L ==> List
  DDRecord ==> Record(factor:EMR,degree:I)
  UDDRecord ==> Record(factor:U,degree:I)
  DDList ==> L DDRecord
  UDDList ==> L UDDRecord


  C == with
    gcd:(U,U,I) -> U
      ++ gcd(f1,f2,p) computes the gcd of the univariate polynomials
      ++ f1 and f2 modulo the integer prime p.
    linears: (U,I) -> U
      ++ linears(f,p) returns the product of all the linear factors
      ++ of f modulo p. Potentially incorrect result if f is not
      ++ square-free modulo p.
    factor:(U,I) -> L U
      ++ factor(f1,p) returns the list of factors of the univariate
      ++ polynomial f1 modulo the integer prime p.
      ++ Error: if f1 is not square-free modulo p.
    ddFact:(U,I) -> UDDList
      ++ ddFact(f,p) computes a distinct degree factorization of the
      ++ polynomial f modulo the prime p, i.e. such that each factor
      ++ is a product of irreducibles of the same degrees. The input
      ++ polynomial f is assumed to be square-free modulo p.
    separateFactors:(UDDList,I) -> L U
      ++ separateFactors(ddl, p) refines the distinct degree factorization
      ++ produced by \spadfunFrom{ddFact}{ModularDistinctDegreeFactorizer}
      ++ to give a complete list of factors.
    exptMod:(U,I,U,I) -> U
      ++ exptMod(f,n,g,p) raises the univariate polynomial f to the nth
      ++ power modulo the polynomial g and the prime p.

  T == add
    reduction(u:U,p:I):U ==
       zero? p => u
       map(positiveRemainder(#1,p),u)
    merge(p:I,q:I):Union(I,"failed") ==
       p = q => p
       p = 0 => q
       q = 0 => p
       "failed"
    modInverse(c:I,p:I):I ==
        (extendedEuclidean(c,p,1)::Record(coef1:I,coef2:I)).coef1
    exactquo(u:U,v:U,p:I):Union(U,"failed") ==
        invlcv:=modInverse(leadingCoefficient v,p)
        r:=monicDivide(u,reduction(invlcv*v,p))
        not zero? reduction(r.remainder,p) => "failed"
        reduction(invlcv*r.quotient,p)
    EMR := EuclideanModularRing(Integer,U,Integer,
                                reduction,merge,exactquo)

    probSplit2:(EMR,EMR,I) -> Union(List EMR,"failed")
    trace:(EMR,I,EMR) -> EMR
    ddfactor:EMR -> L EMR
    ddfact:EMR -> DDList
    sepFact1:DDRecord -> L EMR
    sepfact:DDList -> L EMR
    probSplit:(EMR,EMR,I) -> Union(L EMR,"failed")
    makeMonic:EMR -> EMR
    exptmod:(EMR,I,EMR) -> EMR

    lc(u:EMR):I == leadingCoefficient(u::U)
    degree(u:EMR):I == degree(u::U)
    makeMonic(u) == modInverse(lc(u),modulus(u)) * u

    i:I

    exptmod(u1,i,u2) ==
      negative? i => error("negative exponentiation not allowed for exptMod")
      ans:= 1$EMR
      while positive? i repeat
        if odd?(i) then ans:= (ans * u1) rem u2
        i:= i quo 2
        u1:= (u1 * u1) rem u2
      ans

    exptMod(a,i,b,q) ==
      ans:= exptmod(reduce(a,q),i,reduce(b,q))
      ans::U

    ddfactor(u) ==
      if not one?(c:= lc(u)) then u:= makeMonic(u)
      ans:= sepfact(ddfact(u))
      cons(c::EMR,[makeMonic(f) for f in ans | positive? degree(f)])

    gcd(u,v,q) == gcd(reduce(u,q),reduce(v,q))::U

    factor(u,q) ==
      v:= reduce(u,q)
      dv:= reduce(differentiate(u),q)
      positive? degree gcd(v,dv) =>
        error("Modular factor: polynomial must be squarefree")
      ans:= ddfactor v
      [f::U for f in ans]

    ddfact(u) ==
      p:=modulus u
      w:= reduce(monomial(1,1)$U,p)
      m:= w
      d:I:= 1
      if not one?(c:= lc(u)) then u:= makeMonic u
      ans:DDList:= []
      repeat
        w:= exptmod(w,p,u)
        g:= gcd(w - m,u)
        if positive? degree g then
          g:= makeMonic(g)
          ans:= [[g,d],:ans]
          u:= (u quo g)
        degree(u) = 0 => return [[c::EMR,0$I],:ans]
        d:= d+1
        d > (degree(u):I quo 2) =>
               return [[c::EMR,0$I],[u,degree(u)],:ans]

    ddFact(u,q) ==
      ans:= ddfact(reduce(u,q))
      [[(dd.factor)::U,dd.degree]$UDDRecord for dd in ans]$UDDList

    linears(u,q) == 
       uu:=reduce(u,q)
       m:= reduce(monomial(1,1)$U,q)
       gcd(exptmod(m,q,uu)-m,uu)::U

    sepfact(factList) ==
      "append"/[sepFact1(f) for f in factList]

    separateFactors(uddList,q) ==
      ans:= sepfact [[reduce(udd.factor,q),udd.degree]$DDRecord for
        udd in uddList]$DDList
      [f::U for f in ans]

    decode(s:Integer, p:Integer, x:U):U ==
      s<p => s::U
      qr := divide(s,p)
      qr.remainder :: U + x*decode(qr.quotient, p, x)

    sepFact1(f) ==
      u:= f.factor
      p:=modulus u
      (d := f.degree) = 0 => [u]
      if not one?(c:= lc(u)) then u:= makeMonic(u)
      d = (du := degree(u)) => [u]
      ans:L EMR:= []
      x:U:= monomial(1,1)
      -- for small primes find linear factors by exhaustion
      d=1 and p < 1000  =>
        for i: local in 0.. while positive? du repeat
          if u(i::U) = 0 then
            ans := cons(reduce(x-(i::U),p),ans)
            du := du-1
        ans 
      y:= x
      s:I:= 0
      ss:I := 1
      stack:L EMR:= [u]
      until null stack repeat
        t:= reduce(((s::U)+x),p)
        if not ((flist:= probSplit(first stack,t,d)) case "failed") then
          stack:= rest stack
          for fact in flist repeat
            f1:= makeMonic(fact)
            (df1:= degree(f1)) = 0 => nil
            df1 > d => stack:= [f1,:stack]
            ans:= [f1,:ans]
        p = 2 =>
          ss:= ss + 1
          x := y * decode(ss, p, y)
        s:= s+1
        s = p =>
          s:= 0
          ss := ss + 1
          x:= y * decode(ss, p, y)
          not one? leadingCoefficient(x) =>
              ss := p ** degree x
              x:= y ** (degree(x) + 1)
      [c * first(ans),:rest(ans)]

    probSplit(u,t,d) ==
      (p:=modulus(u)) = 2 => probSplit2(u,t,d)
      f1:= gcd(u,t)
      r:= ((p**(d:NNI)-1) quo 2):NNI
      n:= exptmod(t,r,u)
      f2:= gcd(u,n + 1)
      (g:= f1 * f2) = 1 => "failed"
      g = u => "failed"
      [f1,f2,(u quo g)]

    probSplit2(u,t,d) ==
      f:= gcd(u,trace(t,d,u))
      f = 1 => "failed"
      degree u = degree f => "failed"
      [1,f,u quo f]

    trace(t,d,u) ==
      p:=modulus(t)
      d:= d - 1
      tt:=t
      while positive? d repeat
        tt:= (tt + (t:=exptmod(t,p,u))) rem u
        d:= d - 1
      tt

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