//
// Created by dukantic on 29/09/2026.
//
#ifndef RABBIT_HOLE_H
#define RABBIT_HOLE_H
#include "rabbit_hole.h"

#include <stdio.h>

#include <math.h>

#include "mt19937ar.h"
__uint128_t i,j;
double x_gauss[2];

void init_rabbits(rabbit_hole_t *rabbits) {
    rabbits->female_litter=0;
    for (i = 0; i < N_MATURE; i++) {
        rabbits->male_mature[i] =0;
        rabbits->female_mature[i] =0;
    }

    for (i = 0; i < N_KITTEN; i++) {
        rabbits->kitten_male[i] =0;
        rabbits->kitten_female[i] =0;
    }
}

double sk_sigma_m, sk_mu_m, sk_sigma_f, sk_mu_f;
double sk_draw;
__uint128_t sk_temp_male, sk_temp_female;
void shift_kitten(rabbit_hole_t *rabbits, __uint128_t *new_male, __uint128_t *new_female) {
    (*new_male) = rabbits->kitten_male[N_KITTEN-1];
    (*new_female) = rabbits->kitten_female[N_KITTEN-1];
    for (i = N_KITTEN - 3; i < N_KITTEN-1; i++) {
        sk_temp_male = 0;
        sk_temp_female = 0;
        if (rabbits->kitten_male[i]<100 || rabbits->kitten_female[i]<100) {
            for ( j = 0; j < rabbits->kitten_male[i]; j++) {
                sk_draw = genrand_real1();
                if (sk_draw <= 0.66)
                    sk_temp_male ++;
            }
            for (j = 0; j < rabbits->kitten_female[i]; j++) {
                sk_draw = genrand_real1();
                if (sk_draw <= 0.66)
                    sk_temp_female ++;
            }
            rabbits->kitten_male[i]-= sk_temp_male;
            rabbits->kitten_female[i]-= sk_temp_female;
            (*new_male) += sk_temp_male;
            (*new_female) += sk_temp_female;
        }
        else
        {
            gauss(x_gauss);
            sk_mu_m = (double)rabbits->kitten_male[i] * 0.66;
            sk_sigma_m = sqrt(sk_mu_m * 0.33);
            sk_mu_f = (double)rabbits->kitten_female[i] * 0.66;
            sk_sigma_f = sqrt(sk_mu_f * 0.33);
            x_gauss[0] = CLAMP_BTW(x_gauss[0] * sk_sigma_m + sk_mu_m, 0.0, rabbits->kitten_male[i]);
            x_gauss[1] = CLAMP_BTW( x_gauss[1] * sk_sigma_f + sk_mu_f, 0.0, rabbits->kitten_female[i]);
            (*new_male) += (__uint128_t)roundl(x_gauss[0]);
            (*new_female) += (__uint128_t)roundl(x_gauss[1]);
            rabbits->kitten_male[i] -= (__uint128_t)roundl(x_gauss[0]);
            rabbits->kitten_female[i] -= (__uint128_t)roundl(x_gauss[1]);
        }

    }
    for (i = N_KITTEN-1; i >=1; i--) {
        rabbits->kitten_male[i] = rabbits->kitten_male[i-1];
        rabbits->kitten_female[i] = rabbits->kitten_female[i-1];
    }
    rabbits->kitten_male[0]   = 0;
    rabbits->kitten_female[0] = 0;
}

__uint128_t *males;
__uint128_t *females;
double draw;
void shift_mature(rabbit_hole_t* rabbits) {
    males = rabbits->male_mature;
    females = rabbits->female_mature;

    for (i = N_MATURE-1; i >=1; i--) {
        males[i] = males[i-1];
        females[i] = females[i-1];
    }
    males[0] = 0;
    females[0] = 0;
}

