Revision: 75341
http://sourceforge.net/p/brlcad/code/75341
Author: starseeker
Date: 2020-04-10 19:51:28 +0000 (Fri, 10 Apr 2020)
Log Message:
-----------
Make a stab at swapping in the newer BSD licensed version of mt19937ar. This
isn't the newer faster one, but the version from
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c
Modified Paths:
--------------
brlcad/trunk/src/libbn/CMakeLists.txt
Added Paths:
-----------
brlcad/trunk/src/libbn/mt19937ar.c
Removed Paths:
-------------
brlcad/trunk/src/libbn/randmt.c
Modified: brlcad/trunk/src/libbn/CMakeLists.txt
===================================================================
--- brlcad/trunk/src/libbn/CMakeLists.txt 2020-04-10 18:29:58 UTC (rev
75340)
+++ brlcad/trunk/src/libbn/CMakeLists.txt 2020-04-10 19:51:28 UTC (rev
75341)
@@ -22,6 +22,7 @@
marker.c
mat.c
msr.c
+ mt19937ar.c
multipoly.c
noise.c
plane.cpp
@@ -29,7 +30,6 @@
poly.c
qmath.c
rand.c
- randmt.c
randsph.c
scale.c
sobolseq.c
Added: brlcad/trunk/src/libbn/mt19937ar.c
===================================================================
--- brlcad/trunk/src/libbn/mt19937ar.c (rev 0)
+++ brlcad/trunk/src/libbn/mt19937ar.c 2020-04-10 19:51:28 UTC (rev 75341)
@@ -0,0 +1,215 @@
+/*
+ A C-program for MT19937, with initialization improved 2002/1/26.
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+
+ Before using, initialize the state by using bn_randmt_seed(seed)
+ or init_by_array(init_key, key_length).
+
+ Copyright (C) 1997 - 2002, 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.
+
+
+ Any feedback is very welcome.
+ http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+ email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
+*/
+
+#include "common.h"
+#include <stdio.h>
+
+/* Period parameters */
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL /* constant vector a */
+#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
+#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
+
+static unsigned long mt[N]; /* the array for the state vector */
+static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
+
+/* initializes mt[N] with a seed */
+void bn_randmt_seed(unsigned long s)
+{
+ mt[0]= s & 0xffffffffUL;
+ for (mti=1; mti<N; mti++) {
+ mt[mti] =
+ (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
+ /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+ /* In the previous versions, MSBs of the seed affect */
+ /* only MSBs of the array mt[]. */
+ /* 2002/01/09 modified by Makoto Matsumoto */
+ mt[mti] &= 0xffffffffUL;
+ /* for >32 bit machines */
+ }
+}
+
+
+/* generates a random number on [0,0xffffffff]-interval */
+static unsigned long
+genrand_int32(void)
+{
+ unsigned long y;
+ static unsigned long mag01[2]= {0x0UL, MATRIX_A};
+ /* mag01[x] = x * MATRIX_A for x=0,1 */
+
+ if (mti >= N) { /* generate N words at one time */
+ int kk;
+
+ if (mti == N+1) /* if bn_randmt_seed() has not been called, */
+ bn_randmt_seed(5489UL); /* a default initial seed is used */
+
+ for (kk=0; kk<N-M; kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
+ }
+ for (; kk<N-1; kk++) {
+ y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
+ mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
+ }
+ y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
+ mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
+
+ mti = 0;
+ }
+
+ y = mt[mti++];
+
+ /* Tempering */
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9d2c5680UL;
+ y ^= (y << 15) & 0xefc60000UL;
+ y ^= (y >> 18);
+
+ return y;
+}
+
+/* generates a random number on [0,1]-real-interval */
+double bn_randmt(void)
+{
+ return genrand_int32()*(1.0/4294967295.0);
+ /* divided by 2^32-1 */
+}
+
+#if 0
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+/* slight change for C++, 2004/2/26 */
+void init_by_array(unsigned long init_key[], int key_length)
+{
+ int i, j, k;
+ bn_randmt_seed(19650218UL);
+ i=1;
+ j=0;
+ k = (N>key_length ? N : key_length);
+ for (; k; k--) {
+ mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ + init_key[j] + j; /* non linear */
+ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++;
+ j++;
+ if (i>=N) {
+ mt[0] = mt[N-1];
+ i=1;
+ }
+ if (j>=key_length) j=0;
+ }
+ for (k=N-1; k; k--) {
+ mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
+ - i; /* non linear */
+ mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++;
+ if (i>=N) {
+ mt[0] = mt[N-1];
+ i=1;
+ }
+ }
+
+ mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
+}
+
+/* generates a random number on [0,0x7fffffff]-interval */
+long genrand_int31(void)
+{
+ return (long)(genrand_int32()>>1);
+}
+
+
+/* generates a random number on [0,1)-real-interval */
+double genrand_real2(void)
+{
+ return genrand_int32()*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/* generates a random number on (0,1)-real-interval */
+double genrand_real3(void)
+{
+ return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
+ /* divided by 2^32 */
+}
+
+/* generates a random number on [0,1) with 53-bit resolution*/
+double genrand_res53(void)
+{
+ unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6;
+ return (a*67108864.0+b)*(1.0/9007199254740992.0);
+}
+/* These real versions are due to Isaku Wada, 2002/01/09 added */
+
+int main(void)
+{
+ int i;
+ unsigned long init[4]= {0x123, 0x234, 0x345, 0x456}, length=4;
+ init_by_array(init, length);
+ printf("1000 outputs of genrand_int32()\n");
+ for (i=0; i<1000; i++) {
+ printf("%10lu ", genrand_int32());
+ if (i%5==4) printf("\n");
+ }
+ printf("\n1000 outputs of genrand_real2()\n");
+ for (i=0; i<1000; i++) {
+ printf("%10.8f ", genrand_real2());
+ if (i%5==4) printf("\n");
+ }
+ return 0;
+}
+#endif
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
Property changes on: brlcad/trunk/src/libbn/mt19937ar.c
___________________________________________________________________
Added: svn:eol-style
## -0,0 +1 ##
+native
\ No newline at end of property
Added: svn:mime-type
## -0,0 +1 ##
+text/plain
\ No newline at end of property
Deleted: brlcad/trunk/src/libbn/randmt.c
===================================================================
--- brlcad/trunk/src/libbn/randmt.c 2020-04-10 18:29:58 UTC (rev 75340)
+++ brlcad/trunk/src/libbn/randmt.c 2020-04-10 19:51:28 UTC (rev 75341)
@@ -1,176 +0,0 @@
-/* R A N D M T . C
- * BRL-CAD
- *
- * A C-program for MT19937: Real number version
- * genrand() generates one pseudorandom real number (double)
- * which is uniformly distributed on [0, 1]-interval, for each
- * call. sgenrand(seed) set initial values to the working area
- * of 624 words. Before genrand(), sgenrand(seed) must be
- * called once. (seed is any 32-bit integer except for 0).
- * Integer generator is obtained by modifying two lines.
- * Coded by Takuji Nishimura, considering the suggestions by
- * Topher Cooper and Marc Rieffel in July-Aug. 1997.
- *
- * This library is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public
- * License as published by the Free Software Foundation; either
- * version 2 of the License, or (at your option) any later version.
- *
- * This library is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
- * See the GNU Lesser General Public License for more details.
- * You should have received a copy of the GNU Library General
- * Public License along with this library; if not, write to the
- * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
- * 02111-1307 USA
- *
- * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.
- * Any feedback is very welcome. For any question, comments,
- * see http://www.math.keio.ac.jp/matumoto/emt.html or email
- * [email protected]
- */
-
-/* TODO - note that there is now a BSD licensed version of this
- * file at
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
- *
- * Also, there is a BSD licensed SIMD-oriented Fast Mersenne Twister (SFMT)
- * at the same source.
- *
- * Need to investigate whether we can swap out this LGPL version for
- * one or more of the BSD licensed variations
- */
-
-#include "common.h"
-
-#include "bu/malloc.h"
-#include "bu/exit.h"
-#include "bu/vls.h"
-#include "bn/randmt.h"
-
-
-/* Period parameters */
-#define N 624
-#define M 397
-#define MATRIX_A 0x9908b0df /* constant vector a */
-#define UPPER_MASK 0x80000000 /* most significant w-r bits */
-#define LOWER_MASK 0x7fffffff /* least significant r bits */
-
-
-/* Tempering parameters */
-#define TEMPERING_MASK_B 0x9d2c5680
-#define TEMPERING_MASK_C 0xefc60000
-#define TEMPERING_SHIFT_U(y) (y >> 11)
-#define TEMPERING_SHIFT_S(y) (y << 7)
-#define TEMPERING_SHIFT_T(y) (y << 15)
-#define TEMPERING_SHIFT_L(y) (y >> 18)
-
-#define MERSENNE_MAGIC 0x4D54524E
-
-struct _internal_state_s {
- uint32_t magic;
- int mti; /* state index */
- uint32_t mt[N]; /* state vector */
-};
-
-static struct _internal_state_s global_state = { MERSENNE_MAGIC, N+1, {0} };
-
-void *
-bn_randmt_state_create(void)
-{
- struct _internal_state_s *is;
-
- BU_ALLOC(is, struct _internal_state_s);
- is->magic = MERSENNE_MAGIC;
- is->mti = N+1;
- return (void *)is;
-}
-
-/* initializing the array with a NONZERO seed */
-void
-bn_randmt_state_seed(struct _internal_state_s *is, uint32_t seed)
-{
- /*
- * setting initial seeds to mt[N] using
- * the generator Line 25 of Table 1 in
- * [KNUTH 1981, The Art of Computer Programming
- * Vol. 2 (2nd Ed.), pp102]
- */
- is->mt[0] = seed & 0xffffffff;
- for (is->mti = 1; is->mti<N; is->mti++)
- is->mt[is->mti] = (69069 * is->mt[is->mti-1]) & 0xffffffff;
-}
-
-double
-bn_randmt_state(struct _internal_state_s *is)
-{
- uint32_t y;
- static uint32_t mag01[2]={0x0, MATRIX_A};
-
- /* generate N words at one time */
- if (is->mti >= N) {
- int kk;
-
- /* if sgenrand() has not been called, a default initial seed is used */
- if (is->mti == N+1)
- bn_randmt_seed(4357);
-
- for (kk = 0; kk < N-M; kk++) {
- y = (is->mt[kk]&UPPER_MASK)|(is->mt[kk+1]&LOWER_MASK);
- is->mt[kk] = is->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
- }
-
- for (; kk < N-1; kk++) {
- y = (is->mt[kk]&UPPER_MASK)|(is->mt[kk+1]&LOWER_MASK);
- is->mt[kk] = is->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
- }
-
- y = (is->mt[N-1]&UPPER_MASK)|(is->mt[0]&LOWER_MASK);
- is->mt[N-1] = is->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
-
- is->mti = 0;
- }
-
- y = is->mt[is->mti++];
- y ^= TEMPERING_SHIFT_U(y);
- y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
- y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
- y ^= TEMPERING_SHIFT_L(y);
-
- return (double)y / (uint32_t)0xffffffff;
-}
-
-/* straight binhex would result in a 8+4992 character encoding. We should be
able
- * to compress it quite a bit with better encoding? maybe uuencode? */
-void
-bn_randmt_state_serialize(struct _internal_state_s *UNUSED(is), struct bu_vls
*UNUSED(s))
-{
- bu_bomb("Not implemented yet.\n");
-}
-
-void
-bn_randmt_state_deserialize(struct _internal_state_s *UNUSED(is), struct
bu_vls *UNUSED(s))
-{
- bu_bomb("Not implemented yet.\n");
-}
-
-void
-bn_randmt_seed(unsigned long seed)
-{
- bn_randmt_state_seed(&global_state, (uint32_t)seed);
-}
-
-double bn_randmt(void)
-{
- return bn_randmt_state(&global_state);
-}
-
-/*
- * Local Variables:
- * mode: C
- * tab-width: 8
- * indent-tabs-mode: t
- * c-file-style: "stroustrup"
- * End:
- * ex: shiftwidth=4 tabstop=8
- */
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
_______________________________________________
BRL-CAD Source Commits mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/brlcad-commits
