MagicIceVersion1_StreamCipher.c

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

typedef struct
{
    uint64_t State[16];
    uint64_t UpdatedState[16]; //FSM result
    uint64_t TransformedState[16];
    
    //Registers
    uint64_t R1, R2, R3, R4;
    uint64_t NLFSR, LFSR;
} MagicIceCipher;

uint32_t leftRotate32(uint32_t n, unsigned int d) {
   return (n << d)|(n >> (32 - d));
}

uint32_t rightRotate32(uint32_t n, unsigned int d) {
   return (n >> d)|(n << (32 - d));
}

uint64_t leftRotate64(uint64_t n, unsigned int d) {
   return (n << d)|(n >> (64 - d));
}

uint64_t rightRotate64(uint64_t n, unsigned int d) {
   return (n >> d)|(n << (64 - d));
}

//confusion layer
//Reference ISAAC algorithm mix function
void ComplexMix(uint32_t* a, uint32_t* b, uint32_t* c, uint32_t* d, uint32_t* e, uint32_t* f, uint32_t* g, uint32_t* h)
{
    //Shift Count = (0x243F6A88(pi) * 0x9E3779B9(phi)) mod 31
    //Shift Count 1 = (Shift Count * 0x9E3779B9(phi)) mod 31
    //...
    //{15, 9, 19, 16, 20, 25, 8, 10}

    *h = *h ^ (*g << 15); *f += *a; *h += *g;
    *g = *g ^ (*f << 9); *e += *h; *g += *f;
    *f = *f ^ (*e << 19); *d += *g; *f += *e;
    *e = *e ^ (*d << 16); *c += *f; *e += *d;
    *d = *d ^ (*c << 20); *b += *e; *d += *c;
    *c = *c ^ (*b << 25); *a += *d; *c += *b;
    *b = *b ^ (*a << 8); *h += *c; *b += *a;
    *a = *a ^ (*h << 10); *g += *b; *a += *h;
}

uint64_t LinearTransform0(uint64_t a, uint64_t b, uint64_t c, uint64_t d) {
    uint32_t a1 = a >> 32, a2 = a;
    uint32_t b1 = b >> 32, b2 = b;
    uint32_t c1 = c >> 32, c2 = c;
    uint32_t d1 = d >> 32, d2 = d;
    
    // 32-bit operations
    uint64_t x = ((uint64_t)(leftRotate32(a1, 13) ^ rightRotate32(b1, 7)) << 32) | (leftRotate32(a2, 13) ^ rightRotate32(b2, 7));
    uint64_t y = ((uint64_t)(leftRotate32(c1, 17) ^ rightRotate32(d1, 5)) << 32) | (leftRotate32(c2, 17) ^ rightRotate32(d2, 5));
    
    // 64-bit operations
    return leftRotate64(x, 12) ^ leftRotate64(y, 27);
}

uint64_t LinearTransform1(uint64_t a, uint64_t b, uint64_t c) {
    uint32_t a1 = a >> 32, a2 = a;
    uint32_t b1 = b >> 32, b2 = b;
    uint32_t c1 = c >> 32, c2 = c;
    
    // 32-bit operations
    uint64_t x = ((uint64_t)(rightRotate32(a1, 11) ^ leftRotate32(b1, 19)) << 32) | (rightRotate32(a2, 11) ^ leftRotate32(b2, 19));
    uint64_t y = ((uint64_t)(rightRotate32(c1, 23)) << 32) | rightRotate32(c2, 23);
    
    // 64-bit operations
    return rightRotate64(x, 49) ^ rightRotate64(y, 36);
}

uint64_t LinearTransform2(uint64_t a, uint64_t b) {
    uint32_t a1 = a >> 32, a2 = a;
    uint32_t b1 = b >> 32, b2 = b;
    
    // 32-bit operations
    uint64_t x = ((uint64_t)(leftRotate32(a1, 4)) << 32) | leftRotate32(a2, 4);
    uint64_t y = ((uint64_t)(rightRotate32(b1, 25)) << 32) | rightRotate32(b2, 25);
    
    // 64-bit operations
    return leftRotate64(x, 37) ^ leftRotate64(y, 43);
}

uint64_t LinearTransform3(uint64_t a) {
    uint32_t a1 = a >> 32, a2 = a;
    
    // 32-bit operations
    uint64_t x = ((uint64_t)(leftRotate32(a1, 37)) << 32) | leftRotate32(a2, 37);
    
    // 64-bit operations
    return rightRotate64(x, 9);
}

//https://users.ece.cmu.edu/~koopman/lfsr/index.html
uint64_t LFSR(uint64_t value) {
    uint64_t lsb = value & 1;  // Get least significant bit
    value >>= 1;                // Shift register
    if (lsb) {
        value ^= 0x80000000000019E2;  // Apply a reverse of polynomial if lsb was 1
    }
    return value;
}

//https://people.kth.se/~dubrova/nlfsr.html
//https://eprint.iacr.org/2012/166
uint64_t CombinedNLFSR(uint64_t value)
{
    //0 XOR 7 XOR 8 XOR (6 AND 8 AND 10)
    uint16_t nlfsr1 = value & 0xFFFF;
    uint16_t nlfsr1_feedback = ( nlfsr1 ^ (nlfsr1 >> 7) ^ (nlfsr1 >> 8) ^ ((nlfsr1 >> 6) & (nlfsr1 >> 8) & (nlfsr1 >> 10)) ) & 1;
    if(nlfsr1_feedback)
        nlfsr1 = (nlfsr1 >> 1) | (nlfsr1_feedback << 15);
    
    //0 XOR 7 XOR 8 XOR (5 AND 12 AND 14)
    uint16_t nlfsr2 = (value >> 16) & 0xFFFF;
    uint16_t nlfsr2_feedback = ( nlfsr2 ^ (nlfsr2 >> 7) ^ (nlfsr2 >> 8) ^ ((nlfsr2 >> 5) & (nlfsr2 >> 12) & (nlfsr2 >> 14)) ) & 1;
    if(nlfsr2_feedback)
        nlfsr2 = (nlfsr2 >> 1) | (nlfsr2_feedback << 15);

    //0 XOR 1 XOR 8 XOR (7 AND 8 AND 11)
    uint16_t nlfsr3 = (value >> 32) & 0xFFFF;
    uint16_t nlfsr3_feedback = ( nlfsr3 ^ (nlfsr3 >> 1) ^ (nlfsr3 >> 8) ^ ((nlfsr3 >> 7) & (nlfsr3 >> 8) & (nlfsr3 >> 11)) ) & 1;
    if(nlfsr3_feedback)
        nlfsr3 = (nlfsr3 >> 1) | (nlfsr3_feedback << 15);
    
    //0 XOR 4 XOR 8 XOR 9 XOR 10 (8 AND 12)
    uint16_t nlfsr4 = (value >> 48) & 0xFFFF;
    uint16_t nlfsr4_feedback = ( nlfsr4 ^ (nlfsr4 >> 4) ^ (nlfsr4 >> 8) ^ (nlfsr4 >> 9) ^ (nlfsr4 >> 10) ^ ((nlfsr4 >> 8) & (nlfsr4 >> 12)) ) & 1;
    if(nlfsr4_feedback)
        nlfsr4 = (nlfsr4 >> 1) | (nlfsr4_feedback << 15);

    return ((uint64_t)nlfsr4 << 48) | ((uint64_t)nlfsr3 << 32) | ((uint64_t)nlfsr2 << 16) | nlfsr1;
}

void InitializeState(MagicIceCipher* instance, const uint32_t* KeyAndIV, uint64_t KeyAndIV_Size)
{
    if(KeyAndIV == NULL)
        return;

    uint32_t random_array[8] = {0};
    for(size_t i = 0; i < 8; i++)
        random_array[i] = (uint32_t)0xB7E151628; //e

    //Scramble it
    for ( size_t i = 0; i < 16; i++ )
		ComplexMix(&random_array[0],&random_array[1],&random_array[2],&random_array[3],&random_array[4],&random_array[5],&random_array[6],&random_array[7]);

    if(KeyAndIV_Size > 0 && KeyAndIV_Size <= 32)
    {
        //Mix with Key And IV
        uint32_t TemporaryState[32] = {0};
        for(size_t i = 0; i < 32; i++)
        {
            TemporaryState[i % 16] = random_array[i % 8] ^ KeyAndIV[i];
            if(i % 8 == 0)
                ComplexMix(&random_array[0],&random_array[1],&random_array[2],&random_array[3],&random_array[4],&random_array[5],&random_array[6],&random_array[7]);
        }

        //Scrambled array mix to State
        for(size_t i = 0; i < 16; i++)
            instance->State[i] = ((uint64_t)TemporaryState[i] << 32) | TemporaryState[31 - i];
    }
    else
    {
        for(size_t i = 0; i < (16 * 8); i++)
        {
            instance->State[i % 16] ^= ((uint64_t)random_array[i % 8] << 32) | random_array[(i + 1) % 8];
            if(i % 8 == 0)
                ComplexMix(&random_array[0],&random_array[1],&random_array[2],&random_array[3],&random_array[4],&random_array[5],&random_array[6],&random_array[7]);
        }
    }
}

