/* Sequence.c */
/* David Weenink, 19991022 */

#include <time.h>
#include "Sequence.h"

#include "oo_DESTROY.h"
#include "Sequence_def.h"
#include "oo_COPY.h"
#include "Sequence_def.h"
#include "oo_EQUAL.h"
#include "Sequence_def.h"
#include "oo_WRITE_ASCII.h"
#include "Sequence_def.h"
#include "oo_WRITE_BINARY.h"
#include "Sequence_def.h"
#include "oo_READ_ASCII.h"
#include "Sequence_def.h"
#include "oo_READ_BINARY.h"
#include "Sequence_def.h"
#include "oo_DESCRIPTION.h"
#include "Sequence_def.h"

class_methods (Sequence, Data)
	class_method_local (Sequence, destroy)
	class_method_local (Sequence, copy)
	class_method_local (Sequence, equal)
	class_method_local (Sequence, writeAscii)
	class_method_local (Sequence, writeBinary)
	class_method_local (Sequence, readAscii)
	class_method_local (Sequence, readBinary)
	class_method_local (Sequence, description)
class_methods_end

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 2.3e-16
#define RNMX (1.0-EPS)

static double NUMran1 (Sequence me)
{
    long j, k; double ans;

    /* start here when not initializing */
    k = my seed / IQ;
    /* seed = (IA*seed) % IM without overflows by Schrage's method */
    my seed = IA * (my seed - k * IQ) - IR * k;
    if ( my seed < 0) my seed += IM;
    j = 1 + my iy / NDIV;
    my iy = my tab[j];
    my tab[j] = my seed;
    return ((ans = AM * my iy) > RNMX) ? RNMX : ans;
}

int Sequence_init (Sequence me, long n, int action, int seed)
{
    my n = n; my ntab = NTAB; my iy = 0; my tab = NULL;
    if (((my p = NUMlvector (1, n)) == NULL) ||
    	((my tab = NUMlvector (1, my ntab)) == NULL)) return 0;
    my action = action;
    my seed = seed;
    Sequence_restart (me, seed);
    return 1;
}

Sequence Sequence_create (long n, int action, int seed)
{
    Sequence me = new (Sequence);
    if ( ! me || ! Sequence_init (me, n, action, seed)) forget (me);
    return me;
}

void Sequence_restart (Sequence me, long seed)
{
	long j, k;
    if (seed <= 0) seed = time (NULL);
	for (j = my ntab + 8; j > 0; j--)
	{	/* initialize shuffle table of random generator */
		k = seed / IQ;
		seed = IA * (seed - k * IQ) - IR * k;
		if (seed < 0) seed += IM;
		if (j <= my ntab) my tab[j] = seed;
	}
	my iy = my tab[1];
	my seed = seed;
	(void) Sequence_newSequence (me);
}

void Sequence_newSequence (Sequence me)
{
    long i;
    if (my action == SEQUENCE_RANDOMIZE)
    {
    	for (i=1; i <= my n; i++) my p[i] = 1 + my n * NUMran1 (me);
    	return;
    }
    for (i=1; i <= my n; i++) my p[i] = i;
    if (my action == SEQUENCE_PERMUTE)
    {
		for (i=1; i < my n; i++)
		{	/* offset in range [0, my n - i] */
			long tmp, offset = (my n - i + 1) * NUMran1 (me);
			tmp = my p[i]; my p[i] = my p[i+offset]; my p[i+offset] = tmp;
		}
    }
}

long Sequence_ith (Sequence me, long i)
{
    Melder_assert (i > 0 && i <= my n); 
    return my p[i];
}

#undef NTAB
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NDIV
#undef EPS
#undef RNMX

/* End of Sequence.c */
