From af368351ab761895aa166011986fbbdb8f9e71d8 Mon Sep 17 00:00:00 2001 From: claudiol Date: Fri, 27 Apr 2018 12:56:22 -0400 Subject: [PATCH] Finished FACTORS command --- newrpl/arithmetic.c | 63 +++++++++++++++---- newrpl/arithmetic.h | 1 + newrpl/lib-64-arithmetic.c | 125 ++++++++++++++++++++++++++++++++++++- 3 files changed, 174 insertions(+), 15 deletions(-) diff --git a/newrpl/arithmetic.c b/newrpl/arithmetic.c index 0863314..1b6cda1 100644 --- a/newrpl/arithmetic.c +++ b/newrpl/arithmetic.c @@ -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((startnum0)&&(startnum>ni))); + + + + if(d<0) swapReal(result,&RReg[7]); + return d; } diff --git a/newrpl/arithmetic.h b/newrpl/arithmetic.h index a64d667..0f9423e 100644 --- a/newrpl/arithmetic.h +++ b/newrpl/arithmetic.h @@ -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); diff --git a/newrpl/lib-64-arithmetic.c b/newrpl/lib-64-arithmetic.c index a4c8612..68c2549 100644 --- a/newrpl/lib-64-arithmetic.c +++ b/newrpl/lib-64-arithmetic.c @@ -3615,15 +3615,98 @@ case FACTORS: if(Context.precdigitsMAX_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(k0) { + rplNewBINTPush(k,DECBINT); + rplNewBINTPush(count,DECBINT); + } + k=nextprimeBINT(k); + } + + } else + { + BINT k=2,count; + BINT64 prev; + while((k0) { + 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; }