/*
gestion de grands nombres - juil 97 - Sylvain Huet
*/

#include "scolPrerequisites.h"

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "common/vscol.h"
#include "common/cpu_detect.h"
#if SCOL_PLATFORM == SCOL_PLATFORM_LINUX
#include <sys/time.h>
#include <unistd.h>
#endif

#include "bignum.h"
#include "base.h"




int32 CRtime(int16 *r)
{
	int32 v0,v1; 
	int32 i;
#if SCOL_PLATFORM == SCOL_PLATFORM_LINUX
	struct timeval tm;
#endif
	for(i=0;i<SIZEV;i++) r[i]=0;
#if SCOL_PLATFORM == SCOL_PLATFORM_WINDOWS
	__asm
	{
		rdtsc
			mov	v0,eax
			mov v1,edx
	}
	r[0]=v0;
	v0>>=16;
	r[1]=v0;
	r[2]=v1;
	v1>>=16;
	r[3]=v1;
#elif SCOL_PLATFORM == SCOL_PLATFORM_LINUX
	gettimeofday(&tm,NULL);
	v0=(tm.tv_sec*1000000)+(tm.tv_usec);
	r[0]=v0;
	v0>>=16;
	r[1]=v0;
#endif
	
	return 0;
}





int32 CRrand(int16 *r)
{
	int32 i;
	for(i=0;i<SIZEV;i++) r[i]=rand()&0xffff;
	r[SIZEV-1]&=0x3fff;
	return 0;
}




int32 CRfromAsc(char *buf,int16 *v)
{
	int32 i,c,l,f;
	
	CRnul(v);
	l=strlen(buf);
	if (l>SIZEV*4) return -1;
	
	f=0;
	for(i=l-1;i>=0;i--)
    {
		c=buf[i];
		if ((c>='0')&&(c<='9')) c-='0';
		else if ((c>='a')&&(c<='f')) c+=10-'a';
		else if ((c>='A')&&(c<='F')) c+=10-'A';
		else return -1;
		(*v)|=c<<f;
		f+=4;
		if (f>=16)
		{
			f=0;
			v++;
		}
    }
	return 0;
}





int32 CRfromDec(char *buf,int16 *v)
{
	int32 i,c,j,l,x;
	
	CRnul(v);
	
	l=strlen(buf);
	
	for( i=0; i<l; i++ )
    {
		c = buf[i] ;
		if( (c>='0') && (c<='9') )
		{
			c -= '0' ;
			for( j=0; j<SIZEV; j++ )
			{
				x = ( (v[j] & 0xffff) * 10 ) + c ;
				v[j] = x ;
				c =  (x>>16) & 0xffff ;
			}
		}
		else if ( i || ( c != '-' ) ) return 1 ;
    }
	
	if( buf[0] == '-' ) CRneg(v,v) ;
	
	return 0 ;
}



#define SIZED 64

/* effectue l'operation *4+x */
int32 CRmul4(char *buf,int32 c)
{
	int32 i,x;
	
	for(i=0;i<SIZED;i++)
	{
		x=buf[i]*4+c;
		c=0;
		while(x>9)
		{
			c++;
			x-=10;
		}
		buf[i]=x;
	}
	return 0;
}



int32 CRtoDec(int16 *u,char *dst)
{
	int32 i,l,x;
	char buf[SIZED];
	int16 v[SIZEV];

	if( u[SIZEV-1] & 0x8000 ) 
	{
		*(dst++) = '-' ;
		CRneg(u,v) ;
	}
	else CRcopy(u,v);
	
	for( i=0; i<SIZED; i++ ) buf[i]=0;
	
	for(i=SIZEV-1;i>=0;i--)
    {
		x=v[i];
		CRmul4(buf, 0 ) ;
		CRmul4(buf, (x>>12) & 15 ) ;
		CRmul4(buf, 0 ) ;
		CRmul4(buf, (x>>8) & 15 ) ;
		CRmul4(buf, 0 ) ;
		CRmul4(buf, (x>>4) & 15 ) ;
		CRmul4(buf, 0 ) ;
		CRmul4(buf, x & 15 ) ;
	}
	
	l = SIZED-1 ;
	
	while( (l>0 ) && (buf[l]==0) ) l-- ;
	while( l>=0 ) *(dst++) = buf[l--]+48 ;
	*dst = 0 ;
	return 0 ;
}



