numerical.c

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include <sys/stat.h>
#include <fcntl.h>

#include <ROOT-Sim.h>

#include <core/init.h>
#include <communication/communication.h>
#include <mm/mm.h>
#include <scheduler/process.h>
#include <scheduler/scheduler.h>
#include <statistics/statistics.h> // To have _mkdir helper function

static seed_type master_seed;

static double do_random(void)
{
    uint32_t *seed1;
    uint32_t *seed2;

    seed1 = (uint32_t *) & (current->numerical.seed);
    seed2 =
        (uint32_t *) ((char *)&(current->numerical.seed) +
              (sizeof(uint32_t)));

    *seed1 = 36969u * (*seed1 & 0xFFFFu) + (*seed1 >> 16u);
    *seed2 = 18000u * (*seed2 & 0xFFFFu) + (*seed2 >> 16u);

    // The magic number below is 1/(2^32 + 2).
    // The result is strictly between 0 and 1.
    return (((*seed1 << 16u) + (*seed1 >> 16u) + *seed2) +
        1.0) * 2.328306435454494e-10;

}

double Random(void)
{
    double ret;
    switch_to_platform_mode();

    ret = do_random();

    switch_to_application_mode();
    return ret;
}

int RandomRange(int min, int max)
{
    double ret;
    switch_to_platform_mode();

    ret = (int)floor(do_random() * (max - min + 1)) + min;

    switch_to_application_mode();
    return ret;
}

int RandomRangeNonUniform(int x, int min, int max)
{
    double ret;
    switch_to_platform_mode();

    ret = (((RandomRange(0, x) | RandomRange(min, max))) % (max - min + 1)) + min;

    switch_to_application_mode();
    return ret;
}

double Expent(double mean)
{
    double ret;
    switch_to_platform_mode();

    if (unlikely(mean < 0)) {
        rootsim_error(true,
                  "Expent() has been passed a negative mean value\n");
    }

    ret = (-mean * log(1 - do_random()));

    switch_to_application_mode();
    return ret;
}

double Normal(void)
{
    double fac, rsq, v1, v2;
    bool *iset;
    double *gset;

    iset = &current->numerical.iset;
    gset = &current->numerical.gset;

    if (*iset == false) {
        do {
            v1 = 2.0 * Random() - 1.0;
            v2 = 2.0 * Random() - 1.0;
            rsq = v1 * v1 + v2 * v2;
        } while (rsq >= 1.0 || D_EQUAL_ZERO(rsq));

        fac = sqrt(-2.0 * log(rsq) / rsq);

        // Perform Box-Muller transformation to get two normal deviates. Return one
        // and save the other for next time.
        *gset = v1 * fac;
        *iset = true;
        return v2 * fac;
    } else {
        // A deviate is already available
        *iset = false;
        return *gset;
    }
}

double Gamma(int ia)
{
    int j;
    double am, e, s, v1, v2, x, y;

    if (unlikely(ia < 1)) {
        rootsim_error(false,
                  "Gamma distribution must have a ia value >= 1. Defaulting to 1...");
        ia = 1;
    }

    if (ia < 6) {
        // Use direct method, adding waiting times
        x = 1.0;
        for (j = 1; j <= ia; j++)
            x *= Random();
        x = -log(x);
    } else {
        // Use rejection method
        do {
            do {
                do {
                    v1 = Random();
                    v2 = 2.0 * Random() - 1.0;
                } while (v1 * v1 + v2 * v2 > 1.0);

                y = v2 / v1;
                am = (double)(ia - 1);
                s = sqrt(2.0 * am + 1.0);
                x = s * y + am;
            } while (x < 0.0);

            e = (1.0 + y * y) * exp(am * log(x / am) - s * y);
        } while (Random() > e);
    }

    return x;
}

double Poisson(void)
{
    return Gamma(1);
}

int Zipf(double skew, int limit)
{
    double a = skew;
    double b = pow(2., a - 1.);
    double x, t, u, v;
    do {
        u = Random();
        v = Random();
        x = floor(pow(u, -1. / a - 1.));
        t = pow(1. + 1. / x, a - 1.);
    } while (v * x * (t - 1.) / (b - 1.) > (t / b) || x > limit);
    return (int)x;
}