uint64_t ghash_multiply(uint64_t a, uint64_t b) {
    if(a == 0 || b == 0)
        return 0;

    uint64_t product = 0;
    const uint64_t POLYNOMIAL = 0xE100000000000000; // This is x^64 + x^4 + x^3 + x + 1 in binary

    for(int i = 0; i < 64; i++) {
        if(b & 1) {
            product ^= a;
        }

        //boolean type
        int high_bit_set = a & 0x8000000000000000;

        a <<= 1;

        if(high_bit_set) {
            a ^= POLYNOMIAL;
        }

        b >>= 1;
    }

    return product;
}

void TransformState(MagicIceCipher* instance)
{
    const uint64_t MASK1 = 0xAAAAAAAAAAAAAAAA;
    const uint64_t MASK2 = 0x5555555555555555;
    
    // Shift rows: Rotate cells
    // Row 1 need left rotate 0 cells, Row 2 need left rotate 1 cell, Row 3 need left rotate 2 cells, Row 3 need left rotate 3 cells.
    memcpy(instance->TransformedState, instance->State, 16 * sizeof(uint64_t));
    uint64_t S;
    for (int i = 1; i < 4; i++) {
        int Counter = 0;
        while (Counter < i) {
            S = instance->TransformedState[4 * i];
            
            for (int k = 1; k < 4; k++) {
                instance->TransformedState[4 * i + k - 1] = instance->TransformedState[4 * i + k];
            }

            instance->TransformedState[4 * i + 4 - 1] = S;
            Counter++;
        }
    }

    uint64_t A = 0, B = 0;

    // Mix columns: Finite field multiplication with GHASH part
    // For simplicity, here we implement a kind of mixing operation inspired by AES
    for(size_t round = 0; round < 8; round++) {
        A = instance->State[round] ^ (instance->State[round] & MASK1) + (instance->State[15 - round] & MASK2);
        B = ((~instance->TransformedState[(round + 1) % 16]) >> 32 | instance->TransformedState[(round + 2) % 16] << 32) & instance->TransformedState[(round + 3) % 16];

        // Update transformed state with computed A and B
        instance->TransformedState[round] ^= ghash_multiply(A, B);
    }
    
    instance->NLFSR = CombinedNLFSR(A ^ B);
    instance->LFSR = LFSR(A + B);
}

void ShiftRightBy(uint64_t* array, uint32_t number)
{
    if(number == 0 || number >= 16)
        return;

    uint64_t temp;
    for (int i = 0; i < number; i++) {
        temp = array[15];
        for (int j = 15; j > 0; j--) {
            array[j] = array[j - 1];
        }
        array[0] = temp;
    }
}

//enhanced confusion layer
void FiniteStateAutomaton(MagicIceCipher* instance)
{
    uint32_t random_array[8] = {0};

    //Compute UpdatedState
    for (int i = 0; i < 16; i++)
    {
        uint64_t T1 = instance->State[16 - i - 1] & 0xFFFFFFFF00000000 | (instance->State[16 - i] >> 32);
        uint64_t T2 = ((instance->State[16 - i] & 0x00000000FFFFFFFF) << 32) | instance->State[16 - i - 1] & 0x00000000FFFFFFFF;
        uint64_t T3 = (instance->State[0] + instance->State[1]);
        uint64_t T4 = (instance->State[0] ^ instance->State[1]);

        // Extract high and low parts from T1, T2, T3, T4 and R1, R2, R3, R4
        uint32_t T1_H = (uint32_t)(T1 >> 32);
        uint32_t T1_L = (uint32_t)(T1 & 0xFFFFFFFF);
        uint32_t T2_H = (uint32_t)(T2 >> 32);
        uint32_t T2_L = (uint32_t)(T2 & 0xFFFFFFFF);
        uint32_t T3_H = (uint32_t)(T3 >> 32);
        uint32_t T3_L = (uint32_t)(T3 & 0xFFFFFFFF);
        uint32_t T4_H = (uint32_t)(T4 >> 32);
        uint32_t T4_L = (uint32_t)(T4 & 0xFFFFFFFF);
        uint32_t R1_H = (uint32_t)(instance->R1 >> 32);
        uint32_t R1_L = (uint32_t)(instance->R1 & 0xFFFFFFFF);
        uint32_t R2_H = (uint32_t)(instance->R2 >> 32);
        uint32_t R2_L = (uint32_t)(instance->R2 & 0xFFFFFFFF);
        uint32_t R3_H = (uint32_t)(instance->R3 >> 32);
        uint32_t R3_L = (uint32_t)(instance->R3 & 0xFFFFFFFF);
        uint32_t R4_H = (uint32_t)(instance->R4 >> 32);
        uint32_t R4_L = (uint32_t)(instance->R4 & 0xFFFFFFFF);

        // Call ComplexMix with the required arguments
        ComplexMix(&T1_H, &R1_L, &T2_H, &R2_L, &T3_H, &R3_L, &T4_H, &R4_L);
        ComplexMix(&T3_L, &R3_H, &T4_L, &R4_H, &T1_L, &R1_H, &T2_L, &R2_H);

        // Combine high and low parts back into T1, T2, T3, T4 and R1, R2, R3, R4
        T1 = ((uint64_t)T1_H << 32) | T1_L;
        T2 = ((uint64_t)T2_H << 32) | T2_L;
        T3 = ((uint64_t)T3_H << 32) | T3_L;
        T4 = ((uint64_t)T4_H << 32) | T4_L;
        instance->R1 = ((uint64_t)R1_H << 32) | R1_L;
        instance->R2 = ((uint64_t)R2_H << 32) | R2_L;
        instance->R3 = ((uint64_t)R3_H << 32) | R3_L;
        instance->R4 = ((uint64_t)R4_H << 32) | R4_L;

        //Assgin UpdatedStaete
        instance->UpdatedState[0] = leftRotate64( T1 + (instance->R1 ^ instance->R2), 32);
        instance->UpdatedState[1] = rightRotate64( (T2 ^ instance->R3) + instance->R4, 32);
        instance->UpdatedState[(uint32_t)(15 - i - 1) % 16] ^= (T2 ^ instance->R3) + instance->R4;
        instance->UpdatedState[15 - i] ^= T1 + (instance->R1 ^ instance->R2);
        
        //Update R1, R2, R3, R4
        instance->NLFSR = CombinedNLFSR(instance->NLFSR);
        instance->LFSR = LFSR(instance->LFSR);
        instance->R3 ^= instance->NLFSR;
        instance->R4 ^= instance->LFSR;
        instance->R1 = instance->R2 + (instance->R3 ^ instance->R4);
        instance->R2 = instance->R1 ^ (instance->R4 + instance->R3);

        ShiftRightBy(instance->State, 4);
        ShiftRightBy(instance->UpdatedState, 2);
    }
}

void UpdateState(MagicIceCipher* instance)
{
    TransformState(instance);
    FiniteStateAutomaton(instance);
}

void GenerateKeyStream(MagicIceCipher* instance, uint64_t* KeyStream, uint64_t KeyStreamSize)
{
    uint64_t LinearTransformedState[16] = {0};
    uint64_t UnpredictableState[16] = {0};

    for(int i=0; i<16; i+=4)
    {
        LinearTransformedState[i] = LinearTransform0(instance->State[i], instance->State[(i + 1) % 16], instance->State[(i + 2) % 16], instance->State[(i + 3) % 16]);
        LinearTransformedState[i + 1] = LinearTransform1(instance->State[(i + 5) % 16], instance->State[(i + 6) % 16], instance->State[(i + 7) % 16]);
        LinearTransformedState[i + 2] = LinearTransform2(instance->State[(i + 8) % 16], instance->State[(i + 9) % 16]);
        LinearTransformedState[i + 3] = LinearTransform3(instance->State[(i + 10) % 16]);
    }

    for(int i=0; i<16; i++)
        UnpredictableState[i] = instance->State[i] + (instance->TransformedState[i] ^ instance->UpdatedState[i]);

    size_t ElementSize = KeyStreamSize < 32 ? KeyStreamSize : 32;
    memcpy(KeyStream, UnpredictableState, ElementSize * sizeof(uint64_t));

    UpdateState(instance);

    for(int i=0; i<16; i++)
        instance->State[i] ^= LinearTransformedState[i];
}

void ResetState(MagicIceCipher* instance, const uint32_t* KeyAndIV, uint64_t KeyAndIV_Size)
{
    //Wipe State
    for(size_t i = 0; i < 16; i++)
        instance->State[i] = instance->UpdatedState[i] = instance->TransformedState[i] = 0;
    
    instance->R1 = instance->R2 = instance->R3 = instance->R4 = 0;
    instance->NLFSR = instance->LFSR = 0;

    InitializeState(instance, KeyAndIV, KeyAndIV_Size);
    UpdateState(instance);
}

