update

sin.c

@@ -59,20 +59,20 @@ static double atanlist[100] = {
     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;
     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];
@@ -82,11 +82,16 @@ int trig(double theta, double *xy, double *dxy){
     return i;
 }
 
-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 > PI) return -sin1(theta - PI, it, dsin);
-    if(theta > PI/2) return sin1(PI - theta, it, 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 > PI)
+        return -sin1(theta - PI, it, dsin);
+    if (theta > PI / 2)
+        return sin1(PI - theta, it, dsin);
 
     double xy[2];
     double dxy[2];
@@ -96,16 +101,61 @@ double sin1(double theta, int *it, double *dsin){
     return xy[1];
 }
 
-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 > PI) return -sin2(theta - PI, it, dsin);
-    if(theta > PI/2) return sin2(PI - theta, it, dsin);
+// 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)
+        return -sin2(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 term = theta;
     int i = 0;
-    while (i<100 && fabs(term)>1e-16){
+    while (i < 100 && fabs(term) > 1e-16)
+    {
         result += term;
         i += 1;
         term *= -1;
@@ -118,20 +168,27 @@ double sin2(double theta, int *it, double *dsin){
     return result;
 }
 
-int main(){
+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;
 
+    for (int i = 0; i < 5; i++)
+    {
         start = clock();
-    for (int i=0; i<10000; i++) sin = sin1(theta, &it, &dsin);
+        for (int j = 0; j < 10000; j++)
+            sin = sin11(thetalist[i], &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);
+        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);
+        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", 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;
 }