Linear Algebra: Programing Homework3
2012/12/20
目：定向量，使用向量行Gram Schmidt Process，用以取得互相正交的向量集合
加分件：若可以算出orthonormal的向量集合，即可加10分。
入：名input.txt的案，容：
第一列 ： (1)每向量的elementM
(2)向量的N ( 0<N<=M)
第二 ~ N+1列 ： 每列分是拿做 Gram Schmidt Process的向量
出：名output.txt的案，容：
(1). NGram Schmidt Process之後算出的向量
(2). Northonormal的向量(加分)
注意：入及出每一列的每值皆以一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];
}
nbsp; 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
*/
文章定位：
