//=================================================================================
// Decoding (rank/unrank) algorithms for the bounded-weight Grandaddy DB sequence,
// and its complement. When the there is an upper bound on the weight, the decoding
// algorithms can be applied to decode UCs for both k-subsets and t-multisets of an
// of an n-set using a difference representative. t-multisets can also be decoded
// using a shorthand frequency representation.

// The bounded weight Granddaddy ranking requires O(n^3 k^2) time; theunranking
// requires O(n^4 k^2 log k) time and use O(n^3 k) space.

// Programmed by Joe Sawada Feb, 2026
//=================================================================================
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAX 63      // n=62 is largest feasible for long long integers when k=2
#define MAXW 500    // max value of n*k
#define Max(a,b) (a > b) ? a : b
#define Min(a,b) (a < b) ? a : b

int mu[MAX] = { 0,1,-1,-1,0,-1,1,-1,0,0,1,-1,0,-1,1,1,0,-1,0,-1,0,
                1,1,-1,0,0,1,0,0,-1,-1,-1,0,1,1,1,0,-1,1,1,0,-1,
               -1,-1,0,0,1,-1,0,0,0,1,0,-1,0,1,0,1,1,-1,0,-1,1};

int phi[MAX] = { 0,1,1,2,2,4,2,6,4,6,4,10,4,12,6,8,8,16,6,18,8,12,
               10,22,8,20,12,18,12,28,8,30,16,20,16,24,12,36,18,
               24,16,40,12,42,20,24,22,46,16,42,20,32,24,52,18,
               40,24,36,28,58,16,60,30};

int n,k,w,first[MAX];
int suf[MAX][MAX];
long long int tableB[MAX][MAX][MAXW+1];
long long int tableP[MAX][MAX][MAXW+1];
long long int tableS[MAX][MAXW+1];