char *CRtoAsc(char *buf,int16 *v)
{
	int32 i,j,f,x,y;
	
	f=0;
	for(i=SIZEV-1;i>=0;i--)
    {
		x=v[i]&0xffff;
		for(j=12;j>=0;j-=4)
		{
			y=(x>>j)&15;
			if ((y)||(f))
			{
				if (y<10) y+='0'; else y+='a'-10;
				f=1;
				*(buf++)=y;
			}
		}
    }
	if (!f) *(buf++)='0';
	*buf=0;
	return buf;
}




int32 CRfromString(char *buf,int16 *v,int32 n)
{
	int32 i;
	
	CRnul(v);
	
	if (n>SIZEV*2) return 0;
	
	for( i=0; i<n; i+=2 ) v[i>>1] = (buf[i+1] << 8) + (buf[i] & 0xFF) ;
	if( i>n ) v[(n-1)>>1] &= 0xFF ;
	
	return 0;
}



char* CRtoString(char *buf,int16 *v,int32 *n)
{
	int32 i,x;
	*n=0;
	for(i=0;i<SIZEV;i++)
    {
		x=v[i]&0xffff;
		if (x&255) *n=i+i+1;
		buf[i+i]=x;
		x>>=8;
		if (x&255) *n=i+i+2;
		buf[i+i+1]=x;
    }
	buf[*n]=0;
	return buf;
}



int32 CRfromInt(int32 x,int16 *v)
{
	CRnul(v);
	
	v[0] = x & 0xffff ;
	v[1] = x>>16 ;
	
	return 0;
}



int32 CRtoInt(int16 *v)
{
	int32 a,b;
	
	a=v[0]&0xffff;
	b=v[1]&0xffff;
	return a+(b<<16);
}

char* CRtoStringn(char *buf,int16 *v,int32 n)
{
	int32 i,x;
	for( i=0; i<SIZEV; i++ )
    {
		x = v[i]&0xffff ;
		buf[i+i] = x ;
		x >>= 8 ;
		buf[i+i+1] = x ;
    }
	buf[n] = 0 ;
	return buf ;
}




int32 CRsar(int16 *p, int16 *r)
{
	int32 i,cc,x,z;
	
	cc=p[SIZEV-1]&0x8000;
	for(i=SIZEV-1;i>=0;i--)
    {
		x=p[i]&0xffff;
		z=cc|(x>>1);
		r[i]=z;
		cc=(x&1)<<15;
    }
	return 0;
}




int32 CRg(int16 *p,int16 *q)
{
	int16 v[SIZEV];
	
	CRsub(p,q,v);
	if (v[SIZEV-1]&0x8000) return 0;
	return 1;
}




int32 CRl(int16 *p,int16 *q)
{
	int16 v[SIZEV];
	
	CRsub(q,p,v);
	if (v[SIZEV-1]&0x8000) return 0;
	return 1;
}





int32 CRisone(int16 *v)
{
	int32 i;
	if (v[0]!=1) return 0;
	for(i=1;i<SIZEV;i++) if (v[i]) return 0;
	return 1;
}




int32 CRisminusone(int16 *v)
{
	int32 i,x;
	for(i=0;i<SIZEV;i++)
    {
		x=v[i]&0xffff;
		if (x!=0xffff) return 0;
    }
	return 1;
}




int32 CRestpair(int16 *x)
{
	return !(x[0]&1);
}





