2013-01-10 22:51:54Morris
[線性代數][作業] Gram Schmidt Process
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
| 輸入範例 |
| 檔 名：input.txt |
| 檔案內容： |
| 2 2 |
| 1 1 |
| 2 1 |
| 輸出範例 |
| 檔 名：output.txt |
| 檔案內容： |
| othorgonal vectors: 1 1 |
| 1 /2 -1/2 |
案例2
| 輸入範例 |
| 檔 名：input.txt |
| 檔案內容： |
| 3 3 1 0 1 |
| 2 1 4 3 1 6 |
| 輸出範例 |
| 檔 名：output.txt |
| 檔案內容： |
| othorgonal vectors: 1 0 1 -1 1 1 |
| -1/6 -1/3 1/6 |
案例3
| 輸入範例 |
| 檔 名：input.txt |
| 檔案內容： |
| 4 3 1 0 0 0 0 0 1 1 0 1 1 1 |
| 輸出範例 |
| 檔 名：output.txt |
| 檔案內容： |
| othorgonal vectors:
1 0 0 0 0 0 1 1 0 1 0 0 |
先不考慮大數的輸入，然後為了完成加分題，在此稍做輸出格式的修改。
然後助教沒有規定輸入的型態，特別問過之後居然有可能會有浮點數。
#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
*/