#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

// uncompressor (arithmetic) for a sparse file, which has random 0s and 1s in it.
// This (un)compressor uses a probabilistic model that
// expects the bias of the coin to be from beta distribn with same params as defined in ACdefns.h
// As the data is read in, the model learns.
//
// usage: ACdecode N <  file.ACZ > file.decoded
//   where N is the file length

// (c) David J.C. MacKay
// License: GPL                                  http://www.gnu.org/copyleft/gpl.html

// Originates from:  http://www.inference.phy.cam.ac.uk/mackay/itprnn/code/c/compress/

#include "ACdefns.h"

int main( int argc , char *argv[] ){
  FILE *fp;
  double c0=BETA0, c1=BETA1 , cT , p0  ; // the counts of the data so far,
                                       // and the predictive probability
  long a=0, b=ONE , w ;
  long boundary ; 
  int c ;
  int model_needs_updating =  1;

  int N = 10000 ;
  long u=0, v=ONE;  // the binary interval
  long halfway ;    // midpoint of the binary interval

  // Read N from command line
  if(argc>1) {
    sscanf(argv[argc-1], "%d", &N ) ;
  } else {
    fprintf(stderr,"Assuming " ) ;
  }
  fprintf(stderr,"N=%d\n", N ) ;
  // done reading command line
  // N is number of source characters remaining. N runs down from N to 0.

  fp=stdin;

  while( (N>0) && ((c=getc(fp)) != EOF) ) {
    // (u,v) is the current "encoded alphabet" binary interval, and halfway is its midpoint.
    // (a,b) is the current "source alphabet" interval, and boundary is the "midpoint"

    // Read in enough bits to resolve whether the boundary
    // is above or below "u,v"

    halfway = u + (v-u)/2 ;
    if( c == '1' ) {
      u = halfway ;
    } else if ( c=='0' ) {
      v = halfway ;
    } else {
      // carriage returns or white space, ignore
    } 

    // Read bits until we can decide what the source symbol was.
    // Then emulate the encoder's computations, and tie (u,v) to tag along for the ride.

    do {
      if(model_needs_updating) {
	w = b-a ;   cT = c0 + c1 ;   p0 = c0 / cT ;  
	boundary = a + (long) (p0 * (double)w)  ;
	// Special handling of unusual cases
	if ( b < a + 2*BMIN ) {
	  boundary = a + (b-a)/2 ; 
	} else {
	  if ( boundary < a + BMIN ) {	  boundary = a + BMIN ;	 }
	  if ( boundary > b - BMIN ) {	  boundary = b - BMIN ;	 }
	}
	// End Special handling of unusual cases
	model_needs_updating = 0 ;  
      }

      if  ( boundary <= u ) {
	printf("1\n");	c1 += 1.0 ;
	a = boundary ;
	model_needs_updating = 1 ; 	N-- ;
      } else if ( boundary >= v )  {
	printf("0\n");	c0 += 1.0 ;
	b = boundary ;
	model_needs_updating = 1 ; 	N-- ;
	// every time we discover a source bit, implement exactly the 
	// computations that were done by the encoder (below). 
      } else {
	// not enough bits have yet been read to know the decision.
      }

      // emulate outputting of bits by the encoder, and tie (u,v) to tag along for the ride.
      while ( (a>=HALF) || (b<=HALF) ) {
	if (a>=HALF) {
	  //	printf("1");
	  a = a-HALF ;
	  b = b-HALF ;
	  u = u-HALF ;
	  v = v-HALF ;
	} else {
	  //	printf("0");
	}
	a *= 2 ;      b *= 2 ;      u *= 2 ;      v *= 2 ;
	model_needs_updating = 1 ; 
      }
      if ( ( a < 0 ) || ( b> ONE ) || ( a >= b-2 ) ) {
	fprintf(stderr,"ERROR in a,b range\n") ; 
      }
      if ( ( u < 0 ) || ( v> ONE ) || ( u >=v-1 )  ) {
	fprintf(stderr,"ERROR in u,v range %ld %ld %g %g a=%ld boundary=%ld b=%ld\n" , u , v , c0 , c1 , a , boundary, b ) ; 
	//	exit(-1);
      }
    } while ( (N>0) && model_needs_updating ) ;
  }

  return(0);
}