int main() {
    MagicIceCipher cipher;
    
    // Initialize the KeyAndIV (Array size 32) in the cipher
    uint32_t KeyAndIV[32] = {0};
    for (int i = 0; i < 32; i++) {
        KeyAndIV[i] = i;
    }

    uint64_t Keystream[16] = {0};
    
    ResetState(&cipher, KeyAndIV, 32);
    UpdateState(&cipher);
    GenerateKeyStream(&cipher, Keystream, 16);
    
    // Print out the State and TransformedState, Generated Keystream
    printf("State:\n");
    for (int i = 0; i < 16; i++) {
        printf("%llu\n", cipher.State[i]);
    }
    
    printf("\nTransformedState:\n");
    for (int i = 0; i < 16; i++) {
        printf("%llu\n", cipher.TransformedState[i]);
    }
    
    printf("\nGenerated Keystream:\n");
    for (int i = 0; i < 16; i++) {
        printf("%llu\n", Keystream[i]);
    }

    // Prepare the file for writing
    FILE *file = fopen("MagicIce_Keystream.bin", "wb");
    if (file == NULL) {
        printf("Could not open file for writing.\n");
        return 1;
    }

    ResetState(&cipher, KeyAndIV, 32);
    // Generate and write 128 64-bit values to the file for a total of 128KB
    for (int i = 0; i < 1024; i++) {
        uint64_t Keystream[16] = {0};
       
        GenerateKeyStream(&cipher, Keystream, 16);
        
        fwrite(Keystream, sizeof(uint64_t), 16, file);
    }
    ResetState(&cipher, KeyAndIV, 32);
    
    return 0;
}

MagicIceVersion2_StreamCipher.c

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

typedef struct
{
    uint64_t State[16];
    uint64_t UpdatedState[16]; //FSM result
    uint64_t TransformedState[16];
    
    //Registers
    uint64_t R1, R2, R3, R4;
    uint64_t NLFSR, LFSR;

} MagicIceCipher;

uint32_t leftRotate32(uint32_t n, unsigned int d) {
   return (n << d)|(n >> (32 - d));
}

uint32_t rightRotate32(uint32_t n, unsigned int d) {
   return (n >> d)|(n << (32 - d));
}

uint64_t leftRotate64(uint64_t n, unsigned int d) {
   return (n << d)|(n >> (64 - d));
}

uint64_t rightRotate64(uint64_t n, unsigned int d) {
   return (n >> d)|(n << (64 - d));
}

//confusion layer
//Reference ISAAC algorithm mix function
void ComplexMix(uint32_t* a, uint32_t* b, uint32_t* c, uint32_t* d, uint32_t* e, uint32_t* f, uint32_t* g, uint32_t* h)
{
    //Shift Count = (0x243F6A88(pi) * 0x9E3779B9(phi)) mod 31
    //Shift Count 1 = (Shift Count * 0x9E3779B9(phi)) mod 31
    //...
    //{15, 9, 19, 16, 20, 25, 8, 10}

    *h = *h ^ (*g << 15); *f += *a; *h += *g;
    *g = *g ^ (*f << 9); *e += *h; *g += *f;
    *f = *f ^ (*e << 19); *d += *g; *f += *e;
    *e = *e ^ (*d << 16); *c += *f; *e += *d;
    *d = *d ^ (*c << 20); *b += *e; *d += *c;
    *c = *c ^ (*b << 25); *a += *d; *c += *b;
    *b = *b ^ (*a << 8); *h += *c; *b += *a;
    *a = *a ^ (*h << 10); *g += *b; *a += *h;
}

// Bit-by-bit conversion of 4 of an 8-bit number (Reference Serpent Algorithm)
void AvalancheTransformation(uint8_t* a, uint8_t* b, uint8_t* c, uint8_t* d)
{
    uint8_t tv;
    *d ^= *a;
    tv = *b;
    *b &= *d;
    tv ^= *c;
    *b ^= *a;
    *a |= *d;
    *a ^= tv;
    tv ^= *d;
    *d ^= *c;
    *c |= *b;
    *c ^= tv;
    tv = ~tv;
    tv |= *b;
    *b ^= *d;
    *b ^= tv;
    *d |= *a;
    *b ^= *d;
    tv ^= *d;
    *d = *a;
    *a = *b;
    *b = tv;
    tv = 0;
}

// Break down a 32-bit number into four 8-bit numbers, process each 8-bit number,
// then combine them back into a 32-bit number
void ProcessWithAvalanche(uint32_t* number) {
    uint8_t numbers[4];
    numbers[0] = *number & 0xff;
    numbers[1] = (*number >> 8) & 0xff;
    numbers[2] = (*number >> 16) & 0xff;
    numbers[3] = (*number >> 24) & 0xff;
    
    // Apply the transformation to each 8-bit number
    AvalancheTransformation(&numbers[0], &numbers[1], &numbers[2], &numbers[4]);
    AvalancheTransformation(&numbers[1], &numbers[2], &numbers[3], &numbers[0]);
    AvalancheTransformation(&numbers[2], &numbers[3], &numbers[0], &numbers[1]);
    AvalancheTransformation(&numbers[3], &numbers[0], &numbers[1], &numbers[2]);
    
    // Combine the 8-bit numbers back into a 32-bit number
    *number = numbers[0] | (numbers[1] << 8) | (numbers[2] << 16) | (numbers[3] << 24);
}

//https://users.ece.cmu.edu/~koopman/lfsr/index.html
uint64_t LFSR(uint64_t value) {
    uint64_t lsb = value & 1;  // Get least significant bit
    value >>= 1;                // Shift register
    if (lsb) {
        value ^= 0x80000000000019E2;  // Apply a reverse of polynomial if lsb was 1
    }
    return value;
}

//https://people.kth.se/~dubrova/nlfsr.html
//https://eprint.iacr.org/2012/166
uint64_t CombinedNLFSR(uint64_t value)
{
    //0 XOR 7 XOR 8 XOR (6 AND 8 AND 10)
    uint16_t nlfsr1 = value & 0xFFFF;
    uint16_t nlfsr1_feedback = ( nlfsr1 ^ (nlfsr1 >> 7) ^ (nlfsr1 >> 8) ^ ((nlfsr1 >> 6) & (nlfsr1 >> 8) & (nlfsr1 >> 10)) ) & 1;
    if(nlfsr1_feedback)
        nlfsr1 = (nlfsr1 >> 1) | (nlfsr1_feedback << 15);
    
    //0 XOR 7 XOR 8 XOR (5 AND 12 AND 14)
    uint16_t nlfsr2 = (value >> 16) & 0xFFFF;
    uint16_t nlfsr2_feedback = ( nlfsr2 ^ (nlfsr2 >> 7) ^ (nlfsr2 >> 8) ^ ((nlfsr2 >> 5) & (nlfsr2 >> 12) & (nlfsr2 >> 14)) ) & 1;
    if(nlfsr2_feedback)
        nlfsr2 = (nlfsr2 >> 1) | (nlfsr2_feedback << 15);

    //0 XOR 1 XOR 8 XOR (7 AND 8 AND 11)
    uint16_t nlfsr3 = (value >> 32) & 0xFFFF;
    uint16_t nlfsr3_feedback = ( nlfsr3 ^ (nlfsr3 >> 1) ^ (nlfsr3 >> 8) ^ ((nlfsr3 >> 7) & (nlfsr3 >> 8) & (nlfsr3 >> 11)) ) & 1;
    if(nlfsr3_feedback)
        nlfsr3 = (nlfsr3 >> 1) | (nlfsr3_feedback << 15);
    
    //0 XOR 4 XOR 8 XOR 9 XOR 10 (8 AND 12)
    uint16_t nlfsr4 = (value >> 48) & 0xFFFF;
    uint16_t nlfsr4_feedback = ( nlfsr4 ^ (nlfsr4 >> 4) ^ (nlfsr4 >> 8) ^ (nlfsr4 >> 9) ^ (nlfsr4 >> 10) ^ ((nlfsr4 >> 8) & (nlfsr4 >> 12)) ) & 1;
    if(nlfsr4_feedback)
        nlfsr4 = (nlfsr4 >> 1) | (nlfsr4_feedback << 15);

    return ((uint64_t)nlfsr4 << 48) | ((uint64_t)nlfsr3 << 32) | ((uint64_t)nlfsr2 << 16) | nlfsr1;
}

void InitializeState(MagicIceCipher* instance, const uint32_t* KeyAndIV, uint64_t KeyAndIV_Size)
{
    if(KeyAndIV == NULL)
        return;

    uint32_t random_array[8] = {0};
    for(size_t i = 0; i < 8; i++)
        random_array[i] = (uint32_t)0xB7E151628; //e

    //Scramble it
    for ( size_t i = 0; i < 16; i++ )
	{
        ComplexMix(&random_array[0],&random_array[1],&random_array[2],&random_array[3],&random_array[4],&random_array[5],&random_array[6],&random_array[7]);
    }

    if(KeyAndIV_Size > 0 && KeyAndIV_Size <= 32)
    {
        uint32_t TemporaryState[32] = {0};

        memcpy(TemporaryState, KeyAndIV, KeyAndIV_Size * sizeof(uint32_t));

        //Mix with Key And IV
        for(size_t i = 0; i < 32; i += 8)
        {
            for (size_t j = 0; j < 8; ++j)
                random_array[j] += TemporaryState[i + j];
            
            ComplexMix(&random_array[0],&random_array[1],&random_array[2],&random_array[3],&random_array[4],&random_array[5],&random_array[6],&random_array[7]);

            //Scrambled array mix to State
            for(size_t j = 0; j < 8; ++j)
                instance->State[(i + j) % 16] += ((uint64_t)random_array[j] << 32) | random_array[7 - j];
        }
    }
    else
    {
        for(size_t counter = 0; counter < (16 * 8); counter++)
        {
            for(size_t i = 0; i < 16; i += 8)
            {
                for (size_t j = 0; j < 8; ++j)
                    instance->State[i + j] += ((uint64_t)random_array[j] << 32) | random_array[7 - j];

                ComplexMix(&random_array[0],&random_array[1],&random_array[2],&random_array[3],&random_array[4],&random_array[5],&random_array[6],&random_array[7]);
            }
        }
    }
}

