// OptimalArrayPermute.cpp : Defines the entry point for the console application.
// Copyright (C) April 30, 2011 -  George W. Taylor - All rights reserved
// See http://www.tropicalcoder.com/APermutationOnCombinatorialAlgorithms.htm

// This console app takes as the command line argument the number of digits in
// the desired set. The digits of the permutations will held in a byte array. 
// This version employs an table of indexes for the value 'k' so it will not
// need to be calculated in Narayana's Algorithm. The starting set and final
// permutation will be displayed.

// All permutations of the desired set will be generated using the algorithm discussed
// in "A Permutation on Combinatorial Algorithms" Each permutation can be verified
// against another set permuted independently via conventional means if desired.

// This was compiled as a 64 bit application, but should work as a 32 bit app.
// A 13 digit set (0,1,2,3,4,5,6,7,8,9,a,b,c) requires only 2.7 Gb
// of memory.

// On my Win 7 64 bit platform with this compiled as a 64 bit release
// version, the optimized table will generate the permutations of a 13
// digit set in 4.7 minutes, as compared to 6.3 minutes without the help
// of the meta-Algorithm. Of course this code is far from optimal, having been
// developed with the sole objective of illustrating the algorithm.

// See http://www.tropicalcoder.com/APermutationOnCombinatorialAlgorithms.htm


#include <Windows.h>
#include <stdio.h>
#include <conio.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <new>
using namespace std;

#define MAX_ELEMENTS 16
#define MIN_ELEMENTS 4
#define MAX_GIGABYTES 6
#define GIGABYTE (1024 * 1024 * 1024)

#define MAX_DIGITS 16
#define MIN_DIGITS 4
#define USE_TIMER
#define MAX_STACK (1024 * 40) // * 40 for 13 digits
// #define VALIDATE // uncomment this to enable validation
// #define DEBUG // uncomment this to enable validation
// #define REMOVE_FOR_TEST // uncomment this to speed test Narayana's algorithm unassisted

UINT64 *Factorials;
UINT64 *Offsets;
UINT64 kStackSize;
LPBYTE subsetStack;
LPBYTE kStack;
UINT64 countPermutations, maxPermutations;
UINT64 countGenerated;
UINT64 countSingletons;
int stackSize;

BYTE ValidationSet[MAX_ELEMENTS];
BYTE CurrentPermutation[MAX_ELEMENTS];

void ShowPermutation(LPBYTE Permutation, BYTE setSize);

#ifdef USE_TIMER
// #define NANOSECS // Timer delivers nanosecs with this defined - #undef this if you want time in microsecs
LARGE_INTEGER prevPerformanceCount, PerformanceCount, Frequency;


void TimerInit(void)
{ 
 prevPerformanceCount.LowPart = 0; 
 prevPerformanceCount.HighPart = 0; 
 PerformanceCount.LowPart = 0; 
 PerformanceCount.HighPart = 0;
 Frequency.LowPart = 0; 
 Frequency.HighPart = 0; 
 QueryPerformanceFrequency(&Frequency);
}

void TimerStart(void)
{
 QueryPerformanceCounter(&prevPerformanceCount);
}

void TimerEnd(void)
{
 QueryPerformanceCounter(&PerformanceCount);
}

UINT64 TimerGetInterval(void)
{
 UINT64 ticks;
 long double time;

 ticks = PerformanceCount.LowPart - prevPerformanceCount.LowPart;
#ifdef NANOSECS
 time =  (long double)ticks / (long double)Frequency.LowPart) * 1000000000;
#else
 time = ((long double)ticks / (long double)Frequency.LowPart) * 1000000;
#endif
 return (UINT64)floor(time);
}
#endif

// Get the factorial of NumDigits >= 3
UINT64 Factorial(unsigned NumDigits)
{
 UINT64 factorial = 2;
 for(unsigned n = 3; n <= NumDigits; n++)
 {
     factorial *= n;
 }
 return factorial;
}

void GetFactorials(unsigned setSize)
{
 for(unsigned i = 3; i <= setSize; i++)
 {
     Factorials[i] = Factorial(i);
 }
}