static seed_type sanitize_seed(seed_type cur_seed)
{

    uint32_t *seed1 = (uint32_t *) &(cur_seed);
    uint32_t *seed2 = (uint32_t *) ((char *)&(cur_seed) + (sizeof(uint32_t)));

    // Sanitize seed1
    // Any integer multiple of 0x9068FFFF, including 0, is a bad state
    uint32_t state_orig;
    uint32_t temp;

    state_orig = *seed1;
    temp = state_orig;

    // The following is equivalent to % 0x9068FFFF, without using modulo
    // operation which may be expensive on embedded targets. For
    // uint32_t and this divisor, we only need 'if' rather than 'while'.
    if (temp >= UINT32_C(0x9068FFFF))
        temp -= UINT32_C(0x9068FFFF);
    if (temp == 0) {
        // Any integer multiple of 0x9068FFFF, including 0, is a bad state.
        // Use an alternate state value by inverting the original value.
        temp = state_orig ^ UINT32_C(0xFFFFFFFF);
        if (temp >= UINT32_C(0x9068FFFF))
            temp -= UINT32_C(0x9068FFFF);
    }
    *seed1 = temp;

    // Sanitize seed2
    // Any integer multiple of 0x464FFFFF, including 0, is a bad state.
    state_orig = *seed2;
    temp = state_orig;

    // The following is equivalent to % 0x464FFFFF, without using modulo
    // operation which may be expensive on some targets. For
    // uint32_t and this divisor, it may loop up to 3 times.
    while (temp >= UINT32_C(0x464FFFFF))
        temp -= UINT32_C(0x464FFFFF);
    if (temp == 0) {
        // Any integer multiple of 0x464FFFFF, including 0, is a bad state.
        // Use an alternate state value by inverting the original value.
        temp = state_orig ^ UINT32_C(0xFFFFFFFF);
        while (temp >= UINT32_C(0x464FFFFF))
            temp -= UINT32_C(0x464FFFFF);
    }

    *seed2 = temp;

    return cur_seed;
}

static void load_seed(void)
{

    static bool single_print = false; // To print only once a message about manual initialization
    seed_type new_seed;
    char conf_file[512];
    FILE *fp;

    // Get the path to the configuration file
    sprintf(conf_file, "%s/.rootsim/numerical.conf", getenv("HOME"));

    // Check if the file exists. If not, we have to create configuration
    if ((fp = fopen(conf_file, "r+")) == NULL) {

        // Try to build the path to the configuration folder.
        sprintf(conf_file, "%s/.rootsim", getenv("HOME"));
        _mkdir(conf_file);

        // Create and initialize the file
        sprintf(conf_file, "%s/.rootsim/numerical.conf", getenv("HOME"));
        if ((fp = fopen(conf_file, "w")) == NULL) {
            rootsim_error(true, "Unable to create the numerical library configuration file %s. Aborting...", conf_file);
        }

        // We now initialize the first long random number. Thanks Unix!
        // TODO: THIS IS NOT PORTABLE!
        int fd;
        if ((fd = open("/dev/random", O_RDONLY)) == -1) {
            rootsim_error(true, "Unable to initialize the numerical library configuration file %s. Aborting...", conf_file);

        }
        read(fd, &new_seed, sizeof(seed_type));
        close(fd);
        fprintf(fp, "%llu\n", (unsigned long long)new_seed);    // We cast, so that we get an integer representing just a bit sequence
        fclose(fp);

    }

    // Is seed manually specified?
    if (rootsim_config.set_seed > 0) {

        if (!single_print) {
            single_print = true;
            printf("Manually setting master seed to %llu\n", (unsigned long long)rootsim_config.set_seed);
        }
        master_seed = rootsim_config.set_seed;
    }

    // Load the configuration for the numerical library
    if ((fp = fopen(conf_file, "r+")) == NULL) {
        rootsim_error(true, "Unable to load numerical distribution configuration: %s. Aborting...", conf_file);
    }

    // Load the initial seed
    fscanf(fp, "%llu", (unsigned long long *)&master_seed);

    // Replace the initial seed
    if (rootsim_config.deterministic_seed == false) {
        rewind(fp);
        srandom(master_seed);
        new_seed = random();
        fprintf(fp, "%llu\n", (unsigned long long)new_seed);
    }

    fclose(fp);
}


// TODO: con un (non tanto) alto numero di processi logici (> numero di bit di un intero!!), lo shift
// circolare restituisce a due LP differenti lo stesso seme iniziale. Questo fa sì che, se la logica di
// programma è la stessa e basata soltanto su numeri casuali, due o più nodi eseguano sempre la stessa
// sequenza di eventi agli stessi timestamp!!!

#define RS_WORD_LENGTH (8 * sizeof(seed_type))
#define ROR(value, places) (value << (places)) | (value >> (RS_WORD_LENGTH - places))   // Circular shift
void numerical_init(void)
{
    // Initialize the master seed
    load_seed();

    // Initialize the per-LP seed
    foreach_lp(lp) {
        lp->numerical.seed = sanitize_seed(ROR((int64_t) master_seed, lp->gid.to_int % RS_WORD_LENGTH));
    }

}

#undef RS_WORD_LENGTH
#undef ROR


__attribute__((const))
double NeumaierSum(unsigned cnt, double addendums[cnt]) {
    if(!cnt)
        return 0.0;
    double sum = addendums[0];
    double crt = 0.0;
    double tmp;
    while (--cnt) {
        tmp = sum + addendums[cnt];
        if(fabs(sum) >= fabs(addendums[cnt])) {
            crt += (sum - tmp) + addendums[cnt];
        } else {
            crt += (addendums[cnt] - tmp) + sum;
        }
        sum = tmp;
    }
    return sum + crt;
}


__attribute__((const))
struct _sum_helper_t PartialNeumaierSum(struct _sum_helper_t sh, double addendum){
    double tmp = sh.sum + addendum;
    if(fabs(sh.sum) >= fabs(addendum)) {
        sh.crt += (sh.sum - tmp) + addendum;
    } else {
        sh.crt += (addendum - tmp) + sh.sum;
    }
    sh.sum = tmp;
    return sh;
}