collatz24/main.cpp

#include <utility>
#include <map>
#include <ostream>
#include <iostream>
#include <iomanip>
#include <vector>
//#include <gmp.h>
#include <boost/multiprecision/cpp_dec_float.hpp>

typedef long long num;

//typedef double approx;
//typedef long double approx;
typedef boost::multiprecision::cpp_dec_float_50 approx;

num gcd(num a, num b) {
  while (a) {
    // a, b = b%a, a
    num tmp = b%a; b = a; a = tmp;
  }
  return b;
}

template<typename T>
class ord {
  friend bool operator==(T const &a, T const &b) { return cmp(a, b) == 0; }
  friend bool operator!=(T const &a, T const &b) { return cmp(a, b) != 0; }
  friend bool operator<(T const &a, T const &b) { return  cmp(a, b) < 0; }
  friend bool operator<=(T const &a, T const &b) { return cmp(a, b) <= 0; }
  friend bool operator>(T const &a, T const &b) { return  cmp(a, b) > 0; }
  friend bool operator>=(T const &a, T const &b) { return cmp(a, b) >= 0; }
};

struct roll : ord<roll> {
  num r, s;

  constexpr roll(num r_, num s_) : r(r_), s(s_) {}

  static const roll zero;
  static const roll one;
  static const roll inf;  // don't use in algebra

  friend roll operator*(const roll& i, const roll& j) {
    return { i.r + i.s*j.r , i.s*j.s };
  }

  friend roll operator/(const roll& i, const roll& j) {
    return { (i.r-j.r) / j.s , i.s / j.s };
  }

  static roll lcm(roll i, roll j) {
    if (i == inf or j == inf) return inf;
    num g = gcd(i.s, j.s);
    // if (g == 0) return i == j ? i : inf;  // new case
    if ((i.r-j.r) % g) return inf;
    num l = i.s*j.s/g;
    while (i.r != j.r) {
      if (i.r < j.r) i.r += i.s * ((j.r-i.r+i.s-1)/i.s);
      if (i.r > j.r) j.r += j.s * ((i.r-j.r+j.s-1)/j.s);
    }
    return { i.r , l };
  }

  friend std::ostream& operator<<(std::ostream& os, const roll& i) {
    return os << i.r << "(" << i.s << ")";
  }

  friend int cmp(const roll& i, const roll& j) {
    if (i.s < j.s or (i.s == j.s and i.r < j.r)) return -1;
    if (i.s == j.s and i.r == j.r) return 0;
    return 1;
  }

};

constexpr const roll roll::zero = { 0 , 0 };
constexpr const roll roll::one = { 0 , 1 };
constexpr const roll roll::inf = { 0 , -1 };

struct pix : ord<pix> {
  roll i, j;

  constexpr pix(roll i_, roll j_) : i(i_), j(j_) {}

  static const pix zero;
  static const pix one;

  friend pix operator*(pix a, pix b) {
    if (a == zero or b == zero) return zero;
    roll l = roll::lcm(a.j, b.i);
    if (l == roll::inf) return zero;
    return { a.i * (l/a.j) , b.j * (l/b.i) };
  }

  static pix vert(roll i) {
    return { i , roll::zero };
  }

  static pix horz(roll j) {
    return { roll::zero , j };
  }

  friend std::ostream& operator<<(std::ostream& os, const pix& a) {
    return os << a.i << a.j;
  }

  friend int cmp(const pix& a, const pix& b) {
    if (a.i < b.i or (a.i == b.i and a.j < b.j)) return -1;
    if (a.i == b.i and a.j == b.j) return 0;
    return 1;
  }
};

constexpr const pix pix::zero = { roll::inf , roll::inf };
constexpr const pix pix::one = { roll::one , roll::one };

struct mat {
  std::multimap<pix, approx> m;

  mat() {}
  mat(std::multimap<pix, approx> m_) : m(m_) {}
  mat(approx x) : m{ {pix::one, x} } {}

