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,34 +59,44 @@ static double atanlist[100] = {
     3.1554436208840472216469142611311e-30, 1.5777218104420236108234571305656e-30
 };
 
-int trig(double theta, double *xy, double *dxy){
+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;
     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){
+    for (i = 0; i < 56; i++)
+    {
         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];
-	i++;
         t /= 2.0;
     }
     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 dxy[2];
@@ -96,16 +106,60 @@ double sin1(double theta, int *it, double *dsin){
     return xy[1];
 }
 
-double sin2(double theta, int *it, double *dsin){
+// just do sin, instead of trig()
-    if(theta < 0) return -sin2(theta, it, dsin);
+double sin11(double theta, int *it, double *dsin)
-    if(theta > 2*PI) return sin2(theta - 2*PI, it, dsin);
+{
-    if(theta > PI) return -sin2(theta - PI, it, dsin);
+    if (theta < 0)
-    if(theta > PI/2) return sin2(PI - theta, it, dsin);
+        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;
+    for (i = 0; i < 55; i++)
+    {
+        // 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];
+        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 +172,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 < 10000000; 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 %f 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 < 10000000; 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 %f s.\n", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
+    }
+
     return 0;
 }