uint64_t ghash_multiply(uint64_t a, uint64_t b) {
    if(a == 0 || b == 0)
        return 0;

    uint64_t product = 0;
    const uint64_t POLYNOMIAL = 0xE100000000000000; // This is x^64 + x^4 + x^3 + x + 1 in binary

    for(int i = 0; i < 64; i++) {
        if(b & 1) {
            product ^= a;
        }

        //boolean type
        int high_bit_set = a & 0x8000000000000000;

        a <<= 1;

        if(high_bit_set) {
            a ^= POLYNOMIAL;
        }

        b >>= 1;
    }

    return product;
}

void TransformState(MagicIceCipher* instance)
{
    const uint64_t MASK1 = 0xAAAAAAAAAAAAAAAA;
    const uint64_t MASK2 = 0x5555555555555555;
    
    // Shift rows: Rotate cells
    // Row 1 need left rotate 0 cells, Row 2 need left rotate 1 cell, Row 3 need left rotate 2 cells, Row 3 need left rotate 3 cells.
    memcpy(instance->TransformedState, instance->State, 16 * sizeof(uint64_t));
    uint64_t S;
    for (uint32_t i = 1; i < 4; i++) {
        uint32_t counter = 0;
        while (counter < i) {
            S = instance->TransformedState[4 * i];
            
            for (uint32_t k = 1; k < 4; k++) {
                instance->TransformedState[4 * i + k - 1] = instance->TransformedState[4 * i + k];
            }

            instance->TransformedState[4 * i + 4 - 1] = S;
            counter++;
        }
    }

    uint64_t A = 0, B = 0;

    // Mix columns: Finite field multiplication with GHASH part
    // For simplicity, here we implement a kind of mixing operation inspired by AES
    for(size_t round = 0; round < 16; round++) {
        A ^= instance->State[round] ^ (instance->State[round] & MASK1) + (instance->State[15 - round] & MASK2);
        B ^= ((~instance->TransformedState[(round + 1) % 16]) >> 32 | instance->TransformedState[(round + 2) % 16] << 32) & instance->TransformedState[(round + 3) % 16];

        // Update transformed state with computed A and B
        instance->TransformedState[round] ^= ghash_multiply(A, B);
    }

    instance->NLFSR = CombinedNLFSR(A ^ B);
    instance->LFSR = LFSR(A + B);
}

void ShiftRightBy(uint64_t* array, uint32_t number)
{
    if(number == 0 || number >= 16)
        return;

    uint64_t temp;
    for (int i = 0; i < number; i++) {
        temp = array[15];
        for (int j = 15; j > 0; j--) {
            array[j] = array[j - 1];
        }
        array[0] = temp;
    }
}

//enhanced confusion layer
void FiniteStateAutomaton(MagicIceCipher* instance)
{
    uint32_t random_array[8] = {0};

    //From Memory Random
    uint64_t T1 = 0, T2 = 0, T3 = 0, T4 = 0;

    //Compute UpdatedState
    for(size_t round = 4; round < 16; round++)
    {
        for (int i = 0; i < 16; i++)
        {
            T1 = instance->State[(uint32_t)(15 - i - 1) % 16] & 0xFFFFFFFF00000000 | (instance->State[15 - i] >> 32);
            T2 = ((instance->State[15 - i] & 0x00000000FFFFFFFF) << 32) | instance->State[(uint32_t)(15 - i - 1) % 16] & 0x00000000FFFFFFFF;
            T3 = instance->NLFSR;
            T4 = instance->LFSR;

            // Extract high and low parts from T1, T2, T3, T4 and R1, R2, R3, R4
            uint32_t T1_H = (uint32_t)(T1 >> 32);
            uint32_t T1_L = (uint32_t)(T1 & 0xFFFFFFFF);
            uint32_t T2_H = (uint32_t)(T2 >> 32);
            uint32_t T2_L = (uint32_t)(T2 & 0xFFFFFFFF);
            uint32_t T3_H = (uint32_t)(T3 >> 32);
            uint32_t T3_L = (uint32_t)(T3 & 0xFFFFFFFF);
            uint32_t T4_H = (uint32_t)(T4 >> 32);
            uint32_t T4_L = (uint32_t)(T4 & 0xFFFFFFFF);
            uint32_t R1_H = (uint32_t)(instance->R1 >> 32);
            uint32_t R1_L = (uint32_t)(instance->R1 & 0xFFFFFFFF);
            uint32_t R2_H = (uint32_t)(instance->R2 >> 32);
            uint32_t R2_L = (uint32_t)(instance->R2 & 0xFFFFFFFF);
            uint32_t R3_H = (uint32_t)(instance->R3 >> 32);
            uint32_t R3_L = (uint32_t)(instance->R3 & 0xFFFFFFFF);
            uint32_t R4_H = (uint32_t)(instance->R4 >> 32);
            uint32_t R4_L = (uint32_t)(instance->R4 & 0xFFFFFFFF);

            ProcessWithAvalanche(&T1_H); ProcessWithAvalanche(&T1_L); ProcessWithAvalanche(&T2_H); ProcessWithAvalanche(&T2_L);
            ProcessWithAvalanche(&T3_H); ProcessWithAvalanche(&T3_L); ProcessWithAvalanche(&T4_H); ProcessWithAvalanche(&T4_L);
            ProcessWithAvalanche(&R1_H); ProcessWithAvalanche(&R1_L); ProcessWithAvalanche(&R2_H); ProcessWithAvalanche(&R2_L);
            ProcessWithAvalanche(&R3_H); ProcessWithAvalanche(&R3_L); ProcessWithAvalanche(&R4_H); ProcessWithAvalanche(&R4_L);

            // Mix T1, T2, T3, T4 and R1, R2, R3, R4 high and low parts
            ComplexMix(&T1_H, &R1_L, &T2_H, &R2_L, &T3_H, &R3_L, &T4_H, &R4_L);
            ComplexMix(&T3_L, &R3_H, &T4_L, &R4_H, &T1_L, &R1_H, &T2_L, &R2_H);
            ComplexMix(&R4_H, &T3_L, &R2_H, &T1_L, &R1_H, &T2_L, &R3_H, &T4_L);
            ComplexMix(&R1_L, &T2_H, &R3_L, &T4_H, &R4_L, &T3_H, &R2_L, &T1_H);

            // ZUC LinearTransform 1 and 2 based MDS
            T1_H = T1_H ^ leftRotate32(T1_H, 2) ^ leftRotate32(T1_H, 10) ^ leftRotate32(T1_H, 18) ^ leftRotate32(T1_H, 24);
            T1_L = T1_L ^ leftRotate32(T1_L, 8) ^ leftRotate32(T1_L, 14) ^ leftRotate32(T1_L, 22) ^ leftRotate32(T1_L, 30);
            T2_H = T2_H ^ leftRotate32(T2_H, 2) ^ leftRotate32(T2_H, 10) ^ leftRotate32(T2_H, 18) ^ leftRotate32(T2_H, 24);
            T2_L = T2_L ^ leftRotate32(T2_L, 8) ^ leftRotate32(T2_L, 14) ^ leftRotate32(T2_L, 22) ^ leftRotate32(T2_L, 30);
            T3_H = T3_H ^ leftRotate32(T3_H, 2) ^ leftRotate32(T3_H, 10) ^ leftRotate32(T3_H, 18) ^ leftRotate32(T3_H, 24);
            T3_L = T3_L ^ leftRotate32(T3_L, 8) ^ leftRotate32(T3_L, 14) ^ leftRotate32(T3_L, 22) ^ leftRotate32(T3_L, 30);
            T4_H = T4_H ^ leftRotate32(T4_H, 2) ^ leftRotate32(T4_H, 10) ^ leftRotate32(T4_H, 18) ^ leftRotate32(T4_H, 24);
            T4_L = T4_L ^ leftRotate32(T4_L, 8) ^ leftRotate32(T4_L, 14) ^ leftRotate32(T4_L, 22) ^ leftRotate32(T4_L, 30);

            //Custom LinearTransform 1 and 2
            R1_H = R1_H ^ rightRotate32(R1_H, 26) ^ rightRotate32(R1_H, 18) ^ rightRotate32(R1_H, 10) ^ rightRotate32(R1_H, 8);
            R1_L = R1_L ^ rightRotate32(R1_L, 22) ^ rightRotate32(R1_L, 20) ^ rightRotate32(R1_L, 12) ^ rightRotate32(R1_L, 2);
            R2_H = R2_H ^ rightRotate32(R2_H, 26) ^ rightRotate32(R2_H, 18) ^ rightRotate32(R2_H, 10) ^ rightRotate32(R2_H, 8);
            R2_L = R2_L ^ rightRotate32(R2_L, 22) ^ rightRotate32(R2_L, 20) ^ rightRotate32(R2_L, 12) ^ rightRotate32(R2_L, 2);
            R3_H = R3_H ^ rightRotate32(R3_H, 26) ^ rightRotate32(R3_H, 18) ^ rightRotate32(R3_H, 10) ^ rightRotate32(R3_H, 8);
            R3_L = R3_L ^ rightRotate32(R3_L, 22) ^ rightRotate32(R3_L, 20) ^ rightRotate32(R3_L, 12) ^ rightRotate32(R3_L, 2);
            R4_H = R4_H ^ rightRotate32(R4_H, 26) ^ rightRotate32(R4_H, 18) ^ rightRotate32(R4_H, 10) ^ rightRotate32(R4_H, 8);
            R4_L = R4_L ^ rightRotate32(R4_L, 22) ^ rightRotate32(R4_L, 20) ^ rightRotate32(R4_L, 12) ^ rightRotate32(R4_L, 2);

            // Combine high and low parts back into T1, T2, T3, T4 and R1, R2, R3, R4
            T1 = ((uint64_t)T1_H << 32) | R1_L;
            T2 = ((uint64_t)T2_H << 32) | R2_L;
            T3 = ((uint64_t)T3_H << 32) | R3_L;
            T4 = ((uint64_t)T4_H << 32) | R4_L;
            instance->R1 = ((uint64_t)R1_L << 32) | T1_H;
            instance->R2 = ((uint64_t)R2_L << 32) | T2_H;
            instance->R3 = ((uint64_t)R3_L << 32) | T3_H;
            instance->R4 = ((uint64_t)R4_L << 32) | T4_H;

            //Mix T1, T2, T3, T4 and R1, R2, R3, R4 with ghash multiply
            instance->UpdatedState[i] = (T1 + instance->R4) ^ ghash_multiply(T3, instance->R2);
            instance->UpdatedState[(i + 1) % 16] = (T2 ^ instance->R3) + ghash_multiply(T1, instance->R4);
            instance->UpdatedState[(uint32_t)(15 - i - 1) % 16] = ghash_multiply(instance->R4, T2);
            instance->UpdatedState[15 - i] = ghash_multiply(instance->R2, T4);

            //Update R1, R2, R3, R4
            instance->NLFSR = CombinedNLFSR(instance->NLFSR);
            instance->LFSR = LFSR(instance->LFSR);
            instance->R3 ^= instance->NLFSR;
            instance->R4 ^= instance->LFSR;
            instance->R1 = instance->R2 + (instance->R3 ^ instance->R4);
            instance->R2 = instance->R1 ^ (instance->R4 + instance->R3);

            ShiftRightBy(instance->State, 4);
            ShiftRightBy(instance->UpdatedState, 2);
        }
    }
}

