Solution for Live Archive 6495 - Probability Paradox

The fraction class is not that good, (very very bad implementation sorry )

/*
 * 6495.cpp
 *
 *  Created on: Jul 26, 2014
 *      Author: ghooo
 *       Email: ghooo.ghooo@gmail.com
 */

#include <cstring>
#include <map>
#include <deque>
#include <queue>
#include <stack>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <algorithm>
#include <vector>
#include <set>
#include <complex>
#include <list>
#include <climits>
#include <cctype>
#include <bitset>
//#include <windows.h>

using namespace std;

#define pb push_back
#define all(v) v.begin(),v.end()
#define rall(v) v.rbegin(),v.rend()
#define sz(v) ((ll)v.size())
#define rep(i,m) for(int i=0;i<(int)(m);++i)
#define rep2(i,n,m) for(int i=n;i<(int)(m);i++)
#define For(it,c) for(__typeof(c.begin()) it=c.begin();it!=c.end();++it)
#define reset(a,b) memset(a,b,sizeof(a))
#define mp make_pair
#define dot(a,b) ((conj(a)*(b)).X)
#define X real()
#define Y imag()
#define length(V) (hypot((V).X,(V).Y))
#define vect(a,b) ((b)-(a))
#define cross(a,b) ((conj(a)*(b)).imag())
#define normalize(v) ((v)/length(v))
#define rotate(p,about,theta) ((p-about)*exp(point(0,theta))+about)
#define pointEqu(a,b) (comp(a.X,b.X)==0 && comp(a.Y,b.Y)==0)
#define clrq(x) while(!x.empty()) x.pop();
#define clrvv(v) rep(i,sz(v))v[i].clear();

typedef stringstream ss;
typedef long long ll;
typedef pair<int, int> pii;
typedef vector<pii> vpii;
typedef vector<string> vs;
typedef vector<int> vi;
typedef vector<double> vd;
typedef vector<vector<int> > vii;
typedef long double ld;
typedef complex<double> point;
typedef pair<point, point> segment;
typedef pair<double, point> circle;
typedef vector<point> polygon;
typedef unsigned long long ull;
typedef long long ll;

//==============================================================
// handling triples
typedef pair<ll, pair<ll, ll> > triple;
#define tfirst first
#define tsecond second.first
#define tthird second.second
#define mt(x,y,z) mp(x,mp(y,z))
//---------------------------------------------------------------

const int oo = (int) 1e9;
const double PI = 2 * acos(0.0);
const double EPSILON = 1e-9;
struct fraction{
 ll nom, denom;
 fraction(ll n, ll d):nom(n),denom(d){fix();};
 fraction(ll n):nom(n),denom(1){};
 fraction():nom(0),denom(1){};
 void div(int val){
  denom*=val;
  //fix();
 }
 fraction fix(){
  ll g = gcd(nom, denom);
  nom/=g,denom/=g;
//  if(denom<0) denom*=-1, nom*=-1;
  return *this;
 }
 ll gcd(ll a, ll b){
  a = llabs(a), b = llabs(b);
  if(a == 0 || b == 0) return max(a,b);
  while(b!=0){
   ll t = b;
   b = a%b;
   a = t;
  }
  return a;
 }
 inline fraction operator-(){
  fraction ret =  fraction(-nom,denom);
  //ret.fix();
  return ret;
 }
 inline fraction operator+(const fraction &f){
  fraction ret = fraction(nom*f.denom+f.nom*denom,denom*f.denom);
  ret.fix();
  return ret;
 }
 inline fraction operator+=(const fraction&f){
  ll d = denom*f.denom, n = nom*f.denom+f.nom*denom;
  nom=n, denom=d;
  return *this;
 }
 inline fraction operator*(const fraction &f){
  fraction ret = fraction(nom*f.nom,denom*f.denom);
  ret.fix();
  return ret;
 }
 inline fraction operator*=(const fraction&f){
  nom *=f.nom,denom*=f.denom;
  return *this;
 }
 inline fraction operator/(const fraction &f){
  fraction ret = fraction(nom*f.denom,denom*f.nom);
  ret.fix();
  return ret;
 }
 inline fraction operator/=(const fraction&f){
  nom*=f.denom,denom*=f.nom;
  return *this;
 // fix();
 }
 void print(){
  cout << nom<<'/'<<denom;
 }
};
enum SOL {
 NOSOL, UNIQUE, INFINIT
};

inline bool zeros(vector<fraction> & row) {

// foreach(i, row)
 For(i,row)
 {
//  if (fabs(*i - 0) > EPSILON)
  if((*i).nom)
   return false;
 }
 return true;
}

inline bool swapRow(int i, vector<vector<fraction> > &mat) {
 rep2(j,i+1,mat.size()-1) {
//  if (fabs(mat[j][i] - 0) > EPSILON) {
  if(mat[j][i].nom){
   swap(mat[i], mat[j]);
   return true;
  }
 }
 return false;
}

inline void devideRow(int i, fraction x, vector<vector<fraction> > &mat) {
// if(x.nom/x.denom == 1) return;
 For(j, mat[i]){
  if((*j).nom == 0) continue;
  (*j) = (*j) / x;
//  (*j) /= x;

 }
}

