// ARENC.cpp : Defines the entry point for the console application.
//
//This is DEMO version of range encoder written for research purposes by 
//Andrew Polar (www.ezcodesample.com, Dec. 24, 2006). The program reproduce
//algorithm of G.N.N. Martin, published in 1979. It was tested, however, 
//you can use it on your own risk without guaranteed successful encoding 
//and decoding for every possible data sample. Author does not accept any 
//responsibility for the use or misuse of this program.  
//The permission to copy and modify code is hereby granted as long as this 
//disclaimer stays. In case this code is chosen for production the user is 
//responsible for purchasing of necessary licenses to avoid possible patent 
//infringement. Last correction Jan.24, 2007.
//No warranty, however, any data sample that crash program or cause wrong
//result will be carefully investigated when sent to author. Find contact 
//information by WHOIS query for www.ezcodesample.com

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <time.h>

typedef struct {
	int* value;
	int size;
}INT_DATA;

double entropy24(int* data, int size) {
	
	int max_size = 1 << 24;
	int min, counter;
	int* buffer;
	double entropy;
	double log2 = log(2.0);
	double prob;

	min = data[0];
	for (counter=0; counter<size; counter++) {
		if (data[counter] < min)
			min = data[counter];
	}

	for (counter=0; counter<size; counter++) {
		data[counter] -= min;
	}

	buffer = (int*)malloc(max_size*sizeof(int));
	memset(buffer, 0x00, max_size*sizeof(int));
	for (counter=0; counter<size; counter++) {
		buffer[data[counter]]++;
	}

	entropy = 0.0;
	for (counter=0; counter<max_size; counter++) {
		if (buffer[counter] > 0) {
			prob = (double)buffer[counter];
			prob /= (double)size;
			entropy += log(prob)/log2*prob;
		}
	}
	entropy *= -1.0;

	for (counter=0; counter<size; counter++) {
		data[counter] += min;
	}

	if (buffer)
		free(buffer);

	return  entropy;
}

double round(double x) {
	double y  = floor(x);
	double yy = x - y;
	if (yy >= 0.5)
		return ceil(x);
	else
		return floor(x);
}

double MakeData(INT_DATA* int_data, int size, int MIN, int MAX, int redundancy_factor) {

	int counter, cnt, min, max;
	double buf;

	if (redundancy_factor <= 1)
		redundancy_factor = 1;

	if (MAX <= MIN)
		MAX = MIN + 1;

	if (int_data->size > 0) {
		if (int_data->value)
			free(int_data->value);
		int_data->size = 0;
	}

	int_data->size = size;
	int_data->value = (int*)malloc(int_data->size*sizeof(int));

	srand((unsigned)time(0)); 
	for (counter=0; counter<int_data->size; counter++) {
		buf = 0;
		for (cnt=0; cnt<redundancy_factor; cnt++) {
			buf += (double)rand();
		}
		int_data->value[counter] = ((int)buf)/redundancy_factor;
	}

	min = int_data->value[0];
	max = int_data->value[0];
	for (counter=0; counter<int_data->size; counter++) {
		if (int_data->value[counter] > max)
			max = int_data->value[counter];
		if (int_data->value[counter] < min)
			min = int_data->value[counter];
	}

	for (counter=0; counter<int_data->size; counter++) {
		buf = (double)(int_data->value[counter] - min);
		buf /= (double)(max - min);
		buf *= (double)(MAX - MIN);
		buf = round(buf);
		int_data->value[counter] = (int)buf + MIN;
	}

	return entropy24(int_data->value, int_data->size);
}

////////////////////////////////////Encoder////////////////////////////////

unsigned current_byte;
unsigned char bit_mask;
bool skip_zeros = true; 
unsigned char* g_bin_data = 0;

long long SafeProduct(long long x, int y, int mask, int& shift) {
	long long p = x * (long long)y;
	shift = 0;
	while (p > mask) {
		p >>= 1;
		shift++;
	}
	return p;
}

long long getBits(long long msg, int bits) {
	for (int i=0; i<bits; i++) {
		msg <<= 1;
		if ((g_bin_data[current_byte] & bit_mask) == bit_mask)
			msg |= 1;
		bit_mask >>= 1;
		if (bit_mask == 0) {
			current_byte++;
			bit_mask = 128;
		}
	}
	return msg;
}

void setBits(long long msg, int bits) {
	unsigned char shift = 63;
	for (int i=0; i<bits; i++) {
		if (((msg >> shift) & 1) > 0) {
			g_bin_data[current_byte] |= bit_mask;
			skip_zeros = false;
		}
		shift--;
		if (skip_zeros == false) {
			bit_mask >>= 1;
			if (bit_mask == 0) {
				current_byte++;
				bit_mask = 128;
			}
		}
	}
}

////////////////////////////////////////////////////////////////////////////