//===================================================================================
void Print(int a[], int j) {
    
    for (int i=1; i<=j; i++) printf("%d", a[i]); printf("  ");
}
//===================================================================================
// Returns the sum of the digits in a[1..j]
//===================================================================================
int Weight(int a[], int j) {
    int w = 0;
    
    for (int i=1; i <=j; i++)  w += a[i];
    return w;
}
//===================================================================================
// Returns: S_k(n,W) = the number of k-ary strings of length n with weight >= w
//===================================================================================
long long int S(int n, int w, int k) {
    long long int total = 0;
    
    if (n == 0 && w == 0) return 1;
    if ((n==0 && w>0) || w >k*n) return 0;
    if (w<=n) return pow(k,n);
    
    if (tableS[n][w] == -1) {
        for (int i = 1; i <= k; i++) total += S(n - 1, w - i,k);
        tableS[n][w] = total;
    }
    return tableS[n][w];
}
//===================================================================================
// Returns: B(t,j,w)
//===================================================================================
long long int B(int t, int j, int w, int k, int a[]) {
    long long int total=0;
    
    if (w < 0) w = 0;   // TO AVOID OUT OF BOUNDS
    if (t == 0 && j == 0 && w <= 0) return 1;
    if (t == j)  return 0;

    if (tableB[t][j][w] == -1){
        total = B(t,j+1,w,k,a);
        for (int x = a[j+1]+1; x <= k; x++) {
            total += B(t-j-1,0, w-x-Weight(a,j), k,a);
        }
        tableB[t][j][w] = total;
    }
    return tableB[t][j][w];
}
//===================================================================================
// Returns: P(t,j,w)
//===================================================================================
long long int P(int t, int j, int w, int k, int a[]) {
    long long int total=0;
    
    if (t==0 && j==0 && w==0) return 1;
    if (t==j || w < 0) return 0;

    if (tableP[t][j][w] == -1){
        total = P(t,j+1,w,k,a);
        for (int x=a[j+1]+1; x<=k; x++)  total += P(t-j-1,0, w-x-Weight(a,j),k,a);
        tableP[t][j][w] = total;
    }
    return tableP[t][j][w];
}
//===================================================================================
// Returns: A(t,j)
//===================================================================================
long long int A(int n, int t, int j, int w, int k, int a[]){
    long long int total = 0;
    int z,x,w2;
    
    if (n==j) return 0;
                                
    // Case 1
    if (t+j <= n){
        for (x=1; x<= a[j+1]-1; x++) {
            for (w2 = t-1; w2 <= k*(t-1); w2++) total += P(t-1,0,w2,k,a)*S(n-t-j,Max(0,w-w2-Weight(a,j)-x),k);
        }
    }
    // Case 2
    else {
        if (j < n-t+2) z = 0;
        else z = suf[n-t+2][j];
        if (a[z+1] < a[j+1]) {
            total = B(n-j+z,z+1,w-Weight(a,j-z),k,a);
            for (x=a[z+1]+1;x<=a[j+1]-1;x++) total += B(n-j-1,0,w-Weight(a,j)-x,k,a);
        }
    }
    return total;
}
//===================================================================================
// Find the longest prefix of a[1..n] that is a Lyndon word
//===================================================================================
int Lyn(int a[]) {
    int i,p=1;
    
    for (i=2; i<=n; i++) {
        if (a[i] < a[i-p]) return p;
        if (a[i] > a[i-p]) p = i;
    }
    return p;
}
//===================================================================================
// Return whether or not a[1..n] is a necklace
//===================================================================================
int IsNecklace(int a[]) {
    int i,p=1;

    for (i=2; i<=n; i++) {
        if (a[i] < a[i-p]) return 0;
        if (a[i] > a[i-p]) p=i;
    }
    if (n%p == 0) return 1;
    return 0;
}
//===================================================================================
// Compute the largest necklace neck[1..n] <= a[1..n]
//===================================================================================
void LargestNecklace(int n, int a[], int neck[]) {
    int i,p;
    
    for (i=1; i<=n; i++) neck[i] = a[i];
    while (!IsNecklace(neck)) {
        p = Lyn(neck);
        neck[p]--;
        for (i=p+1; i<=n; i++) neck[i] = k;
    }
}
//============================================================
// Return the number of strings whose necklace is <= str[1..n]
//============================================================
long long int T(int n, int str[]) {
    int i,j,l,p,t,neck[MAX];
    long long int tot=0;
    
    // Sets neck[1..n] to the largest necklace less than or equal to str[1..n]
    LargestNecklace(n, str, neck);
        
    // Initialize tables to be filled later via dynamic programming
    for (i=0;i<=n;i++){
        for (j=0;j<=n;j++){
            for (l=0;l<=k*n;l++) tableP[i][j][l] = tableB[i][j][l] = -1;
        }
    }
    
    // Compute suf[i][j] = longest suffix of neck[i..j] that is a prefix of neck[1..n]
    for (i=2; i<=n; i++) {
        p = i-1;
        for (j=i; j<=n; j++) {
            if (neck[j] > neck[j-p]) p = j;
            suf[i][j] = j-p;
        }
    }
    
    // Compute T(...)
    if (Weight(neck,n) >= w) tot = Lyn(neck);
    for (t=1; t<=n; t++) {
        for (j=0; j<n; j++) tot += A(n,t,j,w,k,neck);
    }
    
    return(tot);
}
//=========================================================================================
// Return the rank of str[1..n] in the bounded weight Granddaddy DB sequence (lex smallest)
//=========================================================================================
long long int RankDB(int str[]) {
    int i,t,s=0,neck[MAX],prev[MAX];
    
    //------------------------------------------------------------------------
    // WRAPAROUND a[1..n] = k^t x where t > 0 and x is a prefix of first[1..n]
    //------------------------------------------------------------------------
    t = 0;
    while (str[t+1] == k && t < n) t++;
    if (t == n && w == n*k) return 1;  // Special case
    if (t >= 1) {
        for (i=1; i<=n-t; i++) if (str[t+i] != first[i]) break;
        if (i > n-t) return(S(n,w,k)-t+1);
    }
    //-----------------------------------------------------------------------------------------------
    // str[1..n] is a necklace
    if (IsNecklace(str)) return( T(n,str) - Lyn(str) + 1 );  // MATCHES PAPER SINCE T(n,w) INCLUDES w
    
    // Find the necklace neck[1..n] of str[1..n] and its offset
    for (i=1; i<=n; i++) neck[i] = str[i];
    while (!IsNecklace(neck)) {
        s++;
        for (i=1; i<=n; i++) {
            if (i+s <= n) neck[i] = str[i+s];
            else neck[i] =str[i+s-n];
        }
    }
    if (s != t)  return (T(n,neck) - s + 1);  // HERE neck IS THE NECKLACE beta_1 IN THE PAPER
    if (Lyn(neck) < n)  return (T(n,neck) - Lyn(neck) - s + 1);  // MATCHES PAPER SINCE T(n,neck) INCLUDES neck
    
    for (i=n-s+1; i<=n; i++) neck[i] = 1;
    LargestNecklace(n,neck,prev);
    return (T(n,prev) - s + 1);  // HERE prev IS THE NECKLACE beta_1 IN THE PAPER
}
//------------------------------------------------------------------------------
// Find the smallest necklace neck[1..n] whose rank is >= to r via binary search
//------------------------------------------------------------------------------
void Smallest(int r, int neck[]) {
    int j,min,max,w2,t,last;
    
    for (j=1; j<=n; j++) neck[j] = k;
    w2 = n*k;
    
    // Determine character neck[j] from left to right using a binary search
    for (j=1; j<=n; j++) {
        min = 1;  max = k;
        while (min < max) {
            last = neck[j];
            t = (min+max)/2;
            neck[j] = t;
            w2 = w2 - last + t;
            if (w2 >= w && IsNecklace(neck) && RankDB(neck) >= r) max = t;
            else {
                neck[j] = last;
                w2 = w2 - t + last;
                min = t+1;
            }
        }
    }
}
//===================================================================================
// Return the string starting at index r
//===================================================================================
void UnRankDB(long long int r, int a[]) {
    int i,j,neck[MAX],prev[MAX];
    long long int d;
    
    // Special case for strings in the wraparound
    d =  (long long int) S(n,w,k) - r+1;
    if (d < n) {
        for (j=1; j<= d; j++)  a[j] = k;
        for (j=1; j<=n-d; j++) a[j+d] = first[j];
        return;
    }
        
    Smallest(r,neck);           // neck[1..n] = smallest necklace with rank >=r
    d = RankDB(neck);
    Smallest(Max(1,d-n), prev); // prev[1..n] = largest necklace with rank < r

    d = d - r;      // Take last d symbols from prev[1..n], and first n-d symbols from neck[1..n]
    j = 1;
    for (i=1; i<=d; i++)   a[j++] = prev[n-d+i];
    for (i=1; i<=n-d; i++) a[j++] = neck[i];
}
//===================================================================================
int main() {
    int i,j,count,option,k2,n2,v,a[MAX],s[MAX];
    long long int r;
 
    //-------------------- INPUT OPTION -----------------------------------------
    printf("1 = Rank Grandaddy \n");
    printf("2 = Rank UC with lower bound on weight \n");
    printf("3 = Rank UC with upper bound on weight  \n");
    printf("4 = Rank UC for k-subsets of n-set (difference reps) \n");
    printf("5 = Rank UC for k-multisets of n-set (difference reps) \n");
    printf("6 = Rank UC for k-multisets of n-set (frequency reps) \n\n");
    
    printf("7 = Unrank Grandaddy \n");
    printf("8 = Unrank UC with lower bound on weight \n");
    printf("9 = Unrank UC with upper bound on weight  \n");
    printf("10 = Unrank UC for k-subsets of n-set (difference reps)\n");
    printf("11 = Unrank UC for k-multisets of n-set (difference reps)\n");
    printf("12 = Unrank UC k-multisets of n-set (frequency reps) \n");

    printf("Select option #: ");  scanf("%d", &option);
    //---------------------------------------------------------------------------

    if (option >=1 && option <=12) {
        
        // Get n k [w]
        if (option == 2 || option == 3 || option == 8 || option == 9) { printf("Enter n k w: "); scanf("%d %d %d", &n, &k, &w); }
        else  { printf("Enter n k: "); scanf("%d %d", &n, &k);}
        
        // Map to bounded weight values n,k,w for subsets/multisets
        k2 = k; n2 = n;
        if (option == 4 || option == 10) { k = n2-k2+1; n = k2;     w = n2;  }
        else if (option == 5 || option == 11) { k = n2;      n = k2;     w = n2+k2-1;     }
        else if (option == 6 || option == 12) { k = k2+1;    n = n2-1;   w = n2+k2-1;   }
        else if (k <= 1) { printf("k must be greater than 1\n"); return 0; }
                
        // Map the upper bound to lower bound, and output the complement
        if ((option >= 3 && option <= 6) || (option >= 9 &&  option <= 12))  w = n*k - w + n;
        
        // Initialize S(n,w,k) to be filled later via dynamic programming
        for (i=0;i<=n;i++){
            for (j=0;j<=k*n;j++) tableS[i][j] = -1;
        }
        
        // Smallest necklace first[1..n] with weight w is of form 1^ixk^{n-i-1} with weight w, careful of when k=1
        j = 1;
        while (j + (n-j)*k >= w && k > 1) first[j++] = 1;
        first[j] = w - (n-j)*k - (j-1);
        for (i=j+1; i<=n; i++) first[i] = k;
                
        //--------------------------------------
        // RANKING THE BOUNDED WEIGHT GRANDDADDY
        //--------------------------------------
        if (option <=6) {
            
            if (option <=3) {
                for (i=1; i<=n; i++) {
                    printf("Enter s[%d]: ", i);
                    scanf("%d", &a[i]);
                    if (a[i] < 1 || a[i] > k) { printf("Invalid symbol, must be in the range [1..%d]\n", k); return 0; }
                }
            }
            if (option == 4) {
                printf("Enter subset in increasing order of elements over {1,2, ...., %d}\n", n2);
                for (i=1; i<=n; i++) {
                    printf("  Enter s[%d]: ", i);
                    scanf("%d", &s[i]);
                    if (s[i] < 1 || s[i] > n2) { printf("Invalid symbol, must be in the range [1..%d]\n", n2); return 0; }
                    if (i > 1 && s[i] <= s[i-1]) { printf("Invalid symbol, must be greater than the previous element\n"); return 0; }
                }
                // Convert to difference representation and complemented
                a[1] = k-s[1]+1;
                for (i=2; i<=k2; i++) a[i] = k-(s[i]-s[i-1])+1;
            }
            if (option == 5 || option == 6) {
                printf("Enter multiset in non-decreasing order of elements over {0,1, ..., %d}\n", n2-1);
                for (i=1; i<=k2; i++) {
                    printf("  Enter s[%d]: ", i);
                    scanf("%d", &s[i]);
                    s[i]++;
                    if (s[i] < 1 || s[i] > n2) { printf("Invalid symbol, must be in the range [0..%d]\n", n2-1); return 0; }
                    if (i > 1 && s[i] < s[i-1]) { printf("Invalid symbol, must be greater than or equal to the previous element\n"); return 0; }
                }
                if (option == 5) { // Convert to difference representation
                    a[1] = k-s[1]+1;
                    for (i=2; i<=k2; i++) a[i] = k-(s[i]-s[i-1]+1)+1;
                }
                if (option == 6) { // Convert to shorthand frequency representation
                    j = 1;
                    for (i=1; i<=n; i++) {
                        count = 1; // Adjusting to alphabet starting at 1 instead of 0
                        while (s[j] == i && j <= k2) { j++; count++;}
                        a[i] = k-count+1;
                    }
                }
            }

            printf("\nRank = %lld\n", RankDB(a));
        
        }
        //----------------------------------------
        // UNRANKING THE BOUNDED WEIGHT GRANDDADDY
        //----------------------------------------
        else {
            printf("Enter rank: ");   scanf("%lld", &r);
            if (r < 1 || r > S(n,w,k)) { printf("Invalid rank, must be in range [1..%lld]\n", S(n,w,k)); return 0;}

            UnRankDB(r,a);
            
            // CONVERT OUTPUT DEPENDING ON THE OPTION
            if (option >=9) for (i=1; i<=n; i++) a[i] = k - a[i] +1;    // Convert to complement
            if (option == 11 || option == 12) for (i=1; i<=n; i++) a[i]--; // Convert to alphabet {0, 1, ... , k-1}
            
            // OUTPUT
            printf("String at rank %lld: ", r);
            for (i=1; i<=n; i++) printf("%d", a[i]);
            
            // Add subset (difference rep)
            if (option == 10) {
                v = a[1];
                printf(" and the corresponding subset is { ");
                for (i=2; i<=n; i++) {
                    printf("%d, ", v);
                    v += a[i];
                }
                printf("%d }", v);
            }
            // Add multiset (difference rep)
            if (option == 11) {
                v = a[1]+1;
                printf(" and the corresponding multiset is { ");
                for (i=2; i<=n; i++) {
                    printf("%d, ", v-1);
                    v += a[i];
                }
                printf("%d }",v-1);
            }
            // Add multiset (shorthand frequency rep)
            if (option == 12) {
                count = 0;
                printf(" and the corresponding multiset is { ");
                for (i=1; i<=n; i++) {
                    for (j=1; j<=a[i]; j++) {
                        count++;
                        if (count < k2) printf("%d, ", i-1);
                        else printf("%d ", i-1);
                    }
                }
                while (count < k2) {
                    count++;
                    if (count < k2) printf("%d, ", n);
                    else printf("%d ", n);
                }
                printf("}");
            }
            
            printf("\n");
        }
    }
    else { printf("Option must be between 1 and 12\n"); return 0; }
}