#include <stdio.h>
#include <sys/random.h>
#include <stdlib.h>
#include <gmp.h>

/**
* This program simulates a 2D quantum field inspired by ideas from quantum gravity
* and topological quantum field theory (TQFT). Think of it like a grid where each
* point represents a little quantum state with its complex field, energy, and
* some topological charge (basically a fancy number between -1 and 1). The fun
* part is watching how these fields evolve, interact, and distribute energy
* across the whole grid.
*
* We're using a grid to represent space-time. Each point (or cell) evolves based
* on its own quantum state and the states of its neighbors. This is similar to
* lattice gauge theory, where you break up continuous quantum fields into little
* chunks so we can simulate them on a computer.
*
* Each point has a quantum field, which is just a complex number (real and
* imaginary parts). These fields control the energy at each point. The fields
* evolve recursively, which means they change based on nearby points and
* sometimes even points far away due to quantum entanglement.
*
* Every point on the grid has a topological charge, which is a number between -1
* and 1. This charge acts as an invariant and affects how the quantum fields
* evolve over time. Topological charges are key in theories like TQFT and help
* stabilize the quantum states.
*
* The energy at each point comes from the quantum field's magnitude (|φ|²). The
* field evolves through a recursive, fractal-like process. This makes the system
* behave in complex, sometimes chaotic ways. There's also a small chance that
* points far away from each other interact, kind of like how quantum entanglement
* works.
*
* To prevent the energy from becoming uncontrollable, we use a technique called
* renormalization. This process redistributes the energy across the grid,
* ensuring a balanced and stable energy distribution.
*
* There's also an optional energy normalization step that makes sure the total
* energy stays within a reasonable range as the system evolves.
*
* Probably not scientifically accurate, but it's a fun way to play with quantum
* field theory concepts and fractal recursion.
**/

// Quantum lattice size
// Note: The larger the grid size, the slower the simulation will be. You better have a good CPU for this...
// Reduce the grid size if you want to run this on a potato.
#define GRID_SIZE 120
// Number of simulation steps. This is arbitrary though, the important part is the journey, not the destination.
#define TIME_STEPS 10000
// Probability of non-local entanglement
#define NONLOCAL_PROB 0.001
// Renormalization scale
#define RENORM_SCALE 5
// Precision for high-precision arithmetic
#define PRECISION 512
// Small epsilon value to avoid zero division or underflow
#define EPSILON 1e-10
// Enable energy normalization
#define ENABLE_ENERGY_NORMALIZATION 0
// Enable fireworks
#define ENABLE_FIREWORKS 1
// Print state every N time steps
#define REPORT_EVERY 100
// Print quantum field state every N time steps (must be more and multiple of REPORT_EVERY)
#define VISUALS_EVERY 1000


typedef struct {
    mpf_t real; // Real part of quantum field (high precision)
    mpf_t imag; // Imaginary part of quantum field (high precision)
    mpf_t energy_density; // Energy density at this point (high precision)
    mpf_t topological_charge; // Topological charge (continuous values)
} QuantumField;

void initialize_field(QuantumField field[GRID_SIZE][GRID_SIZE]) {
    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            mpf_init2(field[i][j].real, PRECISION);
            mpf_init2(field[i][j].imag, PRECISION);
            mpf_init2(field[i][j].energy_density, PRECISION);
            mpf_init2(field[i][j].topological_charge, PRECISION);

            // Set random real and imaginary values for the quantum field
            // in the range [-10, 10]. This is the initial state of the field and is random.
            mpf_set_d(field[i][j].real, rand() / (double) RAND_MAX * 20.0 - 10.0);
            mpf_set_d(field[i][j].imag, rand() / (double) RAND_MAX * 20.0 - 10.0);

            // Set energy density to a non-zero value
            mpf_set_d(field[i][j].energy_density, 10.0);

            // Set random continuous topological charge in range [-1, 1]
            mpf_set_d(field[i][j].topological_charge, rand() / (double) RAND_MAX * 2.0 - 1.0);
        }
    }
}

void compute_scaling_factor(const QuantumField *field, mpf_t scaling_factor) {
    mpf_t temp;
    mpf_init2(temp, PRECISION);

    // scaling_factor = 10.0 / (1 + energy_density)
    mpf_set_d(temp, 1.0);
    mpf_add(temp, temp, field->energy_density);
    mpf_set_d(scaling_factor, 10.0);
    mpf_div(scaling_factor, scaling_factor, temp);

    mpf_mul(scaling_factor, scaling_factor, field->topological_charge);

    mpf_clear(temp);
}

