Morris' Blog 我在摸索我存在的意義、生命的意義、內涵及價值。 沒有真正神手般的厲害，套用模式是大部分的情況。 沒那麼厲害的我存在的意義為何，襯托出神的存在？ 有時候不能怪別人不給我機會，我已經開始省思， 為什麼我沒有給別人機會的原因，而是全要靠別人。 有時候我認為我與眾不同，那是因為我做不到很多人能辦到的事情。 我聽不到、玩不起來、說得不夠明瞭、理解得不夠多， 我缺少的語文能力，就是一切一切無能的判定。 為此，我不曉得能做些對別人有幫助的事情。 我知道我自己是個不服輸的性格，但現在我不得不認輸。 我沒有能力，因此我必須認輸。 我到底來這世界作什麼？處處充滿不確定性， 而我總是無法出類拔萃，想當個一般人似乎也回不了頭了。 越強大的挑戰需要越強力的支持，而我的支持是什麼？ 最近嘗試與人談話，但都是失敗的結局，沒人能接受， 不需要言語，不需要溝通，一切所想都表達不出來。

49愛的鼓勵 6訂閱站台

Linear Algebra: Programing Homework3

2012/12/20

 

 

 

 

題　　目：給定幾個向量，使用這幾個向量進行Gram Schmidt Process，用以取得互相正交的向量集合

 

 

 

加分條件：若可以算出orthonormal的向量集合，即可加10分。

 

輸　　入：檔名為input.txt的檔案，內容為：

                    第一列       ：   (1)每個向量的element個數M

(2)向量的個數N ( 0<N<=M)

                    第二 ~ N+1列 ： 每列分別是拿來做 Gram Schmidt Process的向量

 

輸　　出：檔名為output.txt的檔案，內容為：

(1).   N個經過Gram Schmidt Process之後計算出來的向量

(2).   N個orthonormal的向量(加分題)

 

注意：輸入及輸出每一列的每個數值皆以一個TAB鍵做分隔

 

 

收件時間：1/10 24:00以前（未避免網路壅塞，請儘早上傳，不接受遲交）

 

 

評分說明：

1.    一個測試案例(下述案例)執行正確，給60分

 

2.    兩個測試案例(未給)執行正確，給80分

 

3.    三個測試案例(未給)執行正確，給100分

 

4.    避免因浮點數運算造成的誤差，測試數據的計算過程中不會出現無限小數，四捨五入至小數點第二位。

若轉換過程中以分數計算並表示結果，助教將斟酌加分。

抄襲者雙方則以零分計算。

案例1

 

 

 

案例2   

 

 

 

案例3   

 

 

---

先不考慮大數的輸入，然後為了完成加分題，在此稍做輸出格式的修改。

