///////////////////////////////  ARBINT.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 when needed for rational arithmetic

#include "arbint.h"



const ARBPERLONG = sizeof(long) / sizeof(ARB_TYPE) ;
const ARBPERINT  = 1;        // = size(int) / sizeof(ARB_TYPE);
const int CHUNKSIZE = 8 ;    // acceptable space to waste in an Arbint
/////////////////////////  auxilary functions //////////////////////////////



int hexVALUE( char c) {  
	if (isdigit(c)) return c-'0';
	else if (c >= 'a' && c <= 'f') return c - 'a' + 10;
	else return c - 'A' + 10;
}

inline int isodigit(int c) {  
	return isdigit(c) && c < '8';
}


void Arbint :: addonedigit( unsigned int multiplier, int addend){ 
	ARB_TYPE * i_s = & p[2] ;
	ARB_LTYPE k = addend;
	/*  assumes  p[1] = LENGTH is positive   */
	for (int l = 0; l < LENGTH; l++){  
		k += (ARB_LTYPE) i_s[l] * multiplier;
		i_s[l] = ARB_TYPE (k);
		k = k >> ARB_SIZE;
	}
	if ( k != 0 )  {           // extend field 
		i_s = new ARB_TYPE[LENGTH+3];
		memcpy((char *)i_s, (char*)p, (LENGTH+2)*sizeof(ARB_TYPE));
		i_s[LENGTH+2] = ARB_TYPE(k);
		delete [] p;
		p = i_s ;
		(LENGTH)++ ;
	}
}

Arbint copy( const Arbint & u ) {  
	Arbint w;									// null pointer
	int ulen = abs( ARB_SIGNED(u.LENGTH) );
	w.p = new ARB_TYPE[ ulen + 2 ];
	w.p[0] = 1;									// no other references to this new VALUE
	memcpy ( &w.p[1], &u.p[1], (ulen+1) * sizeof(ARB_TYPE) );
	return w;
}


/////////////////////////   constructors and destructors  //////////////////


Arbint::Arbint()   {          // Arbint x;  or new Arbint
	// not initialized
    p = NULL;
}

Arbint::Arbint(const Arbint &x)  {  // Arbint x = y;
    p=x.p;
    REFCNT++;
}

Arbint::~Arbint() {  
	if ((p != NULL) && (--REFCNT == 0))   delete [] p;
}


Arbint::Arbint(int n)    {    // Arbint  x = -123;
	p = new ARB_TYPE[ARBPERINT+2];            // allocate data area
	REFCNT = 1;
	if ( n>=0 ) {
		LENGTH = ARBPERINT;
		VALUE = n;
	}
	else  {
		LENGTH = -1 ;
		VALUE = - n;
	}
}


Arbint::Arbint(unsigned int n)  {             // Arbint  x = 123
	p = new ARB_TYPE[ARBPERINT+2];            // allocate data area
	REFCNT = 1;
	LENGTH = ARBPERINT;
	VALUE = n;
}


Arbint::Arbint(long n) {                       // Arbint  x = -12334445
	ARB_LTYPE k = ( n < 0 ) ? -n : n;         // ignore overflow
	int nlen = ( k > ARB_MAX) ? 2 : 1;
	int j ;
	p = new ARB_TYPE[nlen+2];                 // allocate data area
	REFCNT = 1;
	LENGTH = ( n < 0 ) ? -nlen : nlen;
	for  ( j=1; j <= nlen;){  
		p[++j] = ARB_TYPE(k);
		k = k >> ARB_SIZE;
	}
}


Arbint::Arbint(unsigned long n)   {           // Arbint  x = 12334445
	int nlen =  ( n > ARB_MAX ) ? 2 : 1;
	int j ;
	p = new ARB_TYPE[nlen+2];                 // allocate data area
	REFCNT = 1;
	LENGTH = nlen ;
	for  ( j = 1; j <= nlen ; )   {             // no leading zeroes
		p[++j] = ARB_TYPE(n);
		n = n >> ARB_SIZE;
	}
}