// Narayana's Algorithm
// Generates the next permutation in the series in place
// This is a referrence implimentation for providing
// independent validation.
void GeneratePermutation(LPBYTE Digits, BYTE setSize)
{
 // Find the largest index k such that a[k] < a[k + 1].
 short k, l, largestIndex = -1;
 for(k = (setSize - 2); k >= 0; k--)
 {
     if(Digits[k] < Digits[k+1])
     {
        if(k > largestIndex)largestIndex = k; 
     }
 }

 // If no such index exists, the permutation is the last permutation.
 if(largestIndex == -1)return;

 // Find the largest index l such that a[k] < a[l].
 // Since k + 1 is such an index, l is well defined and satisfies k < l.
 k = largestIndex;
 largestIndex = -1;
 for(l = (setSize - 1); l >= 0; l--) // <-- without the optimization here given the others
 {
     if(Digits[k] < Digits[l])
     {
        if(l > largestIndex)largestIndex = l; 
     }
 }

 // Swap a[k] with a[l].
 l = largestIndex;
 BYTE tmp = Digits[k];
 Digits[k] = Digits[l];
 Digits[l] = tmp;

 // Reverse the sequence from a[k + 1] up to and including the final element a[n].
 ++k; l = setSize - 1;
 while(l > k)
 {
       tmp = Digits[k];
       Digits[k] = Digits[l];
       Digits[l] = tmp;
       ++k; --l;
 }
}

// This version of Narayana's algorithm does not
// interate to find index 'k', the largest index
// such that a[k] < a[k + 1]. Instead, k is passed in,
// thereby eliminating the time that would have required.
void GeneratePermutation(short k, BYTE setSize)
{
 short l, largestIndex = -1;

 largestIndex = -1;
 // Note an optimization here. Was l >= 0, now l >= 1
 // Up to set of 14 digits at least, never gets below 1
 for(l = (setSize - 1); l >= 1; l--) 
 {
     if(CurrentPermutation[k] < CurrentPermutation[l])
     {
        if(l > largestIndex)largestIndex = l; 
     }
 }

 // Swap a[k] with a[l].
 l = largestIndex;
 BYTE tmp = CurrentPermutation[k];
 CurrentPermutation[k] = CurrentPermutation[l];
 CurrentPermutation[l] = tmp;

 // Reverse the sequence from a[k + 1] up to and including the final element a[n].
 ++k; l = setSize - 1;
 while(l > k)
 {
       tmp = CurrentPermutation[k];
       CurrentPermutation[k] = CurrentPermutation[l];
       CurrentPermutation[l] = tmp;
       ++k; --l;
 }
}


// This function will verify each permutation generated
// by checking it against an independently generated permutation.
// ValidationSet must initialized with the same set before we begin.
// This function must be called every single time a new permutation
// is created, to stay in sync.
bool ValidatePermutation(LPBYTE pCurrentPermutation, BYTE setSize)
{
 GeneratePermutation(ValidationSet, setSize);
 for(BYTE tmp = 0; tmp < setSize; tmp++)
 {
     if(pCurrentPermutation[tmp] != ValidationSet[tmp])
     {
        printf("Error: Permuation is invalid!\n");
        ShowPermutation(ValidationSet, setSize);
        return false;
     }
 }
 return true;
}


// A table is employed to set amd get the offset
// to the start of each subset.
UINT64 getOffset(BYTE digits, UINT64 indexTable)
{
 if(Offsets[digits] == 0)
 {
    Offsets[digits] = indexTable - 2;
    // printf("Created table for %d digits\n", digits);
    if(digits == 3)
    {
       printf("Added to table %d\n", 3);
    }else
    {
       printf("Added to table %I64d\n", Offsets[digits] - Offsets[digits - 1]);
    }
 }
 return Offsets[digits];
}

// We get k indexes reading backwards through the table
// to get the reflection of the upper half of each subset.
void getReflection(UINT64 count, UINT64 offset, BYTE setSize)
{
 for(UINT64 n = 0; n < count; n++)
 {
#ifdef DEBUG
     if((offset == 0) || (offset == 0xFFFFFFFFFFFFFFFF))
     {
        printf("Decremented offset too far!\n");
        break;
     }
#endif
	 GeneratePermutation(kStack[offset], setSize);
#ifdef VALIDATE
     if(!ValidatePermutation(CurrentPermutation, setSize))
     {
        ShowPermutation(CurrentPermutation, setSize);
        break;
     }	 
#endif
     ++countPermutations;
     --offset;
 }
}

