---

SimSeqs.h

Go to the documentation of this file.

/***************************************************************************
                                                SimSeqs.h  -  description
                                                -----------------
                begin                           : Sun May 29 2000
                copyright                       : (C) 2001 by Reinhold Kainhofer
                email                           : reinhold@kainhofer.com
 ***************************************************************************/

/***************************************************************************
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 ***************************************************************************/


#ifndef SIMSEQS_H
#define SIMSEQS_H
//#include <stdlib.h>
//#include <stdio.h>
//#include <math.h>
#include "BaseDefinitions.h"

#define USEMC 1
#define USEHALTON 2
#define USENET 4
#define USENIEDER 8
#define USESOBOL 16
#define USEFAURE 32
#define USEATANASSOV 64
#define USEHAMMERSLEY 128
#define USENALPHA 256
#define USEMAX USENALPHA


        // constants for the recursive construction of the Halton sequence
#define GENAU           1E-9
#define BGENAU   29
#define MAXLONG 2147483647
#define MAXPREC 30
#define primefile "primes.dat"

#define NAlphaOnes .6180339887498948482045969343656391177203091798058
#define NAlphaTwos .4142135623730950488016887242096980785696718753769
#define NAlphaThrees .3027756377319946465596106337352479731256482869226
#define NAlphaFours .2360679774997896964091736687312762354406183596115

/* *********************************************************** */
/* Konstanten, Typen                                                                                                                                                             */
/* *********************************************************** */
//TODO
#define maxanz 310
#define max_nieder 900000000 // Max. Anz.der Niederr. Vektoren
#define maxdim 310      // max. dimension for the sequence. set this to a high value to allow 400-dim. sequences as we need them.
#define MAX_I_NORMAL 1025

typedef double double_vektor[maxanz];
typedef int             int_vektor[maxanz];
typedef struct
                                {
                                        int                             deg;                                            // degree of the polynomial
                                        int_vektor coeff;                                       // coefficients
                                                                                                                                                 // p(x)= coeff[0]+coeff[1]*x+ ...
                                } polynom;
typedef polynom polarr[maxdim+4];
//endTODO */

class longpoly {
protected:
        long deg;
        long base;
public:
        longpoly(long cflen=15, long bs=0);
        ~longpoly();
        long GetDegree();
        long SetDegree(long deg);
        long SetCoeff(long k, long val);
        long GetCoeff(long k);

        longpoly* times(longpoly& poly1);
//      longpoly* operator *(longpoly& poly1);

protected:
        long*coeff;
        long coefflen;
};


//c lass LDSqBase;
long ReadPrimes(long* primes, long m);
long NextPrime(long current);
long*InitGenericPrimes(long dm);
long*InitGenericInts(long dm);




#define NTAB 32

class LDSqBase {
protected:
        char* name, *ext;
        long dim;
        long len;
        long* bases;
        double* numbers;
        bool created;
                        /* needed for the Pseudo random generators. Stores state and some other coefficients for shuffling
                        */
        long state, iy, iv[NTAB], type;
public:
        long lastnr;

public:
        /*const*/ bool qmc;
        LDSqBase(long*b=NULL, long dm=0,long iterations=0,
                        long genau=0, double genau1=0., char* ex="",char* nm="");
        virtual ~LDSqBase();

        virtual long GetDimension() {return dim;};
        virtual long GetLength() {return len;};
        virtual char* GetName() {return name;};
        virtual char* GetExtension() {return ext;};
        virtual long SetBases(long*b, long dm, long genau, double genau1);
        virtual long GetBases(long*b, long dm) {
                        for (long i=0;i<min(dm, dim);i++) b[i]=bases[i]; return min(dm, dim);}
        virtual long NextElement(double*buffer, long bufflen, long nr=-1);
        virtual bool CreateNumbers();
        virtual void SetNames(char* ex="",char* nm="");
        virtual void SetSeed(long int status) {state=status;}
        double ran1(long*idum);
        double ran2(long*idum);
        double ran3(long*idum);

protected:
        virtual void InitMethod(long*b=NULL, long dm=0,long iterations=0,
                        long genau=0, double genau1=0., char* ex="",char* nm="");
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau, double genau1);
        virtual long ExitData() { return 0;}
        virtual long*InitGenericBases(long dm) {return NULL;};
};



#define MC1 1   // "minimal" random number generator of Park and Miller with Bays-Durham shuffle and added safeguards
#define MC2 2   // long period rng of L'Ecuyer with Bays-Durham shuffle and added safeguards
#define MC3 3   // Knuth's subtractive method rng

class LDSqMonteCarlo : public LDSqBase  {
public:
        LDSqMonteCarlo(long tp=MC1, long status=-44891368, char* ex="mc",char* nm="Monte Carlo");
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen) {
                switch (type) {
                        case MC2:       for (long i=0;i<bufflen;i++) { buffer[i]=ran2(&state);}; break;
                        case MC3:       for (long i=0;i<bufflen;i++) { buffer[i]=ran3(&state);}; break;
                        case MC1:default: for (long i=0;i<bufflen;i++) { buffer[i]=ran1(&state);}; break;
                }
                return dim;
        }
        virtual long InitData(long genau, double genau1){return dim;};
        virtual long ExitData(){return 0;};
};


class LDSqAfflerbach : public LDSqMonteCarlo {
        long altx;
public:
        LDSqAfflerbach(char* ex="aff", char*nm="Afflerbach"):LDSqMonteCarlo(MC1,0,ex,nm) {altx=12785;}
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen) {
                for (long i=0;i<bufflen;i++) buffer[i]=createNumber();
                return dim;
        }
