wyj/sin
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update: format

Browse Source
This commit is contained in:
Yingjie Wang 2024-06-29 15:35:48 -04:00
parent c8874275c3
commit ad173aa8ae
1 changed files with 104 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
sin.c
View File

@ -59,7 +59,8 @@ static double atanlist[100] = {
3.1554436208840472216469142611311e-30, 1.5777218104420236108234571305656e-30 3.1554436208840472216469142611311e-30, 1.5777218104420236108234571305656e-30
}; };
int trig(double theta, double *xy, double *dxy){ int trig(double theta, double *xy, double *dxy)
{
int i = 0; int i = 0;
double t = 1.0; // t=2^-i double t = 1.0; // t=2^-i
int d = 0; int d = 0;
@ -68,7 +69,8 @@ int trig(double theta, double *xy, double *dxy){
xy[0] = x0; xy[0] = x0;
xy[1] = y0; xy[1] = y0;
while(i<100 && dnorm>1e-32){ while (i < 100 && dnorm > 1e-32)
{
d = alpha > 0 ? 1 : (alpha < 0 ? -1 : 0); d = alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
dxy[0] = -d * xy[1] * t; dxy[0] = -d * xy[1] * t;
dxy[1] = d * xy[0] * t; dxy[1] = d * xy[0] * t;
@ -82,11 +84,16 @@ int trig(double theta, double *xy, double *dxy){
return i; return i;
} }
double sin1(double theta, int *it, double *dsin){ double sin1(double theta, int *it, double *dsin)
if(theta < 0) return -sin1(theta, it, dsin); {
if(theta > 2*PI) return sin1(theta - 2*PI, it, dsin); if (theta < 0)
if(theta > PI) return -sin1(theta - PI, it, dsin); return -sin1(theta, it, dsin);
if(theta > PI/2) return sin1(PI - theta, it, dsin); if (theta > 2 * PI)
return sin1(theta - 2 * PI, it, dsin);
if (theta > PI)
return -sin1(theta - PI, it, dsin);
if (theta > PI / 2)
return sin1(PI - theta, it, dsin);
double xy[2]; double xy[2];
double dxy[2]; double dxy[2];
@ -96,16 +103,22 @@ double sin1(double theta, int *it, double *dsin){
return xy[1]; return xy[1];
} }
double sin2(double theta, int *it, double *dsin){ double sin2(double theta, int *it, double *dsin)
if(theta < 0) return -sin2(theta, it, dsin); {
if(theta > 2*PI) return sin2(theta - 2*PI, it, dsin); if (theta < 0)
if(theta > PI) return -sin2(theta - PI, it, dsin); return -sin2(theta, it, dsin);
if(theta > PI/2) return sin2(PI - theta, it, dsin); if (theta > 2 * PI)
return sin2(theta - 2 * PI, it, dsin);
if (theta > PI)
return -sin2(theta - PI, it, dsin);
if (theta > PI / 2)
return sin2(PI - theta, it, dsin);
double result = 0; double result = 0;
double term = theta; double term = theta;
int i = 0; int i = 0;
while (i<100 && fabs(term)>1e-16){ while (i < 100 && fabs(term) > 1e-16)
{
result += term; result += term;
i += 1; i += 1;
term *= -1; term *= -1;
@ -118,20 +131,24 @@ double sin2(double theta, int *it, double *dsin){
return result; return result;
} }
int main(){ int main()
{
int it; int it;
double theta = 1.0; double theta = 1.0;
double sin, dsin; double sin, dsin;
clock_t start, end; clock_t start, end;
start = clock(); start = clock();
for (int i=0; i<10000; i++) sin = sin1(theta, &it, &dsin); for (int i = 0; i < 10000; i++)
sin = sin1(theta, &it, &dsin);
end = clock(); end = clock();
printf("sin1(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC); printf("sin1(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
start = clock(); start = clock();
for (int i=0; i<10000; i++) sin = sin2(theta, &it, &dsin); for (int i = 0; i < 10000; i++)
sin = sin2(theta, &it, &dsin);
end = clock(); end = clock();
printf("sin2(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC); printf("sin2(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
return 0; return 0;
} }