然後助教沒有規定輸入的型態，特別問過之後居然有可能會有浮點數。


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
#include <iostream>
#define LL long long
using namespace std;
class Frac {
    public:
        LL a, b;
        Frac() {
            a = 0, b = 1;
        }
        Frac(LL x, LL y) {
            a = x, b = y;
            reduce();
        }
        Frac operator+(const Frac &y) {
            LL ta, tb;
            tb = this->b/gcd(this->b, y.b)*y.b;
            ta = this->a*(tb/this->b) + y.a*(tb/y.b);
            Frac z(ta, tb);
            return z;
        }
        Frac operator-(const Frac &y) {
            LL ta, tb;
            tb = this->b/gcd(this->b, y.b)*y.b;
            ta = this->a*(tb/this->b) - y.a*(tb/y.b);
            Frac z(ta, tb);
            return z;
        }
        Frac operator*(const Frac &y) {
            LL tx, ty, tz, tw, g;
            tx = this->a, ty = y.b;
            g = gcd(tx, ty), tx /= g, ty /= g;
            tz = this->b, tw = y.a;
            g = gcd(tz, tw), tz /= g, tw /= g;
            Frac z(tx*tw, ty*tz);
            return z;
        }
        Frac operator/(const Frac &y) {
            LL tx, ty, tz, tw, g;
            tx = this->a, ty = y.a;
            g = gcd(tx, ty), tx /= g, ty /= g;
            tz = this->b, tw = y.b;
            g = gcd(tz, tw), tz /= g, tw /= g;
            Frac z(tx*tw, ty*tz);
            return z;
        }
        void print() {
            printf("%lld", a);
            if(b != 1)
                printf("/%lld", b);
        }
    private:
        static LL gcd(LL x, LL y) {
            if(!y)  return x;
            if(x < 0)   x *= -1;
            if(y < 0)   y *= -1;
            LL t;
            while(x%y)
                t = x, x = y, y = t%y;
            return y;
        }
        void reduce() {
            LL g = gcd(a, b);
            a /= g, b /= g;
            if(b < 0)   a *= -1, b *= -1;
        }
};
Frac norm_square(vector<Frac>& v) {
    Frac sum(0,1);
    for(vector<Frac>::iterator it = v.begin();
        it != v.end(); it++) {
        sum = sum + (*it)*(*it);
    }
    return sum;
}
Frac inner_product(vector<Frac>& a, vector<Frac>& b) {
    if(a.size() != b.size()) {
        puts("error");
        return Frac(0,1);
    }
    Frac res(0,1);
    for(int i = 0; i < a.size(); i++)
        res = res + a[i]*b[i]; // standard inner product
    return res;
}
void floatSol(vector<double> v[], int M, int N) {
    int i, j, k;
    puts("<float expression>");
    puts("othorgonal vectors:");
    for(i = 0; i < N; i++) {
        double v_2, inner_res, K;
        for(j = 0; j < i; j++) {
            for(k = 0, v_2 = 0; k < M; k++)
                v_2 += v[j][k]*v[j][k];
            for(k = 0, inner_res = 0; k < M; k++)
                inner_res += v[i][k]*v[j][k];
            K = inner_res / v_2; // ||v[i]|| length of vi
            for(k = 0; k < M; k++)
                v[i][k] = v[i][k] - K*v[j][k];
        }
        for(k = 0, v_2 = 0; k < M; k++)
            v_2 += v[j][k]*v[j][k];
        if(fabs(v_2) < 1e-6)  continue;
        for(j = 0; j < M; j++) {
            if(j)   putchar(' ');
            cout << v[i][j];
        }
        puts("");
    }
    // orthonormal vector process
    puts("\northonormal vectors:");
    for(i = 0; i < N; i++) {
        double v_2;
        for(j = 0, v_2 = 0; j < M; j++)
            v_2 += v[i][j]*v[i][j];
        v_2 = sqrt(v_2);
        if(fabs(v_2) < 1e-6)
            continue;
        for(j = 0; j < M; j++) {
            if(j)   putchar(' ');
            cout << v[i][j]/v_2;
        }
        puts("");
    }
}
void FracSol(vector<Frac> v[], int M, int N) {
    int i, j, k;
    puts("<fraction expression>");
    puts("othorgonal vectors:");
        for(i = 0; i < N; i++) {
        Frac v_2, inner_res, K;
        for(j = 0; j < i; j++) {
            v_2 = norm_square(v[j]);
            inner_res = inner_product(v[i], v[j]);
            K = inner_res / v_2; // ||v[i]|| length of vi
            for(k = 0; k < M; k++)
                v[i][k] = v[i][k] - K*v[j][k];
        }
        v_2 = norm_square(v[i]);
        if(v_2.a == 0)  continue;
        for(j = 0; j < M; j++) {
            if(j)   putchar(' ');
            v[i][j].print();
        }
        puts("");
    }
    // orthonormal vector process
    puts("\northonormal vectors:");
    for(i = 0; i < N; i++) {
        Frac sq;
        sq = norm_square(v[i]);
        if(sq.a == 0)   continue;
        for(j = 0; j < M; j++) {
            Frac res(v[i][j].a, v[i][j].b*sq.a), K(1,1);
            long long tmp = sq.a*sq.b, sqt = sqrt(tmp); // _/tmp
            for(k = 2; k <= sqt && k*k <= tmp; k++) {
                while(tmp%(k*k) == 0) {
                    tmp /= k*k;
                    K = K * Frac(k, 1);
                }
            }
            res = res * K;
            if(j)   putchar(' ');
            if(tmp != 1 && (res.a == 1 || res.a == -1)) {
                if(res.a < 0)   putchar('-');
            } else    printf("%lld", res.a);
            if(tmp != 1 && res.a)    printf("_/%lld", tmp);
            if(res.b != 1)  printf("/%lld", res.b);
        }
        puts("");
    }
}
double parseDouble(char s[]) {
    double d;
    sscanf(s, "%lf", &d);
    return d;
}
Frac parseFrac(char s[]) {
    LL i, j, cnt = 0;
    for(i = 0; s[i]; i++)
        if(s[i] == '.') break;
    if(s[i] == '.') {
        for(i++; s[i]; i++) {
            s[i-1] = s[i];
            cnt++;
        }
        s[i-1] = '\0';
        sscanf(s, "%lld", &i);
        j = 1;
        while(cnt) {
            j *= 10;
            cnt--;
        }
    } else {
        sscanf(s, "%lld", &i);
        j = 1;
    }
    return Frac(i, j);
}
int main() {
    int M, N, i, j, k, x;
    char cmd[128];
    freopen("input.txt", "r+t", stdin);
    freopen("output.txt", "w+t", stdout);
    while(scanf("%d %d", &M, &N) == 2) {
        vector<Frac> v1[N];
        vector<double> v2[N];
        for(i = 0; i < N; i++) {
            for(j = 0; j < M; j++) {
                scanf("%s", cmd);
                v2[i].push_back(parseDouble(cmd));
                v1[i].push_back(parseFrac(cmd));
            }
        }
        floatSol(v2, M, N);
        puts("");
        FracSol(v1, M, N);
        puts("");
    }
    return 0;
}
/*
4    3
1    0    0    0
0    0    1    1
0    1    1    1

3    3
1    0    1
2    1    4
3    1    6

*/