Initial population.

1 files changed, 282 insertions, 0 deletions

diff --git a/src/algebra/moddfact.spad.pamphlet b/src/algebra/moddfact.spad.pamphlet
new file mode 100644
index 00000000..19692d42
--- /dev/null
+++ b/src/algebra/moddfact.spad.pamphlet
@@ -0,0 +1,282 @@
+\documentclass{article}
+\usepackage{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))
+        reduction(r.remainder,p) ^=0 => "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) ==
+      i < 0 => error("negative exponentiation not allowed for exptMod")
+      ans:= 1$EMR
+      while i > 0 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 (c:= lc(u)) ^= 1$I then u:= makeMonic(u)
+      ans:= sepfact(ddfact(u))
+      cons(c::EMR,[makeMonic(f) for f in ans | degree(f) > 0])
+
+    gcd(u,v,q) == gcd(reduce(u,q),reduce(v,q))::U
+
+    factor(u,q) ==
+      v:= reduce(u,q)
+      dv:= reduce(differentiate(u),q)
+      degree gcd(v,dv) > 0 =>
+        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 (c:= lc(u)) ^= 1$I then u:= makeMonic u
+      ans:DDList:= []
+      repeat
+        w:= exptmod(w,p,u)
+        g:= gcd(w - m,u)
+        if degree g > 0 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 (c:= lc(u)) ^= 1$I 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 in 0.. while du > 0 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) =>
+          not (leadingCoefficient(x) = 1) =>
+              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 d > 0 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}