int32 CRinv(int16 *p, int16 *q, int16 *r)
{
	int16 x[SIZEV];
	int16 y[SIZEV];
	int16 g[SIZEV];
	int16 a0[SIZEV];
	int16 a1[SIZEV];
	int16 a2[SIZEV];
	
	CRcopy(q,x);
	CRcopy(p,y);
	CRcopy(p,g);
	CRnul(a0); a0[0]=1;
	CRnul(a1);
	while(CRnotnull(x))
    {
		CRdiv(y,x,a2);
		CRmuln(a2,a1,q,a2);
		CRsub(a0,a2,a0); 
		CRmodn(a0,q,a0);
		CRcopy(a0,a2); CRcopy(a1,a0); CRcopy(a2,a1);
		
		CRcopy(x,g);
		CRmodn(y,x,x);
		CRcopy(g,y);
    }
	CRcopy(a0,r);
	return 0;
}




int32 CRrandPrimal(int16 *r,int32 i0)
{
	int32 i,m,k;
	
	i0--;
	if ((i0<=0)||(i0>SIZEV*16-2))
    {
		CRnul(r);
		return 0;
    }
	do
    {
		CRrand(r);
		r[0]|=1;
    }
	while(CRisone(r));
	
	k=i0>>4;
	m=1<<(i0&15);
	
	r[k]&=m-1; r[k]|=m;
	for(i=k+1;i<SIZEV;i++) r[i]=0;
	return 0;
}




int32 CRmillerRabin(int16 *n)
{
	int16 a[SIZEV];
	int16 nm[SIZEV];
	int16 t[SIZEV];
	int32 s;
	
	if (CRestpair(n)) return 0;
	CRcopy(n,nm);
	CRnul(a); a[0]=1;
	CRsub(nm,a,nm);
	
	do
    {
		CRrand(a);
		CRmodn(a,n,a);
    }
	while(!CRnotnull(a));
	
	CRcopy(nm,t);
	s=0;
	while(CRestpair(t))
    {
		CRshr(t,t);
		s++;
    }
	
	CRexpn(a,t,n,a);
	if (CRisone(a)) return 1;
	if (CRe(a,nm)) return 1;
	
	while(s)
    {
		s--;
		CRmuln(a,a,n,a);
		if (CRisone(a)) return 0;
		if (CRe(a,nm)) return 1;
    }
	return 0;
}




int32 CRfindPrimal(int16 *u,int32 i0)
{
	int32 i;
	int16 trois[SIZEV];
	int16 un[SIZEV];
	int16 v[SIZEV];

	if( i0 <= 0 ) return 1 ;
	else if( i0 == 1 ) 
	{
		CRnul(u) ; u[0]=1 ;
		return 0 ;
	}
	else
	{
		CRnul(un); un[0]=1;
		CRnul(trois); trois[0]=3;
		
		while(1)
		{
			CRrandPrimal(u,i0);
			CRsub(u,un,v);
			CRpgcd(v,trois,v);
			if (CRisone(v))
			{
				i=16;
				while((i)&&(CRmillerRabin(u))) i--;
				if (i==0) return 0;
			}
		}
		return 0;
	}
}




int32 CRnotnull(int16 *x)
{
	int16 i ;
	for( i=0; i<SIZEV; i++ ) if(x[i]) return 1;
	return 0;
}







int32 CRadd(int16 *p, int16 *q, int16 *r)
{
#ifdef OPTI_P4
	if( infoP4 )
	{
		_asm
		{
			pushad

			xor ecx, ecx
			pxor mm0, mm0
			mov eax, r
			mov esi, p
			mov edi, q

		boucle: 		
			cmp ecx, 4
			je end

			movd mm1, dword ptr [esi+ecx*4]
			movd mm2, dword ptr [edi+ecx*4]
			
			psrlq mm0, 32
			paddq mm0, mm1
			paddq mm0, mm2

			movd dword ptr [eax+ecx*4], mm0
			
			inc ecx
			jmp boucle

		end:	
			popad
		}
	}
	else
#endif
	{
		int16 i ;
		int32 cc ;
		
		cc = 0 ;
		
		for( i=0; i<SIZEV; i++ )
		{
			cc += ( (p[i] & 0xFFFF) + (q[i] & 0xFFFF) ) ;
			r[i] = cc ;
			cc >>= 16 ;
		}
	}
	
	return 0 ;
}



