Miroirs/newrpl
1
0
Fork 0
mirror of https://git.code.sf.net/p/newrpl/sources synced 2024-11-16 19:51:25 +01:00
Code Issues Projects Releases Packages Wiki Activity Actions

Added atan argument reduction.

Browse source
This commit is contained in:
claudio 2017-06-19 14:31:27 -04:00
parent ddfd119b55
commit 6edd0834a7
1 changed files with 75 additions and 29 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

104
tools/decimal-optimization/main.c
View file

@ -1290,7 +1290,7 @@ void trig_tanpower(REAL *angle, BINT angmode)
void atanpower()
{
int neg,invert;
int neg,invert,pass;
int savedprec=Context.precdigits;
int k;
// ARGUMENT REDUCTION
@ -1303,26 +1303,62 @@ void atanpower()
// MAKE ANGLE SMALLER TO IMPROVE CONVERGENCE SPEED
// USE THE IDENTITY ATAN(X) = 2 TAN(X/2)/(1-TAN(X/2)^2)
// USE THE IDENTITY ATAN(X) = 1/2 ATAN(2X/(1-X^2))
// APPLY TWICE AND WORK WITH X/4 OR 4/X
if(gtReal(&RReg[0],&RReg[1])) {
addReal(&RReg[2],&RReg[0],&RReg[0]);
addReal(&RReg[0],&RReg[2],&RReg[2]); // 1/4X
invert=1;
divReal(&RReg[0],&RReg[1],&RReg[0]); // INVERT THE ARGUMENT
} else {
invert=0;
RReg[1].data[0]=25;
RReg[1].exp=-2;
mulReal(&RReg[2],&RReg[0],&RReg[1]);
swapReal(&RReg[2],&RReg[0]);
}
// HERE 0<=X<=0.25
// HERE 0<=X<=1
// USE SERIES
// y=2x/(1-x^2)
RReg[2].data[0]=57;
RReg[2].exp=-2;
RReg[2].len=1;
RReg[2].flags=0;
if(gtReal(&RReg[0],&RReg[2])) {
// When x>0.57 use: atan(x)=pi/4-atan(1/y)/2 with only 1 pass
pass=1;
mulReal(&RReg[3],&RReg[0],&RReg[0]); // X^2
subReal(&RReg[2],&RReg[1],&RReg[3]); // 1-X^2
addReal(&RReg[3],&RReg[0],&RReg[0]); // 2X
divReal(&RReg[0],&RReg[2],&RReg[3]); // (1-X^2)/2X = 1/Y
} else {
RReg[2].data[0]=42;
if(gtReal(&RReg[0],&RReg[2])) {
// When x>0.42 and x<0.57 use: atan(x)=pi/8+atan(-1/y)/4 with 2 passes
pass=2;
mulReal(&RReg[3],&RReg[0],&RReg[0]); // X^2
subReal(&RReg[2],&RReg[1],&RReg[3]); // 1-X^2
addReal(&RReg[3],&RReg[0],&RReg[0]); // 2X
divReal(&RReg[0],&RReg[3],&RReg[2]); // 2X/(1-x^2) = Y
mulReal(&RReg[3],&RReg[0],&RReg[0]); // X^2
subReal(&RReg[2],&RReg[1],&RReg[3]); // 1-X^2
addReal(&RReg[3],&RReg[0],&RReg[0]); // 2X
divReal(&RReg[0],&RReg[2],&RReg[3]); // (1-X^2)/2X = 1/Y
RReg[0].flags^=F_NEGATIVE; // -1/Y
} else {
// LESS THAN 0.42, JUST USE THE NUMBER AS-IS
pass=0;
}
}
RReg[2].flags=F_NEGATIVE;
RReg[2].exp=0;
@ -1336,7 +1372,7 @@ void atanpower()
copyReal(&RReg[4],&RReg[0]); // ACCUMULATOR
// DO AS MANY TERMS AS NEEDED
for(k=3;k<=REAL_PRECISION_MAX/2;k+=2)
for(k=3;1;k+=2)
{
mulReal(&RReg[1],&RReg[0],&RReg[3]); // TERM*X^2
@ -1355,28 +1391,38 @@ void atanpower()
printf("atan iters=%d\n",(k-3)/2);
// NOW APPLY THE IDENTITY BACKWARDS
// ATAN(X) = 2 ATAN(X/2)/(1-ATAN(X/2)^2)
MACROOneToRReg(3);
mulReal(&RReg[1],&RReg[4],&RReg[4]); // ATAN(X/4)^2
subReal(&RReg[2],&RReg[3],&RReg[1]); // 1-ATAN(X/4)^2
divReal(&RReg[1],&RReg[4],&RReg[2]); // ATAN(X/4)/(1-ATAN(X/4)^2
addReal(&RReg[4],&RReg[1],&RReg[1]); // ATAN(X/2)=2*...
REAL pi_2;
decconst_PI_2(&pi_2);
// AND ONE MORE TIME...
mulReal(&RReg[1],&RReg[4],&RReg[4]); // ATAN(X/2)^2
subReal(&RReg[2],&RReg[3],&RReg[1]); // 1-ATAN(X/2)^2
divReal(&RReg[1],&RReg[4],&RReg[2]); // ATAN(X/2)/(1-ATAN(X/2)^2
addReal(&RReg[4],&RReg[1],&RReg[1]); // ATAN(X)=2*...
switch(pass)
{
default:
case 0:
break;
case 1:
// When x>0.57 use: atan(x)=pi/4-atan(1/y)/2 with only 1 pass
subReal(&RReg[2],&pi_2,&RReg[4]);
RReg[1].data[0]=5;
RReg[1].exp=-1;
RReg[1].flags=0;
RReg[1].len=1;
mulReal(&RReg[4],&RReg[2],&RReg[1]);
break;
case 2:
// When x>0.42 and x<0.57 use: atan(x)=pi/8+atan(-1/y)/4 with 2 passes
addReal(&RReg[2],&pi_2,&RReg[4]);
RReg[1].data[0]=25;
RReg[1].exp=-2;
RReg[1].flags=0;
RReg[1].len=1;
mulReal(&RReg[4],&RReg[2],&RReg[1]);
break;
}
if(invert) {
REAL pi_2;
decconst_PI_2(&pi_2);
subReal(&RReg[0],&pi_2,&RReg[4]);
}
else swapReal(&RReg[4],&RReg[0]);
@ -1456,7 +1502,7 @@ int main()
#define TEST_DIGITS 32
#define TEST_DIGITS 2000
Context.precdigits=TEST_DIGITS;