// QF evolution function (fractal recursion with hyperloops (the other kind, go away, Elon))
void recursive_evolve( // NOLINT(*-no-recursion)
    QuantumField field[GRID_SIZE][GRID_SIZE],
    const int x,
    const int y,
    const int depth,
    mpf_t real_res,
    mpf_t imag_res
) {
    if (depth <= 0) {
        mpf_set(real_res, field[x][y].real);
        mpf_set(imag_res, field[x][y].imag);
        return;
    }

    const int nx = (x + depth + GRID_SIZE) % GRID_SIZE;
    const int ny = (y + depth + GRID_SIZE) % GRID_SIZE;

    mpf_t temp_real, temp_imag;
    mpf_init2(temp_real, PRECISION);
    mpf_init2(temp_imag, PRECISION);

    recursive_evolve(field, nx, ny, depth - 1, temp_real, temp_imag);

    mpf_t scaling_factor;
    mpf_init2(scaling_factor, PRECISION);
    compute_scaling_factor(&field[x][y], scaling_factor);

    mpf_mul(temp_real, temp_real, scaling_factor);
    mpf_mul(temp_imag, temp_imag, scaling_factor);

    mpf_add(real_res, field[x][y].real, temp_real);
    mpf_add(imag_res, field[x][y].imag, temp_imag);

    mpf_clear(temp_real);
    mpf_clear(temp_imag);
    mpf_clear(scaling_factor);
}

// Quantum gravity-inspired field interaction (non-local + anisotropic effects)
void evolve_field(QuantumField field[GRID_SIZE][GRID_SIZE]) {
    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            // Random recursion depth (min 3, max 5)
            const int depth = rand() % 3 + 3;

            mpf_t real_res, imag_res;
            mpf_init2(real_res, PRECISION);
            mpf_init2(imag_res, PRECISION);

            recursive_evolve(field, i, j, depth, real_res, imag_res);

            mpf_set(field[i][j].real, real_res);
            mpf_set(field[i][j].imag, imag_res);

            mpf_mul(field[i][j].energy_density, real_res, real_res); // real^2
            mpf_t imag_sq;
            mpf_init2(imag_sq, PRECISION);
            mpf_mul(imag_sq, imag_res, imag_res); // imag^2
            mpf_add(field[i][j].energy_density, field[i][j].energy_density, imag_sq);

            mpf_clear(real_res);
            mpf_clear(imag_res);
            mpf_clear(imag_sq);
        }
    }
}

// Renormalization step to avoid infinities in the quantum field (numerical stabilization)
void renormalize_field(QuantumField field[GRID_SIZE][GRID_SIZE]) {
    // This is going to be so fun to decode in a few years...
    for (int scale = 1; scale <= RENORM_SCALE; scale++) {
        for (int i = 0; i < GRID_SIZE; i += scale) {
            for (int j = 0; j < GRID_SIZE; j += scale) {
                mpf_t total_real, total_imag, total_energy;
                mpf_init2(total_real, PRECISION);
                mpf_init2(total_imag, PRECISION);
                mpf_init2(total_energy, PRECISION);

                // Sum over local fields to renormalize on different scales
                for (int dx = 0; dx < scale; dx++) {
                    for (int dy = 0; dy < scale; dy++) {
                        const int nx = (i + dx) % GRID_SIZE;
                        const int ny = (j + dy) % GRID_SIZE;
                        mpf_add(total_real, total_real, field[nx][ny].real);
                        mpf_add(total_imag, total_imag, field[nx][ny].imag);
                        mpf_add(total_energy, total_energy, field[nx][ny].energy_density);
                    }
                }

                // Renormalize by averaging the quantum states
                mpf_t scale_factor;
                mpf_init2(scale_factor, PRECISION);
                mpf_set_d(scale_factor, scale * scale);
                mpf_div(total_real, total_real, scale_factor);
                mpf_div(total_imag, total_imag, scale_factor);
                mpf_div(total_energy, total_energy, scale_factor);

                // Redistribute the averaged state to all lattice points in the block
                for (int dx = 0; dx < scale; dx++) {
                    for (int dy = 0; dy < scale; dy++) {
                        const int nx = (i + dx) % GRID_SIZE;
                        const int ny = (j + dy) % GRID_SIZE;
                        mpf_set(field[nx][ny].real, total_real);
                        mpf_set(field[nx][ny].imag, total_imag);
                        mpf_set(field[nx][ny].energy_density, total_energy);
                    }
                }

                mpf_clear(total_real);
                mpf_clear(total_imag);
                mpf_clear(total_energy);
                mpf_clear(scale_factor);
            }
        }
    }
}