int32 CRsub(int16 *p, int16 *q, int16 *r)
{
#ifdef OPTI_P4
	if( infoP4 )
	{
		uint64 mask = 0xFFFFFFFF00000000 ;

		_asm
		{
			pushad

			xor ecx, ecx
			pxor mm0, mm0
			mov esi, p
			mov edi, q

		boucle: 		
			cmp ecx, 4
			je end

			movd mm1, dword ptr [esi+ecx*4]
			movd mm2, dword ptr [edi+ecx*4]
			movq mm3, qword ptr mask
			
			movq mm4, mm0
			psrlq mm0, 32
			psrlq mm4, 63
			
			movd eax, mm4
			cmp eax, 0x00000000
			je positif
			
			por mm0, mm3

		positif:
			paddq mm0, mm1
			psubq mm0, mm2

			mov eax, r
			movd dword ptr [eax+ecx*4], mm0
			
			inc ecx
			jmp boucle

		end:	
			popad
		}
	}
	else
#endif
	{
		int32 i,cc;
		
		cc = 0 ;
		
		for( i=0; i<SIZEV; i++ )
		{
			cc += ( (p[i] & 0xFFFF) - (q[i] & 0xFFFF) ) ;
			r[i] = cc ;
			cc >>= 16 ;
		}
	}
	
	return 0 ;
}






int32 CRneg(int16 *p,int16 *r)
{
#ifdef OPTI_P4
	if( infoP4 )
	{
		__my128i mask = { 0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF } ;
		__my128i un   = { 0x0000000000000001, 0x0000000000000000 } ;	
		
		_asm
		{
			pushad

			xor ecx, ecx
			mov edi, r
			mov esi, p

			movdqu xmm0, xmmword ptr [esi]
			movdqa xmm1, xmmword ptr mask
			movdqa xmm2, xmmword ptr un
			andnpd xmm0, xmm1
			
			movdqa xmm3, xmm0
			
		boucle:
			cmp ecx, 4
			je end

			movd eax, xmm3
			pshufd xmm3, xmm3, _MM_SHUFFLE(0,3,2,1)
			inc ecx

			cmp eax, 0xFFFFFFFF
			je decalage
			jmp end
			
		decalage:
			pslldq xmm1, 4
			pslldq xmm2, 4
			jmp boucle

		end:	
			
			paddq xmm0, xmm2
			pand xmm0, xmm1
			
			movdqu xmmword ptr [edi], xmm0

			popad
		}
	}
	else
#endif
	{
		int16 n[SIZEV];
	
		CRnul(n);
		CRsub(n,p,r);
	}

	return 0;
}





int32 CRnul(int16 *r)
{
#ifdef OPTI_P4
	if( infoP4 )
	{
		_asm
		{
			pushad

			pxor xmm0, xmm0
			mov eax, r
			movdqu xmmword ptr [eax], xmm0	

			popad
		}
	}
	else
#endif
	{
		int16 i;
		for(i=0;i<SIZEV;i++) r[i]=0;
	}

	return 0 ;
}





int32 CRxor(int16 *p, int16 *q, int16 *r)
{
#ifdef OPTI_P4
	if( infoP4 )
	{
		_asm
		{
			pushad

			mov eax, r
			mov esi, p
			mov edi, q

			movdqu xmm0, xmmword ptr [esi]
			movdqu xmm1, xmmword ptr [edi]
			pxor xmm0, xmm1
			movdqu xmmword ptr [eax], xmm0

			popad
		}
	}
	else
#endif
	{
		int16 i;
		for(i=0;i<SIZEV;i++) r[i]=p[i]^q[i] ;
	}

	return 0 ;
}





int32 CRcopy(int16 *p, int16 *r)
{
	int16 i ;
	for( i=0; i<SIZEV; i++ ) r[i] = p[i] ;

	return 0 ;
}