double dp_threshold_death, dp_draw_death;
double sigma_m,sigma_f, mu_m, mu_f;
__uint128_t temp_male, temp_female;
void death_process(rabbit_hole_t *rabbits) {
    for (i =0; i < N_MATURE; i++) {
        dp_threshold_death = 1 - (0.4 + ((i >= 9)? (double)(i-9) * 0.1: 0)) / 12;
        if (rabbits->male_mature[i] < 100 || rabbits->female_mature[i] < 100) {
            temp_male = 0;
            temp_female = 0;
            for (j = 0; j < rabbits->male_mature[i]; j++) {
                dp_draw_death = genrand_real1();
                if (dp_draw_death <= dp_threshold_death)
                    temp_male ++;
            }

            for (j = 0; j < rabbits->female_mature[i]; j++) {
                dp_draw_death = genrand_real1();
                if (dp_draw_death <= dp_threshold_death)
                    temp_female ++;
            }
            rabbits->male_mature[i] = temp_male;
            rabbits->female_mature[i] = temp_female;
        }else {
            mu_f = (double)rabbits->female_mature[i] * dp_threshold_death;
            sigma_f = sqrt(mu_f * (1.0-dp_threshold_death));
            mu_m = (double)rabbits->male_mature[i] * dp_threshold_death;
            sigma_m = sqrt(mu_m * (1.0-dp_threshold_death));

            gauss(x_gauss);
            x_gauss[0] = CLAMP_BTW(x_gauss[0] * sigma_f + mu_f, 0.0, (double)rabbits->female_mature[i]);
            x_gauss[1] = CLAMP_BTW(x_gauss[1] * sigma_m + mu_m, 0.0, (double)rabbits->male_mature[i]);
            rabbits->female_mature[i] = (__uint128_t)roundl(x_gauss[0]);
            rabbits->male_mature[i] = (__uint128_t)roundl(x_gauss[1]);
        }
    }

    for (i = 0; i < N_KITTEN; i++) {
        if (rabbits->kitten_male[i] < 100 || rabbits->kitten_female[i] < 100) {
            temp_male = 0;
            temp_female = 0;
            for (j = 0; j < rabbits->kitten_male[i]; j++) {
                dp_draw_death = genrand_real1();
                if (dp_draw_death > 0.05) //0.35/7
                    temp_male ++;
            }
            for (j = 0; j < rabbits->kitten_female[i]; j++) {
                dp_draw_death = genrand_real1();
                if (dp_draw_death > 0.05) //0.35/7
                    temp_female ++;
            }
            rabbits->kitten_male[i] = temp_male;
            rabbits->kitten_female[i] = temp_female;
        }else {
            mu_f = (double)rabbits->kitten_female[i] * 0.95;
            sigma_f = sqrt(mu_f * 0.05);
            mu_m = (double)rabbits->kitten_male[i] * 0.95;
            sigma_m = sqrt(mu_m * 0.05);

            gauss(x_gauss);
            x_gauss[0] = x_gauss[0] * sigma_f + mu_f;
            x_gauss[1] = x_gauss[1] * sigma_m + mu_m;
            rabbits->kitten_female[i] = (__uint128_t)roundl(x_gauss[0]);
            rabbits->kitten_male[i] = (__uint128_t)roundl(x_gauss[1]);
        }
    }
}

__uint128_t nb_females;
int taille1;
__uint128_t new_litter;
double draw_litter;
double ul_mu_k,ul_sigma_k, ul_mu_sex, ul_sigma_sex;
double ul_mu_litter, ul_sigma_litter;
void update_litter(rabbit_hole_t *rabbits, __uint128_t *new_kitten_male, __uint128_t *new_kitten_female) {
    (*new_kitten_female) = 0;
    (*new_kitten_male) = 0;
if (rabbits->female_litter < 100) {
    for (i = 0; i < rabbits->female_litter; i++) {
        taille1 = 3 + (int)(4.0 * genrand_real2());
        for (int k = 0; k < taille1; k++) {
            if (genrand_real1() > 0.5) (*new_kitten_male)++;
            else                       (*new_kitten_female)++;
        }

    }
}else {
    ul_mu_k = 4.5*(double)rabbits->female_litter; // 4.5 = (b+a)/2
    ul_sigma_k = sqrt(1.118 * (double)rabbits->female_litter); //(k²−1)/12 = 1.118
    gauss(x_gauss);
    x_gauss[0] = CLAMP_BTW(x_gauss[0] * ul_sigma_k+ ul_mu_k, 0.0 , rabbits->female_litter * 9.0);

    ul_mu_sex = x_gauss[0] * 0.5;
    ul_sigma_sex = sqrt(ul_mu_sex * 0.5);
    x_gauss[1] = CLAMP_BTW(x_gauss[1] * ul_sigma_sex + ul_mu_sex , 0.0, rabbits->female_litter * 9.0);

    (*new_kitten_female) = (__uint128_t)roundl(x_gauss[1]) ;
    (*new_kitten_male) = (__uint128_t)roundl(x_gauss[0] - x_gauss[1]);
}

    nb_females = 0;
    for (i = 0; i < N_MATURE; i++) {
        nb_females += rabbits->female_mature[i];
    }

    if (nb_females < 100) {
        new_litter = 0;
        for (i = 0; i < nb_females; i++) {
            draw_litter = genrand_real1();
            if (draw_litter > 0.5)
                new_litter++;
        }
        rabbits->female_litter = new_litter;
    }else {
        ul_mu_litter = (double)nb_females * 0.5;
        ul_sigma_litter = sqrt(ul_mu_litter * 0.5);
        gauss(x_gauss);
        x_gauss[0] = CLAMP_BTW(x_gauss[0]* ul_sigma_litter + ul_mu_litter ,0.0,nb_females);
        rabbits->female_litter = (__uint128_t)roundl(x_gauss[0]);

    }
}


