mirror
/
libdspl-2.0
kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

added Mersene twister random generator

pull/6/merge
Dsplib 2019-10-20 19:49:09 +03:00
rodzic 7a66faf4e0
commit 89a3a34e42
10 zmienionych plików z 416 dodań i 88 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

62
dspl.project.win.geany
Wyświetl plik

@ -28,7 +28,7 @@ long_line_behaviour=1
long_line_column=72
[files]
current_page=43
current_page=29
FILE_NAME_0=463;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ccheby_poly1_test.c;0;2
FILE_NAME_1=507;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Ciir_test.c;0;2
FILE_NAME_2=80;None;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Cbin%5Cgnuplot%5Ciir_test.plt;0;2
@ -37,42 +37,28 @@ FILE_NAME_4=2672;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_ser
FILE_NAME_5=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Cru%5Ccontent%5Cfourier_series%5Cc%5Cfourier_series_rec.c;0;2
FILE_NAME_6=714;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Carray.c;0;2
FILE_NAME_7=163;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cblas.h;0;2
FILE_NAME_8=56628;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.h;0;2
FILE_NAME_9=14081;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.c;0;2
FILE_NAME_10=323;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Carray_test.c;0;2
FILE_NAME_11=921;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cmatrix.c;0;2
FILE_NAME_12=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_print.c;0;2
FILE_NAME_13=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_transpose.c;0;2
FILE_NAME_14=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_mul.c;0;2
FILE_NAME_15=523;Make;0;EUTF-8;1;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5CMakefile.dspl;0;2
FILE_NAME_16=630;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_eig.c;0;2
FILE_NAME_17=47;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Ccheby.c;0;2
FILE_NAME_18=5878;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Ccomplex.c;0;2
FILE_NAME_19=17821;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cconv.c;0;2
FILE_NAME_20=8338;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cdft.c;0;2
FILE_NAME_21=8117;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cellipj.c;0;2
FILE_NAME_22=9969;None;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_dox%5Cru%5Cellipj.dox;0;2
FILE_NAME_23=909;None;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_dox%5Cru%5Carray.dox;0;2
FILE_NAME_24=1074;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cdspl_internal.h;0;2
FILE_NAME_25=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfft.c;0;2
FILE_NAME_26=11081;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfft_subkernel.c;0;2
FILE_NAME_27=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfillarray.c;0;2
FILE_NAME_28=11043;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_an.c;0;2
FILE_NAME_29=872;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_ap.c;0;2
FILE_NAME_30=4563;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_fir.c;0;2
FILE_NAME_31=888;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_ft.c;0;2
FILE_NAME_32=5440;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_iir.c;0;2
FILE_NAME_33=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfourier_series.c;0;2
FILE_NAME_34=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cgoertzel.c;0;2
FILE_NAME_35=1042;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cinout.c;0;2
FILE_NAME_36=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cmath.c;0;2
FILE_NAME_37=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cpolyval.c;0;2
FILE_NAME_38=3428;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Crandgen.c;0;2
FILE_NAME_39=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cresampling.c;0;2
FILE_NAME_40=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Csignals.c;0;2
FILE_NAME_41=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cstatistic.c;0;2
FILE_NAME_42=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Ctrapint.c;0;2
FILE_NAME_43=9115;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cwin.c;0;2
FILE_NAME_8=14081;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.c;0;2
FILE_NAME_9=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_print.c;0;2
FILE_NAME_10=11081;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfft_subkernel.c;0;2
FILE_NAME_11=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfillarray.c;0;2
FILE_NAME_12=10893;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_an.c;0;2
FILE_NAME_13=872;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_ap.c;0;2
FILE_NAME_14=4563;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_fir.c;0;2
FILE_NAME_15=888;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_ft.c;0;2
FILE_NAME_16=5440;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfilter_iir.c;0;2
FILE_NAME_17=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cfourier_series.c;0;2
FILE_NAME_18=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cgoertzel.c;0;2
FILE_NAME_19=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cmath.c;0;2
FILE_NAME_20=0;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cpolyval.c;0;2
FILE_NAME_21=3739;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Crandgen.c;0;2
FILE_NAME_22=5368;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cmt19937.c;0;2
FILE_NAME_23=48054;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cinclude%5Cdspl.h;0;2
FILE_NAME_24=2887;C++;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cmt19937.h;0;2
FILE_NAME_25=623;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cmatrix_eig.c;0;2
FILE_NAME_26=490;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cverif%5Csrc%5Creadbin_verif.c;0;2
FILE_NAME_27=361;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cfilter_iir_test.c;0;2
FILE_NAME_28=313;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cexamples%5Csrc%5Cdspl_info_test.c;0;2
FILE_NAME_29=829;C;0;EUTF-8;0;1;0;F%3A%5Cdsplib.org%5Clibdspl-2.0%5Cdspl%5Cdspl_src%5Cinout.c;0;2
[build-menu]
NF_00_LB=_Собрать
@ -91,7 +77,7 @@ EX_00_WD=%p/examples/bin
[prjorg]
source_patterns=*.c;*.C;*.cpp;*.cxx;*.c++;*.cc;*.m;
header_patterns=*.h;*.H;*.hpp;*.hxx;*.h++;*.hh;
ignored_dirs_patterns=.*;CVS;
ignored_dirs_patterns=.*;CVS;libblas;liblapack;
ignored_file_patterns=*.o;*.obj;*.a;*.lib;*.so;*.dll;*.lo;*.la;*.class;*.jar;*.pyc;*.mo;*.gmo;
generate_tag_prefs=0
external_dirs=