void UpdateState(MagicIceCipher* instance)
{
    for(size_t counter = 0; counter < 16; counter++)
    {
        TransformState(instance);
        FiniteStateAutomaton(instance);
    }
}

void GenerateKeyStream(MagicIceCipher* instance, uint64_t* KeyStream, uint64_t KeyStreamSize)
{
    uint64_t UnpredictableState[16] = {0};

    for(int i=0; i<16; i++)
        UnpredictableState[i] = instance->State[i] ^ (instance->TransformedState[i] + instance->UpdatedState[i]);

    size_t ElementSize = KeyStreamSize < 32 ? KeyStreamSize : 32;
    memcpy(KeyStream, UnpredictableState, ElementSize * sizeof(uint64_t));

    UpdateState(instance);

    for(int i=0; i<16; i++)
        instance->State[i] ^= UnpredictableState[i];
}

void ResetState(MagicIceCipher* instance, const uint32_t* KeyAndIV, uint64_t KeyAndIV_Size)
{
    //Wipe State
    for(size_t i = 0; i < 16; i++)
        instance->State[i] = instance->UpdatedState[i] = instance->TransformedState[i] = 0;
    
    instance->R1 = instance->R2 = instance->R3 = instance->R4 = 0;
    instance->NLFSR = instance->LFSR = 0;

    InitializeState(instance, KeyAndIV, KeyAndIV_Size);
    UpdateState(instance);
}

int main() {
    MagicIceCipher cipher;

    // Initialize the KeyAndIV (Array size 32) in the cipher
    uint32_t KeyAndIV[32] = {0};
    
    KeyAndIV[0] = 1;
    
    // srand(1);
    // for (int i = 0; i < 32; i++) {
    //     KeyAndIV[i] = rand();
    // }
    
    // for (int i = 0; i < 32; i++) {
    //     KeyAndIV[i] = i;
    // }

    uint64_t Keystream[16] = {0};
    
    #if 0
    ResetState(&cipher, KeyAndIV, 32);
    GenerateKeyStream(&cipher, Keystream, 16);
    
    // Print out the State and TransformedState, Generated Keystream
    printf("State:\n");
    for (int i = 0; i < 16; i++) {
        printf("%llu\n", cipher.State[i]);
    }
    
    printf("\nTransformedState:\n");
    for (int i = 0; i < 16; i++) {
        printf("%llu\n", cipher.TransformedState[i]);
    }
    
    printf("\nGenerated Keystream:\n");
    for (int i = 0; i < 16; i++) {
        printf("%llu\n", Keystream[i]);
    }
    #endif
    
    // Prepare the file for writing
    FILE *file = fopen("MagicIce_Keystream_OldVersion.bin", "wb");
    if (file == NULL) {
        printf("Could not open file for writing.\n");
        return 1;
    }
    
    ResetState(&cipher, KeyAndIV, 32);
    // Generate and write 128 64-bit values to the file for a total of 128KB
    for (int i = 0; i < 1024; i++) {
        
        //printf("Use cipher round: %llu\n", i);
        GenerateKeyStream(&cipher, Keystream, 16);
        
        fwrite(Keystream, sizeof(uint64_t), 16, file);
    }
    
    ResetState(&cipher, KeyAndIV, 32);
    
    return 0;
}

MagicIceVersion3_StreamCipher.c

#define _CRT_SECURE_NO_WARNINGS	

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <time.h>

uint32_t leftRotate32(uint32_t n, unsigned int d) {
	return (n << d) | (n >> (32 - d));
}

uint32_t rightRotate32(uint32_t n, unsigned int d) {
	return (n >> d) | (n << (32 - d));
}

uint64_t leftRotate64(uint64_t n, unsigned int d) {
	return (n << d) | (n >> (64 - d));
}

uint64_t rightRotate64(uint64_t n, unsigned int d) {
	return (n >> d) | (n << (64 - d));
}

// CHACHA20 context structure
typedef struct {
	uint32_t state[16];
	uint32_t pool[16];
	int pool_index;
}
Chacha20Context;

// CHACHA20 quarter round operation
void chacha20_quarter_round(uint32_t* a, uint32_t* b, uint32_t* c, uint32_t* d) {
	*a += *b;
	*d ^= *a;
	*d = (*d << 16) | (*d >> 16);

	*c += *d;
	*b ^= *c;
	*b = (*b << 12) | (*b >> 20);

	*a += *b;
	*d ^= *a;
	*d = (*d << 8) | (*d >> 24);

	*c += *d;
	*b ^= *c;
	*b = (*b << 7) | (*b >> 25);
}

// Initialize CHACHA20 context with key, nonce, and counter
void chacha20_seed(Chacha20Context* context,
	const uint32_t* key,
	const uint32_t* nonce) {
	// Initialize state with constants and key
	const uint32_t initial_state[16] = {
		0x61707865,
		0x3320646e,
		0x79622d32,
		0x6b206574,
		key[0],
		key[1],
		key[2],
		key[3],
		key[4],
		key[5],
		key[6],
		key[7],
		0,
		nonce[0],
		nonce[1],
		nonce[2]
	};

	// Set the initial state
	memcpy(context->state, initial_state, 16 * sizeof(uint32_t));

	// Reset pool index
	context->pool_index = 0;
}

// Generate random number from CHACHA20 context
uint32_t chacha20_csprng(Chacha20Context* context) {
	// Check if pool is empty, refill it
	if (context->pool_index == 0) {
		// Generate state from current context
		uint32_t state[16];
		memcpy(state, context->state, 16 * sizeof(uint32_t));

		// Perform CHACHA20 rounds
		for (int round = 0; round < 16; round += 2) {
			// Column rounds
			chacha20_quarter_round(&state[0], &state[4], &state[8], &state[12]);
			chacha20_quarter_round(&state[1], &state[5], &state[9], &state[13]);
			chacha20_quarter_round(&state[2], &state[6], &state[10], &state[14]);
			chacha20_quarter_round(&state[3], &state[7], &state[11], &state[15]);

			// Diagonal rounds
			chacha20_quarter_round(&state[0], &state[5], &state[10], &state[15]);
			chacha20_quarter_round(&state[1], &state[6], &state[11], &state[12]);
			chacha20_quarter_round(&state[2], &state[7], &state[8], &state[13]);
			chacha20_quarter_round(&state[3], &state[4], &state[9], &state[14]);
		}

		// Add state to pool
		memcpy(context->pool, state, 16 * sizeof(uint32_t));

		// Increment counter
		context->state[12]++;
	}

	// Get random number from pool and update index
	uint32_t randomNumber = context->pool[context->pool_index];
	context->pool_index = (context->pool_index + 1) % 16;

	return randomNumber;
}

/**************************************************************************************************************/

typedef struct {
	uint64_t state[312];
	int state_index;
}
MT19937_64;

