//----------------------------------------------------------------------
// Stefan Pochmann
// www.stefan-pochmann.info
// pochmann@gmx.de
// April 15, 2007
//----------------------------------------------------------------------
//
// This program analyzes the standard 3x3x3 scrambles with 25 random
// moves. Well, actually all lengths up to MOVES moves.
//
// It was inspired by this post in the cubing yahoo group:
//   http://games.groups.yahoo.com/group/speedsolvingrubikscube/message/34938
//
// My first post mentioning my results is here:
//   http://games.groups.yahoo.com/group/speedsolvingrubikscube/message/34941
//
// This version only analyzes the edge orientation, using the orientation
// definition that's used by some blindfolders and that can be described
// as using the centers as fixed reference, and only F and B quarter turns
// changing the orientation of all four involved edges, no other turns
// changing any orientations.
//
// The edge-flip-state of the cube is a 12 bit vector encoded as an
// integer from 0 to 2047, with the highest bit cut off because it can
// be computed from the other 11.
//
// The core idea of the below algorithm is to use dynamic programming.
// Imagine this table:
// 
//     probability[S][N] := probability of arriving in state S
//                          after N moves
//
// Computing that table is the goal of the program. Though actually the
// main part walks over the lengths and then calls the "analyze" method
// for each length n, giving it the table probability[s]. That method
// can then do some statistical analysis.
// 
// Well, two scrambles ending up in the same state with the same number
// of moves aren't necessarily equivalent, as they can end with different
// moves and thus allow different further moves. So imagine this table
// instead:
// 
//     probability[S][N][A] := probability of arriving in state S
//                             after N moves, so that for the next move
//                             only the sides in set A are allowed.
// 
// This table gets initialised with
//
//     probability[0][0][63] = 1
//
// meaning that with probability 1 you'll be in state 0 (no edges flipped)
// after 0 moves and all sides allowed to be turned next (63=111111 binary).
//
// Then the probability gets spread, filling the table.
//
//     probability[*][0][*] represents the scrambles with 0 moves and
//       spreads its probabilites to probability[*][1][*].
//
//     probability[*][1][*] represents the scrambles with 1 move and
//       spreads its probabilites to probability[*][2][*].
//
//     etc.
//
// Each time we've computed the probabilities for a certain scramble
// length, we can sum up over the sets of allowed sides to get the
// table we're actually more interested in, probability[S][N].
//
// To save memory, the table isn't stored for all N. Only for the
// current scramble length and the next scramble length.
//
// So that's it, fairly straightforward actually. I've arranged the
// program so that adaption to analyze the corner orientation and
// the corner permutation should be pretty easy. Edge permutation
// though is a different beast, as its 12! states are quite a lot.
// In any case, for these adaption you only need to change the
// configuration section, telling the program
//
//   - up to how many moves you want to analyze
//   - how many states there are
//   - how to analyze a state probability distribution
//   - how to turn one state into another with a given move
//
//----------------------------------------------------------------------

//----------------------------------------------------------------------
// Initial stuff section.
//----------------------------------------------------------------------

#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;

typedef long double probability;

int bitcount ( int i ) {
  return i ? bitcount( i & (i-1) ) + 1 : 0;
}

//----------------------------------------------------------------------
// This is the configuration section.
//----------------------------------------------------------------------

const int MOVES = 100;          // maximum number of moves
const int STATES = 2048;        // number of possible states

void analyze ( int n, probability probabilities[STATES] ) {
  
  //--- Announce the scramble length.
  cout << endl << "scramble length: " << n << endl;
  
  //--- Sum up the probabilities for each number of flipped edges (0 to 12) ...
  probability flipCtrProbs[13] = {};
  for( int s=0; s<STATES; s++ ){
    int flipped = bitcount( s );
    flipped += flipped & 1;
    flipCtrProbs[flipped] += probabilities[s];
  }
  
  //--- ... and then tell them.
  for( int i=0; i<=12; i++ )
    if( flipCtrProbs[i] )
      printf( "  probability for %2d flipped edges: %.20Lf\n", i, flipCtrProbs[i] );
  
  //--- Find and show the least/most probable state.
  int least = min_element( probabilities, probabilities+STATES ) - probabilities;
  int most = max_element( probabilities, probabilities+STATES ) - probabilities;
  printf( "  state %4d is least probable with %.20Lf\n", least, probabilities[least] );
  printf( "  state %4d is most  probable with %.20Lf\n", most, probabilities[most] );
  

}