// We get k indexes reading forwards through the table from
// the start offset to get the upper half of each subset.
void topHalf(UINT64 count, UINT64 offset, BYTE setSize)
{
 UINT64 n;

 for(n = 0; n < count; n++)
 {
#ifdef DEBUG
     if((offset == 0) || (offset == 0xFFFFFFFFFFFFFFFF))
     {
        printf("Decremented offset too far!\n");
        break;
     }
#endif
	 GeneratePermutation(kStack[offset], setSize);
#ifdef VALIDATE
     if(!ValidatePermutation(CurrentPermutation, setSize))
     {
        ShowPermutation(CurrentPermutation, setSize);
        break;
     }	 
#endif
     ++countPermutations;
     ++offset;
 }
}

// Calculate the offset for bottom of subset in the table
// and fetch the k indexes via topHalf()
void GetTopHalfFromTable(UINT64 count, BYTE digits, BYTE setSize)
{
 UINT64 offset;

 // We already have this in the table
 if(digits == 3)
 {
    // count = 3;
    offset = 1;
 }else
 {
#ifdef DEBUG
    if(count != (getOffset(digits, kStackSize) - getOffset(digits - 1, kStackSize)))
    {
       // We never enter here. This code is more for illustrative purposes,
       // showing another way of calculating num differences stored for each subset.
       printf("Discrepancy in count!\n"); 
    }
#endif

    offset = 2 + getOffset(digits, kStackSize) - count;
 }
 topHalf(count, offset, setSize);
}

// This version of Narayana's algorithm
// generates 'count' permutations in a loop.
// It also writes each k index found to the kStack
// for future referrence.
void GeneratePermutations(UINT64 count, BYTE setSize)
{
 for(UINT64 i = 0; i < count; i++)
 {
     // Find the largest index k such that a[k] < a[k + 1].
     short k, l, largestIndex = -1;
     for(k = (setSize - 2); k >= 0; k--)
     {
         if(CurrentPermutation[k] < CurrentPermutation[k+1])
         {
            if(k > largestIndex)largestIndex = k; 
         }
     }
     // If no such index exists, the permutation is the last permutation.
     if(largestIndex == -1)break;

     // Find the largest index l such that a[k] < a[l].
     k = largestIndex;
#ifndef REMOVE_FOR_TEST
	 // Save 'k' for later...
     kStack[kStackSize++] = (BYTE)k;
#endif
     largestIndex = -1;
     // Note the optimization here. Was l >= 0, now l >= 1
     // Up to set of 13 digits at least, never gets below 2
     for(l = (setSize - 1); l >= 1; l--)
     {
         if(CurrentPermutation[k] < CurrentPermutation[l])
         {
            if(l > largestIndex)largestIndex = l; 
         }
     }
     // Since k + 1 is such an index, l is well defined and satisfies k < l.

     // Swap a[k] with a[l].
     l = largestIndex;
     BYTE tmp = CurrentPermutation[k];
     CurrentPermutation[k] = CurrentPermutation[l];
     CurrentPermutation[l] = tmp;
     // Reverse the sequence from a[k + 1] up to and including the final element a[n].
     ++k; l = setSize - 1;
     while(l > k)
     {
           tmp = CurrentPermutation[k];
           CurrentPermutation[k] = CurrentPermutation[l];
           CurrentPermutation[l] = tmp;
           ++k; --l;
     }
#ifdef VALIDATE
     if(!ValidatePermutation(CurrentPermutation, setSize))
     {
        break; // <-- debug breakpoint
     }	 
#endif
     ++countGenerated;
 }
}


// Calculate how many permutations are in the main body
// of a set, excluding its subsets.
UINT64 getUpperBodyCount(BYTE digits)
{
 // How many permutations required for current subset?
 UINT64 count = Factorials[digits]/2;
 if(digits > 3)++count;

 // From this subtract what we already have from previous subsets...
 int numDigits = digits;
 while(--numDigits >= 3)
 {
       count -= Factorials[numDigits];
 }
 return count;
}