7
dspl/dspl_src/inout.c
Wyświetl plik

@ -31,11 +31,12 @@ Print DSPL info
void DSPL_API dspl_info()
{
printf("\n\n D S P L - 2.0\n");
printf(" version 2.19.10.13\n");
printf(" version 2.19.10.20\n");
printf("\n Copyright (C) 2015-2019\n");
printf(" Sergey Bakhurin www.dsplib.org\n");
printf(" ---------------------------------------------
printf(" BLAS and LAPACK ver.: 3.8.0 www.netlib.org\n")
printf(" ---------------------------------------------\n");
printf(" BLAS and LAPACK ver.: 3.8.0 www.netlib.org\n");
}

181
dspl/dspl_src/mt19937.c 100644
Wyświetl plik

@ -0,0 +1,181 @@
/*
A C-program for MT19937-64 (2004/9/29 version).
Coded by Takuji Nishimura and Makoto Matsumoto.
This is a 64-bit version of Mersenne Twister pseudorandom number
generator.
Before using, initialize the state by using init_genrand64(seed)
or init_by_array64(init_key, key_length).
Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
References:
T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
ACM Transactions on Modeling and
Computer Simulation 10. (2000) 348--357.
M. Matsumoto and T. Nishimura,
``Mersenne Twister: a 623-dimensionally equidistributed
uniform pseudorandom number generator''
ACM Transactions on Modeling and
Computer Simulation 8. (Jan. 1998) 3--30.
Any feedback is very welcome.
http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html
email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)
*/
#include <stdio.h>
#include "dspl.h"
#include "mt19937.h"
#define NN RAND_MT19937_NN
#define MM 156
#define MATRIX_A 0xB5026F5AA96619E9ULL
#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
#define LM 0x7FFFFFFFULL /* Least significant 31 bits */
/* initializes mt[NN] with a seed */
void mt19937_init_genrand64(unsigned long long seed, random_t* prnd)
{
unsigned long long *mt = prnd->mt19937_mt;
int mti = prnd->mt19937_mti;
mt[0] = seed;
for (mti=1; mti<NN; mti++)
mt[mti] = (6364136223846793005ULL * (mt[mti-1] ^ (mt[mti-1] >> 62)) + mti);
}
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
void mt19937_init_by_array64(unsigned long long init_key[],
unsigned long long key_length,
random_t* prnd)
{
unsigned long long i, j, k;
unsigned long long *mt = prnd->mt19937_mt;
/* int mti = prnd->mt19937_mti; */
mt19937_init_genrand64(19650218ULL, prnd);
i=1; j=0;
k = (NN>key_length ? NN : key_length);
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 3935559000370003845ULL))
+ init_key[j] + j; /* non linear */
i++; j++;
if (i>=NN) { mt[0] = mt[NN-1]; i=1; }
if (j>=key_length) j=0;
}
for (k=NN-1; k; k--) {
/* non linear */
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 62)) * 2862933555777941757ULL)) - i;
i++;
if (i>=NN)
{
mt[0] = mt[NN-1];
i=1;
}
}
mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */
}
/* generates a random number on [0, 2^64-1]-interval */
unsigned long long mt19937_genrand64_int64(random_t* prnd)
{
int i;
unsigned long long x;
static unsigned long long mag01[2]={0ULL, MATRIX_A};
unsigned long long *mt = prnd->mt19937_mt;
int mti = prnd->mt19937_mti;
if (mti >= NN) { /* generate NN words at one time */
/* if init_genrand64() has not been called, */
/* a default initial seed is used */
if (mti == NN+1)
mt19937_init_genrand64(5489ULL, prnd);
for (i=0;i<NN-MM;i++) {
x = (mt[i]&UM)|(mt[i+1]&LM);
mt[i] = mt[i+MM] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
}
for (;i<NN-1;i++) {
x = (mt[i]&UM)|(mt[i+1]&LM);
mt[i] = mt[i+(MM-NN)] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
}
x = (mt[NN-1]&UM)|(mt[0]&LM);
mt[NN-1] = mt[MM-1] ^ (x>>1) ^ mag01[(int)(x&1ULL)];
mti = 0;
}
x = mt[mti++];
x ^= (x >> 29) & 0x5555555555555555ULL;
x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
x ^= (x << 37) & 0xFFF7EEE000000000ULL;
x ^= (x >> 43);
return x;
}
/* generates a random number on [0, 2^63-1]-interval */
long long mt19937_genrand64_int63(random_t* prnd)
{
return (long long)(mt19937_genrand64_int64(prnd) >> 1);
}
/* generates a random number on [0,1]-real-interval */
double mt19937_genrand64_real1(random_t* prnd)
{
return (mt19937_genrand64_int64(prnd) >> 11) * (1.0/9007199254740991.0);
}
/* generates a random number on [0,1)-real-interval */
double mt19937_genrand64_real2(random_t* prnd)
{
return (mt19937_genrand64_int64(prnd) >> 11) * (1.0/9007199254740992.0);
}
/* generates a random number on (0,1)-real-interval */
double mt19937_init_genrand64_real3(random_t* prnd)
{
return ((mt19937_genrand64_int64(prnd) >> 12) + 0.5) * (1.0/4503599627370496.0);
}