ARB_TYPE * dassign( double num ){  
	const int dlen = DBL_MAX_EXP /(CHAR_BIT * (int) sizeof(ARB_TYPE));
	ARB_TYPE s[dlen+1];
	
	int nn (1);              // convert double to positive first
	if ( num < 0 ){
		num = - num;
		nn = -1;
	}
	
	ARB_TYPE *s1 = s ;   // go thru digits
	do {  
		*s1++ = (ARB_TYPE) fmod(num, ARB_BASE);
		num /= ARB_BASE;
	} while (num >= 1.0);
	
	int len = ARB_TYPE(s1 - s);    // number of digits needed
	ARB_TYPE * q = new ARB_TYPE[len+2];
	q[0] = 1;
	q[1] = nn * len;
	for ( nn=0; nn<len; nn++) q[2+nn] = s[nn] ;    // reverse digits
	
	return q;
}

Arbint::Arbint( double n )  {          // Arbint  x = 35.0;
	p = dassign( n );
}

Arbint& Arbint::operator= ( double n ) // x = 346.0;
{
    if (( p != NULL ) && (--(REFCNT) == 0))   delete [] p;
    p = dassign( n );
    return *this;
}



Arbint::Arbint(ARB_TYPE *w, int wlen, int offset) {
    // w not normalized, wlen can be positive or negative
	// offset can be 2 to indicate array for Arbint
	// with data for number starting at  w+offset
	ARB_TYPE * wv = w + offset ;
	int len = abs( wlen );
	while ( (--len > 0) && (wv[len] == 0 ) );
	len++ ;
	p = new ARB_TYPE[len+2];
	REFCNT = 1;
	LENGTH = len;
	if (( len == 1 )  && (wv[0] == 0) ) LENGTH = 1;          // generate unique 0: not -0!
	else if ( wlen < 0  ) LENGTH = - LENGTH;
	memcpy( &VALUE, wv, len * sizeof(ARB_TYPE));
	delete [] w;
}

Arbint::Arbint(int wlen, ARB_TYPE *w) {  
	// array w can be used, at most the high order digit is zero
	p = w;     // w also includes reference count and length
	
}


////////////////////////    Assignment operators    ////////////////////////

Arbint& Arbint::operator=(const Arbint  &x)  { //  y = x
	x.REFCNT++;             // increase reference count for new VALUE
	if ( p &&  (--(REFCNT) ==0))
		delete [] p;            // decrease reference count for old VALUE
	p = x.p;
	return( *this );
}



Arbint& Arbint::operator=(int n)  {     //  x = -1234;
	if ( p ) {  
		if (--REFCNT > 0)   {             // get rid of old data
			p = new ARB_TYPE[3];
		}
		else if (abs((ARB_SIGNED) (LENGTH)) != 1)  {      // reallocate old VALUE
			delete [] p;
			p = new ARB_TYPE[3];
		}
	}
	else    {                             // number did not yet exist
		p = new ARB_TYPE[3];
	}
	REFCNT = 1;
	if ( n < 0 ) {  
		LENGTH = - 1;
		VALUE = - n;
	}
	else {  
		LENGTH = 1;
		VALUE = n;
	}
	return( *this );
}


Arbint& Arbint::operator=(unsigned int n)  {          //  x = 655353
	if ( p ) {  
		if (--REFCNT > 0)   {             // get rid of old data
			p = new ARB_TYPE[3];
		}
		else if (abs((ARB_SIGNED) (LENGTH)) != 1)  {     // reallocate old VALUE
			delete [] p;
			p = new ARB_TYPE[3];
		}
	}
	else   {                            // number did not yet exist
		p = new ARB_TYPE[3];
	}
	REFCNT = LENGTH = 1;
	VALUE = n;                          // store new VALUE
	return( *this );
}