inline void makeZero(int i, int j, vector<vector<fraction> > &mat) {
 fraction d = -mat[j][i];
 rep(k,mat[i].size()) {
  if(mat[i][k].nom == 0) continue;
  mat[j][k] = mat[j][k] + d * mat[i][k];
//  mat[j][k] += d*mat[i][k];
 }
}
void printMat(vector<vector<fraction> > &mat){
 rep(i,sz(mat)){
  rep(k,sz(mat[i])){
   mat[i][k].print(); cout << " ";
  }
  cout << endl;
 }
 cout << endl;
}

SOL solveLinearEquation(vector<vector<fraction> > &mat) {
 rep(i,mat.size()) {
  if (zeros(mat[i])) {
   mat.erase(mat.begin() + i);
   i--;
   continue;
  }
  if (i == (int) mat[i].size())
   break;
  //if (fabs(mat[i][i] - 0) < EPSILON)
  if(mat[i][i].nom == 0)
   if (!swapRow(i, mat))
    return NOSOL;
  devideRow(i, mat[i][i], mat);
  //printMat(mat);
  rep(j,mat.size()) {
   if (i >= j)
    continue;
   if(mat[j][i].nom == 0) continue;
   makeZero(i, j, mat);
  // printMat(mat);
  }
 // printMat(mat);
 // cout << endl;
 }
 //printMat(mat);
 for(int i = sz(mat)-1; i >= 0; i--){
  fraction f = -mat[i][sz(mat)];
  for(int k = i+1; k < sz(mat); k++){
   f = f+(mat[k][sz(mat)] * mat[i][k]);
//   f.fix()
  }
  mat[i][sz(mat)] = -f;
 }
 //printMat(mat);
 if (mat.size() + 1 == mat[0].size())
  return UNIQUE;
 if (mat.size() + 1 < mat[0].size())
  return INFINIT;
 return NOSOL;
}
void test(){
 vector<fraction> v;v.pb(fraction(3,5));v.pb(fraction(2,7));v.pb(fraction(2));
 vector<vector<fraction> > mat;
 mat.pb(v);
 v.clear();v.pb(fraction(3));v.pb(fraction(2));v.pb(fraction(12));//v.pb(5);v.pb(3);v.pb(10);
 mat.pb(v);
 solveLinearEquation(mat);
}

inline int getBin(string s){
 reverse(all(s));
 int ret = 0;
 rep(i,sz(s)){
  if(s[i]=='H')
   ret |= (1<<i);
 }
 return ret;
}
fraction mem[1<<9];
void build(string s, string t){
 int a = getBin(s), b = getBin(t);
 vector<vector<fraction> >mat;
 int n = 1<<sz(s);
 mat.resize(n);
 rep(i,n){
  mat[i].resize(n+1, fraction(0));
  if(i == a) mat[i][i] = fraction(1), mat[i][n] = fraction(-1);
  else if(i == b) mat[i][i] = 1, mat[i][n] = 0;
  else {
   int ih = ((i<<1)&(~n)), it = ih|1;
   mat[i][i] = fraction(2), mat[i][ih] = mat[i][ih]+fraction(-1), mat[i][it] = mat[i][it]+fraction(-1);
  }
 }
 solveLinearEquation(mat);
// printMat(mat);
 rep(i,n){
  rep(k,n){
   if(mat[i][k].nom != 0){
    mem[i] = -mat[i][n];
   }
  }
 }
// printMat(mat);
}
int vis[1<<9], ID, target, n;
pair<fraction,fraction> rec(int i, bool flag =false){
 if(vis[i] == ID) return mp(mem[i],fraction(0));
 if(i == target && flag) return mp(fraction(0),1);
 pair<fraction,fraction>p1,p2;
 int ih = ((i<<1)&(~n)), it = ih|1;
 p1 = rec(ih,true),p1.first = p1.first * fraction(1,2), p1.second = p1.second * fraction(1,2);
 p2 = rec(it,true),p2.first = p2.first * fraction(1,2), p2.second = p2.second * fraction(1,2);
 return mp(p1.first + p2.first, p1.second+p2.second);
}
void build2(string s, string t){
 int a = getBin(s), b = getBin(t);
 mem[a] = fraction(1), mem[b] = fraction(0);
 ID++;
 vis[a]=vis[b]=ID;
 int len = sz(s);
 n = (1<<len);
 rep(i,(1<<len)){
  if(i == a || i == b) continue;
  target = i;
  pair<fraction, fraction> p = rec(i);
  vis[i] = ID;
  mem[i] = p.first / (fraction(1) + -p.second);
 }
}
int len;
fraction rec(int i, int l){
 if(l == len) return mem[i];
 int ih = ((i<<1)&(~n)), it = ih|1;
 fraction ret = rec(ih,l+1) + rec(it,l+1);
 ret.div(2);
 return ret;
}
int main() {
 ios_base::sync_with_stdio(false);
 cin.tie(NULL);
 cout.tie(NULL);
#ifndef ONLINE_JUDGE
 freopen("in.txt", "r", stdin);
 //freopen("out.txt", "w", stdout);
#endif
 string s, t;
 while(cin >> s,s!="$"){
  cin >> t;
  build(s,t);
  len = sz(t);
  n = (1<<len);
  rec(0,0).fix().print();
  cout << endl;
 }
 return 0;
}