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