#ifndef GAUSSIANH
#define GAUSSIANH

/* ************************************************************************ *
 *                                                                          *
 *              Gaussian.h  (C) 1999    Michael R. Toth                     *
 *                                                                          *
 *    Gaussian integer class:                                               *
 *    The class parameter "Intg" should implement  +, -, *, /, +=, -=, *=,  *
 *  =, ==, !=, <, << and >> to and from streams, and  abs .  It should be   *
 *  constructable from an "int".                                            *
 *                                                                          *
 * ************************************************************************ */

#include <stdlib.h>
#include <iostream.h>

template<class Intg =int> class  Gaussian
{
  friend Gaussian operator- ( const Gaussian& );
  friend Gaussian operator~ ( const Gaussian& );   // Complex conjugate.
  friend Intg abs( const Gaussian& );   // Square of the absolute value!
  friend Intg& Re( Gaussian& );
  friend Intg& Im( Gaussian& );
  friend ostream& operator<< ( ostream&, const Gaussian& );
  friend istream& operator>> ( istream&, Gaussian& );
private:
  Intg  real, imag;
  void Idiv( Intg&, const Intg& );
  bool Invert( );    // Only  1, -1, i, and -i  are invertable.
public:
  Gaussian ( )  { }
  Gaussian ( const Gaussian& that ): real(that.real), imag(that.imag)  { }
  Gaussian ( int x ): real(Intg(x)), imag(0)  { }
  Gaussian ( const Intg& x, const Intg& y =0 ): real(x), imag(y)  { }
  Gaussian& operator= ( const Gaussian& that )
                      {
                        if ( this!=&that )  real=that.real, imag=that.imag;
                        return *this;
                      }
  bool operator== ( const Gaussian& that ) const
                  { return  real==that.real && imag==that.imag; }
  bool operator!= ( const Gaussian& that ) const
                  { return  real!=that.real || imag!=that.imag; }
  Gaussian& operator+= ( const Gaussian& );
  Gaussian& operator+= ( const Intg& x )  { real+=x;  return *this; }
  Gaussian& operator-= ( const Gaussian& );
  Gaussian& operator-= ( const Intg& x )  { real-=x;  return *this; }
  Gaussian& operator*= ( const Gaussian& );
  Gaussian& operator*= ( const Intg& x )  { real*=x, imag*=x;  return *this; }
  Gaussian& operator/= ( const Gaussian& );
  Gaussian& operator%= ( const Gaussian& );
  Gaussian operator[] ( int ) const;   // Small integer power.
  bool Associate( const Gaussian& ) const;
};

template<class Intg> inline
Intg& Re( Gaussian<Intg>& z )  { return z.real; }

template<class Intg> inline
Intg& Im( Gaussian<Intg>& z )  { return z.imag; }

template<class Intg> inline  Intg abs( const Gaussian<Intg>& z )
                                     { return  z.real*z.real+z.imag*z.imag; }

template<class Intg> inline
Gaussian<Intg> operator~ ( const Gaussian<Intg>& z )
                         { return  Gaussian<Intg>(z.real,-z.imag); }

template<class Intg> inline
Gaussian<Intg> operator- ( const Gaussian<Intg>& z )
                         { return  Gaussian<Intg>(-z.real,-z.imag); }

template<class Intg> inline
Gaussian<Intg> operator- ( Gaussian<Intg> x, const Gaussian<Intg>& y )
                         { return  x-=y; }

template<class Intg> inline
Gaussian<Intg> operator+ ( Gaussian<Intg> x, const Gaussian<Intg>& y )
                         { return  x+=y; }

template<class Intg> inline
Gaussian<Intg> operator* ( Gaussian<Intg> x, const Gaussian<Intg>& y )
                         { return  x*=y; }

template<class Intg> inline
Gaussian<Intg> operator/ ( Gaussian<Intg> x, const Gaussian<Intg>& y )
                         { return  x/=y; }

template<class Intg> inline
Gaussian<Intg> operator% ( Gaussian<Intg> x, const Gaussian<Intg>& y )
                         { return  x%=y; }

template<class Intg> inline
Gaussian<Intg>& Gaussian<Intg>::operator+= ( const Gaussian& that )
                { real+=that.real, imag+=that.imag;  return *this; }
          
template<class Intg> inline
Gaussian<Intg>& Gaussian<Intg>::operator-= ( const Gaussian& that )
                { real-=that.real, imag-=that.imag;  return *this; }

template<class Intg>
Gaussian<Intg>& Gaussian<Intg>::operator*= ( const Gaussian& that )
                {
                  Intg r=real*that.real-imag*that.imag;
                  imag=imag*that.real+real*that.imag;
                  real=r;
                  return *this;
                }

template<class Intg> inline
Gaussian<Intg>& Gaussian<Intg>::operator/= ( const Gaussian& that )
                {
                  Intg d=abs(that);
                  operator*=(~that);
                  Idiv(real,d), Idiv(imag,d);
                  return *this;
                }

template<class Intg>  void Gaussian<Intg>::Idiv( Intg& q, const Intg& d )
                      // Set  q  to the closest integer to the rational  q/d .
                           {
                             if ( d==0 || q==0 )  return;
                             Intg n=q;
                             q=q/d;
                             Intg r=n-d*q;
                             if ( abs(r+r)>d )  q+= r>0? 1: -1;
                           }

template<class Intg> inline
Gaussian<Intg>& Gaussian<Intg>::operator%= ( const Gaussian& that )
                { *this-=*this/that*that;  return *this; }


template<class Intg>
bool Gaussian<Intg>::Associate( const Gaussian& z ) const
     {
       if ( real==z.real && imag==z.imag )  return true;
       if ( real==z.imag && imag==-z.real )  return true;
       if ( real==-z.real && imag==-z.imag )  return true;
       if ( real==-z.imag && imag==z.real )  return true;
       return false;
     }

template<class Intg> inline bool Gaussian<Intg>::Invert( )
                                 {
                                   if ( !Associate(1) )  return false;
                                   imag=-imag;  return true;
                                 }

template<class Intg>
Gaussian<Intg> Gaussian<Intg>::operator[] ( int p ) const
               {
                 Gaussian x=*this;
                 if ( p<0 )  if ( x.Invert( ) )  p=-p;
                                  else return 0;
                 Gaussian y= p&1? x: Gaussian(1);
                 while ( p>>=1 )  { x*=x;  if ( p&1 )  y*=x; }
                 return y;
               }

template<class Intg>
ostream& operator<< ( ostream& out, const Gaussian<Intg>& z )
                    { 
                      if ( z.imag==0 )  { out<<z.real;  return out; }
                      if ( z.real!=0 )
                         { out<<z.real;  if ( z.imag>0 )  out<<"+"; }
                      if ( z.imag!=-1 && z.imag!=1 )  out<<z.imag;
                      if ( z.imag==-1 )  out<<"-";
                      out<<"i";
                      return out;
                    }

template<class Intg>
istream& operator>> ( istream& in, Gaussian<Intg>& z )
                    { in>>z.real>>z.imag;  return in; }

inline  void Gprompt( ostream& out )  { out<<" Re   Im ?    "; }

template<class Intg> inline
Gaussian<Intg> Ldivide( Gaussian<Intg>& n, const Gaussian<Intg>& d )
               // Return the result of long division of  n  by  d .
               // Put the remainder in  n .
                      { Gaussian<Intg> q=n/d;  n-=q*d;  return q; }

#endif GAUSSIANH