//--- Specification of the 4-edges cycles the six sides cause.
int cycles[6][4] = {
  { 0, 1, 2, 3 },
  { 4, 5, 6, 7 },
  { 3, 11, 7, 10 },
  { 1, 8, 5, 9 },
  { 0, 8, 4, 11 },
  { 2, 9, 6, 10 }
};

//--- The naive way to compute the successor state of state s after applying move m.
int nextStateNaive ( int s, int m ) {
  
  //--- Turn the number (0-2047) into an array with 12 bools.
  bool flipped[12] = {};
  for( int i=0; i<11; i++ )
    flipped[11] ^= flipped[i] = (s >> i) & 1;
  
  //--- Apply the cycle permutation for the moved side.
  for( int i=0; i<3; i++ ){
    int a = cycles[m][i];
    int b = cycles[m][i+1];
    swap( flipped[a], flipped[b] );
  }
  
  //--- For F and B, flip the four moved edges.
  if( m >= 4 ){
    for( int i=0; i<4; i++ ){
      int a = cycles[m][i];
      flipped[a] ^= 1;
    }
  }
  
  //--- Translate back to a number (0-2047).
  s = 0;
  for( int i=0; i<11; i++ )
    s |= (flipped[i] << i);
  
  //--- Answer.
  return s;
}

//----------------------------------------------------------------------
// The transition tables.
//----------------------------------------------------------------------

int nextAllow[64][6];
int nextState[STATES][6];

int nextAllowNaive ( int a, int m ) {
  a -= (1 << m);
  for( int i=0; i<6; i++ )
    if( i/2 != m/2 )
      a |= (1 << i);
  return a;
}

void precomputeTransitionTables () {
  
  //--- Precompute the state transition table.
  for( int s=0; s<STATES; s++ )
    for( int m=0; m<6; m++ )
      nextState[s][m] = nextStateNaive( s, m );
  
  //--- Precompute the allowed-moves transition table.
  for( int a=0; a<64; a++ )
    for( int m=0; m<6; m++ )
      nextAllow[a][m] = nextAllowNaive( a, m );
}

//----------------------------------------------------------------------

probability currentProbs[STATES][64];
probability nextProbs[STATES][64];
probability probabilities[STATES];

int main () {
  
  //--- Initialize.
  precomputeTransitionTables();
  currentProbs[0][63] = 1;
  
  //--- Dance...
  for( int n=1; n<=MOVES; n++ ){
    
    //--- Compute nextProbs by walking over all states and all allowed-sides-sets.
    for( int s=0; s<STATES; s++ ){
      for( int a=0; a<64; a++ ){
	
	//--- Nothing to share? Skip!
	if( ! currentProbs[s][a] )
	  continue;
	
	//--- The probability will be spread how many ways? Three ways per allowed side.
	probability spreadProb = currentProbs[s][a] / (3 * bitcount( a ));
	
	//--- Apply the allowed moves, each three times.
	for( int m=0; m<6; m++ ){
	  if( a & (1<<m) ){
	    int s2 = s;
	    for( int x=0; x<3; x++ ){
	      
	      //--- Compute the successor state and allowed-sides-set.
	      s2 = nextState[s2][m];
	      int a2 = nextAllow[a][m];
	      
	      //--- Spread the love... uh, probability, I mean.
	      nextProbs[s2][a2] += spreadProb;
	    }
	  }
	}
      }
    }
    
    //--- Finished the next round, now call it current and clean it again.
    for( int s=0; s<STATES; s++ ){
      for( int a=0; a<64; a++ ){
	currentProbs[s][a] = nextProbs[s][a];
	nextProbs[s][a] = 0;
      }
    }
    
    
    //--- Summarize over the allowed-sides-sets.
    for( int s=0; s<STATES; s++ ){
      probabilities[s] = 0;
      for( int a=0; a<64; a++ )
	probabilities[s] += currentProbs[s][a];
    }
    
    //--- Now let the analyzer do whatever it wants to.
    analyze( n, probabilities );
  }
}