有理数を使って、ベクトルの正規化

水曜日に困っていたところ。

• 割合の合計して1にならない - Seeking for my unique color.

コメントで有理数を使うとよいと教えてもらった。boostのほうは分数の形で渡さないといけないので、今回のにはちょっと使えなさそう。ということで、こちらのを使わせてもらうことに。

• Spaghetti Source - 有理数

こいつで実験してみると

/Users/syou6162/cpp% ./a.out
The sum is 1.0
9.9999999999991e-0119

となり、通った。ただ、vectorの要素数を多くするとオーバーフローが起こってしまうので、うまくいかないときはあるようだ。以下コード。

#include <iostream>
#include <limits>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <boost/foreach.hpp>
// ref http://www.prefield.com/algorithm/math/rational.html
typedef long long Integer;
Integer gcd(Integer a, Integer b) { return a > 0 ? gcd(b % a, a) : b; }
struct rational {
  Integer p, q;
  void normalize() { // keep q positive
    if (q < 0) p *= -1, q *= -1;
    Integer d = gcd(p < 0 ? -p : p, q);
    if (d == 0) p = 0,  q = 1;
    else        p /= d, q /= d;
  }
  rational(Integer p, Integer q = 1) : p(p), q(q) {
    normalize();
  }
  rational &operator += (const rational &a) {
    p = a.q * p + a.p * q; q = a.q * q; normalize();
    return *this;
  }
  rational &operator -= (const rational &a) {
    p = a.q * p - a.p * q; q = a.q * q; normalize();
    return *this;
  }
  rational &operator *= (const rational &a) {
    p *= a.p; q *= a.q; normalize();
    return *this;
  }
  rational &operator /= (const rational &a) {
    p *= a.q; q *= a.p; normalize();
    return *this;
  }
  rational &operator - () {
    p *= -1;
    return *this;
  }
};
rational operator + (const rational &a, const rational &b) {
  return rational(a) += b;
}
rational operator * (const rational &a, const rational &b) {
  return rational(a) *= b;
}
rational operator - (const rational &a, const rational &b) {
  return rational(a) -= b;
}
rational operator / (const rational &a, const rational &b) {
  return rational(a) /= b;
}
bool operator < (const rational &a, const rational &b) { // avoid overflow
  return (long double) a.p * b.q < (long double) a.q * b.p;
}
bool operator <= (const rational &a, const rational &b) {
  return !(b < a);
}
bool operator > (const rational &a, const rational &b) {
  return b < a;
}
bool operator >= (const rational &a, const rational &b) {
  return !(a < b);
}
bool operator == (const rational &a, const rational &b) {
  return !(a < b) && !(b < a);
}
bool operator != (const rational &a, const rational &b) {
  return (a < b) || (b < a);
}

std::vector<rational> normalize_vector(const std::vector<rational>& v) {
  rational sum = 0.0;
  std::vector<rational> result;
  BOOST_FOREACH(rational d, v) {
	sum += d;
  }
  BOOST_FOREACH(rational d, v) {
	result.push_back(d / sum);
  }
  return result;
}

int main(int argc, char *argv[]) {
  std::vector<rational> tmp;
  for (int i = 0; i < 10000; i++) {
	tmp.push_back(rational(2.49148741798853211e-06 * 1000000.0));
  }
  
  std::vector<rational> v = normalize_vector(tmp);

  double sum = 0.0;
  BOOST_FOREACH(rational r, v) {
	sum += (double) r.p / (double) r.q;
  }
  
  if(-1.0e-13 > 1.0 - sum || 1.0e-13 < 1.0 - sum) {
	fprintf(stderr, "The sum is not 1.0.\n");
  } else {
	fprintf(stderr, "The sum is 1.0.\n");
  }
  std::cout << printf("%15.13e", sum) << std::endl;

  return 0;
}