int32 CRe(int16 *p, int16 *q)
{
#ifdef OPTI_P4
	if( infoP4 )
	{
		int res = 1 ;

		_asm
		{
			pushad

			xor ecx, ecx
			mov esi, p
			mov edi, q

			movdqu xmm0, xmmword ptr [esi]
			movdqu xmm1, xmmword ptr [edi]
			movdqu xmm2, xmm0

			pcmpeqd xmm2, xmm1
			
		boucle:
			cmp ecx, 4
			je end
			
			movd eax, xmm2
			cmp eax, 0
			je zero

			pshufd xmm2, xmm2, _MM_SHUFFLE(0,3,2,1)
			inc ecx
			jmp boucle

		zero:
			mov res, 0

		end:
			popad
		}

		return res ;
	}
	else
#endif
	{
		int16 i ;
		
		for( i=0; i<SIZEV; i++ ) if( p[i]!=q[i] ) return 0 ;

		return 1 ;
	}
}	





int CRgu( uint16 *p, uint16 *q )
{
	int16 i ;
	for(i=SIZEV-1;i>=0;i--)
    {
		if (p[i]>q[i]) return 1;
		if (p[i]<q[i]) return 0;
    }
	return 1;
}




int32 CRcmp( uint16 *p, uint16 *q )
{
	if( CRpositif(p) )
	{
		if( CRpositif(q) )
		{
			int16 i ;
			for( i=SIZEV-1; i>=0; i-- )
			{
				if (p[i]>q[i]) return 1;
				if (p[i]<q[i]) return -1;
			}
			return 0;
		}
		else return 1 ;
	}
	else
	{
		if( CRpositif(q) ) return -1 ;
		else
		{
			int16 i ;
			for( i=SIZEV-1; i>=0; i-- )
			{
				if (p[i]>q[i]) return 1;
				if (p[i]<q[i]) return -1;
			}
			return 0;
		}
	}
}
		

		



int32 CRpositif(int16 *v)
{
	if (v[SIZEV-1]&0x8000) return 0;
	return 1;
}






int32 CRshr(int16 *p, int16 *r)
{
#ifdef OPTI_P4
	__my128i mask1 = { 0x0000000000000000, 0x0000000000000001 } ;
	__my128i mask2 = { 0x0000000000000000, 0x8000000000000000 } ;

	if( infoP4 )
	{
		_asm
		{
			pushad

			mov eax, r
			mov esi, p
			movdqu xmm0, xmmword ptr [esi]
			movdqa xmm1, xmmword ptr mask1
			movdqa xmm2, xmmword ptr mask2
			
			movdqu xmm3, xmm0
			pand xmm3, xmm1
			psrldq xmm3, 1
			psllq xmm3, 7

			movdqu xmm4, xmm0
			pand xmm4, xmm2
			por xmm3, xmm4

			psrlq xmm0, 1
			por xmm0, xmm3

			movdqu xmmword ptr [eax], xmm0

			popad
		}
	}
	else
#endif
	{	
		int16 i ;
		int32 cc,x,z ;
		
		cc = 0 ;
		for( i=SIZEV-1; i>=0; i-- )
		{
			x = p[i]&0xffff ;
			z = cc|(x>>1) ;
			r[i] = z ;
			cc = (x&1)<<15 ;
		}
	}

	return 0;
}








int32 CRshl(int16 *p, int16 *r)
{
/*#ifdef OPTI_P4
	__my128i mask = { 0x8000000000000000, 0x0000000000000000 } ;

	if( infoP4 )
	{
		_asm
		{
			pushad

			mov eax, r
			mov esi, p
			movdqu xmm0, xmmword ptr [esi]
			movdqu xmm1, xmmword ptr mask
			
			movdqu xmm2, xmm0
			pand xmm2, xmm1
			pslldq xmm2, 1
			psrlq xmm2, 7
			psllq xmm0, 1
			por xmm0, xmm2

			movdqu xmmword ptr [eax], xmm0

			popad
		}
	}
	else
#endif*/
	{
		int32 i,x,cc ;
	
		cc=0;
		for(i=0;i<SIZEV;i++)
		{
			x=p[i]&0xffff;
			cc|=x<<1;
			r[i]=cc;
			cc>>=16;
		}
	}

	return 0;
}