// Optionally normalize total energy to conserve energy (see ENABLE_ENERGY_NORMALIZATION)
void normalize_total_energy(QuantumField field[GRID_SIZE][GRID_SIZE], mpf_t initial_total_energy) {
    mpf_t current_total_energy, scaling_factor, allowed_fluctuation;
    mpf_init2(current_total_energy, PRECISION);
    mpf_init2(scaling_factor, PRECISION);
    mpf_init2(allowed_fluctuation, PRECISION);

    mpf_set_d(current_total_energy, 0.0);
    mpf_set_d(allowed_fluctuation, 10); // Allow 10% energy fluctuation. Param?

    // Calculate current total energy
    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            mpf_add(current_total_energy, current_total_energy, field[i][j].energy_density);
        }
    }

    // Check if current total energy is within the allowed fluctuation range
    mpf_t upper_limit, lower_limit;
    mpf_init2(upper_limit, PRECISION);
    mpf_init2(lower_limit, PRECISION);

    mpf_mul(upper_limit, initial_total_energy, allowed_fluctuation);
    mpf_add(upper_limit, upper_limit, initial_total_energy);

    mpf_mul(lower_limit, initial_total_energy, allowed_fluctuation);
    mpf_sub(lower_limit, initial_total_energy, lower_limit);

    if (mpf_cmp(current_total_energy, lower_limit) > 0 && mpf_cmp(current_total_energy, upper_limit) < 0) {
        // Total energy is within the allowed fluctuation range, no need to normalize
        mpf_clear(current_total_energy);
        mpf_clear(scaling_factor);
        mpf_clear(allowed_fluctuation);
        mpf_clear(upper_limit);
        mpf_clear(lower_limit);
        return;
    }

    // scaling_factor = initial_total_energy / current_total_energy
    mpf_div(scaling_factor, initial_total_energy, current_total_energy);

    // Adjust fields to conserve energy
    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            mpf_mul(field[i][j].energy_density, field[i][j].energy_density, scaling_factor);
            mpf_mul(field[i][j].real, field[i][j].real, scaling_factor);
            mpf_mul(field[i][j].imag, field[i][j].imag, scaling_factor);
        }
    }

    mpf_clear(current_total_energy);
    mpf_clear(scaling_factor);
    mpf_clear(allowed_fluctuation);
    mpf_clear(upper_limit);
    mpf_clear(lower_limit);
}

void visualize_field(QuantumField field[GRID_SIZE][GRID_SIZE]) {
    mpf_t max_energy, min_energy, range;
    mpf_init2(max_energy, PRECISION);
    mpf_init2(min_energy, PRECISION);
    mpf_init2(range, PRECISION);
    mpf_set_d(max_energy, 0.0);
    mpf_set_d(min_energy, 1e10);

    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            if (mpf_cmp(field[i][j].energy_density, max_energy) > 0) {
                mpf_set(max_energy, field[i][j].energy_density);
            }
            if (mpf_cmp(field[i][j].energy_density, min_energy) < 0) {
                mpf_set(min_energy, field[i][j].energy_density);
            }
        }
    }

    mpf_sub(range, max_energy, min_energy);

    char symbols[] = " .:-=+*#%@";

    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            mpf_t normalized_energy;
            mpf_init2(normalized_energy, PRECISION);

            // Normalize energy density between 0 and 1
            mpf_sub(normalized_energy, field[i][j].energy_density, min_energy);
            mpf_div(normalized_energy, normalized_energy, range);

            int symbol_index = (int) (mpf_get_d(normalized_energy) * (sizeof(symbols) - 2));
            symbol_index = symbol_index < 0 ? 0 : symbol_index;

            printf("%c", symbols[symbol_index]);

            mpf_clear(normalized_energy);
        }
        printf("\n");
    }

    mpf_clear(max_energy);
    mpf_clear(min_energy);
    mpf_clear(range);
}