// Get the permutations for the central portion of
// the set, and add to the the reflection of the entire subset,
// who's makeup is determined from the stack.
void GetBodyCurrentSubset(BYTE digits, BYTE setSize)
{
 UINT64 offset, count = getUpperBodyCount(digits);


 if(Offsets[digits] != 0)
 {
    // We already have this in the table
    GetTopHalfFromTable(count, digits, setSize);
 }else
 {
    // Generate the permutations for upper half of main body
    // printf("begin top of body\n");
    GeneratePermutations(count, setSize);
 }

 // Each subset size is saved on the stack so that
 // we might know how to reconstruct it as a reflection.
 subsetStack[stackSize++] = digits;
 

 // We just completed the upper half of a set or subset.
 // We need the reflection of the entire thing,
 // for which count = Factorials[digits]/2-1.
 // In previous versions it would all be in the table,
 // but now we have to create it piece by piece
 // using the subset fragments stored in table.

 // Doing bottom half of body
 count -= 1;
 offset = getOffset(digits, kStackSize); // Was: offset = indexTable - 2;

 // We have the full table of differences
 // printf("begin remainder of body\n");
 getReflection(count, offset, setSize);
 // printf("Completed body\n");

 // Now do the remainder, the reflection of all the subsets that make up this set
 if(digits > 3)
 {
    int i;
    // Go back to the first instance of current set
    for(i = 0; i < stackSize; i++)
    {
        if(subsetStack[i] == digits)break;
    }
    // We begin with the set before that, 
    // since we just did the body of current set above
    while(--i >= 0) // Go on back through the subsets in reverse order
    {
          BYTE subset = subsetStack[i];

          if(digits < setSize)
          {
             // We add these to the stack as well,
             // such that the entire structure is elaborated
             // for future reference.
             subsetStack[stackSize++] = subset;
             // We must ensure that the stack is big enough.
             // Make stack bigger if necessary.
#ifdef DEBUG
             if(stackSize == MAX_STACK)return; 
#endif
          }

          if(subset == 1)
          {
             // Emit the reflected "singleton" between each subset > 3 
             // We use the version that generates a single permutation.
             // We do not save the k index on the kStack for these singletons.
             GeneratePermutation(CurrentPermutation, setSize);
#ifdef VALIDATE
             GeneratePermutation(ValidationSet, setSize);
#endif
             ++countGenerated;
          }else
          {
             // regenerate the body of this subset
             count = getUpperBodyCount(subset);
             GetTopHalfFromTable(count, subset, setSize);
             count -= 1;
             offset = getOffset(subset, kStackSize);
             getReflection(count, offset, setSize);
          }
    }
 }


 if(digits == setSize)
 {
    return; // Done!
 }

 // Did we need bigger arrays?
#ifdef DEBUG
 if(stackSize >= MAX_STACK)return;
#endif

 if(digits > 3)
 {
    // I call this a "singleton" - an unique
    // permutation that is emitted between subsets.
    // Each signals a swap of all digits from near or at
    // start of the set to the very end.
    GeneratePermutation(CurrentPermutation, setSize);
#ifdef VALIDATE
    GeneratePermutation(ValidationSet, setSize);
#endif
    subsetStack[stackSize++] = 1;
    ++countGenerated;
    ++countSingletons;
 }
}

// The meta-Algorithm
void recurseTree(BYTE digits, BYTE setSize)
{
 if(digits > 3)
 {
    recurseTree(digits-1, setSize);
 }


 GetBodyCurrentSubset(digits, setSize);


#ifdef DEBUG
 if(stackSize >= MAX_STACK)
 {
    printf("\nExceeded Stack size!\n");
    return;
 }
#endif

 if(digits > 3)
 {
#ifdef DEBUG
    if(countPermutations >= maxPermutations)return; 
    if(stackSize >= MAX_STACK)return;
#endif
    // Are we done yet?
    if(digits == setSize)
    {
       return; // Done!
    }

    recurseTree(digits-1, setSize);
 }
}

void ShowPermutation(LPBYTE Permutation, BYTE setSize)
{
 BYTE j;
 
 for(j = 0; j < setSize; j++)
 {
     if(Permutation[j] < 10)
     {
        printf("%d", Permutation[j]);
     }else
     {
        printf("%c", (char)((Permutation[j] - 10) + 'a'));
     }
 }
 printf("\n");
}