86
dspl/dspl_src/mt19937.h 100644
Wyświetl plik

@ -0,0 +1,86 @@
/*
A C-program for MT19937-64 (2004/9/29 version).
Coded by Takuji Nishimura and Makoto Matsumoto.
This is a 64-bit version of Mersenne Twister pseudorandom number
generator.
Before using, initialize the state by using init_genrand64(seed)
or init_by_array64(init_key, key_length).
Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. The names of its contributors may not be used to endorse or promote
products derived from this software without specific prior written
permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
References:
T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
ACM Transactions on Modeling and
Computer Simulation 10. (2000) 348--357.
M. Matsumoto and T. Nishimura,
``Mersenne Twister: a 623-dimensionally equidistributed
uniform pseudorandom number generator''
ACM Transactions on Modeling and
Computer Simulation 8. (Jan. 1998) 3--30.
Any feedback is very welcome.
http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html
email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)
*/
#ifndef MT19937_H
#define MT19937_H
#include "dspl.h"
/* initializes mt[NN] with a seed */
void mt19937_init_genrand64(unsigned long long seed, random_t* prnd);
/* initialize by an array with array-length */
/* init_key is the array for initializing keys */
/* key_length is its length */
void mt19937_init_by_array64(unsigned long long init_key[],
unsigned long long key_length,
random_t* prnd);
/* generates a random number on [0, 2^64-1]-interval */
unsigned long long mt19937_genrand64_int64(random_t* prnd);
/* generates a random number on [0, 2^63-1]-interval */
long long mt19937_genrand64_int63(random_t* prnd);
/* generates a random number on [0,1]-real-interval */
double mt19937_genrand64_real1(random_t* prnd);
/* generates a random number on [0,1)-real-interval */
double mt19937_genrand64_real2(random_t* prnd);
/* generates a random number on (0,1)-real-interval */
double mt19937_genrand64_real3(random_t* prnd);
#endif

