2021SC@SDUSC
mrxgcd.c实现扩展最大公约数模块。
static mr_small qdiv(mr_large u,mr_large v)
{ /* fast division - small quotient expected. */
mr_large lq,x=u;
#ifdef MR_FP
mr_small dres;
#endif
x-=v;
if (x<v) return 1;
x-=v;
if (x<v) return 2;
x-=v;
if (x<v) return 3;
x-=v;
if (x<v) return 4;
x-=v;
if (x<v) return 5;
x-=v;
if (x<v) return 6;
x-=v;
if (x<v) return 7;
x-=v;
if (x<v) return 8;
/* do it the hard way! */
lq=8+MR_DIV(x,v);
if (lq>=MAXBASE) return 0;
return (mr_small)lq;
}
#else
实现快除法,小商预期。
int xgcd(_MIPD_ big x,big y,big xd,big yd,big z)
{ /* greatest common divisor by Euclids method *
* extended to also calculate xd and yd where *
* z = x.xd + y.yd = gcd(x,y) *
* if xd, yd not distinct, only xd calculated *
* z only returned if distinct from xd and yd *
* xd will always be positive, yd negative */
int s,n,iter;
mr_small r,a,b,c,d;
mr_small q,m,sr;
#ifdef MR_FP
mr_small dres;
#endif
#ifdef mr_dltype
union doubleword uu,vv;
mr_large u,v,lr;
#else
mr_small u,v,lr;
#endif
BOOL last,dplus=TRUE;
big t;
#ifdef MR_OS_THREADS
miracl *mr_mip=get_mip();
#endif
if (mr_mip->ERNUM) return 0;
MR_IN(30)
#ifdef MR_COUNT_OPS
fpx++;
#endif
copy(x,mr_mip->w1);
copy(y,mr_mip->w2);
s=exsign(mr_mip->w1);
insign(PLUS,mr_mip->w1);
insign(PLUS,mr_mip->w2);
convert(_MIPP_ 1,mr_mip->w3);
zero(mr_mip->w4);
last=FALSE;
a=b=c=d=0;
iter=0;
while (size(mr_mip->w2)!=0)
{
if (b==0)
{ /* update mr_mip->w1 and mr_mip->w2 */
divide(_MIPP_ mr_mip->w1,mr_mip->w2,mr_mip->w5);
t=mr_mip->w1,mr_mip->w1=mr_mip->w2,mr_mip->w2=t; /* swap(mr_mip->w1,mr_mip->w2) */
multiply(_MIPP_ mr_mip->w4,mr_mip->w5,mr_mip->w0);
add(_MIPP_ mr_mip->w3,mr_mip->w0,mr_mip->w3);
t=mr_mip->w3,mr_mip->w3=mr_mip->w4,mr_mip->w4=t; /* swap(xd,yd) */
iter++;
}
else
{
/* printf("a= %I64u b= %I64u c= %I64u d= %I64u \n",a,b,c,d); */
mr_pmul(_MIPP_ mr_mip->w1,c,mr_mip->w5); /* c*w1 */
mr_pmul(_MIPP_ mr_mip->w1,a,mr_mip->w1); /* a*w1 */
mr_pmul(_MIPP_ mr_mip->w2,b,mr_mip->w0); /* b*w2 */
mr_pmul(_MIPP_ mr_mip->w2,d,mr_mip->w2); /* d*w2 */
if (!dplus)
{
mr_psub(_MIPP_ mr_mip->w0,mr_mip->w1,mr_mip->w1); /* b*w2-a*w1 */
mr_psub(_MIPP_ mr_mip->w5,mr_mip->w2,mr_mip->w2); /* c*w1-d*w2 */
}
else
{
mr_psub(_MIPP_ mr_mip->w1,mr_mip->w0,mr_mip->w1); /* a*w1-b*w2 */
mr_psub(_MIPP_ mr_mip->w2,mr_mip->w5,mr_mip->w2); /* d*w2-c*w1 */
}
mr_pmul(_MIPP_ mr_mip->w3,c,mr_mip->w5);
mr_pmul(_MIPP_ mr_mip->w3,a,mr_mip->w3);
mr_pmul(_MIPP_ mr_mip->w4,b,mr_mip->w0);
mr_pmul(_MIPP_ mr_mip->w4,d,mr_mip->w4);
if (a==0) copy(mr_mip->w0,mr_mip->w3);
else mr_padd(_MIPP_ mr_mip->w3,mr_mip->w0,mr_mip->w3);
mr_padd(_MIPP_ mr_mip->w4,mr_mip->w5,mr_mip->w4);
}
if (mr_mip->ERNUM || size(mr_mip->w2)==0) break;
n=(int)mr_mip->w1->len;
if (n==1)
{
last=TRUE;
u=mr_mip->w1->w[0];
v=mr_mip->w2->w[0];
}
else
{
m=mr_mip->w1->w[n-1]+1;
#ifndef MR_SIMPLE_BASE
if (mr_mip->base==0)
{
#endif
#ifndef MR_NOFULLWIDTH
#ifdef mr_dltype
/* use double length type if available */
if (n>2 && m!=0)
{ /* squeeze out as much significance as possible */
uu.h[MR_TOP]=muldvm(mr_mip->w1->w[n-1],mr_mip->w1->w[n-2],m,&sr);
uu.h[MR_BOT]=muldvm(sr,mr_mip->w1->w[n-3],m,&sr);
vv.h[MR_TOP]=muldvm(mr_mip->w2->w[n-1],mr_mip->w2->w[n-2],m,&sr);
vv.h[MR_BOT]=muldvm(sr,mr_mip->w2->w[n-3],m,&sr);
}
else
{
uu.h[MR_TOP]=mr_mip->w1->w[n-1];
uu.h[MR_BOT]=mr_mip->w1->w[n-2];
vv.h[MR_TOP]=mr_mip->w2->w[n-1];
vv.h[MR_BOT]=mr_mip->w2->w[n-2];
if (n==2) last=TRUE;
}
u=uu.d;
v=vv.d;
#else
if (m==0)
{
u=mr_mip->w1->w[n-1];
v=mr_mip->w2->w[n-1];
}
else
{
u=muldvm(mr_mip->w1->w[n-1],mr_mip->w1->w[n-2],m,&sr);
v=muldvm(mr_mip->w2->w[n-1],mr_mip->w2->w[n-2],m,&sr);
}
#endif
#endif
#ifndef MR_SIMPLE_BASE
}
else
{
#ifdef mr_dltype
if (n>2)
{ /* squeeze out as much significance as possible */
u=muldiv(mr_mip->w1->w[n-1],mr_mip->base,mr_mip->w1->w[n-2],m,&sr);
u=u*mr_mip->base+muldiv(sr,mr_mip->base,mr_mip->w1->w[n-3],m,&sr);
v=muldiv(mr_mip->w2->w[n-1],mr_mip->base,mr_mip->w2->w[n-2],m,&sr);
v=v*mr_mip->base+muldiv(sr,mr_mip->base,mr_mip->w2->w[n-3],m,&sr);
}
else
{
u=(mr_large)mr_mip->base*mr_mip->w1->w[n-1]+mr_mip->w1->w[n-2];
v=(mr_large)mr_mip->base*mr_mip->w2->w[n-1]+mr_mip->w2->w[n-2];
last=TRUE;
}
#else
u=muldiv(mr_mip->w1->w[n-1],mr_mip->base,mr_mip->w1->w[n-2],m,&sr);
v=muldiv(mr_mip->w2->w[n-1],mr_mip->base,mr_mip->w2->w[n-2],m,&sr);
#endif
}
#endif
}
最大公约数欧几里得算法