393 lines
12 KiB
C
393 lines
12 KiB
C
//
|
|
// 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;
|
|
}
|
|
|
|
int csv_enabled = 0;
|
|
|
|
void set_csv_enabled(int enabled) {
|
|
csv_enabled = enabled;
|
|
}
|
|
|
|
void csv_header() {
|
|
if (!csv_enabled)
|
|
return;
|
|
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 (!csv_enabled)
|
|
return 1;
|
|
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
|
|
|