update: update readme

Readme.md

@@ -0,0 +1,21 @@
+对比 https://www.bilibili.com/video/BV1Ys3TeFEgi/ 中提到的迭代计算 $\sin$ 函数的算法和使用泰勒展开，意外地发现还是泰勒展开赢麻了……
+
+## 当前结果
+在 i5 10400ES 上，通过
+```bash
+gcc -march=native -O2 -o sin sin.c
+```
+得到结果为（用时是指循环 10000000 次总耗时，否则都太快了计时误差太大）
+```bash
+./sin
+sin11(0.0001) = 9.999999983330736e-05, it = 55, dsin=-5.551115095370206e-17, took 1.013875 s.
+sin2(0.0001) = 9.999999983333334e-05, it = 2, dsin=8.333333333333334e-23, took 0.064005 s.
+sin11(0.5) = 0.4794255386042029, it = 55, dsin=-4.871561831101115e-17, took 0.977128 s.
+sin2(0.5) = 0.479425538604203, it = 7, dsin=-2.333729166204778e-17, took 0.284716 s.
+sin11(1) = 0.8414709848078967, it = 55, dsin=2.999280301164362e-17, took 0.981563 s.
+sin2(1) = 0.8414709848078965, it = 9, dsin=-8.220635246624331e-18, took 0.421509 s.
+sin11(1.3) = 0.9635581854171931, it = 55, dsin=-1.484916792996379e-17, took 0.981924 s.
+sin2(1.3) = 0.963558185417193, it = 10, dsin=4.835779466409166e-18, took 0.471155 s.
+sin11(1.7) = 0.9916648104524686, it = 55, dsin=-7.152306208153797e-18, took 0.991410 s.
+sin2(1.7) = 0.9916648104524686, it = 10, dsin=4.239843809646965e-17, took 0.472351 s.
+```

sin.c

@@ -59,6 +59,12 @@ static double atanlist[100] = {
     3.1554436208840472216469142611311e-30, 1.5777218104420236108234571305656e-30
 };
 
+int sign(double alpha)
+{
+    return alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
+    // return (int)copysign(1.0, alpha);
+}
+
 int trig(double theta, double *xy, double *dxy)
 {
     int i = 0;
@@ -68,7 +74,7 @@ int trig(double theta, double *xy, double *dxy)
 
     xy[0] = x0;
     xy[1] = y0;
-    while (i < 100 && fabs(alpha) > 1e-16)
+    for (i = 0; i < 56; i++)
     {
         d = alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
         dxy[0] = -d * xy[1] * t;
@@ -76,7 +82,6 @@ int trig(double theta, double *xy, double *dxy)
         xy[0] += dxy[0];
         xy[1] += dxy[1];
         alpha -= d * atanlist[i];
-        i++;
         t /= 2.0;
     }
     return i;
@@ -124,15 +129,14 @@ double sin11(double theta, int *it, double *dsin)
     xy[0] = x0;
     xy[1] = y0;
     dxy[1] = 1;
-    while (i < 100 && fabs(dxy[1]) > 1e-16)
+    for (i = 0; i < 55; i++)
     {
-        d = alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
-        dxy[0] = -d * xy[1] * t;
-        dxy[1] = d * xy[0] * t;
+        // d = sign(alpha);
+        dxy[0] = -xy[1] * copysign(t,alpha);
+        dxy[1] = xy[0] * copysign(t,alpha);
         xy[0] += dxy[0];
         xy[1] += dxy[1];
         alpha -= d * atanlist[i];
-        i++;
         t /= 2.0;
     }
     *dsin = dxy[1];
@@ -178,16 +182,16 @@ int main()
     for (int i = 0; i < 5; i++)
     {
         start = clock();
-        for (int j = 0; j < 10000; j++)
+        for (int j = 0; j < 10000000; 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);
+        printf("sin11(%.16g) = %.16g, it = %d, dsin=%.16g, took %f s.\n", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
 
         start = clock();
-        for (int j = 0; j < 10000; j++)
+        for (int j = 0; j < 10000000; 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);
+        printf("sin2(%.16g) = %.16g, it = %d, dsin=%.16g, took %f s.\n", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
     }
 
     return 0;