private:
        virtual double createNumber(){
                long xx=(532393*altx+1);
                while (xx>=268435456) xx-=268435456;
                altx=xx;
                return ((double)xx/(double)268435456);
        }
};


class LDSqNAlpha : public LDSqBase {
public:
        LDSqNAlpha(Zahl*alph=NULL, long dm=0, long iterations=0, char*ex="nalpha", char*nm="[n Alpha]"):
                        LDSqBase(NULL, dm, iterations, 0,0,ex,nm) {
                        alpha=alph;
                        InitMethod(NULL, dm, iterations, 0,0,ex,nm); }
protected:
        virtual long*InitGenericBases(long dm) {return NULL;/*TODO!!!*/};
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen) {
                long d=min(bufflen,dim);
                for (long i=0; i<d;i++) buffer[i]=(nr*alpha[i])-floor((nr*alpha[i]));
                return d; }
        Zahl*alpha;
};

class LDSqHalton : public LDSqBase      {
public:
        LDSqHalton(long*b=NULL, long dm=0,long iterations=0, long genau=BGENAU,
                        double genau1=GENAU, char* ex="hal",char* nm="Halton"):
                        LDSqBase(b, dm, iterations, genau, genau1, ex, nm){
                        p=NULL; q=NULL; lastnr=0; InitMethod(b, dm, iterations, genau, genau1,  ex, nm);}
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau, double genau1);
        virtual long ExitData();
        virtual long*InitGenericBases(long dm) {return InitGenericPrimes(dm);};
        long prec;
        double *q,*p;
};



class LDSqAtanassov : public LDSqBase   {
protected: // Protected attributes
        long* factors;
        long* proots;
        long* modifiers;
public:
        LDSqAtanassov(long*b=NULL, long dm=NULL,long iterations=0, long genau=BGENAU,
                        double genau1=GENAU, char* ex="ata",char* nm="Atanassov", long*prts=NULL,long*mdfs=NULL):
                        LDSqBase(b, dm, iterations, genau, genau1, ex, nm){
                        proots=prts; modifiers=mdfs; InitMethod(b, dm, iterations, genau, genau1,       ex, nm);}
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau=0, double genau1=0);
        virtual long ExitData();
        virtual long*InitGenericBases(long dm) {return InitGenericPrimes(dm);};
private:
        bool GeneratePrimitiveRoots(long*bases, long*prts, long dim);
        bool GenerateModifiers(long*bases, long*prts, long*modfs, long dim);
};

class LDSqHammersley : public LDSqHalton        {
protected:
        int Npos;
public:
        LDSqHammersley(long*b=NULL, long dm=0,long iterations=0, long genau=0,
                        double genau1=0., int npos=1, char* ex="ham",char* nm="Hammersley"):
                        LDSqHalton(b,dm-1, iterations, genau, genau1,ex,nm) {
                        Npos=((npos%dm)+dm)%dm;dim++;};
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
};


#define konst 1073741824

class LDSqSobol : public LDSqBase       {
protected:
        long*V;
public:
        LDSqSobol(long*b=NULL, long dm=0,long iterations=0, long genau=0,
                        double genau1=0., char* ex="sob",char* nm="Sobol"):
                        LDSqBase(b, dm, iterations, genau, genau1, ex,nm){
                        V=NULL;InitMethod(b, dm, iterations, genau, genau1,ex,nm);}
                                // write the next element(s) of the sequence in a buffer.
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau, double genau1);
        virtual long ExitData();
private:
        void dirnum_erz();
};



class LDSqFaure : public LDSqBase       {
protected:
        long r, base;
        long*c;
public:
        // fstep contains base^fstp
        long fstep;
public:
        LDSqFaure(long*b=NULL, long dm=0,long iterations=0, long genau=0,
                        double genau1=0., long fstp=0, char* ex="fau",char* nm="Faure");
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau, double genau1);
        virtual long ExitData();
private:
        void cyMultiply(long*cc, long*yold, long*ynew, long dm, long coefflen, long bs);
protected:
        virtual void CreateDirMx(long*cc, long len, long bs);
};

class LDSqNetz : public LDSqBase        {
protected:
        long base, nu, lambda;
        long nprime, ncurr, ncmax;
        long*ncoeff, *ypji;
        long ypjilen;
public:
        LDSqNetz(long*b=NULL, long bas=5, long dm=0,long iterations=0,
                        long lamb=5, char* ex="ts" ,char* nm="(t,s)-Sequence");
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau, double genau1);
        virtual long ExitData();
        virtual long*InitGenericBases(long dm) {return InitGenericInts(dm);};
private:
        long getypji(long p, long j, long i);
        long getypji(long index);
        long setypji(long p, long j, long i, long value);
        long setypji(long index, long value);
        long getindex(long p, long j, long i);
};

class LDSqNiederreiter : public LDSqBase        {
protected:
        long basis;
        long maxnum,log_N;
        int_vektor q;
        long max_num;
        long c[32][32][maxanz]; //      c_ij acc. B-F-N Niederr
        polarr monic_pols;
//      longpoly**MonicPolys;
//      long nu, lambda;
//      long nprime, ncurr, ncmax;
//      long*ncoeff, *ypji;
//      long ypjilen;
public:
        LDSqNiederreiter(long*b=NULL, long bas=5, long dm=0,long iterations=0,
                        char* ex="nie" ,char* nm="Niederreiter (t,s)-Sequence");
protected:
        virtual long CalculateNextElement(long nr, double*buffer, long bufflen);
        virtual long InitData(long genau, double genau1);
        virtual long ExitData();
private:
        void pol_mult(polynom ap, polynom bp, polynom* cp);
        int entw_in_basis (int n, int *a);
        void gen_poly(polynom* pol );
};



#endif

---