int main()
{
	INT_DATA int_data;
	int SIZE = 1000000;
	int MIN  = 1;
	int MAX  = 200;
	int redundancy_factor = 10;	

	double entropy = MakeData(&int_data, SIZE, MIN, MAX, redundancy_factor);
	int Bytes = (int)(ceil((double)(SIZE) * entropy)/8.0);
	printf("Entropy %f, estimated compressed size %d bytes\n", entropy, Bytes);

	g_bin_data = (unsigned char*)malloc(Bytes * 2);
	memset(g_bin_data, 0x00, Bytes * 2);

	//prepare data for encoding
	int max = int_data.value[0];
	int min = int_data.value[0];
	for (int i=0; i<int_data.size; i++) {
		if (max < int_data.value[i]) max = int_data.value[i];
		if (min > int_data.value[i]) min = int_data.value[i];
	}

	for (int i=0; i<int_data.size; i++)
		int_data.value[i] -= min;
	max -= min;

	int range = max + 1;
	int* f = (int*)malloc(range*sizeof(int));
	memset(f, 0x00, range * sizeof(int));
	for (int i=0; i<int_data.size; i++)
		f[int_data.value[i]]++;

	/////////////////////////////////////////////////////
	//correction of frequencies for stability
	int Denominator = SIZE;
	int base = 0;
	while (Denominator > 0) {
		Denominator >>= 1;
		base++;
	}
	Denominator = (1 << base);
	if (Denominator < SIZE) {
		Denominator <<= 1;
	}
	unsigned SIZE_MASK = Denominator - 1;
	if (SIZE_MASK > 0xFFFF) {
		SIZE_MASK = 0xFFFF;
		Denominator = SIZE_MASK + 1;
		base = 16;
	}

	double coeff = (double)(SIZE_MASK)/((double)(int_data.size));
	for (int i=0; i<range; i++) {
		if (f[i] > 0) {
			double ff = (double)(f[i]);
			ff *= coeff;
			f[i] = (int)(ff);
			if (f[i] == 0)
				f[i] = 1;
		}
	}
	int sm = 0;
	for (int i=0; i<range; i++) {
		sm += f[i];
	}
	int diff = SIZE_MASK - sm;
	if (diff > 0) {
		int mx = f[0];
		int ps = 0;
		for (int i=0; i<range; i++) {
			if (f[i] > mx) {
				mx = f[i];
				ps = i;
			}
		}
		f[ps] += diff;
	}
	if (diff < 0) {
		i = 0;
		while (diff < 0) {
			if (f[i] > 1) {
				f[i]--;
				diff++;
			}
			i++;
			if (i == (range-1))
				i = 0;
		}
	}
	//end of correction
	//////////////////////////////////////////////////////

	int* c = (int*)malloc((range+1)*sizeof(int));
	memset(c, 0x00, (range+1)*sizeof(int));
	c[0] = 0;
	for (i=1; i<range; i++) {
		c[i] = c[i-1] + f[i-1];
	}
	int denominator = 1 << base;
	c[range] = denominator;
	//data is ready

	//encoding
	printf("Encoding...\n");
	int Mask = 0xFFFFFFFF >> (32 - base);
	long long Message = 1;  //We start message from first positive bit
	//it must be removed from the stream in decoding.  The reason for 
	//that is that first cumulative frequency may have leading zeros
	//like follows 00010101. By setting first positive bit 100010101
	//we mark the beginning of codes.
	long long Freqs = 1;
	int shift = 0;
	current_byte = 0;
	bit_mask = 128;
	for (int i=0; i<int_data.size; i++) {
		setBits(Message, base-shift);
		Message <<= (base - shift);
		Message += ((long long)(c[int_data.value[i]]) * Freqs);
		Freqs = SafeProduct(Freqs, f[int_data.value[i]], Mask, shift);
	}
	setBits(Message, 64);
	printf("End encoding, size %d bytes\n\n", current_byte+1);
	//end encoding

	//prepare data for decoding
	int* symbol = (int*)malloc(denominator*sizeof(int));
	memset(symbol, 0x00, denominator*sizeof(int));
	for (int k=0; k<range; k++) {
		for (int i=c[k]; i<c[k+1]; i++) {
			symbol[i] = k;
		}
	}
	//data is ready

	//decoding
	printf("Decoding...\n");
	current_byte = 0;
	bit_mask = 128;
	long long ID;
	Freqs = 1;
	shift = 0;
	long long MSG = 0;
	int print_counter = 0;
	bool error_occur = false;
	//this is how we skip first positive bit
	getBits(MSG, 1);
	MSG = 0;
	////////////////////////////////////////
	for (int i=0; i<int_data.size; i++) {
		MSG = getBits(MSG, base-shift);
		ID = MSG/Freqs;
		if ((symbol[ID]) != (int_data.value[i])) {
			printf("Error decoding %d %d %d %d %d %d\n", i, 
				symbol[ID], int_data.value[i], 
				f[symbol[ID]], f[int_data.value[i]], shift);
			error_occur = true;
			break;
		}
		//else {
		//	print_counter++;
		//	if (print_counter == 8) {
		//		printf("%7d\n", symbol[ID] + min);
		//		print_counter = 0;
		//	}
		//	else {
		//		printf("%7d", symbol[ID] + min);
		//	}
		//}
		MSG -= Freqs*((long long)(c[symbol[ID]]));
		Freqs = SafeProduct(Freqs, f[symbol[ID]], Mask, shift);
	}
	if (error_occur == false) {
		printf("Decoding is correct, 100 percent match\n");
	}
	printf("\n");
	//end decoding

	if (f)
		free(f);
	if (c)
		free(c);
	if (symbol)
		free(symbol);
	if (g_bin_data)
		free(g_bin_data);
	if (int_data.value)
		free(int_data.value);

	return 0;
}