  static mat zero() { return mat(); }
  static mat one() { return mat({ {pix::one, 1.0} }); }
  static mat single(pix a) {
    mat X;
    X.m.emplace(a, 1.0);
    X.compress();
    return X;
  }

  friend mat operator+(const mat& X, const mat& Y) {
    mat Z;
    Z.m.insert(X.m.begin(), X.m.end());
    Z.m.insert(Y.m.begin(), Y.m.end());
    Z.compress();
    return Z;
  }

  friend mat operator*(const mat& X, const mat& Y) {
    //std::cout << "input: " << X.m.size() << "x" << Y.m.size() << std::endl;
    mat Z;
    for (auto i = X.m.begin(); i != X.m.end(); i++)
      for (auto j = Y.m.begin(); j != Y.m.end(); j++)
        Z.m.emplace(i->first * j->first, i->second * j->second);
    //std::cout << "before: {" << Z.m.size() << "} " << Z << std::endl;
    Z.compress();
    //std::cout << "after: {" << Z.m.size() << "} " << Z << std::endl;
    return Z;
  }

  static mat pow(mat X, num n) {
    mat Y = one();
    for ( ; n; n /= 2) {
      if (n % 2) Y = X*Y;
      if (n != 1) X = X*X;
    }
    return Y;
  }

  approx technical() {
    approx s = 0.0;
    for (auto i = m.begin(); i != m.end(); i++)
      s += i->second;
    return s;
  }

  approx quadranc1() {
    approx s = 0.0;
    for (auto i = m.begin(); i != m.end(); i++) {
      s += i->second / i->first.j.s;
    }
    return s;
  }

  approx inner(roll u, roll v) {
    mat uX = (single(pix::horz(u)) * (*this));
    approx uXv = 0.0;
    for (auto i = uX.m.begin(); i != uX.m.end(); i++) {
      num lcm = roll::lcm(i->first.j, v).s;
      if (lcm != -1)
        uXv += i->second / lcm;
    }
    return uXv;
  }

  friend std::ostream& operator<<(std::ostream& os, const mat& X) {
    if (X.m.empty()) return os << "0";
    for (auto i = X.m.begin(); i != X.m.end(); i++)
      os << i->second << "[" << i->first << "]";
    return os;
  }

private:
  void compress() {
    std::multimap<pix, approx> w;
    for (auto i = m.begin(); i != m.end(); ) {
      std::pair<pix, approx> p = {pix::zero, 0};
      for (p.first = i->first; p.first == i->first; i++)
        p.second += i->second;
      if (p.first != pix::zero and p.second != 0)
        w.insert(p);
    }
    m = w;
  }
};

mat A = mat::single({{2 , 2} , {1 , 1}});
mat B = mat::single({{1 , 2} , {4 , 6}});

mat C_plain = A + B;

mat C(approx p, approx q) { return p*A + q*B; }

void quadranc1_experiment(const approx p, const approx q) {
  mat X = mat::one();
  approx t = 0;
  std::cout << std::setprecision(30);
  for (int i = 0; i < 40; i++) {
    std::cout << std::setw(3) << i << ", "
              << std::setw(30) << X.technical() << ", "
              << (X.quadranc1() / t) << std::endl;
    t = X.quadranc1();
    X = X*C(p,q);
  }
}

void inner_experiment(const roll u, const roll v) {
  mat X = mat::one();
  std::cout << u << v << std::endl;
  for (int i = 0; i < 25; i++) {
    if (i >= 20)
      std::cout << std::setw(3) << i << ", "
                << std::setprecision(3) << std::scientific
                << X.technical() << ", "
                << std::setprecision(30) << std::fixed
                << X.inner(u, v) << std::endl;
    X = X*C(.75,1);  // max eig = 1.0
  }
}