int CRmuln(int16 *p, int16 *q, int16 *n, int16 *res)
{
	int16 m[SIZEV];
	int16 x[SIZEV];
	int16 r[SIZEV];
	int b,i;
	
	CRmodn(p,n,x);
	CRmodn(q,n,m);
	CRnul(r);
	
	for(i=0;i<SIZEV;i++)
    {
		b=1;
		while(b<0x10000)
		{
			if (x[i]&b)
			{
				CRadd(r,m,r);
				if (CRgu(r,n)) CRsub(r,n,r);
			}
			CRshl(m,m);
			if (CRgu(m,n)) CRsub(m,n,m);
			b<<=1;
		}
    }
	CRcopy(r,res);
	
	return 0;
}





int32 CRmul(int16 *p, int16 *q, int16 *r )
{
	int signe = 1 ;

	if( p[SIZEV-1]&0x8000)
	{
		signe = -signe ;
		CRneg( p, p ) ;
	}

	if( q[SIZEV-1]&0x8000)
	{
		signe = -signe ;
		CRneg( q, q) ;
	}

#ifdef OPTI_P4
	if( infoP4 )
	{
		__my128i t[4] ;

		_asm
		{
			pushad

			xor ecx, ecx
			mov esi, p
			mov edi, q
			lea eax, t
			
			pxor xmm2, xmm2								// xmm2 = 0 0 0 0
			movq xmm2, qword ptr [edi]					// xmm2 = 0 0 x y
			pshufd xmm2, xmm2, _MM_SHUFFLE(3,1,2,0)     // xmm2 = 0 x 0 y

			pxor xmm3, xmm3								// xmm3 = 0 0 0 0	
			movq xmm3, qword ptr [edi+8]				// xmm3 = 0 0 x y
			pshufd xmm3, xmm3, _MM_SHUFFLE(3,1,2,0)		// xmm3 = 0 x 0 y

			pxor xmm4, xmm4
			pxor xmm5, xmm5


		mult:
			cmp ecx, 4
			je end

			movd xmm0, dword ptr [esi+ecx*4]			// xmm0 = 0 0 0 x
			pshufd xmm0, xmm0, _MM_SHUFFLE(1,0,1,0)		// xmm0 = 0 x 0 x
			movdqa xmm1, xmm0							// xmm1 = 0 x 0 x

			pmuludq xmm0, xmm2
			pmuludq xmm1, xmm3

			paddq xmm0, xmm4
			paddq xmm1, xmm5
		
			movd dword ptr [eax+ecx*4], xmm0			// xmm0 = a b c d
			pshufd xmm0, xmm0, _MM_SHUFFLE(1,0,3,2)		// xmm0 = c d a b
			movd dword ptr [eax+ecx*4+16], xmm0
			pshufd xmm0, xmm0, _MM_SHUFFLE(1,0,3,2)		// xmm0 = a b c d

			psrlq xmm0, 32								// xmm0 = 0 b 0 c
			movdqa xmm4, xmm0							// xmm4 = 0 b 0 c

			
			movd dword ptr [eax+ecx*4+32], xmm1			// xmm1 = a b c d
			pshufd xmm1, xmm1, _MM_SHUFFLE(1,0,3,2)		// xmm1 = c d a b
			movd dword ptr [eax+ecx*4+48], xmm1
			pshufd xmm1, xmm1, _MM_SHUFFLE(1,0,3,2)		// xmm1 = a b c d

			psrlq xmm1, 32								// xmm1 = 0 b 0 c
			movdqa xmm5, xmm1							// xmm5 = 0 b 0 c
		
			inc ecx
			jmp mult

			
		end:		
			movdqa xmm1, xmmword ptr [eax+16]
			pslldq xmm1, 4
			movdqa xmmword ptr [eax+16], xmm1

			movdqa xmm2, xmmword ptr [eax+32]
			pslldq xmm2, 8
			movdqa xmmword ptr [eax+32], xmm2

			movdqa xmm3, xmmword ptr [eax+48]
			pslldq xmm3, 12
			movdqa xmmword ptr [eax+48], xmm3

			popad
		}

		CRadd( (int16*)t, (int16*)(t+1), r ) ;
		CRadd( (int16*)(t+2), r, r ) ;
		CRadd( (int16*)(t+3), r, r ) ;
	}
	else
#endif
	{
		int16 t[SIZEV] ;
		int16 u[SIZEV] ;
		int16 i,j ;
		uint32 cc ;
		
		CRnul( u ) ;

		for( i=0; i<SIZEV; i++ )
		{
			cc = 0 ;
			CRnul( t ) ;
			for( j=0; j<SIZEV-i; j++ )
			{
				cc += ( (p[i] & 0xFFFF) * (q[j] & 0xFFFF) ) ;
				t[j+i] = cc & 0xFFFF ;
				cc >>= 16 ;
			}
			CRadd( t, u, u ) ;
		}
		CRcopy( u, r ) ;
	}

	if( signe == -1 ) CRneg( r, r ) ;

	return 0 ;
}