void MT19937_64_Seed(MT19937_64* instance, uint64_t seed) {
	if (seed == 0)
		seed = 1;
	instance->state[0] = seed;
	for (int i = 1; i < 312; i++) {
		instance->state[i] = (uint64_t)6364136223846793005 * (instance->state[i - 1] ^ (instance->state[i - 1] >> 62)) + i;
	}
	instance->state_index = 0;
}

uint64_t MT19937_64_Generate(MT19937_64* instance) {
	const uint64_t lower_mask = ((uint64_t)1 << 31) - 1;
	const uint64_t upper_mask = ~lower_mask;

	if (instance->state_index == 0) {
		for (int i = 0; i < 312; i++) {
			uint64_t y = (instance->state[i] & upper_mask) | (instance->state[(i + 1) % 312] & lower_mask);
			instance->state[i] = instance->state[(i + 156) % 312] ^ (y >> 1);
			if (y % 2 != 0) {
				instance->state[i] ^= (uint64_t)0xB5026F5AA96619E9;
			}
		}
	}

	uint64_t y = instance->state[instance->state_index];
	y ^= (y >> 29) & (uint64_t)0x5555555555555555;
	y ^= (y << 17) & (uint64_t)0x71D67FFFEDA60000;
	y ^= (y << 37) & (uint64_t)0xFFF7EEE000000000;
	y ^= (y >> 43);

	instance->state_index = (instance->state_index + 1) % 312;
	return y;
}

/**************************************************************************************************************/

typedef enum {
	AND_STATE,
	OR_STATE,
	NOT_STATE,
	XOR_STATE,
	XNOR_STATE,
	NAND_STATE,
	NUM_STATES
}
PRBG_StateRule;

typedef struct {
	uint64_t automaton_state[64];
}
CellularAutomaton;

void cellular_automaton_init(CellularAutomaton* automaton, MT19937_64* mt19937_64) {
	for (int i = 0; i < 64; i++)
		automaton->automaton_state[i] = MT19937_64_Generate(mt19937_64);
}

//https://rosettacode.org/wiki/Elementary_cellular_automaton/Random_number_generator
uint64_t cellular_automaton_step(CellularAutomaton* automaton, int rule) {
	// An array to store the new state of each cell in the automaton.
	uint64_t NewStateValues[64] = { 0 };

	// Variable to store the binary representation of the state for each step.
	uint64_t binaryRepresentation = 0;

	// Loop through each cell in the automaton.
	for (int i = 0; i < 64; i++)
	{
		// Get the current state of the cell at position i.
		uint64_t currentState = automaton->automaton_state[i];

		// Reset binaryRepresentation for each new cell state.
		binaryRepresentation = 0;

		// Calculate the binary representation of the current cell's state.
		for (int bitIndex = 63; bitIndex >= 0; bitIndex--)
		{
			binaryRepresentation |= (currentState & 1) << bitIndex;
			currentState = (currentState >> 1) | (currentState << 63);
		}

		// Initialize the new state for the current cell.
		uint64_t NewStateValue = 0;

		// Update the state of each cell according to the rule.
		for (int j = 0; j < 64; j++)
		{
			int shift = (7 & (binaryRepresentation>>(j == 0 ? 0 : j-1) | binaryRepresentation<<(64+1-j)));
			if (rule & ((uint64_t)1 << shift))
				NewStateValue |= ((uint64_t)1 << j);
		}

		// Store the new state of the current cell.
		NewStateValues[i] = NewStateValue;
	}

	// Update the state of all cells in the automaton after applying the rule.
	memcpy(automaton->automaton_state, NewStateValues, 64 * sizeof(uint64_t));

	// Return the binary representation of the last cell's state.
	return binaryRepresentation;
}

typedef struct {
	PRBG_StateRule state;
	CellularAutomaton* automaton;
	int bit; // Changed the type to int (boolean)
	uint64_t CA_RandomNumber;
	PRBG_StateRule stateOrder[NUM_STATES];
	int selectStateOrder[NUM_STATES];
	Chacha20Context* chacha20Context; //CHACHA20 context
}
PseudorandomBitGenerator;

void generator_init(PseudorandomBitGenerator* prbg, CellularAutomaton* automaton, MT19937_64* mt19937_64) {
	prbg->stateOrder[0] = AND_STATE;
	prbg->stateOrder[1] = OR_STATE;
	prbg->stateOrder[2] = NOT_STATE;
	prbg->stateOrder[3] = XOR_STATE;
	prbg->stateOrder[4] = XNOR_STATE;
	prbg->stateOrder[5] = NAND_STATE;

	prbg->selectStateOrder[0] = 0;
	prbg->selectStateOrder[1] = 1;
	prbg->selectStateOrder[2] = 2;
	prbg->selectStateOrder[3] = 3;
	prbg->selectStateOrder[4] = 4;
	prbg->selectStateOrder[5] = 5;

	prbg->state = MT19937_64_Generate(mt19937_64) % NUM_STATES;
	prbg->automaton = automaton;
	prbg->bit = automaton->automaton_state[0] % 2;
	prbg->CA_RandomNumber = 0;
}

int generate_bit(PseudorandomBitGenerator* prbg) {
	// Shuffle stateOrder using Fisher-Yates algorithm
	for (int i = NUM_STATES - 1; i > 0; i--)
	{
		int j = chacha20_csprng(prbg->chacha20Context) % (i + 1);
		PRBG_StateRule tempState = prbg->stateOrder[i];
		prbg->stateOrder[i] = prbg->stateOrder[j];
		prbg->stateOrder[j] = tempState;
	}

	// Randomly update prbg->state using Knuth Shuffle algorithm
	for (int i = NUM_STATES - 1; i > 0; i--)
	{
		int j = chacha20_csprng(prbg->chacha20Context) % (i + 1);
		int tempIndex = prbg->selectStateOrder[i];
		prbg->selectStateOrder[i] = prbg->selectStateOrder[j];
		prbg->selectStateOrder[j] = tempIndex;
	}

	// Execute the state cases in the shuffled order
	switch (prbg->stateOrder[prbg->selectStateOrder[prbg->CA_RandomNumber % NUM_STATES]])
	{
		case AND_STATE:
			prbg->bit = prbg->bit & (chacha20_csprng(prbg->chacha20Context) % 2);
			break;
		case OR_STATE:
			prbg->bit = prbg->bit | (chacha20_csprng(prbg->chacha20Context) % 2);
			break;
		case NOT_STATE:
			prbg->bit = !prbg->bit;
			break;
		case XOR_STATE:
			prbg->bit = prbg->bit ^ (chacha20_csprng(prbg->chacha20Context) % 2);
			break;
		case XNOR_STATE:
			prbg->bit = !(prbg->bit ^ (chacha20_csprng(prbg->chacha20Context) % 2));
			break;
		case NAND_STATE:
			prbg->bit = !(prbg->bit & (chacha20_csprng(prbg->chacha20Context) % 2));
			break;
	}
	
	int CA_Rule = chacha20_csprng(prbg->chacha20Context) % 256;
	
	if (CA_Rule)
		prbg->CA_RandomNumber = cellular_automaton_step(prbg->automaton, CA_Rule);
	else
		prbg->CA_RandomNumber = cellular_automaton_step(prbg->automaton, 30);
		
	return prbg->bit;
}

uint64_t generate_bit_64(PseudorandomBitGenerator* prbg) {
	uint64_t RandomNumber = 0;
	for (size_t bit_index = 0; bit_index < 64; bit_index++) {
		RandomNumber |= ((uint64_t)generate_bit(prbg) << bit_index);
	}

	return RandomNumber;
}

/**************************************************************************************************************/

typedef struct {
	uint64_t State[16];
	uint64_t UpdatedState[16]; //FSM result
	uint64_t TransformedState[16];

	//Registers
	uint64_t R1, R2, R3, R4;
	uint64_t NLFSR, LFSR;

	//Pseudo Random Bit Generator
	PseudorandomBitGenerator* prbg;

	//Pseudo Random Number Generator
	MT19937_64* mt19937_64;

}
MagicIceCipher;

//https://users.ece.cmu.edu/~koopman/lfsr/index.html
uint64_t LFSR(uint64_t value) {
	uint64_t lsb = value & 1; // Get least significant bit
	value >>= 1; // Shift register
	if (lsb) {
		value ^= 0x80000000000019E2; // Apply a reverse of polynomial if lsb was 1
	}
	return value;
}