Arbint& Arbint::operator=(long n)  {     //  x = -12345678
	ARB_LTYPE k = ( n < 0 ) ? -n : n;   // check this if it works
	int nlen = ( k > ARB_MAX) ? 2 : 1;
	if ( p != NULL ){  
		if (--REFCNT > 0)   {             // get rid of old data
			p = new ARB_TYPE[nlen+2] ;
		}
		else if (abs((ARB_SIGNED) (LENGTH)) != nlen)  {  // reallocate old VALUE
			delete [] p;
			p = new ARB_TYPE[nlen+2];
		}
	}
	else {  
		p = new ARB_TYPE[nlen+2];
	}
	REFCNT = 1;
	LENGTH = (n>=0) ? nlen : - nlen ;
	
	for ( int j = 1; j <= nlen ; ) {   // assign the new VALUE
		p[++j] = ARB_TYPE(k) ;
		k = k >> ARB_SIZE;     // need to be checked !!!!!!!!!!!!!!
	}
	return * this;
}


Arbint& Arbint::operator=(unsigned long n) {    //  x = 12345678
	int nlen = (n > ARB_MAX) ? 2 : 1;
	if ( p != NULL )  {  
		if (REFCNT > 1)   {                       // get rid of old data
			REFCNT--;
			p = new ARB_TYPE[nlen+2];
		}
		else if (abs((ARB_SIGNED) (LENGTH)) != nlen)   {         // reallocate old VALUE
			delete [] p;
			p = new ARB_TYPE[nlen+2];
		}
	}
	else {  
		p = new ARB_TYPE[nlen+2];
	}
	REFCNT = 1;
	LENGTH = nlen++;
	for ( int j = 1; j <= nlen ; )  {                 // assign the new VALUE
		p[++j] = ARB_TYPE(n);
		n = n >> ARB_SIZE;
	}
	return *this;
}


////////////////////////// arithmetic operators /////////////////////////

