建立數學向量 (Vector)

    數學向量的抽象資料結構

    此處的向量 (vector) 是指在數學上的向量運算。現在已經有許多程式語言,像是 MATLAB 或 R,支援這類運算。實務上,當我們要用 C 或 Fortran 進行向量或矩陣運算時,我們會使用 BLAS 系列函式庫,因為使用 BLAS 會比手刻來得有效率。此處的重點是了解其原理,而非重造輪子來用。

    向量的抽象資料結構如下:

    T, U, V are vectors.
    s is a scalar.
    
    sub Size(V): sz, sz is an integer and sz >= 0
    sub At(V, index): data
    sub SetAt(V, index, value): void
    sub IsEqual(U, V): bool
    sub Add(U, V): T
    sub Add(s, V): T
    sub Sub(U, V): T
    sub Sub(s, V): T
    sub Sub(V, s): T
    sub Mul(U, V): T
    sub Mul(s, V): T
    sub Div(U, V): T
    sub Div(s, V): T
    sub Div(V, s): T
    sub Dot(U, V): S
    sub Cross(U, V): T
    

    以下運算和向量的狀態有關:

    • Size(V):得知向量的長度 (或大小)
    • At(V, index):取得向量中任意位置的值
    • SetAt(V, index, value):設置向量中位意任置的值
    • IsEqual(U, V):檢查兩向量是否相等

    根據向量和純量的相對位置,有以下三種運算:

    • Add(U, V):向量和向量相加
    • Add(s, V):純量和向量相加,純量在前
    • Add(V, s):向量和純量相加,純量在後

    加法和乘法有交換律,因此情境二和情境三等效;但在減法及除法兩者不相等。

    以下是其他的向量運算:

    • Dot(U, V):內積 (product)
    • Cross(U, V):外積 (cross product)

    宣告數學向量類別

    以下是向量類別的宣告:

    class Vector:
        size <- n, n belongs to integer and n > 0
        elements <- a zero-based array
    

    在我們這個範例程式中,向量不會擴展容量,故只用單一 size 來記錄向量的大小。由於向量會頻繁用到隨機存取,故內部以陣列而非串列來實作。以下是等效的 C 語言程式:

    typedef struct vector vector_t;
    
    struct vector {
        size_t size;
        double *elements;
    };
    

    數學向量的建構函式

    以下是 C 語言的數學向量的建構函式:

    vector_t * vector_new(size_t size)               /*  1 */
    {                                                /*  2 */
        assert(size > 0);                            /*  3 */
    
        vector_t *v = malloc(sizeof(vector_t));      /*  4 */
        if (!v) {                                    /*  5 */
            return v;                                /*  6 */
        }                                            /*  7 */
        
        v->size = size;                              /*  8 */
        v->elements = calloc(size, sizeof(double));  /*  9 */
        if (!(v->elements)) {                        /* 10 */
            vector_delete(v);                        /* 11 */
            v = NULL;                                /* 12 */
            return v;                                /* 13 */
        }                                            /* 14 */
    
        return v;                                    /* 15 */
    }                                                /* 16 */
    

    我們在第 4 行建立 vector_t 物件 v,並在第 5 行至第 7 行確認物件 v 不為空。

    接著,我們在第 8 行存入物件 v 的大小。

    我們在第 9 行建立並初始化 v 的內部陣列 elements。緊接著在第 10 行至第 14 行檢查 elements 不為空,若 elements 為空,則釋放掉物件 v

    最後在第 15 行回傳物件 v

    上述建構函式會建立一個大小為 size,初始值皆為 0.0 的向量,這算是比較傳統的建構函式。

    我們也可以在建立向量時一併初始化向量,其使用方式如下:

    vector_t *v = vector_init(4, 2.0, 3.0, 4.0, 5.0);
    

    在此處的 vector_init 建構函式中,第一個參數是向量的大小,第二個以後的參數是實際要填入向量的值。其建構函式如下:

    vector_t * vector_init(size_t size, double value, ...)  /*  1 */
    {                                                       /*  2 */
        assert(size > 0);                                   /*  3 */
    
        vector_t *v = malloc(sizeof(vector_t));             /*  4 */
        if (!v) {                                           /*  5 */
            return v;                                       /*  6 */
        }                                                   /*  7 */
    
        v->size = size;                                     /*  8 */
        v->elements = calloc(size, sizeof(double));         /*  9 */
        if (!(v->elements)) {                               /* 10 */
            vector_delete(v);                               /* 11 */
            v = NULL;                                       /* 12 */
            return v;                                       /* 13 */
        }                                                   /* 14 */
    
        v->elements[0] = value;                             /* 15 */
        va_list args;                                       /* 16 */
        va_start(args, value);                              /* 17 */
    
        for (size_t i = 1; i < size; i++) {                 /* 18 */
            v->elements[i] = va_arg(args, double);          /* 19 */
        }                                                   /* 20 */
    
        va_end(args);                                       /* 21 */
    
        return v;                                           /* 22 */
    }                                                       /* 23 */
    

    第 4 行至第 14 行為初始化 vector_t 物件 v 的動件,和前一個建構函式雷同,不重覆說明。

    從第 15 行到第 17 行,我們初始化傳入的不定參數 arg,並填入第一個值 value

    從第 18 行到第 20 行,我們依序填入其他參數。

    在第 21 行,我們釋放不定參數 arg

    最後,在第 22 行,回傳填好值的物件 v

    數學向量的解構函式

    以下是 C 語言的解構函式:

    void vector_delete(void *self)                         /*  1 */
    {                                                      /*  2 */
        if (!self) {                                       /*  3 */
            return;                                        /*  4 */
        }                                                  /*  5 */
    
        double *elements = ((vector_t *) self)->elements;  /*  6 */
        if (elements) {                                    /*  7 */
            free(elements);                                /*  8 */
        }                                                  /*  9 */
        
        free(self);                                        /* 10 */
    }                                                      /* 11 */
    

    釋放記憶體時,要由內向外逐一釋放。我們先在第 6 行至第 9 行間釋放內部的陣列 elements,最後在第 10 行釋放物件 self

    向量的大小

    以下是取得向量大小的虛擬碼:

    V is a vector.
    
    sub Size(V): sz
        return V.size
    

    由於我們在建構向量時即存有其大小,直接取值即可。以下是等效的 C 語言程式碼:

    size_t vector_size(const vector_t *self)
    {
        assert(self);
    
        return self->size;
    }
    

    隨機存取向量中任意位置的值

    以下是存取向量中任意位置的值的虛擬碼:

    V is a vector.
    
    sub At(V, index): data
        assert 0 <= index < V.size
        
        return V.elements[index]
    

    確認索引值合法後,以索引值從陣列取值即可。以下是等效的 C 語言實作:

    double vector_at(const vector_t *self, size_t index)
    {
        assert(index < vector_size(self));
    
        return self->elements[index];
    }
    

    由於 size_t 一定大於等於 0,故只需檢查索引值上限。

    設置向量中任意位置的值

    以下是設置向量中任意位置的值的虛擬碼:

    V is a vector.
    
    sub SetAt(V, index, value): void
        assert 0 =< index < V.size
    
        V.elements[index] <- value
    

    確認索引值合法後,藉由索引值取得陣列中相對應的位置,再修改其值即可。以下是等效的 C 語言程式碼:

    void vector_set_at(vector_t *self, size_t index, double value)
    {
        assert(index < vector_size(self));
    
        self->elements[index] = value;
    }
    

    確認兩向量相等

    以下是確認兩向量相等的虛擬碼:

    U, V are vectors.
    
    sub IsEqual(U, V): bool
        if Size(U) != Size(V) then
            return false
        end if
    
        for (i from 0 to Size(U) - 1) do
            if At(U, i) != At(V, i) then
                return false
            end if
        end for
        
        return true
    

    先確認兩向量長度相等,再確認兩向量中每個值皆相等即可得知兩向量是相等的。以下是等效的 C 語言實作:

    bool vector_is_equal(vector_t* u, vector_t* v)     /*  1 */
    {                                                  /*  2 */
        if (!(vector_size(u) == vector_size(v))) {     /*  3 */
            return false;                              /*  4 */
        }                                              /*  5 */
    
        for (size_t i = 0; i < vector_size(u); i++) {  /*  6 */
            double a = vector_at(u, i);                /*  7 */
            double b = vector_at(v, i);                /*  8 */
    
            if (fabs(a - b) > 0.000001) {              /*  9 */
                return false;                          /* 10 */
            }                                          /* 11 */
        }                                              /* 12 */
    
        return true;                                   /* 13 */
    }                                                  /* 14 */
    

    由於本範例程式的向量儲存的是浮點數,在兩數間有一些微小誤差,因此我們在第 9 行至第 11 行檢查兩元素間的誤差值。在本實作的認定中,只要誤差小於 10^-6 即代表兩者相等。讀者可視自身需求修改兩者相等的定義。

    兩向量相加

    以下是兩數學向量相加的虛擬碼:

    U, V, T are vectors.
    
    sub Add(U, V): T
        assert Size(U) == Size(V)
        
        T <- New(Size(U))
        for i (0 to Size(U) - 1) do
            SetAt(T, i, At(U, i) + At(V, i))
        end for
        
        return T
    

    要先確認兩向量的長度相等,之後再進行下一步動作。按照其數學定義,將兩向量的元素逐一相加即可。

    我們在用 C 語言實作時,利用 C11 的巨集泛型敘述宣告具有多型的向量加法,其程式碼如下:

    #if __STDC_VERSION__ >= 201112L
    #define vector_add(u, v) \
        _Generic((u), \
            vector_t *: _Generic((v), \
                vector_t *: vector_add_vv, \
                default: vector_add_vs), \
            default: vector_add_sv)((u), (v))
    #else
    vector_t * vector_add(const vector_t *u, const vector_t *v);
    #endif
    

    當編譯器支援 C11 時,宣告具有多型的 vector_add;反之,則將 vector_add 改為一般的向量加法,做為退路 (fallback)。以下是相對應的 C 語言程式碼:

    #if __STDC_VERSION__ < 201112L
    vector_t * vector_add(const vector_t *u, const vector_t *v)
    {
        return vector_add_vv(u, v);
    }
    #endif
    
    vector_t * vector_add_vv(const vector_t *u, const vector_t *v)
    {
        assert(vector_size(u) == vector_size(v));
    
        vector_t *out = vector_new(vector_size(u));
        if (!out) {
            return out;
        }
    
        for (size_t i = 0; i < vector_size(u); i++) {
            vector_set_at(out, i, vector_at(u, i) + vector_at(v, i));
        }
    
        return out;
    }
    

    向量和純量相加

    以下是向量和純量相加的虛擬碼:

    V, T are vectors.
    S is a scalar.
    
    sub Add(V, S): T
        T <- New(sz)
    
        for (i from 0 to Size(V) - 1) do
            SetAt(T, i, At(V, i) + S)
        end for
        
        return T
    

    根據其數學上的定義來實作即可。以下是等效的 C 語言程式碼:

    vector_t * vector_add_vs(const vector_t *v, double s)
    {
        vector_t *out = vector_new(vector_size(v));
        if (!out) {
            return out;
        }
    
        for (size_t i = 0; i < vector_size(v); i++) {
            vector_set_at(out, i, vector_at(v, i) + s);
        }
    
        return out;
    }
    
    vector_t * vector_add_sv(double s, const vector_t *v)
    {
        return vector_add_vs(v, s);
    }
    

    向量內積 (dot product)

    以下是向量內積的虛擬碼:

    U, V are vectors
    
    sub Dot(U, V): result
        assert Size(U) == Size(V)
        
        result <- 0
    
        for i (0 to U.size - 1) do
            result += U.elements[i] * V.elements[i]
        end for
        
        return result
    

    確認兩向量的長度相等後,根據其數學定義實作即可。以下是等效的 C 語言程式碼:

    double vector_dot(const vector_t *u, const vector_t *v)
    {
        assert(vector_size(u) == vector_size(v));
    
        double result = 0.0;
    
        for (size_t i = 0; i < vector_size(u); i++) {
            result += vector_at(u, i) * vector_at(v, i);
        }
    
        return result;
    }
    

    向量外積 (cross product)

    註:本節的實作僅限於三維向量。

    以下是其虛擬碼:

    U, V, T are vectors
    
    sub Cross(U, V): T
        assert Size(U) == 3
        assert Size(V) == 3
        
        T <- a new vector which size == 3
        
        T.elements[0] <- U.elements[1] * V.elements[2] - U.elements[2] * V.elements[1]
        T.elements[1] <- U.elements[2] * V.elements[0] - U.elements[0] * V.elements[2]
        T.elements[2] <- U.elements[0] * V.elements[1] - U.elements[1] * V.elements[0]
    
        return T
    

    我們使用狹義的向量外積,只能在三維向量下進行運算。根據其數學定義實作即可。以下是等效的 C 語言實作:

    vector_t * vector_cross(const vector_t *u, const vector_t *v)
    {
        assert(vector_size(u) == 3);
        assert(vector_size(v) == 3);
    
        vector_t *out = vector_new(vector_size(u));
        if (!out) {
            return out;
        }
    
        vector_set_at(out, 0, 
            vector_at(u, 1) * vector_at(v, 2) - vector_at(u, 2) * vector_at(v, 1));
        vector_set_at(out, 1,
            vector_at(u, 2) * vector_at(v, 0) - vector_at(u, 0) * vector_at(v, 2));
        vector_set_at(out, 2,
            vector_at(u, 0) * vector_at(v, 1) - vector_at(u, 1) * vector_at(v, 0));
    
        return out;
    }
    

    在演算法上的效率

    依照本範例程式的實作,在演算法上的效率如下:

    • Size(V)O(1)
    • At(V, index)O(1)
    • SetAt(V, index, value)O(1)
    • IsEqual(U, V)O(n)
    • Add(U, V)O(n)
    • Add(V, S)O(n)
    • Dot(U, V)O(n)
    • Cross(U, V)O(n)
    【分享本文】
    Facebook Twitter LinkedIn LINE Skype EverNote GMail Yahoo Yahoo
    【追蹤本站】
    Facebook Facebook Twitter Plurk