struct rat {
  num a, b;
  rat (num a_, num b_) : a(a_), b(b_) {}
  rat (num a_) : a(a_), b(1) {}
  friend rat operator+(rat x, rat y) { return { x.a*y.b+x.b*y.a , x.b*y.b }; }
  friend rat operator*(rat x, rat y) { return { x.a*y.a , x.b*y.b }; }
  friend rat operator/(rat x, rat y) { return { x.a*y.b , x.b*y.a }; }
  void compress() { num g = gcd(a, b); a /= g; b /= g; }
  friend std::ostream& operator<<(std::ostream& os, const rat& x) {
    return os << x.a << "/" << x.b;
  }
};

rat sb(approx x, approx tol) {
  if (x > 1) return sb(1/x, tol);
  if (x < tol) return rat(0, 1);
  return 1 / (rat(num(1/x), 1) + sb(1/x-num(1/x), tol));
}

void inner_matrix(num u_max, num v_max) {
  for (num u_num = 1; u_num <= u_max; u_num++) {
    for (num v_num = 1; v_num <= v_max; v_num++) {
      roll u = { 0 , u_num };
      roll v = { 0 , v_num };
      mat X = mat::one();
      int i = 0;
      approx t = -1;
      do {
        t = X.inner(u, v);
        X = X*C(.75,1);
        i++;
      } while (i < 20 or (i < 30 and t > 1e-3 and
                          abs(X.inner(u, v)-t) > 1e-8));
      std::cout << std::setfill(' ')
                << "u:(0," << std::setw(2) << u_num
                << "),v:(0," << std::setw(2) << v_num << ") "
                << std::setprecision(3) << std::scientific
                << abs(X.inner(u, v)-t);
      rat r = sb(X.inner(u, v), 1e-4);
      r.compress();
      if (r.b > 100000)
        r = 0;
      std::cout << " " << r << std::endl;
    }
  }
}

void inner_left_triangle(num u_max) {
  for (num u_num = 1; u_num <= u_max; u_num++) {
    for (num u_shift = 1; u_shift <= u_num; u_shift++) {
      roll u = { u_shift , u_num };
      roll v = roll::one;
      mat X = mat::one();
      int i = 0;
      approx t = -1;
      do {
        t = X.inner(u, v);
        X = X*C(.75,1);
        i++;
      } while (i < 20 or (i < 30 and t > 1e-3 and
                          abs(X.inner(u, v)-t) > 1e-8));
      std::cout << std::setfill(' ')
                << u << " "
                << std::setprecision(3) << std::scientific
                << abs(X.inner(u, v)-t);
      rat r = sb(X.inner(u, v), 1e-4);
      r.compress();
      if (r.b > 100000)
        r = 0;
      std::cout << " " << r << std::endl;
    }
  }
}

void inner_top_triangle(num v_max) {
  for (num v_num = 1; v_num <= v_max; v_num++) {
    for (num v_shift = 1; v_shift <= v_num; v_shift++) {
      roll u = roll::one;
      roll v = { v_shift , v_num };
      mat X = mat::one();
      int i = 0;
      approx t = -1;
      do {
        t = X.inner(u, v);
        X = X*C(.75,1);
        i++;
      } while (i < 20 or (i < 30 and t > 1e-3 and
                          abs(X.inner(u, v)-t) > 1e-8));
      std::cout << std::setfill(' ')
                << v << " "
                << std::setprecision(3) << std::scientific
                << abs(X.inner(u, v)-t);
      rat r = sb(X.inner(u, v), 1e-4);
      r.compress();
      if (r.b > 100000)
        r = 0;
      std::cout << " " << r << std::endl;
    }
  }
}