//https://people.kth.se/~dubrova/nlfsr.html
//https://eprint.iacr.org/2012/166
uint64_t CombinedNLFSR(uint64_t value) {
	//0 XOR 7 XOR 8 XOR (6 AND 8 AND 10)
	uint16_t nlfsr1 = value & 0xFFFF;
	uint16_t nlfsr1_feedback = (nlfsr1 ^ (nlfsr1 >> 7) ^ (nlfsr1 >> 8) ^ ((nlfsr1 >> 6) & (nlfsr1 >> 8) & (nlfsr1 >> 10))) & 1;
	if (nlfsr1_feedback)
		nlfsr1 = (nlfsr1 >> 1) | (nlfsr1_feedback << 15);

	//0 XOR 7 XOR 8 XOR (5 AND 12 AND 14)
	uint16_t nlfsr2 = (value >> 16) & 0xFFFF;
	uint16_t nlfsr2_feedback = (nlfsr2 ^ (nlfsr2 >> 7) ^ (nlfsr2 >> 8) ^ ((nlfsr2 >> 5) & (nlfsr2 >> 12) & (nlfsr2 >> 14))) & 1;
	if (nlfsr2_feedback)
		nlfsr2 = (nlfsr2 >> 1) | (nlfsr2_feedback << 15);

	//0 XOR 1 XOR 8 XOR (7 AND 8 AND 11)
	uint16_t nlfsr3 = (value >> 32) & 0xFFFF;
	uint16_t nlfsr3_feedback = (nlfsr3 ^ (nlfsr3 >> 1) ^ (nlfsr3 >> 8) ^ ((nlfsr3 >> 7) & (nlfsr3 >> 8) & (nlfsr3 >> 11))) & 1;
	if (nlfsr3_feedback)
		nlfsr3 = (nlfsr3 >> 1) | (nlfsr3_feedback << 15);

	//0 XOR 4 XOR 8 XOR 9 XOR 10 (8 AND 12)
	uint16_t nlfsr4 = (value >> 48) & 0xFFFF;
	uint16_t nlfsr4_feedback = (nlfsr4 ^ (nlfsr4 >> 4) ^ (nlfsr4 >> 8) ^ (nlfsr4 >> 9) ^ (nlfsr4 >> 10) ^ ((nlfsr4 >> 8) & (nlfsr4 >> 12))) & 1;
	if (nlfsr4_feedback)
		nlfsr4 = (nlfsr4 >> 1) | (nlfsr4_feedback << 15);

	return ((uint64_t)nlfsr4 << 48) | ((uint64_t)nlfsr3 << 32) | ((uint64_t)nlfsr2 << 16) | nlfsr1;
}

uint64_t ghash_multiply(uint64_t a, uint64_t b) {
	uint64_t product = 0;
	const uint64_t POLYNOMIAL = 0xE100000000000000; // This is x^64 + x^4 + x^3 + x + 1 in binary

	for (int i = 0; i < 64; i++) {
		if (b & 1) {
			product ^= a;
		}

		//boolean type
		int high_bit_set = a & 0x8000000000000000;

		a <<= 1;

		if (high_bit_set) {
			a ^= POLYNOMIAL;
		}

		b >>= 1;
	}

	return product;
}

void TransformState(MagicIceCipher* instance) {
	const uint64_t MASK1 = 0xAAAAAAAAAAAAAAAA;
	const uint64_t MASK2 = 0x5555555555555555;

	// Shift rows: Rotate cells
	// Row 1 need left rotate 0 cells, Row 2 need left rotate 1 cell, Row 3 need left rotate 2 cells, Row 3 need left rotate 3 cells.
	memcpy(instance->TransformedState, instance->State, 16 * sizeof(uint64_t));
	uint64_t S;
	for (uint32_t i = 1; i < 4; i++)
	{
		uint32_t counter = 0;
		while (counter < i) {
			S = instance->TransformedState[4 * i];

			for (uint32_t k = 1; k < 4; k++)
				instance->TransformedState[4 * i + k - 1] = instance->TransformedState[4 * i + k];

			instance->TransformedState[4 * i + 4 - 1] = S;
			counter++;
		}
	}

	uint64_t A = 0, B = 0;

	// Mix columns: Finite field multiplication with GHASH part
	// For simplicity, here we implement a kind of mixing operation inspired by AES
	for (size_t round = 0; round < 16; round++) {
		A ^= instance->State[round] ^ (instance->State[round] & MASK1) + (instance->State[15 - round] & MASK2);
		B ^= ((~instance->TransformedState[(round + 1) % 16]) >> 32 | instance->TransformedState[(round + 2) % 16] << 32) & instance->TransformedState[(round + 3) % 16];

		// Update transformed state with computed A and B
		instance->TransformedState[round] ^= ghash_multiply(A, B);
	}

	instance->NLFSR = CombinedNLFSR(A ^ B);
	instance->LFSR = LFSR(A + B);
}

void ShiftRightBy(uint64_t* array, uint32_t number) {
	if (number == 0 || number >= 16)
		return;

	uint64_t temp;
	for (int i = 0; i < number; i++) 
	{
		temp = array[15];
		for (int j = 15; j > 0; j--)
			array[j] = array[j - 1];
		
		array[0] = temp;
	}
}

//enhanced confusion layer
void FiniteStateAutomaton(MagicIceCipher* instance) {
	uint64_t T1 = 0, T2 = 0, T3 = 0, T4 = 0;

	//Compute UpdatedState
	for (size_t round = 4; round < 16; round++) {
		T1 = instance->State[(uint32_t)(15 - round - 1) % 16] & 0xFFFFFFFF00000000 | (instance->State[15 - round] >> 32);
		T2 = ((instance->State[15 - round] & 0x00000000FFFFFFFF) << 32) | instance->State[(uint32_t)(15 - round - 1) % 16] & 0x00000000FFFFFFFF;
		T3 = instance->NLFSR;
		T4 = instance->LFSR;

		// Extract high and low parts from T1, T2, T3, T4 and R1, R2, R3, R4
		uint32_t T1_H = (uint32_t)(T1 >> 32);
		uint32_t T1_L = (uint32_t)(T1 & 0xFFFFFFFF);
		uint32_t T2_H = (uint32_t)(T2 >> 32);
		uint32_t T2_L = (uint32_t)(T2 & 0xFFFFFFFF);
		uint32_t T3_H = (uint32_t)(T3 >> 32);
		uint32_t T3_L = (uint32_t)(T3 & 0xFFFFFFFF);
		uint32_t T4_H = (uint32_t)(T4 >> 32);
		uint32_t T4_L = (uint32_t)(T4 & 0xFFFFFFFF);
		uint32_t R1_H = (uint32_t)(instance->R1 >> 32);
		uint32_t R1_L = (uint32_t)(instance->R1 & 0xFFFFFFFF);
		uint32_t R2_H = (uint32_t)(instance->R2 >> 32);
		uint32_t R2_L = (uint32_t)(instance->R2 & 0xFFFFFFFF);
		uint32_t R3_H = (uint32_t)(instance->R3 >> 32);
		uint32_t R3_L = (uint32_t)(instance->R3 & 0xFFFFFFFF);
		uint32_t R4_H = (uint32_t)(instance->R4 >> 32);
		uint32_t R4_L = (uint32_t)(instance->R4 & 0xFFFFFFFF);

		for (int i = 0; i < 16; i++) {
			// ZUC LinearTransform 1 and 2 based MDS
			T1_H = T1_H ^ leftRotate32(T1_H, 2) ^ leftRotate32(T1_H, 10) ^ leftRotate32(T1_H, 18) ^ leftRotate32(T1_H, 24);
			T1_L = T1_L ^ leftRotate32(T1_L, 8) ^ leftRotate32(T1_L, 14) ^ leftRotate32(T1_L, 22) ^ leftRotate32(T1_L, 30);
			T2_H = T2_H ^ leftRotate32(T2_H, 2) ^ leftRotate32(T2_H, 10) ^ leftRotate32(T2_H, 18) ^ leftRotate32(T2_H, 24);
			T2_L = T2_L ^ leftRotate32(T2_L, 8) ^ leftRotate32(T2_L, 14) ^ leftRotate32(T2_L, 22) ^ leftRotate32(T2_L, 30);
			T3_H = T3_H ^ leftRotate32(T3_H, 2) ^ leftRotate32(T3_H, 10) ^ leftRotate32(T3_H, 18) ^ leftRotate32(T3_H, 24);
			T3_L = T3_L ^ leftRotate32(T3_L, 8) ^ leftRotate32(T3_L, 14) ^ leftRotate32(T3_L, 22) ^ leftRotate32(T3_L, 30);
			T4_H = T4_H ^ leftRotate32(T4_H, 2) ^ leftRotate32(T4_H, 10) ^ leftRotate32(T4_H, 18) ^ leftRotate32(T4_H, 24);
			T4_L = T4_L ^ leftRotate32(T4_L, 8) ^ leftRotate32(T4_L, 14) ^ leftRotate32(T4_L, 22) ^ leftRotate32(T4_L, 30);

			//Custom LinearTransform 1 and 2
			R1_H = R1_H ^ rightRotate32(R1_H, 26) ^ rightRotate32(R1_H, 18) ^ rightRotate32(R1_H, 10) ^ rightRotate32(R1_H, 8);
			R1_L = R1_L ^ rightRotate32(R1_L, 22) ^ rightRotate32(R1_L, 20) ^ rightRotate32(R1_L, 12) ^ rightRotate32(R1_L, 2);
			R2_H = R2_H ^ rightRotate32(R2_H, 26) ^ rightRotate32(R2_H, 18) ^ rightRotate32(R2_H, 10) ^ rightRotate32(R2_H, 8);
			R2_L = R2_L ^ rightRotate32(R2_L, 22) ^ rightRotate32(R2_L, 20) ^ rightRotate32(R2_L, 12) ^ rightRotate32(R2_L, 2);
			R3_H = R3_H ^ rightRotate32(R3_H, 26) ^ rightRotate32(R3_H, 18) ^ rightRotate32(R3_H, 10) ^ rightRotate32(R3_H, 8);
			R3_L = R3_L ^ rightRotate32(R3_L, 22) ^ rightRotate32(R3_L, 20) ^ rightRotate32(R3_L, 12) ^ rightRotate32(R3_L, 2);
			R4_H = R4_H ^ rightRotate32(R4_H, 26) ^ rightRotate32(R4_H, 18) ^ rightRotate32(R4_H, 10) ^ rightRotate32(R4_H, 8);
			R4_L = R4_L ^ rightRotate32(R4_L, 22) ^ rightRotate32(R4_L, 20) ^ rightRotate32(R4_L, 12) ^ rightRotate32(R4_L, 2);

			//TwilightDreamCryptHash: Compress
			T1 = ((uint64_t)(R4_L ^ T3_H) << 32) | (R2_L ^ T1_H);
			T2 = ((uint64_t)(R4_H ^ T3_L) << 32) | (R2_H ^ T1_L);
			T3 = ((uint64_t)(R1_L ^ T2_H) << 32) | (R3_L ^ T4_H);
			T4 = ((uint64_t)(R1_H ^ T2_L) << 32) | (R3_H ^ T4_L);
			instance->R1 = ((uint64_t)(R4_L ^ R4_H) << 32) | (T3_H ^ R2_L);
			instance->R2 = ((uint64_t)(T1_H ^ T1_L) << 32) | (T3_L ^ R2_H);
			instance->R3 = ((uint64_t)(R1_L ^ R1_H) << 32) | (T2_H ^ R3_L);
			instance->R4 = ((uint64_t)(T4_H ^ T4_L) << 32) | (T2_L ^ R3_H);

			//Mix T1, T2, T3, T4 and R1, R2, R3, R4 with ghash multiply
			//instance->UpdatedState[i] = (T1 + instance->R4) ^ ghash_multiply(T3, instance->R2);
			//instance->UpdatedState[(i + 1) % 16] = (T2 ^ instance->R3) + ghash_multiply(T1, instance->R4);
			//instance->UpdatedState[(uint32_t)(15 - i - 1) % 16] = ghash_multiply(instance->R4, T2);
			//instance->UpdatedState[15 - i] = ghash_multiply(instance->R2, T4);

			T2 ^= leftRotate64((T1 & instance->R1), 1);
			T1 ^= (T2 | instance->R2);
			T3 ^= (T4 | instance->R4);
			T4 ^= leftRotate64((T3 & instance->R3), 1);

			instance->UpdatedState[i] ^= T1 + generate_bit_64(instance->prbg);
			instance->UpdatedState[(i + 1) % 16] ^= T2 + generate_bit_64(instance->prbg);
			instance->UpdatedState[(uint32_t)(15 - i - 1) % 16] ^= T3 + generate_bit_64(instance->prbg);
			instance->UpdatedState[15 - i] ^= T4 + generate_bit_64(instance->prbg);

			//Update R1, R2, R3, R4
			instance->NLFSR = CombinedNLFSR(instance->NLFSR);
			instance->LFSR = LFSR(instance->LFSR);
			instance->R1 += instance->NLFSR;
			instance->R2 += instance->LFSR;
			instance->NLFSR = CombinedNLFSR(instance->NLFSR);
			instance->LFSR = LFSR(instance->LFSR);
			instance->R3 += instance->NLFSR;
			instance->R4 += instance->LFSR;

			ShiftRightBy(instance->UpdatedState, 2);
		}

		ShiftRightBy(instance->State, 1);
	}
}