int32 CRmodn(int16 *p, int16 *q, int16 *r)
{
	int16 v[SIZEV], w[SIZEV], qs[SIZEV] ;
	int16 n, f ;

	CRcopy(q,qs);
	CRcopy(q,v);
	CRcopy(p,w);
	
	f=0;

	if (!CRpositif(w))
	{
		f=1;
		CRneg(w,w);
	}
	
	n=1;

	while(CRgu(w,v))
	{
		CRshl(v,v);
		n++;
	}

	CRcopy(w,r);

	while(n)
	{
		if (CRgu(r,v)) CRsub(r,v,r);
		CRshr(v,v);
		n--;
	}

	if(f)
	{
		CRneg(r,r);
		CRadd(r,qs,r);
	}

	return 0;
}




int32 CRpgcd(int16 *p, int16 *q, int16 *r)
{
	int16 x[SIZEV];
	int16 y[SIZEV];
	int16 g[SIZEV];
	
	CRcopy(p,x);
	CRcopy(q,y);
	CRcopy(q,g);

	while(CRnotnull(x))
    {
		CRcopy(x,g);
		CRmodn(y,x,x);
		CRcopy(g,y);
    }
	CRcopy(g,r);

	return 0;
}



int32 CRexpn(int16 *p, int16 *e, int16 *n, int16 *res)
{
	int16 m[SIZEV];
	int16 r[SIZEV];
	int32 b,i;
	
	CRmodn(p,n,m);
	CRnul(r); r[0]=1;
	for(i=0;i<SIZEV;i++)
    {
		b=1;
		while(b<0x10000)
		{
			if( e[i]&b ) CRmuln(r,m,n,r);
			CRmuln(m,m,n,m) ;
	
			b<<=1;
		}
    }
	CRcopy(r,res);

	return 0;
}



int32 CRdiv(int16 *p, int16 *q, int16 *r)
{
	int16 v[SIZEV];
	int16 m[SIZEV];
	int32 n;

	if( CRnotnull( q ) )
	{
		int signe = 1 ;

		if( p[SIZEV-1]&0x8000)
		{
			signe = -signe ;
			CRneg( p, p ) ;
		}

		if( q[SIZEV-1]&0x8000)
		{
			signe = -signe ;
			CRneg( q, q) ;
		}

		CRcopy(q,v);
		n=1;

		while(CRgu(p,v))
		{
			CRshl(v,v);
			n++;
		}
		CRcopy(p,m);
		CRnul(r);
		
		while(n)
		{
			CRshl(r,r);
			if (CRgu(m,v))
			{
				CRsub(m,v,m);
				r[0]+=1;
			}
			CRshr(v,v);
			n--;
		}

		if( signe == -1 ) CRneg( r, r ) ;

		return 0 ;
	}
	else return 1 ;
}