void inner_left_zero(num u_max) {
  for (num u_num = 1; u_num <= u_max; u_num++) {
    {
      num u_shift = 0;
      roll u = { u_shift , u_num };
      roll v = roll::one;
      mat X = mat::one();
      int i = 0;
      approx t = -1;
      do {
        t = X.inner(u, v);
        X = X*C(.75,1);
        i++;
      } while (i < 20 or (i < 30 and t > 1e-3 and
                          abs(X.inner(u, v)-t) > 1e-8));
      std::cout << std::setfill(' ')
                << u << " "
                << std::setprecision(3) << std::scientific
                << abs(X.inner(u, v)-t);
      rat r = sb(X.inner(u, v), 1e-4);
      r.compress();
      if (r.b > 100000)
        r = 0;
      std::cout << " " << r << std::endl;
    }
  }
}

void inner_top_zero(num v_max) {
  {
    for (num v_num = 1; v_num <= v_max; v_num++) {
      roll u = roll::one;
      num v_shift = 0;
      roll v = { v_shift , v_num };
      mat X = mat::one();
      int i = 0;
      approx t = -1;
      do {
        t = X.inner(u, v);
        X = X*C(.75,1);
        i++;
      } while (i < 20 or (i < 30 and t > 1e-3 and
                          abs(X.inner(u, v)-t) > 1e-8));
      std::cout << std::setfill(' ')
                << v << " "
                << std::setprecision(3) << std::scientific
                << abs(X.inner(u, v)-t);
      rat r = sb(X.inner(u, v), 1e-4);
      r.compress();
      if (r.b > 100000)
        r = 0;
      std::cout << " " << r << std::endl;
    }
  }
}

void quadranc1_2nd_eig() {
  const approx p = approx(3)/4;
  const approx q = 1;
  mat X = mat::one();
  approx t = 0;
  std::cout << std::setprecision(30);
  for (int i = 0; i < 20; i++) {
    approx t2 = X.quadranc1() - approx(14)/15;
    std::cout << std::setw(3) << i << ", "
              << std::setprecision(3) << std::scientific
              << X.technical() << ", "
              << std::setprecision(30) << std::fixed
              << (t2 / t) << std::endl;
    t = t2;
    X = X*C(p,q);
  }
}

void fine_experiment() {
  const approx p = approx(3)/4;
  const approx q = 1;
  mat X = mat::single(pix::horz(roll::one));
  std::cout << std::setprecision(30);
  for (int i = 0; i < 40; i++) {

    approx t = 0;
    for (auto j = X.m.begin(); j != X.m.end(); j++)
      if (j->first.j.r == 1)
        t += j->second;

    std::cout << std::setw(3) << i << ", "
              << std::setprecision(3) << std::scientific
              << X.technical() << ", "
              << std::setprecision(30) << std::fixed
              << t << std::endl;
    X = X*C(p,q);
  }
}

void fine_experiment_special() {
  std::vector<std::pair<num, approx>> oui = { { 1 , 1.0 } };

  for (num i = 0; i < 60; i++) {
    std::vector<std::pair<num, approx>> non;
    for (auto j = oui.begin(); j != oui.end(); j++) {
      if (j->first % 6 == 4)
        non.push_back({ (j->first-1)/3 , j->second });
      non.push_back({ j->first*2 , j->second*3/4 });
    }
    oui = non;
    approx t = 0;
    for (auto j = oui.begin(); j != oui.end(); j++)
      t += j->second;
    std::cout << i << " " << std::setprecision(10) << t << std::endl;
  }
}


int main() {

  // quadranc1_experiment(.75,1);  // <-- you will want to vary this

  /*
  for (int m = 1; m < 8; m++) {
    for (int n = 1; n < 8; n++) {
      inner_experiment({ 0 , m }, { 0 , n });
    }
    }*/

  // inner_left_matrix(30);
  // inner_left_triangle(50);
  // inner_left_zero(100);

  // inner_top_triangle(30);
  // inner_top_zero(100);

  // quadranc1_2nd_eig();  // -1/4 pure
  // fine_experiment();
  fine_experiment_special();  // 27/22?

}