char total_rabbits(rabbit_hole_t *rabbits, __uint128_t *total) {
    *total = 0;
    __uint128_t total_temp;

    for (i = 0; i < N_MATURE; i++) {
        total_temp = *total;
        *total += rabbits->male_mature[i];
        *total += rabbits->female_mature[i];
        if (*total < total_temp) //overflow
        {
            return 0;
        }
    }
    for (i = 0; i < N_KITTEN; i++) {
        total_temp = *total;
        *total += rabbits->kitten_male[i];
        *total += rabbits->kitten_female[i];
        if (*total < total_temp) //overflow
        {
            return 0;
        }
    }
    return 1;
}

static void u128_to_str(__uint128_t v, char *buf)
{
    char tmp[45];
    int i = 0, j = 0;

    if (v == 0) {
        buf[0] = '0';
        buf[1] = '\0';
        return;
    }
    while (v > 0) {
        tmp[i++] = '0' + (char)(v % 10);
        v /= 10;
    }
    while (i > 0)
        buf[j++] = tmp[--i];
    buf[j] = '\0';
}

static void print_u128(__uint128_t v)
{
    char buf[45];
    u128_to_str(v, buf);
    printf("%s", buf);
}

static void print_row_u128(const char *label, const __uint128_t *arr, size_t n)
{
    fputs(label, stdout);
    fputs("| ", stdout);
    for (size_t i = 0; i < n; i++) {
        print_u128(arr[i]);
        fputs(" | ", stdout);
    }
    putchar('\n');
}


char print_rabbits(rabbit_hole_t *rabbits, __uint128_t *month)
{
    printf("During time : ");
    print_u128(*month / 12);
    printf(" years, ");
    print_u128(*month % 12);
    printf(" months\n");

    print_row_u128("Male   : ",       rabbits->male_mature,   N_MATURE);
    print_row_u128("Female : ",       rabbits->female_mature, N_MATURE);
    print_row_u128("Male kitten   : ", rabbits->kitten_male,   N_KITTEN);
    print_row_u128("Female kitten : ", rabbits->kitten_female, N_KITTEN);

    printf("\nLitter : | ");
    print_u128(rabbits->female_litter);
    printf("\n=> ");
    __uint128_t total;
    if ( ! total_rabbits(rabbits,&total)) {
        puts("======================OVERFLOW===================");
        return 0;
    }

    print_u128(total);
    printf("\n");

    puts("########################################");
    return 1;
}

void csv_header() {
    fprintf(stderr, "Month, Total Rabbit, Litters,");
    for ( i = 0; i < N_MATURE; i++) {
        fprintf(stderr,"Male (age %d),",i+1);
    }
    for (i = 0; i < N_MATURE; i++) {
        fprintf(stderr,"Female (age %d),",i+1);
    }
    for (i = 0; i < N_KITTEN; i++) {
        fprintf(stderr,"Kitten Male (age %d),",i+1);
    }
    for (i = 0; i < N_KITTEN; i++) {
        fprintf(stderr,"Kitten Female (age %d),",i+1);
    }
    fprintf(stderr, "\n");
}

char csv_rabbits(rabbit_hole_t *rabbits, const __uint128_t *month) {
    __uint128_t total;
    if (!total_rabbits(rabbits,&total))
        return 0;
    char number_string[512];
    u128_to_str(*month, number_string);
    fprintf(stderr, "%s,", number_string);
    u128_to_str(total,number_string);
    fprintf(stderr, "%s,", number_string);
    u128_to_str(rabbits->female_litter,number_string);
    fprintf(stderr, "%s,", number_string);
    for ( i = 0; i < N_MATURE; i++) {
        u128_to_str(rabbits->male_mature[i],number_string);
        fprintf(stderr,"%s,", number_string);
    }
    for (i = 0; i < N_MATURE; i++) {
        u128_to_str(rabbits->female_mature[i],number_string);
        fprintf(stderr,"%s,", number_string);
    }

    for (i = 0; i < N_KITTEN; i++) {
        u128_to_str(rabbits->kitten_male[i],number_string);
        fprintf(stderr,"%s,", number_string);
    }
    for (i = 0; i < N_KITTEN; i++) {
        u128_to_str(rabbits->kitten_female[i],number_string);
        fprintf(stderr,"%s,", number_string);
    }
    fprintf(stderr, "\n");
    return 1;
}

__uint128_t new_male;
__uint128_t new_female;
__uint128_t new_kitten_male;
__uint128_t new_kitten_female;
char process_month(rabbit_hole_t *rabbits,__uint128_t *month) {
    shift_kitten(rabbits, &new_male, &new_female);
    if (*month % 12 == 0) {
        shift_mature(rabbits);
    }
    rabbits->male_mature[0] += new_male;
    rabbits->female_mature[0] += new_female;
    death_process(rabbits);
    update_litter(rabbits, &new_kitten_male, &new_kitten_female);
    rabbits->kitten_male[0] += new_kitten_male;
    rabbits->kitten_female[0] += new_kitten_female;
    (*month)++;
    if (!csv_rabbits(rabbits,month))
        return 0;
    return 1;
}




#endif