Initial population.

1 files changed, 233 insertions, 0 deletions

diff --git a/src/algebra/gaussfac.spad.pamphlet b/src/algebra/gaussfac.spad.pamphlet
new file mode 100644
index 00000000..1b1e3197
--- /dev/null
+++ b/src/algebra/gaussfac.spad.pamphlet
@@ -0,0 +1,233 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra gaussfac.spad}
+\author{Patrizia Gianni}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package GAUSSFAC GaussianFactorizationPackage}
+<<package GAUSSFAC GaussianFactorizationPackage>>=
+)abbrev package GAUSSFAC GaussianFactorizationPackage
+++ Author: Patrizia Gianni
+++ Date Created: Summer 1986
+++ Date Last Updated:
+++ Basic Functions:
+++ Related Constructors:
+++ Also See:
+++ AMS Classifications:
+++ Keywords:
+++ References:
+++ Description: Package for the factorization of complex or gaussian
+++ integers.
+GaussianFactorizationPackage() : C == T
+ where
+  NNI  ==  NonNegativeInteger
+  Z      ==> Integer
+  ZI     ==> Complex Z
+  FRZ    ==> Factored ZI
+  fUnion ==> Union("nil", "sqfr", "irred", "prime")
+  FFE    ==> Record(flg:fUnion, fctr:ZI, xpnt:Integer)
+
+  C  == with
+     factor      :     ZI     ->     FRZ
+       ++ factor(zi) produces the complete factorization of the complex
+       ++ integer zi.
+     sumSquares  :     Z      ->    List Z
+       ++ sumSquares(p) construct \spad{a} and b such that \spad{a**2+b**2}
+       ++ is equal to
+       ++ the integer prime p, and otherwise returns an error.
+       ++ It will succeed if the prime number p is 2 or congruent to 1
+       ++ mod 4.
+     prime?      :     ZI     ->    Boolean
+       ++ prime?(zi) tests if the complex integer zi is prime.
+
+  T  == add
+     import IntegerFactorizationPackage Z
+
+     reduction(u:Z,p:Z):Z ==
+       p=0 => u
+       positiveRemainder(u,p)
+
+     merge(p:Z,q:Z):Union(Z,"failed") ==
+       p = q => p
+       p = 0 => q
+       q = 0 => p
+       "failed"
+
+     exactquo(u:Z,v:Z,p:Z):Union(Z,"failed") ==
+        p=0 => u exquo v
+        v rem p = 0 => "failed"
+        positiveRemainder((extendedEuclidean(v,p,u)::Record(coef1:Z,coef2:Z)).coef1,p)
+
+     FMod := ModularRing(Z,Z,reduction,merge,exactquo)
+
+     fact2:ZI:= complex(1,1)
+
+             ----  find the solution of x**2+1 mod q  ----
+     findelt(q:Z) : Z ==
+       q1:=q-1
+       r:=q1
+       r1:=r exquo 4
+       while ^(r1 case "failed") repeat
+         r:=r1::Z
+         r1:=r exquo 2
+       s : FMod := reduce(1,q)
+       qq1:FMod :=reduce(q1,q)
+       for i in 2.. while (s=1 or s=qq1) repeat
+         s:=reduce(i,q)**(r::NNI)
+       t:=s
+       while t^=qq1 repeat
+         s:=t
+         t:=t**2
+       s::Z
+
+
+     ---- write p, congruent to 1 mod 4, as a sum of two squares ----
+     sumsq1(p:Z) : List Z ==
+       s:= findelt(p)
+       u:=p
+       while u**2>p repeat
+         w:=u rem s
+         u:=s
+         s:=w
+       [u,s]
+
+            ---- factorization of an integer  ----
+     intfactor(n:Z) : Factored ZI ==
+       lfn:= factor n
+       r : List FFE :=[]
+       unity:ZI:=complex(unit lfn,0)
+       for term in (factorList lfn) repeat
+         n:=term.fctr
+         exp:=term.xpnt
+         n=2 =>
+           r :=concat(["prime",fact2,2*exp]$FFE,r)
+           unity:=unity*complex(0,-1)**(exp rem 4)::NNI
+
+         (n rem 4) = 3 => r:=concat(["prime",complex(n,0),exp]$FFE,r)
+
+         sz:=sumsq1(n)
+         z:=complex(sz.1,sz.2)
+         r:=concat(["prime",z,exp]$FFE,
+                 concat(["prime",conjugate(z),exp]$FFE,r))
+       makeFR(unity,r)
+
+           ---- factorization of a gaussian number  ----
+     factor(m:ZI) : FRZ ==
+       m=0 => primeFactor(0,1)
+       a:= real m
+
+       (b:= imag m)=0 => intfactor(a) :: FRZ
+
+       a=0 =>
+         ris:=intfactor(b)
+         unity:= unit(ris)*complex(0,1)
+         makeFR(unity,factorList ris)
+
+       d:=gcd(a,b)
+       result : List FFE :=[]
+       unity:ZI:=1$ZI
+
+       if d^=1 then
+         a:=(a exquo d)::Z
+         b:=(b exquo d)::Z
+         r:= intfactor(d)
+         result:=factorList r
+         unity:=unit r
+         m:=complex(a,b)
+
+       n:Z:=a**2+b**2
+       factn:= factorList(factor n)
+       part:FFE:=["prime",0$ZI,0]
+       for term in factn repeat
+         n:=term.fctr
+         exp:=term.xpnt
+         n=2 =>
+           part:= ["prime",fact2,exp]$FFE
+           m:=m quo (fact2**exp:NNI)
+           result:=concat(part,result)
+
+         (n rem 4) = 3 =>
+           g0:=complex(n,0)
+           part:= ["prime",g0,exp quo 2]$FFE
+           m:=m quo g0
+           result:=concat(part,result)
+
+         z:=gcd(m,complex(n,0))
+         part:= ["prime",z,exp]$FFE
+         z:=z**(exp:NNI)
+         m:=m quo z
+         result:=concat(part,result)
+
+       if m^=1 then unity:=unity * m
+       makeFR(unity,result)
+
+           ----  write p prime like sum of two squares  ----
+     sumSquares(p:Z) : List Z ==
+       p=2 => [1,1]
+       p rem 4 ^= 1 => error "no solutions"
+       sumsq1(p)
+
+
+     prime?(a:ZI) : Boolean ==
+        n : Z := norm a
+        n=0 => false            -- zero
+        n=1 => false            -- units
+        prime?(n)$IntegerPrimesPackage(Z)  => true
+        re : Z := real a
+        im : Z := imag a
+        re^=0 and im^=0 => false
+        p : Z := abs(re+im)     -- a is of the form p, -p, %i*p or -%i*p
+        p rem 4 ^= 3 => false
+        -- return-value true, if p is a rational prime,
+        -- and false, otherwise
+        prime?(p)$IntegerPrimesPackage(Z)
+
+@
+\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 GAUSSFAC GaussianFactorizationPackage>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}