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

update

Browse Source
This commit is contained in:
Yingjie Wang 2024-06-29 15:49:52 -04:00
parent ad173aa8ae
commit d03ad61916
1 changed files with 54 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
sin.c
View File

@ -65,16 +65,14 @@ int trig(double theta, double *xy, double *dxy)
double t = 1.0; // t=2^-i double t = 1.0; // t=2^-i
int d = 0; int d = 0;
double alpha = theta; double alpha = theta;
double dnorm = 1.0;
xy[0] = x0; xy[0] = x0;
xy[1] = y0; xy[1] = y0;
while (i < 100 && dnorm > 1e-32) while (i < 100 && fabs(alpha) > 1e-16)
{ {
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;
dnorm = dxy[0] * dxy[0] + dxy[1] * dxy[1];
xy[0] += dxy[0]; xy[0] += dxy[0];
xy[1] += dxy[1]; xy[1] += dxy[1];
alpha -= d * atanlist[i]; alpha -= d * atanlist[i];
@ -103,6 +101,45 @@ double sin1(double theta, int *it, double *dsin)
return xy[1]; return xy[1];
} }
// just do sin, instead of trig()
double sin11(double theta, int *it, double *dsin)
{
if (theta < 0)
return -sin11(theta, it, dsin);
if (theta > 2 * PI)
return sin11(theta - 2 * PI, it, dsin);
if (theta > PI)
return -sin11(theta - PI, it, dsin);
if (theta > PI / 2)
return sin11(PI - theta, it, dsin);
double xy[2];
double dxy[2];
int i = 0;
double t = 1.0; // t=2^-i
int d = 0;
double alpha = theta;
xy[0] = x0;
xy[1] = y0;
dxy[1] = 1;
while (i < 100 && fabs(dxy[1]) > 1e-16)
{
d = alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
dxy[0] = -d * xy[1] * t;
dxy[1] = d * xy[0] * t;
xy[0] += dxy[0];
xy[1] += dxy[1];
alpha -= d * atanlist[i];
i++;
t /= 2.0;
}
*dsin = dxy[1];
*it = i;
return xy[1];
}
double sin2(double theta, int *it, double *dsin) double sin2(double theta, int *it, double *dsin)
{ {
if (theta < 0) if (theta < 0)
@ -134,21 +171,24 @@ double sin2(double theta, int *it, double *dsin)
int main() int main()
{ {
int it; int it;
double theta = 1.0; double thetalist[5] = {1e-4, 0.5, 1, 1.3, 1.7};
double sin, dsin; double sin, dsin;
clock_t start, end; clock_t start, end;
start = clock(); for (int i = 0; i < 5; i++)
for (int i = 0; i < 10000; i++) {
sin = sin1(theta, &it, &dsin); start = clock();
end = clock(); for (int j = 0; j < 10000; j++)
printf("sin1(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC); sin = sin11(thetalist[i], &it, &dsin);
end = clock();
printf("sin11(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
start = clock(); start = clock();
for (int i = 0; i < 10000; i++) for (int j = 0; j < 10000; j++)
sin = sin2(theta, &it, &dsin); sin = sin2(thetalist[i], &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", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
}
return 0; return 0;
} }