98
dspl/dspl_src/randgen.c
Wyświetl plik

@ -25,20 +25,40 @@
#include "dspl.h"
#include "dspl_internal.h"
#include "mt19937.h"
void DSPL_API random_init(random_t* prnd)
/******************************************************************************
random generator initialization
*******************************************************************************/
int DSPL_API random_init(random_t* prnd, int type)
{
srand(time(NULL));
/* MRG32k3a init */
prnd->mrg32k3a_x[0] = prnd->mrg32k3a_x[1] = 1.0;
prnd->mrg32k3a_y[0] = prnd->mrg32k3a_y[1] = prnd->mrg32k3a_y[2] = 1.0;
prnd->mrg32k3a_x[2] = rand();
switch(type)
{
case RAND_TYPE_MRG32K3A:
/* MRG32k3a init */
prnd->mrg32k3a_x[0] = prnd->mrg32k3a_x[1] = 1.0;
prnd->mrg32k3a_y[0] = prnd->mrg32k3a_y[1] = prnd->mrg32k3a_y[2] = 1.0;
prnd->mrg32k3a_x[2] = rand();
break;
case RAND_TYPE_MT19937:
mt19937_init_genrand64((unsigned long long)rand()*rand(), prnd);
break;
default:
return ERROR_RAND_TYPE;
}
prnd->type = type;
return RES_OK;
}
/******************************************************************************
Uniform random generator mrg32k3a
*******************************************************************************/
int randu_mrg32k3a (double* u, int n, random_t* prnd)
{
@ -88,30 +108,51 @@ int randu_mrg32k3a (double* u, int n, random_t* prnd)
/******************************************************************************
Uniform random numbers generator
*******************************************************************************/
int DSPL_API randu(double* x, int n, random_t* prnd)
{
int i;
if(!x || !prnd)
return ERROR_PTR;
if(n < 0)
return ERROR_SIZE;
if(prnd)
return randu_mrg32k3a(x, n, prnd);
{
switch(prnd->type)
{
case RAND_TYPE_MRG32K3A:
return randu_mrg32k3a(x, n, prnd);
case RAND_TYPE_MT19937:
for(i = 0; i < n; i++)
x[i] = mt19937_genrand64_real1(prnd);
return RES_OK;
default:
return ERROR_RAND_TYPE;
}
}
else
{
if(!x)
return ERROR_PTR;
if(n<1)
return ERROR_SIZE;
int i;
for(i = 0; i < n; i++)
x[n] = (double)rand()/RAND_MAX;
}
return RES_OK;
}
/*******************************************************************************
Gaussian random numbers generator
*******************************************************************************/
@ -130,30 +171,37 @@ int DSPL_API randn(double* x, int n, double mu, double sigma, random_t* prnd)
return ERROR_RAND_SIGMA;
k=0;
while(k < n)
while(k+128 < n)
{
res = randu(x1, 128, prnd);
if(res != RES_OK)
if((res = randu(x1, 128, prnd)) != RES_OK)
goto exit_label;
res = randu(x2, 128, prnd);
if(res != RES_OK)
if((res = randu(x2, 128, prnd)) != RES_OK)
goto exit_label;
m = 0 ;
while(k<n && m < 128)
for(m = 0; m < 128; m++)
{
x[k] = sqrt(-2.0*log(x1[m]))*cos(M_2PI*x2[m])*sigma + mu;
k++;
m++;
if(k<n && m < 128)
if(x1[m] != 0.0)
{
x[k] = sqrt(-2.0*log(x1[m]))*sin(M_2PI*x2[m])*sigma + mu;
x[k] = sqrt(-2.0*log(x1[m]))*cos(M_2PI*x2[m])*sigma + mu;
k++;
m++;
}
}
}
while(k < n)
{
if((res = randu(x1, 1, prnd)) != RES_OK)
goto exit_label;
if((res = randu(x2, 1, prnd)) != RES_OK)
goto exit_label;
if(x1[1] != 0.0)
{
x[k] = sqrt(-2.0*log(x1[1]))*cos(M_2PI*x2[1])*sigma + mu;
k++;
}
}
res = RES_OK;
exit_label:
return res;

4
examples/src/filter_iir_test.c
Wyświetl plik

@ -14,7 +14,7 @@ int main()
random_t rnd;
int k;
random_init(&rnd); // random generator init
random_init(&rnd, RAND_TYPE_MT19937); // random generator init
linspace(0, N, N, DSPL_PERIODIC, t); // fill t vector
randn(n, N, 0, 1.0, &rnd); // generate noise
@ -36,4 +36,4 @@ int main()
// run GNUPLOT script
return system("gnuplot -p gnuplot/filter_iir.plt");
}
}