void UpdateState(MagicIceCipher* instance) {
	for (size_t counter = 0; counter < 16; counter++)
	{
		TransformState(instance);
		FiniteStateAutomaton(instance);
	}
}

void GenerateKeyStream(MagicIceCipher* instance, uint64_t* KeyStream, uint64_t KeyStreamSize) {
	UpdateState(instance);
	
	uint64_t UnpredictableState[16] = { 0 };

	for (int i = 0; i < 16; i++)
		UnpredictableState[i] = instance->State[i] ^ (instance->TransformedState[i] + instance->UpdatedState[i]);

	size_t ElementSize = KeyStreamSize < 32 ? KeyStreamSize : 32;
	memcpy(KeyStream, UnpredictableState, ElementSize * sizeof(uint64_t));

	for (int i = 0; i < 16; i++)
		instance->State[i] ^= UnpredictableState[i];
}

void ResetState
(
	MagicIceCipher* instance,
	Chacha20Context* chacha20_instance,
	MT19937_64* mt19937_64,
	CellularAutomaton* automaton,
	PseudorandomBitGenerator* prbg,
	const uint32_t* KeyAndIV,
	uint64_t KeyAndIV_Size
) {
	//Wipe State
	for (size_t i = 0; i < 16; i++)
		instance->State[i] = instance->UpdatedState[i] = instance->TransformedState[i] = 0;

	instance->R1 = instance->R2 = instance->R3 = instance->R4 = 0;
	instance->NLFSR = instance->LFSR = 0;

	for (int i = 0; i < 312; i++)
		mt19937_64->state[i] = 0;
	mt19937_64->state_index = 0;

	if (KeyAndIV == NULL)
		return;

	if (KeyAndIV_Size > 8 && KeyAndIV_Size <= 32)
	{
		uint32_t TemporaryState[32] = {0};

		memcpy(TemporaryState, KeyAndIV, KeyAndIV_Size * sizeof(uint32_t));

		uint64_t Seed = (uint64_t)(KeyAndIV[KeyAndIV_Size - 2]) << 32 | KeyAndIV[KeyAndIV_Size - 1];
		MT19937_64_Seed(mt19937_64, Seed);

		uint32_t Chacha20Nonce[3] = {0};
		
		for (size_t i = 0; i < 3; i++)
			Chacha20Nonce[i] = MT19937_64_Generate(mt19937_64);

		chacha20_seed(chacha20_instance, TemporaryState, Chacha20Nonce);
		cellular_automaton_init(automaton, mt19937_64);
		generator_init(prbg, automaton, mt19937_64);

		instance->prbg = prbg;
		instance->prbg->automaton = automaton;
		instance->prbg->chacha20Context = chacha20_instance;
		instance->mt19937_64 = mt19937_64;

		for (size_t i = 0; i < 32; i += 8)
			TemporaryState[i] = chacha20_csprng(instance->prbg->chacha20Context);

		memcpy(instance->State, &TemporaryState, KeyAndIV_Size * sizeof(uint32_t));
	}
	else
	{
		uint32_t TemporaryState[32] = {0};
	
		uint32_t Chacha20Key[8] = {0};
		
		uint32_t Chacha20Nonce[3] = {0};
		
		MT19937_64_Seed(mt19937_64, 1);
		
		for (size_t i = 0; i < 8; i++)
			Chacha20Key[i] = MT19937_64_Generate(mt19937_64);
		
		for (size_t i = 0; i < 3; i++)
			Chacha20Nonce[i] = MT19937_64_Generate(mt19937_64);
			
		chacha20_seed(chacha20_instance, Chacha20Key, Chacha20Nonce);
		cellular_automaton_init(automaton, mt19937_64);
		generator_init(prbg, automaton, mt19937_64);
		
		instance->prbg = prbg;
		instance->prbg->automaton = automaton;
		instance->prbg->chacha20Context = chacha20_instance;
		instance->mt19937_64 = mt19937_64;

		for (size_t i = 0; i < 32; i += 8)
			TemporaryState[i] = chacha20_csprng(instance->prbg->chacha20Context);

		memcpy(instance->State, &TemporaryState, KeyAndIV_Size * sizeof(uint32_t));
	}

	UpdateState(instance);
}

int main() {
	MagicIceCipher cipher;

	// Initialize the KeyAndIV (Array size 32) in the cipher
	uint32_t KeyAndIV[32] = {0};

	KeyAndIV[0] = 1;

	// srand(1);
	// for (int i = 0; i < 32; i++) {
	//     KeyAndIV[i] = rand();
	// }

	// for (int i = 0; i < 32; i++) {
	//     KeyAndIV[i] = i;
	// }

	uint64_t Keystream[16] = {0};

	// ResetState(&cipher, KeyAndIV, 32);
	// GenerateKeyStream(&cipher, Keystream, 16);

	// // Print out the State and TransformedState, Generated Keystream
	// printf("State:\n");
	// for (int i = 0; i < 16; i++) {
	//     printf("%llu\n", cipher.State[i]);
	// }

	// printf("\nTransformedState:\n");
	// for (int i = 0; i < 16; i++) {
	//     printf("%llu\n", cipher.TransformedState[i]);
	// }

	// printf("\nGenerated Keystream:\n");
	// for (int i = 0; i < 16; i++) {
	//     printf("%llu\n", Keystream[i]);
	// }

	// Prepare the file for writing
	FILE* file = fopen("MagicIce_Keystream_Version3.bin", "wb");
	if (file == NULL) {
		printf("Could not open file for writing.\n");
		return 1;
	}

	Chacha20Context chacha20_instance;
	MT19937_64 mt19937_64;
	CellularAutomaton automaton;
	PseudorandomBitGenerator prbg;
	ResetState(&cipher, &chacha20_instance, &mt19937_64, &automaton, &prbg, KeyAndIV, 32);
	// Generate and write 128 64-bit values to the file for a total of 128KB
	for (int i = 0; i < 1024; i++) {

		printf("Use cipher round: %d\n", i);

		GenerateKeyStream(&cipher, Keystream, 16);

		fwrite(Keystream, sizeof(uint64_t), 16, file);
	}

	return 0;
}