dnl dnl File   : rng_h.m4
dnl Author : Richard A. O'Keefe
dnl SCCS   : "%Z%%E% %M%	%I%"
dnl Purpose: Run this through M4 to generate the implementation for a random
dnl          number generator.
dnl 
dnl The random number generator described here is based on the
dnl PragmARC.Universal_Random generator, which is copyright 2000
dnl by PragmAda Software Engineering and released under the GPL.
dnl It can be used for any of the three floating point types
dnl available in C:
dnl	m4 -DT=float		rng_h.m4 >rngf.h
dnl	m4 -DT=double		rng_h.m4 >rng.h
dnl	m4 -DT='long double'	rng_h.m4 >rngl.h
dnl
ifdef(`T',,`define(T,double)')dnl
ifelse(T,float,
    `define(suffix,f)define(SUFFIX,F)define(SFMT, e)define(PFMT,.7e)',
       T,double,
    `define(suffix, )define(SUFFIX, )define(SFMT,le)define(PFMT,.17e)',
       T,long double,
    `define(suffix,l)define(SUFFIX,L)define(SFMT,Le)define(PFMT,.34Le)')dnl
define(TAG, RNG`'SUFFIX)dnl
`#include' <stdlib.h>
`#include' <time.h>
`#include' `"'rng`'suffix.h`"'

void free_`'TAG`'(TAG *r) {
    if (r != 0) free(r);
}

TAG *new_`'TAG`'_from(unsigned long seed) {
    register TAG *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;
    T 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;
}

TAG *new_`'TAG`'(void) {
    return new_`'TAG`'_from((12ul<<24)|(34ul<<16)|(56ul<<8)|(78ul));
}

TAG *new_`'TAG`'_now(void) {
    return new_`'TAG`'_from(
		(unsigned long)time((time_t *)0) ^ (unsigned long)clock() );
}

TAG *read_`'TAG`'(FILE *stream) {
    TAG *r = malloc(sizeof *r);
    int x;

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

int write_`'TAG`'(FILE *stream, TAG const *r) {
    int v, w, x;

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

T urandom`'suffix`'(register TAG *r) {
    register int i = r->i;
    register int j = r->j;
    register T 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;
}

T brandom`'suffix`'(TAG *r, T l, T u) {
    return urandom`'suffix`'(r)*(u-l) + l;
}

int irandom`'suffix`'(TAG *r, int l, int u) {
    return (int)(urandom`'suffix`'(r)*(T)(u-l+1)) + l;
}