\documentclass{article}
\usepackage{open-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 not (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 not one? d 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 not one? m then unity:=unity * m
       makeFR(unity,result)

           ----  write p prime like sum of two squares  ----
     sumSquares(p:Z) : List Z ==
       p=2 => [1,1]
       not one?(p rem 4) => 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
        not zero? re and not zero? im => 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}