Miroirs/newrpl
1
0
Fork 0
mirror of https://git.code.sf.net/p/newrpl/sources synced 2024-09-29 17:36:57 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Finished FACTORS command

Browse source
This commit is contained in:
claudiol 2018-04-27 12:56:22 -04:00
parent 57afd69b8b
commit af368351ab
3 changed files with 174 additions and 15 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

63
newrpl/arithmetic.c
View file

@ -731,6 +731,26 @@ BINT64 gcdBINT64(BINT64 a,BINT64 b)
}
return a;
}
// COMPUTE INTEGER SQUARE ROOT
BINT64 sqrtBINT64(BINT64 num)
{
int shift=2;
BINT64 nshifted=num>>shift;
while(nshifted)
{
shift+=2;
nshifted=num>>shift;
}
shift-2;
nshifted=0;
while(shift>=0) {
nshifted<<=1;
if( (nshifted+1)*(nshifted+1) <= (num>>shift) ) nshifted+=1;
shift-=2;
}
return nshifted;
}
// RETURN ONE NON-TRIVIAL FACTOR OF n
// RETURNS n IF n IS PRIME
@ -739,19 +759,23 @@ BINT64 gcdBINT64(BINT64 a,BINT64 b)
// RETURNS THE FACTOR AS A BINT64 OR -1 AND THE RESULT IN result
BINT64 factorReal(REAL *result,REAL *n)
{
BINT64 x=2,y=2,ni, d=1;
BINT64 startnum=2,x,y,ni, d=1;
if(inBINT64Range(n)) ni=getBINT64Real(n);
else ni=-1;
do {
x=y=startnum;
do {
// do x=g(x) with g(x)=x^2+1
if(x>=(1LL<<31)) {
// SWITCH TO REALS
x=-1;
newRealFromBINT64(&RReg[8],x,0);
x=-1;
mulReal(&RReg[1],&RReg[8],&RReg[8]);
divmodReal(&RReg[8],&RReg[2],&RReg[1],n);
// HERE RReg[2]=x^2 MOD n
@ -799,8 +823,8 @@ BINT64 factorReal(REAL *result,REAL *n)
if(y>=(1LL<<31)) {
// SWITCH TO REALS
y=-1;
newRealFromBINT64(&RReg[9],y,0);
y=-1;
mulReal(&RReg[1],&RReg[9],&RReg[9]);
divmodReal(&RReg[9],&RReg[2],&RReg[1],n);
// HERE RReg[2]=y^2 MOD n
@ -844,8 +868,8 @@ BINT64 factorReal(REAL *result,REAL *n)
if(y>=(1LL<<31)) {
// SWITCH TO REALS
y=-1;
newRealFromBINT64(&RReg[9],y,0);
y=-1;
mulReal(&RReg[1],&RReg[9],&RReg[9]);
divmodReal(&RReg[9],&RReg[2],&RReg[1],n);
// HERE RReg[2]=y^2 MOD n
@ -891,15 +915,11 @@ BINT64 factorReal(REAL *result,REAL *n)
// DO d=gcd(abs(x-y),n)
if( (x<0)||(y<0)||(ni<0)) {
if(x>=0) newRealFromBINT64(&RReg[0],x,0);
if(y>=0) newRealFromBINT64(&RReg[3],y,0);
if(x>=0) newRealFromBINT64(&RReg[8],x,0);
if(y>=0) newRealFromBINT64(&RReg[9],y,0);
subReal(&RReg[4],&RReg[3],&RReg[0]);
RReg[5].flags&=~F_NEGATIVE;
// PRESERVE x AND y
swapReal(&RReg[0],&RReg[8]);
swapReal(&RReg[3],&RReg[9]);
subReal(&RReg[4],&RReg[8],&RReg[9]);
RReg[4].flags&=~F_NEGATIVE;
gcdReal(&RReg[3],&RReg[4],n);
@ -927,7 +947,26 @@ BINT64 factorReal(REAL *result,REAL *n)
}
} while(d<=1);
if(d<0) {
if(ni>0) newRealFromBINT64(&RReg[0],ni,0);
if(!eqReal(&RReg[7],&RReg[0])) break;
}
else {
if(d!=ni) break;
}
// d==n --> FAILURE, TRY ANOTHER STARTING POINT
startnum=nextprimeBINT(startnum);
} while((startnum<MAX_PRIME) && !((ni>0)&&(startnum>ni)));
if(d<0) swapReal(result,&RReg[7]);
return d;
}

1
newrpl/arithmetic.h
View file

@ -74,6 +74,7 @@ void powmodReal(REAL *result,REAL *a,REAL *b,REAL *mod);
void gcdReal(REAL *result,REAL *a,REAL *b);
BINT64 gcdBINT64(BINT64 a,BINT64 b);
BINT64 sqrtBINT64(BINT64 num);
BINT64 factorReal(REAL *result,REAL *n);

125
newrpl/lib-64-arithmetic.c
View file

@ -3615,15 +3615,98 @@ case FACTORS:
if(Context.precdigits<prec) Context.precdigits=prec;
if(Context.precdigits>MAX_USERPRECISION) Context.precdigits=MAX_USERPRECISION;
#define FACTORS_TRIVIALLIMIT 100
if(inum<0) {
// FIRST, REMOVE FIRST FEW TRIVIAL FACTORS
BINT k=2,count;
RReg[0].exp=0;
RReg[0].flags=0;
RReg[0].len=1;
while(k<FACTORS_TRIVIALLIMIT) {
RReg[0].data[0]=k;
if(gtReal(&RReg[0],&RReg[6])) break;
count=-1;
do {
divmodReal(&RReg[1],&RReg[2],&RReg[6],&RReg[0]);
swapReal(&RReg[6],&RReg[1]);
++count;
} while(iszeroReal(&RReg[2]));
swapReal(&RReg[6],&RReg[1]);
if(count>0) {
rplNewBINTPush(k,DECBINT);
rplNewBINTPush(count,DECBINT);
}
k=nextprimeBINT(k);
}
} else
{
BINT k=2,count;
BINT64 prev;
while((k<FACTORS_TRIVIALLIMIT) && (k<=inum)) {
count=-1;
do {
prev=inum;
inum/=k;
++count;
} while(inum*k==prev);
inum=prev;
if(count>0) {
rplNewBINTPush(k,DECBINT);
rplNewBINTPush(count,DECBINT);
if(Exceptions) {
DSTop=savestk;
Context.precdigits=prec;
return;
}
newRealFromBINT64(&RReg[6],inum,0);
}
k=nextprimeBINT(k);
}
}
// NOW DO POLLARD RHO
do {
// EXIT IF THE NUMBER TO FACTOR IS ONE
RReg[0].data[0]=1;
RReg[0].exp=0;
RReg[0].flags=0;
RReg[0].len=1;
if(!gtReal(&RReg[6],&RReg[0])) break;
onefactor=factorReal(&RReg[7],&RReg[6]);
if(onefactor<0) {
// CHECK FOR PERFECT SQUARES
//hyp_sqrt(&RReg[7]);
//ipReal(&RReg[1],&RReg[0],0);
//mulReal(&RReg[0],&RReg[1],&RReg[1]);
//if(eqReal(&RReg[0],&RReg[7])) swapReal(&RReg[1],&RReg[7]); // USE THE SQUARE ROOT IF IT'S A PERFECT SQUARE
rplNewRealFromRRegPush(7);
if(eqReal(&RReg[7],&RReg[6])) break;
if(eqReal(&RReg[7],&RReg[6])) {
rplPushData((WORDPTR)one_bint);
break;
}
}
else {
//BINT64 tmp=sqrtBINT64(onefactor);
//if(tmp*tmp==onefactor) onefactor=tmp;
rplNewBINTPush(onefactor,DECBINT);
if(onefactor==inum) break;
if(onefactor==inum) {
rplPushData((WORDPTR)one_bint);
break;
}
}
if(Exceptions) {
@ -3635,11 +3718,33 @@ case FACTORS:
if(inum<0) {
if(onefactor>=0) newRealFromBINT64(&RReg[7],onefactor,0);
divReal(&RReg[1],&RReg[6],&RReg[7]);
BINT count=-1;
do {
divmodReal(&RReg[1],&RReg[0],&RReg[6],&RReg[7]);
swapReal(&RReg[1],&RReg[6]);
++count;
} while(iszeroReal(&RReg[0]));
swapReal(&RReg[1],&RReg[6]);
rplNewBINTPush(count,DECBINT);
if(inBINT64Range(&RReg[6])) inum=getBINT64Real(&RReg[6]);
}
else {
BINT64 prev;
BINT count=-1;
do {
prev=inum;
inum/=onefactor;
++count;
} while(inum*onefactor==prev);
inum=prev;
rplNewBINTPush(count,DECBINT);
newRealFromBINT64(&RReg[6],inum,0);
}
@ -3647,6 +3752,20 @@ case FACTORS:
// HERE WE HAVE A LIST OF FACTORS IN THE STACK
Context.precdigits=prec;
if(Exceptions) {
DSTop=savestk;
return;
}
WORDPTR newlist=rplCreateListN(DSTop-savestk,1,1);
if(Exceptions) {
DSTop=savestk;
return;
}
DSTop=savestk;
rplOverwriteData(1,newlist);
return;
}