15
include/dspl.h
Wyświetl plik

@ -52,12 +52,23 @@ typedef struct
} fft_t;
#define RAND_TYPE_MRG32K3A 0x00000001
#define RAND_TYPE_MT19937 0x00000002
#define RAND_MT19937_NN 312
typedef struct
{
double mrg32k3a_seed;
double mrg32k3a_x[3];
double mrg32k3a_y[3];
/* The array for the MT19937 state vector */
unsigned long long mt19937_mt[RAND_MT19937_NN];
int mt19937_mti;
int type;
}random_t;
@ -129,6 +140,7 @@ typedef struct
/* Q 0x17xxxxxx*/
/* R 0x18xxxxxx*/
#define ERROR_RAND_SIGMA 0x18011909
#define ERROR_RAND_TYPE 0x18012009
#define ERROR_RESAMPLE_RATIO 0x18051801
#define ERROR_RESAMPLE_FRAC_DELAY 0x18050604
/* S 0x19xxxxxx*/
@ -781,7 +793,8 @@ DECLARE_FUNC(int, randn, double*
COMMA double
COMMA random_t* prnd);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(void, random_init, random_t* prnd);
DECLARE_FUNC(int, random_init, random_t* prnd
COMMA int type);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, randu, double*
COMMA int

20
release/include/dspl.c
Wyświetl plik