void simulate_quantum_gravity(QuantumField field[GRID_SIZE][GRID_SIZE]) {
    mpf_t initial_total_energy, prev_total_energy;
    mpf_init2(initial_total_energy, PRECISION);
    mpf_init2(prev_total_energy, PRECISION);
    mpf_set_d(initial_total_energy, 0.0);
    mpf_set_d(prev_total_energy, 0.0);

    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            mpf_add(initial_total_energy, initial_total_energy, field[i][j].energy_density);
        }
    }
    mpf_set(prev_total_energy, initial_total_energy);

    mpf_t energy_change_percentage, energy_diff;
    mpf_init2(energy_change_percentage, PRECISION);
    mpf_init2(energy_diff, PRECISION);

    for (int t = 0; t < TIME_STEPS; t++) {
        evolve_field(field);
        renormalize_field(field);

        if (ENABLE_ENERGY_NORMALIZATION == 1) {
            // Normalize total energy to conserve energy
            normalize_total_energy(field, initial_total_energy);
        }

        if (t % 100 == 0) {
            mpf_t total_energy, max_energy, min_energy, sum_energy;
            mpf_init2(total_energy, PRECISION);
            mpf_init2(max_energy, PRECISION);
            mpf_init2(min_energy, PRECISION);
            mpf_init2(sum_energy, PRECISION);

            mpf_set_d(max_energy, 0.0);
            mpf_set_d(min_energy, 1e10);
            mpf_set_d(sum_energy, 0.0);

            int low_count = 0, medium_count = 0, high_count = 0;

            for (int i = 0; i < GRID_SIZE; i++) {
                for (int j = 0; j < GRID_SIZE; j++) {
                    if (mpf_cmp(field[i][j].energy_density, max_energy) > 0) {
                        mpf_set(max_energy, field[i][j].energy_density);
                    }
                    if (mpf_cmp(field[i][j].energy_density, min_energy) < 0) {
                        mpf_set(min_energy, field[i][j].energy_density);
                    }
                    mpf_add(total_energy, total_energy, field[i][j].energy_density);
                }
            }

            mpf_t low_threshold, medium_threshold;
            mpf_init2(low_threshold, PRECISION);
            mpf_init2(medium_threshold, PRECISION);

            mpf_sub(energy_diff, max_energy, min_energy); // energy_diff = max_energy - min_energy
            mpf_div_ui(energy_diff, energy_diff, 3); // energy_diff /= 3
            mpf_add(low_threshold, min_energy, energy_diff); // low_threshold = min_energy + energy_diff

            mpf_mul_ui(energy_diff, energy_diff, 2); // energy_diff *= 2
            mpf_add(medium_threshold, min_energy, energy_diff);

            for (int i = 0; i < GRID_SIZE; i++) {
                for (int j = 0; j < GRID_SIZE; j++) {
                    if (mpf_cmp(field[i][j].energy_density, low_threshold) < 0) {
                        low_count++;
                    } else if (mpf_cmp(field[i][j].energy_density, medium_threshold) < 0) {
                        medium_count++;
                    } else {
                        high_count++;
                    }
                }
            }

            mpf_sub(energy_diff, total_energy, prev_total_energy);
            mpf_div(energy_change_percentage, energy_diff, prev_total_energy);

            mpf_mul_ui(energy_change_percentage, energy_change_percentage, 100);

            mpf_sub(energy_diff, max_energy, min_energy);

            // ================== Print simulation information ==================
            printf("Time step %d:\n", t);

            gmp_printf("  Energy Change from Previous Step: %.5Ff%%\n", energy_change_percentage);
            gmp_printf("  Energy Range (Max - Min): ");
            gmp_printf("%.E\n", energy_diff);

            printf("  Energy Distribution:\n");
            printf("    Low energy cells (<0.1): %d\n", low_count);
            printf("    Medium energy cells (0.1 - 0.5): %d\n", medium_count);
            printf("    High energy cells (>0.5): %d\n", high_count);

            // :D
            if (ENABLE_FIREWORKS == 1 && low_count == 0 && medium_count == 0) {
                printf("Your universe is too excited.....\n");
                printf("============================================\n");
                printf("============== GAME OVER ===================\n");
                printf("============================================\n");
                printf("Insert coin to continue \n");
                exit(0);
            }

            // Update the previous total energy for the next loop
            mpf_set(prev_total_energy, total_energy);

            if (t % VISUALS_EVERY == 0) {
                visualize_field(field);
            }
            // ================== End of simulation information ==================

            mpf_clear(total_energy);
            mpf_clear(max_energy);
            mpf_clear(min_energy);
            mpf_clear(sum_energy);
        }
    }

    mpf_clear(initial_total_energy);
    mpf_clear(prev_total_energy);
    mpf_clear(energy_change_percentage);
    mpf_clear(energy_diff);
}

void clear_field(QuantumField field[GRID_SIZE][GRID_SIZE]) {
    for (int i = 0; i < GRID_SIZE; i++) {
        for (int j = 0; j < GRID_SIZE; j++) {
            mpf_clear(field[i][j].real);
            mpf_clear(field[i][j].imag);
            mpf_clear(field[i][j].energy_density);
            mpf_clear(field[i][j].topological_charge);
        }
    }
}

void seed_random_from_urandom() {
    unsigned int seed;

    if (getrandom(&seed, sizeof(seed), 0) == -1) {
        perror("Failed to get random");
        exit(EXIT_FAILURE);
    }

    srand(seed);
}

int main() {
    seed_random_from_urandom();

    QuantumField field[GRID_SIZE][GRID_SIZE];
    initialize_field(field);

    simulate_quantum_gravity(field);

    clear_field(field);

    return 0;
}