randomWH.c

Go to the documentation of this file.

#include "randomWH.h"


Congruence congruent={
    .x=0,
    .y=0,
    .z=0,
    .t=0,
    .initialised=isfalse,
    .init=Congruence_init,
    .randomWH=randomWH_null,
    .read=Congruence_read,
    .bernoulli=randomWH_bernoulli,
    .geometric=randomWH_geometric,
    .binomial=randomWH_binomial,
    .fairdie=randomWH_fairdie,
}; 
void Congruence_init(Congruence* self, long int x,long int y,long int z,long int t){
    if(x<=0||y<=0||z<=0||t<=0)
        randomWH_error("Seeds must be positive!");
    self->x = x;
    self->y = y;
    self->z = z;
    self->t = t;
    if (self->initialised)
        return;
    self->initialised = istrue;
    switch(sizeof(long int)/8){
        case 1:
                self->randomWH = randomWH_64;
                break;
        case 0:
                self->randomWH = randomWH_32;
                break;
        default:
                randomWH_error("Could not initialise randomWH function pointer");
    }
}

void randomWH_error(const char* message){
    char prefix[] = "Error in randomWH.c: ";
    fprintf(stderr,"%s %s\n",prefix, message);
    exit(EXIT_FAILURE);
}

long int* Congruence_read(Congruence* self){
    return (long int*)&self->x;
}


double randomWH_null(Congruence* self){
        fprintf(stderr,"Could not initialise randomWH function pointer\n");
        exit(EXIT_FAILURE);
}


double randomWH_32(Congruence* self){
    double result;
    self->x = 11600L * (self->x % 185127L) - 10379L* (self->x/185127L);
    self->y = 47003L * (self->y %  45688L) - 10479L* (self->y/ 45688L);
    self->z = 23000L * (self->z %  93368L) - 19423L* (self->z/ 93368L);
    self->t = 33000L * (self->t %  65075L) -  8123L* (self->t/ 65075L);
    if (self->x < 0)
       self->x += 2147483579L;
    if (self->y < 0)
       self->y += 2147483543L;
    if (self->z < 0)
       self->z += 2147483423L;
    if (self->t < 0)
       self->t += 2147483123L;
    result =  (double)(self->x)*0.0000000004656613022697297188506231646486
       + (double)(self->y)*0.0000000004656613100759859932486569933169
       + (double)(self->z)*0.0000000004656613360968421314794009471615
       + (double)(self->t)*0.0000000004656614011489951998100056779817;
    while (result > 1.0)
       result -= 1.0;
    return result;

}

double randomWH_64(Congruence* self){
    double result;
    self->x = (11600L * self->x) % 2147483579L;
    self->y = (47003L * self->y) %  2147483543L;
    self->z = (23000L * self->z) %  2147483423L;
    self->t = (33000L * self->t) %  2147483123L;
    result =  (double)(self->x)*0.0000000004656613022697297188506231646486
       + (double)(self->y)*0.0000000004656613100759859932486569933169
       + (double)(self->z)*0.0000000004656613360968421314794009471615
       + (double)(self->t)*0.0000000004656614011489951998100056779817;
    while(result > 1.0)
       result -= 1.0;
    return result;

}

long int randomWH_geometric(Congruence* self, double prob){
    double unif = self->randomWH(self);
    double p_bar = 1-prob;
    long int result = 1L;
    if(p_bar > 0){
        result += log(unif)/log(p_bar);
    }
    return result;
}

enum BOOLEAN randomWH_bernoulli(Congruence* self, double prob){
    double unif = self->randomWH(self);
    enum BOOLEAN result = isfalse;

    if(unif < prob){
        result = istrue;
    }
    return result;
}

long int randomWH_binomial(Congruence* self, long int trials, double prob){
    long int result = 0;
    long int count;

    for(count=0;count<trials;count++)
        result += self->bernoulli(self,prob);
    return result;
}

long int randomWH_fairdie(Congruence* self, long int lower_inclusive, long int upper_inclusive){
    long int result;
    double randvar = self->randomWH(self);
    if (randvar == 1.0)
        return upper_inclusive;
    randvar = lower_inclusive + randvar * (upper_inclusive - lower_inclusive +1);
    result = randvar; /* compiler will assign result <- floor(randvar) */
    return result;
}