@ -22,11 +22,11 @@
#ifdef WIN_OS
#include <windows.h>
#endif //WIN_OS
#endif /* WIN_OS */
#ifdef LINUX_OS
#include <dlfcn.h>
#endif //LINUX_OS
#endif /* LINUX_OS */
#include <stdio.h>
@ -157,7 +157,7 @@ p_writetxt_cmplx_im writetxt_cmplx_im ;
p_writetxt_cmplx_re writetxt_cmplx_re ;
#endif //BUILD_LIB
#endif /* BUILD_LIB */
@ -192,20 +192,20 @@ void* dspl_load()
printf("libdspl.dll loading ERROR!\n");
return NULL;
}
#endif //WIN_OS
#endif /* WIN_OS */
#ifdef LINUX_OS
char* error;
void *handle;
// open the *.so
/* open the *.so */
handle = dlopen ("./libdspl.so", RTLD_LAZY);
if (!handle)
{
printf("libdspl.so loading ERROR!\n");
return NULL;
}
#endif //LINUX_OS
#endif /* LINUX_OS */
LOAD_FUNC(acos_cmplx);
LOAD_FUNC(asin_cmplx);
@ -332,7 +332,7 @@ void* dspl_load()
if(handle)
FreeLibrary(handle);
return NULL;
#endif //WIN_OS
#endif /* WIN_OS */
#ifdef LINUX_OS
@ -342,7 +342,7 @@ void* dspl_load()
if(handle)
dlclose(handle);
return NULL;
#endif //LINUX_OS
#endif /* LINUX_OS */
}
@ -355,11 +355,11 @@ void dspl_free(void* handle)
{
#ifdef WIN_OS
FreeLibrary((HINSTANCE)handle);
#endif //WIN_OS
#endif /* WIN_OS */
#ifdef LINUX_OS
dlclose(handle);
#endif //LINUX_OS
#endif /* LINUX_OS */
}

27
release/include/dspl.h
Wyświetl plik

@ -52,12 +52,23 @@ typedef struct
} fft_t;
#define RAND_TYPE_MRG32K3A 0x00000001
#define RAND_TYPE_MT19937 0x00000002
#define RAND_MT19937_NN 312
typedef struct
{
double mrg32k3a_seed;
double mrg32k3a_x[3];
double mrg32k3a_y[3];
/* The array for the MT19937 state vector */
unsigned long long mt19937_mt[RAND_MT19937_NN];
int mt19937_mti;
int type;
}random_t;
@ -129,6 +140,7 @@ typedef struct
/* Q 0x17xxxxxx*/
/* R 0x18xxxxxx*/
#define ERROR_RAND_SIGMA 0x18011909
#define ERROR_RAND_TYPE 0x18012009
#define ERROR_RESAMPLE_RATIO 0x18051801
#define ERROR_RESAMPLE_FRAC_DELAY 0x18050604
/* S 0x19xxxxxx*/
@ -215,15 +227,15 @@ typedef struct
#ifdef BUILD_LIB
// Declare DSPL_API for Windows OS
/* Declare DSPL_API for Windows OS */
#ifdef WIN_OS
#define DSPL_API __declspec(dllexport)
#endif // WIN_OS
// Declare DSPL_API for LINUX OS
#endif /* WIN_OS */
/* Declare DSPL_API for LINUX OS */
#ifdef LINUX_OS
#define DSPL_API
#endif //LINUX_OS
#endif //BUILD_DLL
#endif /* LINUX_OS */
#endif /* BUILD_DLL */
#define COMMA ,
@ -781,7 +793,8 @@ DECLARE_FUNC(int, randn, double*
COMMA double
COMMA random_t* prnd);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(void, random_init, random_t* prnd);
DECLARE_FUNC(int, random_init, random_t* prnd
COMMA int type);
/*----------------------------------------------------------------------------*/
DECLARE_FUNC(int, randu, double*
COMMA int
@ -918,5 +931,5 @@ void dspl_free(void* handle);
#endif //DSPL_H
#endif /* DSPL_H */

4
verif/src/readbin_verif.c
Wyświetl plik

@ -18,9 +18,9 @@ int main()
xc = (complex_t*) malloc(N * sizeof(complex_t));
random_t rnd;
random_init(&rnd);
random_init(&rnd, RAND_TYPE_MT19937);
randn(xr, N, 0, 1.0,&rnd);
randn(xr, N, 0, 1.0, &rnd);
randn((double*)xc, 2*N, 0, 1.0,&rnd);