Arbint operator+(const Arbint& u, const Arbint& v)  {    // addition
	int negans = 0;
	int add_sub = 1;     		      // 1  add, -1 sub
	int j;
	if ( u.isneg() )  {     		      // u < 0
		negans = 1;
		if ( v.isneg() == 0 ) add_sub = -1;     // -(u-v)
	}
	else                                       // u > 0
	{
		if ( v.isneg() ) add_sub = -1;             //  u - v
	}
	
	
	int ulen = abs((ARB_SIGNED) (u.LENGTH) );
	int vlen = abs((ARB_SIGNED) (v.LENGTH) );
	int wlen = (ulen < vlen) ? vlen : ulen ;
	
	ARB_TYPE *uv = & u.VALUE;
	ARB_TYPE *vv = & v.VALUE;
	ARB_TYPE tmp;
	ARB_LTYPE k = 0;
	ARB_TYPE *ww = new ARB_TYPE[wlen+2];   // addition can generate new digit
	ARB_TYPE *wv = ww + 2 ;		           // begin of data
	if (add_sub == 1 ){                     // addition
		// no normalization needed
		
		for (j = 0 ; j < ulen && j < vlen; j++ ) {	 
			k += uv[j];
			k += vv[j];
			wv[j] = ARB_TYPE(k) ;
			k =  k >> ARB_SIZE;
		}
		for (; j < ulen ; j++)  {  
			k +=  uv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		for (; 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[wlen+3] ;
			memcpy( & wv[2], &ww[2], (wlen * sizeof(ARB_TYPE)) );
			delete [] ww ;
			ww = wv ;
			ww[2+wlen] = ARB_TYPE(k);        // highest digit left to set
			wlen ++ ;
		}
		ww[0] = 1 ;						// refernce count
		ww[1] = (negans) ? -wlen : wlen ;
		return Arbint( wlen, ww ) ;     // no normalization needed
	}
	else  {                          // subtraction
		k = 1;
		for ( j = 0; j < ulen && j < vlen ; j++ ) {  
			tmp = ~vv[j];
			k +=  tmp;
			k +=  uv[j];
			wv[j] = ARB_TYPE(k);
			k =  k >> ARB_SIZE;
		}
		for (; j < ulen; j++ ) {	 
			k += V_HIGH;
			k += uv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		for (; j < vlen; j++ )  {	 
			tmp = ~vv[j];
			k += tmp;
			wv[j] = ARB_TYPE(k) ;
			k = k >> ARB_SIZE;
		}
		// high digit does not have to be set but tested for negative
		if ( (ARB_SIGNED) (k + V_HIGH) < 0 )  {    // answer is negative
			negans = 1 - negans;
			// convert to positive
			k = 1;
			for ( j = 0; j < wlen; j++ ) {
				tmp = ~wv[j];
				k += tmp;
				wv[j] = ARB_TYPE(k) ;
				k = k >> ARB_SIZE;
			}
		}
		if (negans) wlen = - wlen;
		if ( wv[j-1] != 0 ) {     // no normalization needed
			ww[0] = 1 ;             // reference count
			ww[1] = wlen ;
			return Arbint( wlen, ww ) ; 
		}
		else	
			return Arbint (ww, wlen, 2);  // normalize and return new Arbint
	}
	
}



Arbint operator-(const Arbint& u, const Arbint& v) {	   // subtraction
	int negans = 0;
	int add_sub = -1;		   // 1  add, -1 sub
	int j ;
	if ( u.isneg() ) {		   // u < 0
		negans = 1;
		if ( v.isneg() == 0 ) add_sub = 1;	   // -(u+v)
	}
	else			// u > 0
	{
		if ( v.isneg() ) add_sub = 1;		   //  u + v
	}
	
	int ulen = abs((ARB_SIGNED) (u.LENGTH) );
	int vlen = abs((ARB_SIGNED) (v.LENGTH) );
	int wlen = (ulen < vlen) ? vlen : ulen ;
	
	ARB_TYPE *uv = & u.VALUE;
	ARB_TYPE *vv = & v.VALUE;
	ARB_TYPE tmp;
	ARB_LTYPE k = 0;
	ARB_TYPE *ww = new ARB_TYPE[wlen+2];   // addition can generate new digit
	ARB_TYPE *wv = ww + 2 ;		           // begin of data
	if (add_sub == 1 ){                   // addition
		// no normalization needed
		
		for (j = 0 ; j < ulen && j < vlen; j++ ) {	 
			k += uv[j];
			k += vv[j];
			wv[j] = ARB_TYPE(k) ;
			k =  k >> ARB_SIZE;
		}
		for (; j < ulen ; j++)  {  
			k +=  uv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		for (; 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[wlen+3] ;
			memcpy( & wv[2], &ww[2], (wlen * sizeof(ARB_TYPE)) );
			delete [] ww ;
			ww = wv ;
			ww[2+wlen] = ARB_TYPE(k);        // highest digit left to set
			wlen ++ ;
		}
		ww[0] = 1 ;						// refernce count
		ww[1] = (negans) ? -wlen : wlen ;
		return Arbint( wlen, ww ) ;     // no normalization needed
	}
	else  {                          // subtraction
		k = 1;
		for ( j = 0; j < ulen && j < vlen ; j++ ) {  
			tmp = ~vv[j];
			k +=  tmp;
			k +=  uv[j];
			wv[j] = ARB_TYPE(k);
			k =  k >> ARB_SIZE;
		}
		for (; j < ulen; j++ ) {	 
			k += V_HIGH;
			k += uv[j];
			wv[j] = ARB_TYPE(k);
			k = k >> ARB_SIZE;
		}
		for (; j < vlen; j++ )  {	 
			tmp = ~vv[j];
			k += tmp;
			wv[j] = ARB_TYPE(k) ;
			k = k >> ARB_SIZE;
		}
		// high digit does not have to be set but tested for negative
		if ( (ARB_SIGNED) (k + V_HIGH) < 0 )     // answer is negative
		{
			negans = 1 - negans;
			// convert to positive
			k = 1;
			for ( j = 0; j < wlen; j++ ) {
				tmp = ~wv[j];
				k += tmp;
				wv[j] = ARB_TYPE(k) ;
				k = k >> ARB_SIZE;
			}
		}
		if (negans) wlen = - wlen;
		if ( wv[j-1] != 0 )       // no normalization needed
		{
			ww[0] = 1 ;             // reference count
			ww[1] = wlen ;
			return Arbint( wlen, ww ) ; 
		}
		else	
			return Arbint (ww, wlen, 2);  // normalize and return new Arbint
	}
	
}





Arbint operator*(const Arbint& p_u, const Arbint& p_v) {  
	int n = abs((ARB_SIGNED) (p_u.LENGTH));
	int m = abs((ARB_SIGNED) (p_v.LENGTH));
	/* Check for zero multiplicands or multiplier */
	
    if ( (n == 1 ) && (p_u.VALUE == 0 ) ) return p_u ;
	if ( (m == 1 ) && (p_v.VALUE == 0 ) ) return p_v ;
	
	int iplusj, wlen = n + m;
	ARB_TYPE *uv = & p_u.VALUE;
	ARB_TYPE *vv = & p_v.VALUE;
	
	
	ARB_TYPE *p_w = new ARB_TYPE[wlen+2];
	ARB_TYPE *wv = & p_w[2];
	
	// initialize space for new number
	memset((char*) wv, 0, wlen * sizeof(ARB_TYPE));
	
	
	// numbers are positive, remember sign of result
	ARB_LTYPE k = 1;
	int negans = 0;
	
	if ( (ARB_SIGNED) p_u.LENGTH < 0 )  negans = 1;
	if ( (ARB_SIGNED) p_v.LENGTH < 0 )  negans = 1 - negans;
	
	for (int j = 0; j < m; j++) {  
		ARB_LTYPE vj = vv[j];
		if (vj) {	 
			k = 0;
			for (int i = 0 ; i < n; i++) { 
				iplusj = i + j ;
				k += (ARB_LTYPE) uv[i] * vj + wv[iplusj];
				wv[iplusj] = ARB_TYPE(k) ;
				k = k >> ARB_SIZE;
			}
			wv[j+i] = ARB_TYPE(k) ;
		}
	}
	/*  will not normalize, at most the high order digit is zero
	check for it and reduce wlen if necessary */
	if ( wv[wlen-1] == 0 ) wlen-- ;
	p_w[1] = (negans) ? - wlen : wlen ;  // set length
	p_w[0] = 1 ;                         // reference count
	return Arbint(wlen, p_w) ;           // return new normalized Arbint
}

void divmod( const Arbint &p_u, const Arbint &p_v,
            Arbint &quotient, Arbint &remainder)
{  
	ARB_TYPE * uv = & p_u.p[2] ;
	ARB_TYPE * vv = & p_v.p[2] ;
	int i,j,cmp = 1;                           // set default for no exit
	int negremainder = p_u.isneg();            // sign same as dividend
	int negquotient  = p_v.isneg();            // sign of quotient see below
	negquotient = negquotient-negremainder;
	
	int m_n = abs((ARB_SIGNED) (p_u.LENGTH));             // m_n for size of quotient
	int n   = abs((ARB_SIGNED) (p_v.LENGTH));             // n for length of divisor

	
	// in case first digits are zero, discard them
	// no overflow can occur we are working with unsigned numbers
	
	while ( (m_n > 1 ) && (uv[m_n-1] == 0) )  --m_n ; // m_n is initialized
	while ( (n > 1 )   && (vv[n-1] == 0)   )   --n ;   //   n is initialized
	
	if (n == 1)
	{
		if (vv[0] == 0) 
		{	 
			perror("Division by zero in Arbint");
			quotient = remainder = 0 ;
			// need to return a VALUE in case someone does access results
		}
		else
		{	 
			ARB_TYPE *qq = new ARB_TYPE[m_n+2]; 
			ARB_TYPE *qv = qq + 2 ;       // max size for quotient
			ARB_TYPE prevu = 0;
			ARB_TYPE v1 = vv[0];					// divisor
			ARB_TYPE tmpq;   // quotient has to fit as unsigned
			ARB_LTYPE t;
			
			for (j = m_n-1; j >= 0 ; j--)
			{  
				t = prevu;
				t = t << ARB_SIZE;
				t += uv[j];
				qv[j] = tmpq = ARB_TYPE(t / v1) ;              
				prevu = ARB_TYPE(t - ARB_LTYPE(v1) * tmpq ) ;
			}
			if ( m_n > 1) 
				if ( qv[m_n-1] == 0 )  m_n-- ;
			else 
				if ( qv[0] == 0 )  negquotient = 0 ;

			qq[0] = 1 ;        // reference count 
			qq[1] = ( negquotient ) ? - m_n : m_n;

			Arbint q(m_n,qq);               // normalized quotient
			quotient = q;
			
			remainder = prevu;			// return remainder as Arbint
			if (negremainder)  remainder = - remainder;
		}
		return ;
	}

	if (m_n < n)                     // degenerate cases
	{  
		quotient = 0;
		remainder  = p_u ;
		return ;
	}
	if (m_n == n)
	{  
		cmp = 0;
		for (int j = m_n-1; j >=0 ; j--) 
		{
			if (uv[j] != vv[j]) 
			{  
				cmp = (uv[j] < vv[j]) ? -1 : 1;
				break;
			}
		}
		if ( cmp < 0 ) { 
			quotient = 0;           // u < v     remainder has to be set
			remainder = p_u ;
			return ;
		}
		else if (cmp == 0)         // u = v
		{  
			quotient = (negquotient) ? -1 : 1 ;
			remainder = 0;
			return ;
		}
	}

	// general case
	// can not change u or v so get new work space
	
	ARB_TYPE *tmp_u = new ARB_TYPE[m_n+1] ;
	ARB_TYPE *tmp_v = new ARB_TYPE[n];
	ARB_TYPE tmp;
	ARB_LTYPE K ;
	ARB_TYPE d = 1 ;
	
	// initialize workspace  

	if ( vv[n-1] > ARB_HIBIT )                             //  d=1
	{
		tmp_u[m_n] = 0 ;
		memcpy( tmp_u, uv, m_n * sizeof(ARB_TYPE) );
		memcpy( tmp_v, vv,   n * sizeof(ARB_TYPE) );
	}
	else                                 // multiply u and v  by d
	{
		d = ARB_TYPE(ARB_BASE / (vv[n-1] + 1 ));          
		K = 0;

		for (i =0; i <m_n; i++)
		{	 
			K += (ARB_LTYPE) uv[i] * d ;
			tmp_u[i] = ARB_TYPE(K) ;
			K = K >> ARB_SIZE;
		}
		tmp_u[i] = ARB_TYPE(K);
		
		K = 0;
		for (i= 0; i < n; i++)
		{	 
			K += (ARB_LTYPE) vv[i] * d;
			tmp_v[i] = ARB_TYPE(K) ;
			K = K >> ARB_SIZE;
		}
	}

	uv = tmp_u;          
	vv = tmp_v;

	ARB_TYPE v1 = vv[n-1];
	ARB_TYPE v2 = vv[n-2];
	int m = m_n - n;
	
	ARB_TYPE *qv = new ARB_TYPE[m+1];
	ARB_TYPE *nv = new ARB_TYPE[n+1];
	
	for ( j = m; j >= 0; j--)
	{  
		ARB_LTYPE q_hat;
		ARB_LTYPE u_j = uv[j+n];
		ARB_LTYPE u_j1 = uv[j+n-1];
		ARB_LTYPE u_j2 = uv[j+n-2];
		if (u_j == v1)
			q_hat = ARB_BASE - 1;
		else
		{	 
			q_hat = u_j;
			q_hat = ((q_hat << ARB_SIZE) + u_j1) / v1;
		}
		
		for (;; q_hat--)
		{	 
			ARB_LTYPE u_j_q_hat = u_j;
			u_j_q_hat = u_j_q_hat<< ARB_SIZE;
			u_j_q_hat += u_j1;
			u_j_q_hat -= v1 * q_hat;
			if (u_j_q_hat  >> ARB_SIZE) break;
			u_j_q_hat *= ARB_BASE;
			u_j_q_hat += u_j2;
			if ((v2 * q_hat) <= u_j_q_hat)  break;
		}
		
		// set nv  to q_hat * (v[1..n])
		K = 0;
	
		for (i =0 ; i < n; i++)
		{	 
			K += vv[i] * q_hat;
			nv[i] = ARB_TYPE(K);
			K = K >> ARB_SIZE;
		}
		nv[i] = ARB_TYPE(K);
		// subtract nv[n]...nv[0] from u[j+n]...u[j]
		// little endian notation
		K = 1;
		for (i = 0 ; i <= n;  i++ )
		{	 
			K +=  uv[j+i];
			tmp = ~nv[i];
			K += tmp;
			uv[i+j] = ARB_TYPE(K);
			K = K >> ARB_SIZE;
		}
		
		if (K == 0)         // result got negative so add one v back in
		{                   // notice no carry is overflow !!
			for (i = 0 ; i <= n; i++ )
			{  
				K += uv[j+i];
				K += vv[i];
				uv[j+i] = ARB_TYPE(K) ;
				K = K >> ARB_SIZE;
			}
			uv[j+i] += ARB_TYPE(K);
			q_hat--;
		}
		qv[j] = ARB_TYPE(q_hat) ;
	}
	
	Arbint q(qv, m+1);              // normalize, and delete extra space
	if (negquotient) quotient = - q ;
	else  quotient = q;
	
	// restore remainder by dividing by d

	if (d > 1) 
	{
		ARB_LTYPE t, prevu = 0;
		ARB_LTYPE tmpq;

		for (int i= n; i >= 0 ; i--)     // division of u[0..n] by  d
		{	 
			t = prevu;
			t = t << ARB_SIZE;
			t += uv[i];
		    tmpq = (ARB_LTYPE) t/d;
			uv[i] = (ARB_TYPE) tmpq ;
			prevu = ARB_TYPE( t - (ARB_LTYPE) d * tmpq );
		}
	}
	
	Arbint r(uv,n);    // remainder starts at uv[0]
	// have to give all of uv in order to  delete uv
	if (negremainder) remainder = - r ;
	else   remainder = r;
	
	delete [] tmp_v;
	delete [] nv;
}

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


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

Arbint operator+(Arbint &j)       // unary plus
{                                 // removed const as return  modifies j
    return j;
}

Arbint operator-(const Arbint& u) // unary minus
{  
	if (u.LENGTH == 1 && u.VALUE == 0)
		return u;                   // do not change sign of 0
	Arbint w = copy( u );          // need new copy of Arbint u
	w.LENGTH = - w.LENGTH;         // to change sign
	return w;
}


/////////////////////  Relational Operators   /////////////////////////////

int arb_cmp(ARB_TYPE *a, ARB_TYPE *b, int len) {
	/* len must be positive */
    /* return 1 if a < b else return 0 */
    for (int i=len-1; i>=0; i--)
    {  
		if (a[i] < b[i]) return 1;
		if (a[i] > b[i]) return 0;
    }
    return 0;
}



int operator==(const Arbint& a, const Arbint& b)
{  
	return ((a.LENGTH == b.LENGTH) &&
           (memcmp( & a.VALUE, & b.VALUE,
           abs((ARB_SIGNED) (a.LENGTH))* sizeof(ARB_TYPE)) == 0));
}

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


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


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


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


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



/////////////////////////////  input  output   ////////////////////////////

ostream& operator<< (ostream& out, const Arbint& val){ 
	if (val.p == NULL) {
		out << "NaN " ;
	}
	else {
		int i ;
		int val_len = abs((ARB_SIGNED) (val.LENGTH));
		ARB_TYPE *val_new = new ARB_TYPE[val_len];     // need private copy
		memcpy(val_new, &val.VALUE, val_len*sizeof(ARB_TYPE));
		if (val.isneg())  out << "-";                // print sign if negative
		ARB_TYPE *digs = new ARB_TYPE[val_len*BIN_DEC_RATIO];
		ARB_TYPE *svdigs = digs;
		
		// create digits in reverse order
		val_len -- ;
		while( val_len >= 0 )  {  
			ARB_LTYPE prevu = 0;
			
			for (i = val_len ; i >=0 ; i--) {	
				prevu = prevu << ARB_SIZE;
				prevu += val_new[i];
				ARB_TYPE tmpq = ARB_TYPE(prevu / DEC_BASE);
				val_new[i] = tmpq;
				prevu -= DEC_L_BASE * tmpq;
			}
			*digs++ = ARB_TYPE(prevu) ;
			
			while( val_len >= 0 ) {
				if (val_new[val_len] != 0) break;
				val_len-- ;
			}
		}
		out << *--digs;
		out.fill('0');
		while (digs-- > svdigs) {  
			out.width(4);
			out << *digs;
		}
		out.fill(' ');
		delete [] svdigs;
		delete [] val_new;
	}
	return out ;
}


istream& operator>> (istream& in, Arbint& x) { 
	char c, sign = 0;
	
	// delete x if arbint x already exists 
	// if illegal characters are entered 0 will be returned
	if ((x.p != NULL) && (--(x.REFCNT)==0 ) ) delete [] x.p ;
	
	x.p = new ARB_TYPE[ARBPERINT + 2 ] ;
	x.REFCNT = x.LENGTH = 1 ; 
	x.VALUE = 0 ;
	
	in >> ws;                        // skip white spaces
	if (!in.get(c)) return in;       // return 0 on no more input
	
	if ( (c == '+') || (c == '-')) {       
		sign = c;	                  //  remember the sign
		if (!in.get(c)) return in;    //  return 0 on no more input
	}
	
	switch (c) {  
	case 'o': case 'O':           // octal  letter O
		for (in.get(c); in && isodigit(c); in.get(c))
			x.addonedigit(8, (c-'0'));
		break;
	case 'x': case 'X':           // hexadecimal letter X
		for (in.get(c); in && isdigit(c); in.get(c))
			x.addonedigit(16, hexVALUE(c));
		break;
	default:                      // decimal number
		for (; in && isdigit(c); in.get(c))
			x.addonedigit(10, (c-'0'));
		break;
	}
	if (in) in.putback(c);
	
	if (sign == '-') x.LENGTH = - x.LENGTH ;
	
	return in;
}

size_t fread( Arbint & a, FILE * stream ){  
	size_t len;
	size_t i = fread( &len, sizeof(size_t), 1, stream);
	if (i<=0) return(i);
	i = abs( (ARB_SIGNED) len ) ;
	a.p = new ARB_TYPE[i + 2];
	a.LENGTH = len;
	a.REFCNT = 1;
	return( fread( &a.VALUE, i, sizeof(ARB_TYPE), stream ) );
}

size_t fwrite( const Arbint & a, FILE * stream ){  
	size_t len = abs((ARB_SIGNED) (a.LENGTH));
	size_t i = fwrite( &a.LENGTH, sizeof(size_t), 1, stream);
	if (i != 1 ) cerr << "failed to write out length of arbint" ;
	return( fwrite( &a.VALUE, len, sizeof(ARB_TYPE), stream ) );
}

Arbint::Arbint(char c []) {  
	char sign = 0;
	while ( *c == ' ' ) c++;             // ignore leading blanks
	if ((*c=='+')||(*c=='-'))   sign = *c++;
	while ( *c == ' ' ) c++;             // ignore more blanks
	
	p = new ARB_TYPE[3];                 // basic size
	REFCNT = LENGTH = 1;
	VALUE= 0 ;
	
	switch (* c) {  
	case 'o': case 'O':               // octal  letter O
		for ( c++; *c != 0; *c++)
			this -> addonedigit(8, ( *c-'0'));
		break;
	case 'x': case 'X':               // hexadecimal letter X
		for ( c++; *c != 0; *c++)
			this -> addonedigit(16, hexVALUE(*c));
		break;
	default:                          // decimal number
		for (;  *c != 0; *c++)
			this -> addonedigit(10, (*c-'0'));
		break;
	}
	if (sign == '-')  LENGTH = - LENGTH ;
	
}

#endif /* RAT_COEF  */