/////////////////////   ARBINTMX.CPP   ///////////////////////////

#include <ctype.h>
#include <float.h>
#include <iostream.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

//#include "user.h"
//#include "coef.h"
#define RAT_COEF

#ifdef RAT_COEF
////////////////////  Compile only for rational arithmetic  //////


#include "arbint.h"

///////////////////// mixed arithmetic operators:  Arbint&, int ///////////

Arbint operator+(const Arbint& u, int v)          // addition
{
	//    u + v    with u = +/-U   v = +/-V
	//	 +U + +V     perform   U + V
	//	 +U + -V               U - V
	//	 -U + +V             -(U - V)
	//	 -U + -V             -(U + V)
	int negans = 0;
	int add_sub = 1;     		 // 1  add, -1 sub
	ARB_TYPE tmp;
	ARB_LTYPE k = 1;
	if ( u.isneg() )       		 // u < 0
	{  
		negans = 1;
		tmp = ~v ;               // avoid overflow
		k = tmp + 1;
		add_sub =  ( v >= 0 ) ? -1 : 1 ;
	}
	else                                  // u > 0
	{  
		negans = 0 ;
		k = ARB_TYPE(v) ;    // avoid extension of sign bit
		add_sub = ( v < 0 ) ? -1 : 1 ; 
	}
	int ulen = abs((ARB_SIGNED) ( u.LENGTH ));    // wlen = ulen 
	int j ;
	
	ARB_TYPE *uv = & u.VALUE;
	ARB_TYPE *ww = new ARB_TYPE[ulen+2];
	ARB_TYPE *wv = ww + 2 ;
	
	if (add_sub == 1 )
	{  
		for ( j =0 ; j < ulen ; j++ ) 
		{
			k += uv[j] ;
			wv[j] = ARB_TYPE(k) ;
			k = k >> ARB_SIZE ;
		}
		if ( k != 0 )               // need one more digit
		{
			wv = new ARB_TYPE[ulen+3] ;
			memcpy( & wv[2], &ww[2], (ulen * sizeof(ARB_TYPE)) );
			delete [] ww ;
			ww = wv ;
			ww[2+ulen] = ARB_TYPE(k);        // highest digit left to set
			ulen ++ ;
		}
		ww[0] = 1 ;						// refernce count
		ww[1] = (negans) ? -ulen : ulen ;
		return Arbint( ulen, ww ) ;     // no normalization needed
	}
	else			                 // add_sub == -1
	{  
		k += uv[0] ;
		wv[0] = ARB_TYPE(k) ;
		k = k >> ARB_SIZE ;
		for ( j = 1 ; j < ulen ; j++ )
		{  
			k += V_HIGH;
			k += uv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		
		if ( (ARB_SIGNED) (k + V_HIGH) < 0 )     // answer is negative		
		{	 
			negans = 1 - negans;
			// convert to positive
			k = 1;
			for ( j = 0; j < ulen; j++ )
			{  
				tmp = ~wv[j];
				k += tmp;
				wv[j] = ARB_TYPE(k);
				k = k >> ARB_SIZE;
			}
		}
		
		if (negans) ulen = - ulen;
		if ( wv[j-1] != 0 ) {     // no normalization needed
			ww[0] = 1 ;             // reference count
			ww[1] = ulen ;
			return Arbint( ulen, ww ) ; 
		}
		else	
			return Arbint (ww, ulen, 2);  // normalize and return new Arbint
	}
	
}




Arbint operator-(const Arbint& u, int v) // subtraction
{
	//    u - v    with u = +/-U   v = +/-V
	//	 +U - +V     perform   U - V
	//	 +U - -V               U + V
	//	 -U - +V             -(U + V)
	//	 -U - -V             -(U - V)
	int j,negans =  0;
	int add_sub = -1;     		 // 1  add, -1 subtract
	ARB_TYPE tmp;
	ARB_LTYPE k = 1;
	if ( u.isneg() )
	{  
		negans = 1;
		if ( v >= 0 )
		{	 
			add_sub = 1;                    // -(u+v)
			k = v;
		}
		else			         //  u < 0, v < 0
		{	 
			tmp = ~v;                       // avoid overflow
			k += tmp;
		}
	}
	else                                  // u > 0
	{  
		if ( v < 0 )     add_sub = 1 ;     // u + (-v)  or u - (+v)
		tmp = ~v;
		k += tmp;
	}
	
	int ulen = abs((ARB_SIGNED) ( u.LENGTH ));    // wlen = ulen 
	ARB_TYPE *uv = & u.VALUE;
	ARB_TYPE *ww = new ARB_TYPE[ulen+2];
	ARB_TYPE *wv = ww + 2 ;	
	
	if (add_sub == 1 )
	{  
		
		for ( j =0 ; j < ulen ; j++ ) 
		{
			k += uv[j] ;
			wv[j] = ARB_TYPE(k) ;
			k = k >> ARB_SIZE ;
		}
		if ( k != 0 )               // need one more digit
		{
			wv = new ARB_TYPE[ulen+3] ;
			memcpy( & wv[2], &ww[2], (ulen * sizeof(ARB_TYPE)) );
			delete [] ww ;
			ww = wv ;
			ww[2+ulen] = ARB_TYPE(k);        // highest digit left to set
			ulen ++ ;
		}
		ww[0] = 1 ;						// refernce count
		ww[1] = (negans) ? -ulen : ulen ;
		return Arbint( ulen, ww ) ;     // no normalization needed
	}
	else			                 // add_sub == -1
	{  
		k += uv[0] ;
		wv[0] = ARB_TYPE(k) ;
		k = k >> ARB_SIZE ;
		for ( j = 1 ; j < ulen ; j++ )
		{  
			k += V_HIGH;
			k += uv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		
		if ( (ARB_SIGNED) (k + V_HIGH) < 0 )     // answer is negative		
		{	 
			negans = 1 - negans;
			// convert to positive
			k = 1;
			for ( j = 0; j < ulen; j++ )
			{  
				tmp = ~wv[j];
				k += tmp;
				wv[j] = ARB_TYPE(k);
				k = k >> ARB_SIZE;
			}
		}
		
		
		if (negans) ulen = - ulen;
		if ( wv[j-1] != 0 ) {     // no normalization needed
			ww[0] = 1 ;             // reference count
			ww[1] = ulen ;
			return Arbint( ulen, ww ) ; 
		}
		else	
			return Arbint (ww, ulen, 2);  // normalize and return new Arbint
	}
	
}







Arbint operator*(const Arbint& p_u, int v)
{  
	return operator*( v, p_u ) ;
}


void divmod( const Arbint &p_u, int v,
			Arbint &quotient, int & remainder)
{  
	int negremainder = p_u.isneg();        // sign same as dividend
	int negquotient  = ( v < 0 ) ? 1 : 0;  // sign of quotient see below
	
	int m_n = abs((ARB_SIGNED) (p_u.LENGTH));             // m_n for size of quotient
	
	// n=1 for lenght of divisor
	
	// can not change u so get new work space
	
	ARB_TYPE *tmp_u = new ARB_TYPE[m_n];
	
	ARB_TYPE v1;
	
	// initialize workspace
	memcpy( tmp_u, & p_u.VALUE, m_n * sizeof(ARB_TYPE));
	
	// change to positive if necessary no overflow is possible here
	
	if (negquotient)                       // tmp_v = - tmp_v;
	{  
		v1 = 1;
		v1 += ~v;                           // can ignore carry out
	}
	else v1 = v;
	
	negquotient = negquotient-negremainder;
	
	ARB_TYPE *uv = tmp_u;              // so that we can release tmp_u
	
	// in case first digits are zero, discard them
	// no overflow can occur we are working with unsigned numbers
	
	for ( ; (*uv == 0) && m_n > 1; uv++, m_n-- );  // m_n is initialized
	
	if (v1 == 0)
	{  
		perror("Division by zero in Arbintmx");
		// need to return a VALUE in case someone does access results
		quotient = 0 ;             //  Arbint
		remainder = 0 ;            //  integer
	}
	else
	{  
		ARB_TYPE *qv = new ARB_TYPE[m_n];         // max size for quotient
		ARB_TYPE prevu = 0;
		ARB_TYPE tmpq;   // quotient has to fit as unsigned
		ARB_LTYPE t;
		
		for (int r = m_n-1; r>= 0 ; r--)
		{  
			t = prevu;
			t = t << ARB_SIZE;
			t += uv[r];
			tmpq = ARB_TYPE(t / v1) ;
			qv[r] = tmpq;              // increase index first
			prevu = ARB_TYPE(t - ARB_LTYPE(v1) * tmpq ) ;
		}
		if (negquotient)  m_n = - m_n;
		Arbint q(qv,m_n);               // normalize quotient
		quotient = q;
		
		remainder = prevu;			// return remainder as Arbint
		if (negremainder)  remainder = - remainder;
	}
	
	delete [] tmp_u;
	return;
}

Arbint operator/( const Arbint &u, int v)
{  
	Arbint quo, rem;
	divmod(u,v,quo,rem);
	return quo;
}


int operator%( const Arbint &u, int v)
{  
	Arbint quo ;
	int  rem;
	divmod(u,v,quo,rem);
	return rem;
}




////////////////// mixed relational operators: Arbint&, int  //////////////

int operator==(const Arbint& a, int b)
{  
	if ( (ARB_SIGNED) (a.LENGTH) > 0 )
		return (a.LENGTH == 1) && (a.VALUE == b);
	else
		return (abs((ARB_SIGNED) (a.LENGTH)) == 1) && (a.VALUE == -b);
}

int operator!=(const Arbint& a, int b)
{  
	return !(a==b);
}

int operator<(const Arbint& a, int b)
{  
	int a_neg = a.isneg();
	int b_neg = (b < 0) ? 1 : 0;
	if (a_neg > b_neg) return 1;
	if (a_neg < b_neg) return 0;
	int a_len = abs((ARB_SIGNED) (a.LENGTH));
	if (!a_neg)                           // a and b positive
	{  
		if (a_len > 1) return 0;
		return a.VALUE < b;
	}
	else                                  // a and b negative
	{  
		if (a_len > 1) return 1;
		return a.VALUE > -b;
	}
}


int operator>(const Arbint& a, int b)
{  
	int a_neg = a.isneg();
	int b_neg = ( b < 0 ) ? 1 : 0;
	if (a_neg < b_neg) return 1;
	if (a_neg > b_neg) return 0;
	int a_len = abs((ARB_SIGNED) (a.LENGTH));
	if (!a_neg)                           // a and b positive
	{  
		if (a_len > 1) return 1;
		return a.VALUE > b;
	}
	else                                  // a and b negative
	{  
		if (a_len > 1) return 0;
		return a.VALUE <  - b;
	}
}


int operator>=(const Arbint& a, int b)
{  
	return !(a < b);
}


int operator<=(const Arbint& a, int b)
{  
	return !( a > b );
}



//////////////// mixed arithmetic operators:  int, Arbint& /////////////////

Arbint operator+(int u, const Arbint& v)
{  
	return operator+( v, u );
}


Arbint operator-(int u, const Arbint& v)        // subtraction
{  //    u - v    with u = +/-U   v = +/-V
	//	 +U - +V     perform   U - V
	//	 +U - -V               U + V
	//	 -U - +V             -(U + V)
	//	 -U - -V             -(U - V)
	int negans = 0;
	int add_sub = -1;     	                // 1  add, -1 sub
	ARB_TYPE tmp;
	ARB_LTYPE k = 1;
	if ( u < 0 )       	                        // u < 0
	{  
		negans = 1;
		tmp = ~u;                                 // avoid overflow
		k += tmp;
		if ( v.isneg() == 0 ) add_sub = 1;        // -(u+v)
	}
	else                                         // u > 0
	{  
		k = u;
		if ( v.isneg() )  add_sub = 1;            //  U + V
	}
	
	int vlen = abs((ARB_SIGNED) ( v.LENGTH) );
	int j ;
	ARB_TYPE *vv = & v.VALUE;	
	ARB_TYPE *ww = new ARB_TYPE[vlen+2];
	ARB_TYPE *wv = ww + 2 ;	
	if (add_sub == 1 )
	{  
		
		for ( j = 0; j < vlen; j++)
		{
			k +=  vv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		if ( k != 0 )               // need one more digit
		{
			wv = new ARB_TYPE[vlen+3] ;
			memcpy( & wv[2], &ww[2], (vlen * sizeof(ARB_TYPE)) );
			delete [] ww ;
			ww = wv ;
			ww[2+vlen] = ARB_TYPE(k);        // highest digit left to set
			vlen ++ ;
		}
		ww[0] = 1 ;						// reference count
		ww[1] = (negans) ? -vlen : vlen ;
		return Arbint( vlen, ww ) ;     // no normalization needed
	}
	
	else                            // add_sub == -1
	{								// overflow not possible
		tmp = ~vv[0] ;
		k +=  tmp;
		wv[0] = ARB_TYPE(++k) ;
		k =  k >> ARB_SIZE;
		for (j=1; j < vlen; j++ )
		{  
			tmp = ~vv[j] ;
			k += tmp ;
			wv[j] = ARB_TYPE(k) ;
			k = k >> ARB_SIZE;
		}
		
		if ( (ARB_SIGNED) (k + V_HIGH) < 0 )    // answer is negative	
		{	 
			negans = 1 - negans;
			// convert to positive
			k = 1;
			for ( j = 0; j < vlen; j++ )
			{  
				tmp = ~wv[j];
				k += tmp;
				wv[j] = ARB_TYPE(k);
				k = k >> ARB_SIZE;
			}
		}
		
		if (negans) vlen = - vlen;
		if ( wv[j-1] != 0 ) {     // no normalization needed
			ww[0] = 1 ;             // reference count
			ww[1] = vlen ;
			return Arbint( vlen, ww ) ; 
		}
		else	
			return Arbint (ww, vlen, 2);  // normalize and return new Arbint
	}
	
}


Arbint operator*(int u, const Arbint& p_v)
{  
	int m = abs((ARB_SIGNED) (p_v.LENGTH));
	if ( u == 0 ) return u ;     // will be changed to Arbint
	if ( (m == 1 ) && (p_v.VALUE == 0 ) ) return p_v ;
	int wlen = m + 1;
	ARB_TYPE *vv = &p_v.VALUE;
	ARB_TYPE *ww = new ARB_TYPE[wlen+2];
	ww[0] = 1 ;   // reference count = 1
	ARB_TYPE * wv = ww+2 ;
	
	// initialize workspace
	memset( wv, 0, wlen * sizeof(ARB_TYPE));
	
	
	// change to positive if necessary no overflow is possible here
	ARB_LTYPE k = 1 , uu ;
	
	int negans = 0;
	if ( u < 0 )
	{  
		negans = 1;
		k += ~u;
		uu = ARB_TYPE(k);
	}
	else uu = u ;
	if ( p_v.isneg() )
	{                          // change to positive
		negans = 1 - negans;
	}	
	
	k = 0;
	
	for (int j = 0; j < m; j++)
	{  
		k +=  uu * vv[j];
		wv[j] = ARB_TYPE(k);
		k =  k >> ARB_SIZE;
	}
	if ((wv[j] = ARB_TYPE(k))==0 )  wlen-- ;     //reduce length, but keep zero field

	ww[1] = (negans) ? -wlen : wlen ;
	
	return Arbint(wlen, ww) ;           // return new normalized Arbint
	
}





// the following three are probably never used

void divmod( int u, const Arbint& p_v, Arbint& quotient, Arbint& remainder)
{  
	Arbint p_u = u;
	divmod( p_u, p_v, quotient, remainder );
}

Arbint operator/( int u, const Arbint& p_v)
{  
	Arbint quo, rem;
	Arbint p_u = u;
	divmod ( p_u, p_v, quo, rem );
	return quo;
}

Arbint operator%( int u, const Arbint& p_v)
{  
	Arbint quo, rem;
	Arbint p_u = u;
	divmod ( p_u, p_v, quo, rem );
	return rem;
}


/////////////// mixed relational operators:  int, Arbint& /////////////////

int operator==(int a, const Arbint& b)
{  
	return b == a;
}


int operator!=(int a, const Arbint& b)
{  
	return b != a;
}


int operator<(int a, const Arbint& b)
{  
	return (b > a);
}


int operator<=(int a, const Arbint& b)
{  
	return !(b < a);
}


int operator>(int a, const Arbint& b)
{  
	return (b < a);
}


int operator>=(int a, const Arbint& b)
{  
	return !(b > a);
}


Arbint power( const Arbint & x, int n ) {
	Arbint y = 1 ;
	Arbint z = x ;
	while ( n > 0 ) {
		if ( n & 1 ) y = z * y ; 
		z = z * z;
		n >>= 1 ;
	}
	return y ;
}


#endif   /* RAT_COEF */


/*
#define INT_type Arbint

INT_type gcd(const INT_type & a,const INT_type & b )
{
	//   return egcd( a,b ) ;
	INT_type  r = a ;
	INT_type  g = b ;
	INT_type  d ;
	if ( r < 0 ) r = - r ;
	if ( b < 0 ) g = - g ;
	if ( r > g )
	{
		d = g ;
		g = r ;
		r = d ;
	}
	while ( r != 0 )
	{
		d = g ;
		g = r ;
		r = d % g ;      
	}
	//	 if (g != egcd(a,b))
	//	 { cout << "\nerror a=" << a <<"  b="<< b <<"  g="<<g <<" ECG="<<egcd(a,b);
	//	 exit(1) ;
	//	 } 
	 	 cout <<" result="<<r<<" d="<<d<<" g="<<g <<endl ;
		 return g ;
}


*/