/////////////////////////////  GCD.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>


////////////////////////////   Compile for rational arithmetic  ///

#include "arbint.h"

///////////////// improved Euclidean algorithm  ////////////
///////////////// for greatest common divisor  of arbint /////////////



Arbint gcd( const Arbint & aa,const Arbint & bb )
{
	/*  use private copies of parameters  */

	Arbint  u( copy(aa) )  ;
	Arbint  v( copy(bb) ) ;
	Arbint  tmp, tmp1, tmp2,tmp3, tmp4 ;
	ARB_TYPE u_hat, v_hat ;
	ARB_LTYPE temp ;
	/* make positive */ ;
	u.LENGTH = abs((ARB_SIGNED) (u.LENGTH) ) ;
	v.LENGTH = abs((ARB_SIGNED) (v.LENGTH) ) ;
	if ( v > u )
	{
		tmp = v ;
		v = u ;
		u = tmp ;
	}
	
	// first word can not be zero
	

	while ( v.LENGTH > 1 )
	{                           // v.length can be between 1 and u.length
		
		temp = u.p[u.LENGTH+1] ;
		temp = temp << ARB_SIZE ;
		temp += u.p[u.LENGTH] ;
		short i = 0 ;                 // counter how far we shifted
		while ( (temp & ARB_LHIBIT) == 0 )    // look for leading bit
		{
			i++ ;
			temp = temp << 1 ;
		}
		temp = temp >> ARB_SIZE ;
		i = ARB_SIZE - i ;
		u_hat = ARB_TYPE(temp) ;  // u_hat has now ARB_size significant bits
		temp = 0 ;                // get bits in same position for v
		if ( u.LENGTH - v.LENGTH <= 1 )       // if v has bits in same place
		{
			temp = v.p[v.LENGTH+1] ;           // at least in one word
			if ( u.LENGTH == v.LENGTH )        // bits in two words
			{
				temp = temp << ARB_SIZE ;
				temp += v.p[v.LENGTH] ;
			}
			temp = temp >> i ;               // same bits as in u
		}
		v_hat =  ARB_TYPE(temp) ;
		int a=1, b=0, c=0, d=1, q, b0 = 1, t ;

		// single precision overflow needs to be checked only at beginning 
		// for u_hat and v_hat, if it happens use multiprecision step.  
		if ((v_hat != 0 ) && (u_hat != ARB_MAX) && (v_hat != ARB_MAX))
		{
			q = (u_hat + a) / v_hat ;
			while ( b0  && (q == (u_hat + b) / (v_hat + d)) )
			{
				t = a - q * c ;
				a = c ;
				c = t ;
				t = b - q * d ;
				b = d ;
				d = t ;
				t = u_hat - q * v_hat ;
				u_hat = v_hat ;
				v_hat = t ;
				t = v_hat + c;
				if (t == 0 ) b0 = 0 ;
				else
				{
					q = (u_hat + a) / t ;
					t = v_hat + d ;
					if ( t == 0 ) b0 = 0 ;
				}
			}
		}
		if ( b == 0 )         // multiprecision step, special case
		{
			tmp = u % v ;
			u = v ;
			v = tmp ;
		}
		else                   // multiprecision step combining several
		{                      // preliminary single precision steps;			
			tmp = a * u + b * v ;
			v = c * u + d * v ;
			u = tmp ;
		}
	}
	while ( v.p[2] != 0 )   // v.LENGTH == 1 
	{		
		tmp = u % v ;
		u = v ;
		v = tmp ;
	} 
	return u ;
}