int main(int argc, char* argv[])
{
 UINT64 time, val, tableSize;
 BYTE tmp, setSize;
 bool bBadInput = false;

 Factorials = NULL;
 Offsets = NULL;
 kStack = NULL;
 subsetStack = NULL;

 if(argc < 2)
 {  
    bBadInput = true;
    goto ERR_OUT;
 }

 printf("Set size: %s\n", argv[1]);
 setSize = atoi(argv[1]);

 if((setSize < MIN_DIGITS) ||  (setSize > MAX_DIGITS))
 {  
    bBadInput = true;
    goto ERR_OUT;
 }

 Factorials = new UINT64 [setSize + 1];
 ZeroMemory(Factorials, sizeof(UINT64) * (setSize + 1));
 GetFactorials(setSize);

#ifndef REMOVE_FOR_TEST
 subsetStack = new BYTE [MAX_STACK];
 ZeroMemory(subsetStack, MAX_STACK);
 stackSize = 0;
 Offsets = new UINT64 [(setSize + 1)];
 ZeroMemory(Offsets, sizeof(UINT64) * (setSize + 1));
#endif

 // How big a table do we need?
 tableSize = 0;
 for(tmp = setSize; tmp >= 3; tmp--)
 {
     tableSize += getUpperBodyCount(tmp);
 }
 tableSize += 1;

#ifndef REMOVE_FOR_TEST

 // What is this in Gb?
 val = tableSize;
 if((val % GIGABYTE) != 0)
 {
     val /= GIGABYTE;
     ++val;
 }else val /= GIGABYTE;
 if(val)
 {
    printf("Max Gb of table memory required: %d\n", (int)val);
 }

 try
 {
     kStack = new BYTE [tableSize];
 }
 catch (bad_alloc&)
 {
     printf("Insufficient memory for table for this set.\n", setSize);
     goto ERR_OUT;
 }

 ZeroMemory(kStack, tableSize);
 kStackSize = 1;
#endif
 // Generate the set in ascending order
 for(tmp = 0; tmp < setSize; tmp++)
 {
     ValidationSet[tmp] = tmp;
 }

 memcpy(CurrentPermutation, ValidationSet, setSize);
 ShowPermutation(CurrentPermutation, setSize);
 printf("\n\n");

 maxPermutations = Factorials[setSize];
 countPermutations = 1;
 countGenerated = 0;
 countSingletons = 0;

 //**************************************
#ifdef USE_TIMER
 TimerInit();
 TimerStart();
#endif


// Compare with time required to generate permutations without
// the help of this algorithm by running GeneratePermutations()
#ifndef REMOVE_FOR_TEST
 recurseTree(setSize, setSize);
#else
 GeneratePermutations(maxPermutations, setSize);
#endif

#ifdef USE_TIMER
 TimerEnd();
 time = TimerGetInterval();
#endif

 // Show last permuatation
 ShowPermutation(CurrentPermutation, setSize);

 printf("\nCount permutations expected: %I64d\n", maxPermutations);
 printf("Count permutations created: %I64d\n", countPermutations);
 // Show count of permutations that had to be generated directly
 printf("Count permutations generated: %I64d\n", countGenerated); 
 printf("Count created + generated: %I64d\n", countPermutations + countGenerated); 
 printf("tableSize: %I64d\n", tableSize); 
 printf("kStackSize: %I64d\n", kStackSize);
 printf("stackSize: %i\n", stackSize);
 printf("Count singletons: %I64d\n", countSingletons);

#ifdef USE_TIMER
 printf("\nMinutes x 10: %I64d\n", time/6000000);
 printf("Seconds: %I64d\n", time/1000000);
 printf("Microseconds: %I64d\n", time);
#endif
 printf("\n\n");

ERR_OUT:
 if(bBadInput)
 {
    printf("Usage: Supply number digits in set to permute %d thru %d as command line arugument\n", MIN_DIGITS, MAX_DIGITS);
 }

 if(kStack != NULL)
    delete kStack;
 if(subsetStack != NULL)
    delete [] subsetStack;
 if(Factorials != NULL)
    delete Factorials;
 if(Offsets != NULL)
    delete Offsets;
 printf("\nAny key exits\n");
 while(!_kbhit());
 return 0;
}