圓弧點的數學運算

一個半徑為 r 圓的方程式,圓線上的點(x, y)滿足方程式 x² + y² = r² 來描述
或是用 f(x, y) = x² + y² - r² = 0 上所有 P(x, y) 的點畫出來就會形成一個圓, 他的座標就是:
           x = r cos(θ)              0 <  θ < 2Π
           y = r sin(θ)               0 < θ < 2Π
已知一個點 (x, y) 等於是指定一個角度 θ 所在的圓點上,換句話說, θ1 座標是:
          x1 =  r cos(θ1) - > cos(θ1) = x1 / r 
          y1 =  r sin(θ1) - > sin(θ1) = y1 / r
同樣的另外一個點 (x2, y2)角度是θ2,等於知道該點的座標在:
          x2 = r cos(θ2) - > cos(θ2) = x2 / r
          y2 = r sin(θ2) - > sin(θ2) = y2 / r
以此類推當 θ3 = θ1 + θ2 位置在那呢? 當然就在 θ3 所指的圓線上的一點:
           x3 = r cos(θ3))
           y3 = r sin(θ3)
再將 θ3 替換成 θ1 + θ2 就可得到:
          x3 = r cos(θ3) = r cos(θ1 + θ2) 
          y3 = r sin(θ3)  = r sin(θ1 + θ2)
根據三角公式將 sin(θ), cos(θ) 換成(x1, y1)與(x2, y2)的運算關係就是:
          x3 = r cos(θ1 + θ2) = r [cos(θ1) cos(θ2) - sin(θ2) cos(θ1)] 
              = r [(x1 / r)(x2 / r) - (y1 / r) (y2 / r)] = (x1 x2 - y1 y2) / r   --- (1)
          y3 = r sin(θ1 + θ2) = r [sin(θ1) cos(θ2) + sin(θ2) cos(θ1)]
                 r [(y1 / r)(x2 / r) + (y2 / r)(x1 / r)] = (x1 y2 + y1 x2) / r   --- (2)
這就是圓座標點加法運算式(Adder),兩圓點相加去算出另一個圓點的座標,用數學記號表示就是:
          P(x1, y1) + Q(x2, y2) = R(x3, y3)     簡寫成 而 P +  Q =  R
因此2倍角座標就很容易推導出來, 因為 (x, y) = (x1, y1) = (x2, y2) 帶入(1)(2)式
         x" = r cos(2θ) = (x² - y²) / r   --- (3)
         y" = r sin(2θ) = 2xy / r             --- (4)
這就是圓座標點加倍運算式(Doubler),簡寫成 2P, 再利用式(1)(2)座標點相加運算式, 將 θ 座標點與兩倍角 2θ 座標點相加就可推導出 3 倍角座標:
         x'" = r cos(3θ) = (x(x² - y²) / r - y 2xy / r) / r = (x³ - 3xy²) / r²
         y'" = r sin(3θ) = (x 2xy / r + y(x² - y²) / r) / r = (3x²y- y³) / r²
以此類推, 用 Double and And 疊帶累加法就可算出 n 倍角的座標,弳度的乘法運算 θ + θ        ... 共 n 個 = n * θ 也是我們熟悉的用法,只要知道一個角度的座標點, 即使不用三角函數也可算出其他角度的座標,再來看弳度 0 的座標,用記號 O 來表示,也就是說 O 的座標位在:
         x0 = r cos(0) =r
         y0 = r sin(0) =0
帶入圓點加法公式去看 P(x, y) + O 是否等於 P(x, y)呢?
           (x*r  - y*0) / r = x
           (x*0 + y*r) / r = y 
確實如此,因此座標(r, 0)的確是中性點 O, 一條線的座標點 f(x, y) = 0, 如果將(x, y)座標同時放大 z 倍時,由原先的座標來看這新的點所形成的圖形會長成什麼樣子呢? 既然新增了一個變數 z, 就把他想成是 z 軸,將(x/z, y/z)映射到(X, Y, Z)來表示,得到新的投射座標,換句話說就是分別用(X/Z,Y/Z)去取代方程式中的(x,y),比方一個圓方程式:

            x² + y² = r²      
將 x = X / Z, y = Y / Z 帶入得到:
            (X / Z)² + (X / Z)² = r²
            X² + Y² = r² Z² = (r Z)² = R²      , R = r Z
可以看到它仍舊是圓的方程式,只不過半徑跟著放大成Z倍.

用 C 語言寫個小程式驗證一下:
// cir.cpp use g++ to compile
// g++  cir.cpp
#include  <iostream>
#include  <math.h>
using namespace std;
class  PointC {
    public:
        double x, y;
        PointC( ) { };
        PointC(double x0, double y0) { x=x0; y=y0; }
};
class  curveC {
        PointC p; // private point Generator
    public:             
        curveC(const PointC &P)         { p = P; } //need to implement PointC() constructor
        curveC(const PointC &P, int n)  { p = P; if ( n > 1 ) scale(n); }
        PointC operator + (const PointC &Another)   { return incPointC(Another);  }
        PointC operator + (const curveC &Another)   { return incPointC(Another.p);}
        PointC operator * (int n)                   { return curveC(p, n).p;      }
        PointC incPointC(const PointC &Another)     { return curveC(p).adder(Another)->p;}
        curveC* adder(const PointC &q) { // self add, call by reference,
            double r = sqrt(p.x * p.x + p.y * p.y);
            double px = p.x;    // copy to prevent update by p.x
            p.x = (px * q.x - p.y * q.y) / r;  // update x
            p.y = (px * q.y + p.y * q.x) / r;   // update y  
            // return curveC(p); // if wrap with class, can use . not -> to obtain member
            return this;
        }
        curveC* doubler() { // self double
            double r = sqrt(p.x * p.x + p.y * p.y);
            double px = p.x;    // copy to prevent update by p.x
            p.x = (p.x * p.x - p.y * p.y) / r ; // update x
            p.y = 2 * px * p.y / r;             // update y
            return this;
        }         
        #define  testBit(n, b)   n & (1L << b)
        curveC* scale(int n)   {   // double-and-add   DAA algorithm
            int k = 32;             // 32 bit integer
            PointC G = p;           // copy to prevent update
            while ( k-- )   if (testBit(n,k))   break;
            while ( k-- ) {         // Double And Add
                doubler();          // self double
                if (testBit(n,k))   adder(G);// self adder with G     
            }
            return this;
        }
        #undef testBit
        curveC* showPoint(void) {
            printf("\tx=%f , y=%f\n", p.x, p.y);
            return this;
        }
};
main() {  
    double degree = 10; 
    double theta  = 2 * 3.14159 * degree / 360;
    double r = 5;      // Circle half diameter = 5
    PointC G = PointC(r * cos(theta), r * sin(theta) );
    curveC c = curveC(G);
    for (int i = 1; i <= 10; i ++) {
        curveC(c * i).showPoint();
        printf("equal%2d\tx=%f , y=%f\n\n",i ,r*cos(theta*i), r*sin(theta*i));     
    }
    c.showPoint();
}