#include <math.h>
#include <stdio.h>
#include <stdlib.h>



const int bits_per_item = sizeof(int) * 8;

// bit 0 is the least significant bit
int get_bit(int number, int bit)
{
  int result = number >> bit;
  result = result & 1;
  return result;
}


// Digit 0 is the least significant digit
int get_digit(int number, int bits_per_digit, int digit_position)
{
  int digits_per_int = sizeof(int)*8 / bits_per_digit;
  int left_shift = (digits_per_int - digit_position - 1) * bits_per_digit;
  int right_shift = (digits_per_int - 1) * bits_per_digit;

  unsigned int result = number << left_shift;
  result = result >> right_shift;
  return result;
}

void print_binary(int number)
{
  int bits = sizeof(int) * 8;

  int bit;
  for (bit = bits - 1; bit >= 0; bit--)
  {
    int binary = get_bit(number, bit);
    printf("%d", binary);
  }
}

void print_array(int * A, int length)
{
  int counter;
  for (counter = 0; counter < length; counter++)
  {
    printf ("%2d: %d\n", counter, A[counter]);
  }
  printf("\n");
}

void print_arrayb(int * A, int length)
{
  int counter;
  for (counter = 0; counter < length; counter++)
  {
    printf("%2d: %10d ", counter, A[counter]);
    print_binary(A[counter]);
    printf("\n");    
  }
  
  printf("\n");
}

void radixMSD_help2(int * items, int left, int right, int * scratch, 
                    int radix, int digit_position)
{
  if (left >= right) return;
  if (digit_position < 0) return;
  
  int bits_per_digit = (int) (log(radix) / log(2));
  int * counts = malloc(radix * sizeof(int));
  int * indices = malloc(radix * sizeof(int));

  int i;
  for (i = 0; i < radix; i++)
  {
    counts[i] = 0;
  }
  
  for (i = left; i <= right; i++)
  {
     int digit = get_digit(items[i], bits_per_digit, digit_position);
     counts[digit]++;
  }
  
  indices[0] = left;
  for (i = 1; i < radix; i++)
  {
    indices[i] = indices[i-1] + counts[i-1];
  }
    
  for (i = left; i <= right; i++)
  {
    int digit = get_digit(items[i], bits_per_digit, digit_position);
    int index = indices[digit];
    scratch[index] = items[i];
    indices[digit]++;
  }
  
  for (i = left; i <= right; i++) 
  {
    items[i] = scratch[i];
  }
  
  int temp_left = left;
  for (i = 0; i < radix; i++)
  {
    int temp_right = temp_left + counts[i] - 1;
    radixMSD_help2(items, temp_left, temp_right, scratch, radix, digit_position-1);
    temp_left = temp_right + 1;
  }
}


void radixMSD2(int * items, int length, int radix)
{
  int * scratch = malloc(sizeof(int) * length);
  
  int bits_per_digit = (int) (round((log(radix) / log(2))));
  int digits_per_item = sizeof(int)*8 / bits_per_digit;
  
  if (pow(2, bits_per_digit) != (float) radix)
  {
    printf("error: radix is not a power of 2\n");
    return;
  }

  if (sizeof(int)*8 % bits_per_digit != 0)
  {
    printf("error: bits_per_digit = %d, not a divisor of sizeof(int)*8 = %d:\n",
           bits_per_digit, sizeof(int)*8);
    return;
  }

  radixMSD_help2(items, 0, length -1, scratch, radix, digits_per_item - 1);
  free(scratch);
}


void radixMSD_help(int * items, int left, int right, int * scratch, int bit)
{
  if (left >= right) return;
  if (bit < 0) return;
  
  int counts[2];  
  int indices[2];

  int i;
  counts[0] = 0;
  counts[1] = 0;
  
  for (i = left; i <= right; i++)
  {
     int binary = get_bit(items[i], bit);
     counts[binary]++;
  }
  
  indices[0] = left;
  indices[1] = indices[0] + counts[0];
    
  for (i = left; i <= right; i++)
  {
    int binary = get_bit(items[i], bit);
    int index = indices[binary];
    scratch[index] = items[i];
    indices[binary]++;
  }
  
  for (i = left; i <= right; i++) 
  {
    items[i] = scratch[i];
  }
  
  int temp_left = left;
  for (i = 0; i < 2; i++)
  {
    int temp_right = temp_left + counts[i] - 1;
    radixMSD_help(items, temp_left, temp_right, scratch, bit-1);
    temp_left = temp_right + 1;
  }
}


void radixMSD(int * items, int length)
{
  int * scratch = malloc(sizeof(int) * length);
  
  radixMSD_help(items, 0, length -1, scratch, bits_per_item - 1);
  free(scratch);
}

void radixLSD_help(int * items, int length, int bit)
{
  int bits_per_item = sizeof(int) * 8;  
  if (bit >= bits_per_item) return;
  
  int * scratch = malloc(sizeof(int) * length);

  int counts[2];  
  int indices[2];

  int i;
  counts[0] = 0;
  counts[1] = 0;
  
  for (i = 0; i < length; i++)
  {
     int binary = get_bit(items[i], bit);
     counts[binary]++;
  }
  
  indices[0] = 0;
  indices[1] = indices[0] + counts[0];
    
  for (i = 0; i < length; i++)
  {
    int binary = get_bit(items[i], bit);
    int index = indices[binary];
    scratch[index] = items[i];
    indices[binary]++;
  }
  
  for (i = 0; i < length; i++) 
  {
    items[i] = scratch[i];
  }
  
  free(scratch);
}


void radixLSD(int * items, int length)
{
  int bits_per_item = sizeof(int) * 8;  

  int bit;
  for (bit = 0; bit < bits_per_item; bit++)
  {
    radixLSD_help(items, length, bit);
    printf("done with bit %d\n", bit);
    print_arrayb(items, length);
  }
}



int main()
{
  int A[] = {4, 93, 5, 104, 53, 90, 208};
//  int A[] = {8,7,6,5,4,3,2};
  printf("size of int: %d\n\n", sizeof(int));

  printf("before radix sort:\n");  
  print_array(A, 7);

//  radixMSD2(A, 7, 4);
//  radixMSD(A, 7);
  radixLSD(A, 7);
  
  printf("after radix sort:\n");
  print_array(A, 7);


  return 0;
}
