fork(2) download
copy 
  1. #include <vector>
  2. #include <iostream>
  3. #include <algorithm>
  4. #include <functional>
  5. #include <numeric>
  6.  
  7. using namespace std;
  8.  
  9. long gcd(long a, long b) {
  10.     if (a < 0 || b < 0) return gcd(abs(a), abs(b));
  11.     else if (min(a, b) == 0)
  12.         return max(1l, max(a, b));
  13.     else if (a == b)
  14.         return a;
  15.     else
  16.         return gcd(min(a, b), max(a, b) - min(a, b));
  17. }
  18.  
  19. struct vec {
  20.     vec(vector<long> v) : v_(v) {}
  21.     vector<long>& v() { return v_; }
  22.     const vector<long>& v() const { return v_; }
  23.  
  24.     vec operator+(const vec&) const;
  25.     vec operator*(long) const;
  26.     long operator*(const vec&) const;
  27.  
  28.     void reduceByGCD();
  29.  
  30.     bool isZero() const { return v() == vector<long>(v().size()); }
  31.  
  32. private:
  33.     vector<long> v_;
  34.     long gcd() const;
  35. };
  36.  
  37. long
  38. vec::gcd() const
  39. {
  40.     if (v_.empty()) return 1;
  41.     long res = v_[0];
  42.     for (int i = 1; i < v_.size(); ++i)
  43.         res = ::gcd(res, v_[i]);
  44.     return res;
  45. }
  46.  
  47. void
  48. vec::reduceByGCD()
  49. {
  50.     long div = gcd();
  51.     for (auto &a: v_) a /= div;
  52. }
  53.  
  54. vec
  55. vec::operator+(const vec& r) const
  56. {
  57.     vector<long> res(v().size());
  58.     transform(this->v_.begin(), this->v_.end(), r.v_.begin(), res.begin(), plus<long>());
  59.     return vec(res);
  60. }
  61.  
  62. vec
  63. vec::operator*(long k) const
  64. {
  65.     vector<long> res(v().size());
  66.     transform(this->v_.begin(), this->v_.end(), res.begin(), bind(multiplies<long>(), placeholders::_1, k));
  67.     return vec(res);
  68. }
  69.  
  70. long
  71. vec::operator*(const vec& r) const
  72. {
  73.     return inner_product(v_.begin(), v_.end(), r.v().begin(), 0l);
  74. }
  75.  
  76. vec gschmidt_and_reduce(const vec& what, const vec& other)
  77. {
  78.     long sqNorm = other*other;
  79.     vec res = what*sqNorm  + other*(what*other)*-1;
  80.     res.reduceByGCD();
  81.     return res;
  82. }
  83.  
  84. vector<vec>
  85. gschmidt_and_reduce(const vector<vec>& vs)
  86. {
  87.     vector<vec> res;
  88.     for (size_t i = 0; i < vs.size(); ++i) {
  89.         vec cur = vs[i];
  90.         for (int k = 0; k < i; ++k)
  91.             cur = gschmidt_and_reduce(cur, res[k]);
  92.         res.push_back(cur);
  93.     }
  94.     return res;
  95. }
  96.  
  97. int main() {
  98.     vector<vec> A = {
  99.         vec({ 3, 5, 4, -4 }),
  100.         vec({ -2, -3, -2, 1 }),
  101.         vec({ -5, -9, -6, 8 }),
  102.         vec({ 1, 0, 0, 0 }),
  103.         vec({ 0, 1, 0, 0 }),
  104.         vec({ 0, 0, 1, 0 }),
  105.         vec({ 0, 0, 0, 1 }),
  106.     };
  107.     vector<vec> orthogonilizedA = gschmidt_and_reduce(A);
  108.     cout << "Orthogonalized a:"<< endl;
  109.     for (int i = 0; i < orthogonilizedA.size(); ++i) {
  110.         if (!orthogonilizedA[i].isZero()) {
  111.             for (auto coordinate: orthogonilizedA[i].v())
  112.                 cout << coordinate << ", ";
  113.             cout << endl;
  114.         }
  115.         if (i == 2)
  116.             cout << "Orthogonal complement of a:" << endl;
  117.     }
  118.  
  119.     vector<vec> B = {
  120.         vec({ -4, -5, -3, 3 }),
  121.         vec({ 2, 7, 0, -6 }),
  122.         vec({ -1, 0, -1, 0 }),
  123.         vec({ 1, 0, 0, 0 }),
  124.         vec({ 0, 1, 0, 0 }),
  125.         vec({ 0, 0, 1, 0 }),
  126.         vec({ 0, 0, 0, 1 }),
  127.     };
  128.     vector<vec> orthogonilizedB = gschmidt_and_reduce(B);
  129.     cout << "Orthogonalized b:"<< endl;
  130.     for (int i = 0; i < orthogonilizedB.size(); ++i) {
  131.         if (!orthogonilizedB[i].isZero()) {
  132.             for (auto coordinate: orthogonilizedB[i].v())
  133.                 cout << coordinate << ", ";
  134.             cout << endl;
  135.         }
  136.         if (i == 2)
  137.             cout << "Orthogonal complement of b:" << endl;
  138.     }
  139.     return 0;
  140. }
  141.  
Success #stdin #stdout 0s 4272KB
comments ()
stdin
copy 
Standard input is empty
stdout
copy 
Orthogonalized a:
3, 5, 4, -4, 
-1, -1, 0, -2, 
-1, -9, 17, 5, 
Orthogonal complement of a:
5, -3, -1, -1, 
Orthogonalized b:
-4, -5, -3, 3, 
-42, 36, -61, -57, 
-16, 38, -3, 39, 
Orthogonal complement of b:
3, 0, -3, 1,
https://ideone.com/yg2pFB
language:
C++14 (gcc 8.3)
created:
4 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!