diff --git a/dspl.project.win.geany b/dspl.project.win.geany index 5840281..985c7f9 100644 --- a/dspl.project.win.geany +++ b/dspl.project.win.geany @@ -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= diff --git a/dspl/dspl_src/inout.c b/dspl/dspl_src/inout.c index 07dbf5b..3db6e75 100644 --- a/dspl/dspl_src/inout.c +++ b/dspl/dspl_src/inout.c @@ -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"); + } diff --git a/dspl/dspl_src/mt19937.c b/dspl/dspl_src/mt19937.c new file mode 100644 index 0000000..f7ef687 --- /dev/null +++ b/dspl/dspl_src/mt19937.c @@ -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 +#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> 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>1) ^ mag01[(int)(x&1ULL)]; + } + for (;i>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); +} diff --git a/dspl/dspl_src/mt19937.h b/dspl/dspl_src/mt19937.h new file mode 100644 index 0000000..cdac0fc --- /dev/null +++ b/dspl/dspl_src/mt19937.h @@ -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 diff --git a/dspl/dspl_src/randgen.c b/dspl/dspl_src/randgen.c index e2dbdaa..2fb764f 100644 --- a/dspl/dspl_src/randgen.c +++ b/dspl/dspl_src/randgen.c @@ -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 -#endif //WIN_OS +#endif /* WIN_OS */ #ifdef LINUX_OS #include -#endif //LINUX_OS +#endif /* LINUX_OS */ #include @@ -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 */ } diff --git a/release/include/dspl.h b/release/include/dspl.h index 8777d6f..7f93dfb 100644 --- a/release/include/dspl.h +++ b/release/include/dspl.h @@ -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 */ diff --git a/verif/src/readbin_verif.c b/verif/src/readbin_verif.c index 621dbb3..4c067ec 100644 --- a/verif/src/readbin_verif.c +++ b/verif/src/readbin_verif.c @@ -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);