Improved sqrt algorithm.

This commit is contained in:
claudiol 2018-04-21 08:43:52 -04:00
parent 298fae296d
commit 1811119738
3 changed files with 35 additions and 3 deletions

View file

@ -677,3 +677,30 @@ void powmodReal(REAL *result,REAL *a,REAL *b,REAL *mod)
// RESULT IS IN result // RESULT IS IN result
} }
// COMPUTE GCD OF 2 REALS
// USE DIVISION OR SUBTRACTION, WHATEVER IS FASTER ON EACH LOOP
/*
void gcdReal(REAL *result,REAL *a,REAL *b)
{
copyReal(&RReg[0],a);
copyReal(&RReg[1],b);
RReg[0].flags&=~F_NEGATIVE;
RReg[1].flags&=~F_NEGATIVE;
BINT idxa;
if(gtReal(a,b)) idxa=0; else idxa=1;
#define REG_a RReg[idxa]
#define REG_b RReg[idxa^1]
while( (!iszeroReal(REG_b))&&(!eqReal(REG_a,REG_b)))
{
}
}
*/

View file

@ -84,7 +84,7 @@
CMD(DIGITS,MKTOKENINFO(6,TITYPE_FUNCTION,3,2)), \ CMD(DIGITS,MKTOKENINFO(6,TITYPE_FUNCTION,3,2)), \
CMD(PROOT,MKTOKENINFO(5,TITYPE_FUNCTION,1,2)), \ CMD(PROOT,MKTOKENINFO(5,TITYPE_FUNCTION,1,2)), \
CMD(PREVPRIME,MKTOKENINFO(9,TITYPE_FUNCTION,1,2)), \ CMD(PREVPRIME,MKTOKENINFO(9,TITYPE_FUNCTION,1,2)), \
CMD(FACTORS,MKTOKENINFO(7,TITYPE_FUNCTION,1,2))) CMD(FACTORS,MKTOKENINFO(7,TITYPE_FUNCTION,1,2))

View file

@ -327,9 +327,10 @@ RReg[2].exp=0;
RReg[2].flags=0; RReg[2].flags=0;
RReg[2].len=1; RReg[2].len=1;
int savedprec=Context.precdigits; int savedprec=Context.precdigits;
int maxprec=(Context.precdigits+15)&~7;
// Halley's method // Halley's method
Context.precdigits=(Context.precdigits+15)&~7; Context.precdigits=8;
int iters=0; int iters=0;
int goodexp,gooddigits=0; int goodexp,gooddigits=0;
do { do {
@ -352,10 +353,14 @@ mulReal(&RReg[4],&RReg[2],&RReg[3]); // x(n+1)=0.125*xn*(15-yn*(10-3*yn))
swapReal(&RReg[4],&RReg[1]); swapReal(&RReg[4],&RReg[1]);
subReal(&RReg[3],&RReg[4],&RReg[1]); subReal(&RReg[3],&RReg[4],&RReg[1]);
if(RReg[3].len>1) continue; if(RReg[3].len>1) continue;
if(RReg[3].data[0]==0) break;
if((Context.precdigits==maxprec) && (RReg[3].data[0]==0)) break;
goodexp=RReg[3].exp+sig_digits(RReg[3].data[0]); goodexp=RReg[3].exp+sig_digits(RReg[3].data[0]);
gooddigits=intdigitsReal(&RReg[1])-RReg[1].exp; gooddigits=intdigitsReal(&RReg[1])-RReg[1].exp;
gooddigits-=goodexp-RReg[1].exp; gooddigits-=goodexp-RReg[1].exp;
Context.precdigits+=Context.precdigits; // DOUBLE THE PRECISION AT EACH STEP
if(Context.precdigits<8) Context.precdigits=8;
if(Context.precdigits>maxprec) Context.precdigits=maxprec;
} while(gooddigits<=savedprec); } while(gooddigits<=savedprec);