#include <stdlib.h>
#include <time.h>
#include "rngf.h"

void free_RNGF(RNGF *r) {
    if (r != 0) free(r);
}

RNGF *new_RNGF_from(unsigned long seed) {
    register RNGF *r;
    int i = ((int)((seed >> 24) & 255u) + 175)%177 + 1;
    int j = ((int)((seed >> 16) & 255u) + 175)%177 + 1;
    int k = ((int)((seed >>  8) & 255u) + 175)%177 + 1;
    int l = ((int)((seed >>  0) & 255u)      )%169;
    int x, y;
    float s, t;

    r = malloc(sizeof *r);
    if (r == 0) return r;
    r->i = 96;
    r->j = 32;
    r->c = 362436.0/16777216.0;	/* exact even in float precision */
    for (x = 0; x < 97; x++) {
	s = 0.0;
	t = 0.5;
	for (y = 0; y < 24; y++) {
	    int m = (((j*i)%179)*k)%179;
	    i = j, j = k, k = m, l = (l*53+1)%169;
	    if ((l*m)%64 >= 32) s += t;
	    t = 0.5*t;
	}
	r->u[x] = s;
    }
    return r;
}

RNGF *new_RNGF(void) {
    return new_RNGF_from((12ul<<24)|(34ul<<16)|(56ul<<8)|(78ul));
}

RNGF *new_RNGF_now(void) {
    return new_RNGF_from(
		(unsigned long)time((time_t *)0) ^ (unsigned long)clock() );
}

RNGF *read_RNGF(FILE *stream) {
    RNGF *r = malloc(sizeof *r);
    int x;

    if (r == 0) goto failed;
    if (3 != fscanf(stream, "%d%d%e", &r->i, &r->j, &r->c)) goto failed;
    for (x = 0; x < 97; x++)
	if (1 != fscanf(stream, "%e", &r->u[x])) goto failed;
    return r;
failed:
    if (r != 0) free(r);
    return 0;
}

int write_RNGF(FILE *stream, RNGF const *r) {
    int v, w, x;

    v = fprintf(stream, "%d %d %.7e", r->i, r->j, r->c);
    if (v < 0) return v;
    for (x = 0; x < 97; x++) {
	w = fprintf(stream, " %.7e", r->u[x]);
	if (w < 0) return w;
	v += w;
    }
    w = fprintf(stream, "\n");
    if (w < 0) return w;
    return v+w;
}

float urandomf(register RNGF *r) {
    register int i = r->i;
    register int j = r->j;
    register float c, t;

    t = r->u[i] - r->u[j];
    if (t < 0.0) t += 1.0;
    r->u[i] = t;
    r->i = i == 0 ? 96 : i-1;
    r->j = j == 0 ? 96 : j-1;
    c = r->c - 7654321.0/16777216.0;
    if (c < 0.0) c += 16777213.0/16777216.0;
    r->c = c;
    t -= c;
    if (t < 0.0) t += 1;
    return t;
}

float brandomf(RNGF *r, float l, float u) {
    return urandomf(r)*(u-l) + l;
}

int irandomf(RNGF *r, int l, int u) {
    return (int)(urandomf(r)*(float)(u-l+1)) + l;
}