Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mt_stream.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago
sources for flexwrf v3.1
File size: 23.1 KB

Line
1	!
2	! A Fortran90/95 module program for Multiple Stream Mersenne Twister.
3	! 2010/02/12, Ken-Ichi Ishikawa [ishikawa[at]theo.phys.sci.hiroshima-u.ac.jp]
4	! 2010/07/20, M. S. Briggs, [michael.s.briggs[at]nasa.gov], minor improvements.
5	!
6	! This module provides the Mersenne Twister pseudo random number generator
7	! (MT-PRNG) with multiple stream ability.
8	! A long period MT generaotr (backbone stream) is divided into disjoint streams,
9	! and the each stream has continuous stream from the backbone stream.
10	! The distance between each stream is chosen to have long enough length.
11	! The status and parameters of the each stream are encapsulated in
12	! the f90 type components so that we can use multiple streams simultaneously.
13	! The stream length is fixed to 2^nj.
14	!
15	! This code is converted from original/sample codes located at
16	! Mersenne Twister Home Page:
17	! [http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html]
18	!
19	! For the Original MT PRNG, see also:
20	! M. Matsumoto and T. Nishimura,
21	! "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom
22	! number generator",
23	! ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1,
24	! January pp.3-30 (1998) DOI:10.1145/272991.272995.
25	!
26	! For the jump ahead mechanism, see also:
27	! H. Haramoto, M. Matsumoto, T. Nishimura, F. Panneton, and P. L'Ecuyer,
28	! "Efficient Jump Ahead for F_2-Linear Random Number Generators",
29	! GERAD Report G-2006-62. INFORMS Journal on Computing, 20, 3 (2008), 385-390.
30	!
31	! This routine uses;
32	! Fast arithmetic in GF(2)[x], [http://wwwmaths.anu.edu.au/~brent/software.html]
33	! NTL : A Library for doing Number Theory, [http://www.shoup.net/ntl/index.html]
34	!
35	!
36	! Copyright (c) 2010, Ken-Ichi Ishikawa [ishikawa[at]theo.phys.sci.hiroshima-u.ac.jp]
37	! All rights reserved.
38	!
39	! Redistribution and use in source and binary forms, with or without
40	! modification, are permitted provided that the following conditions are
41	! met:
42	!
43	! * Redistributions of source code must retain the above copyright
44	! notice, this list of conditions and the following disclaimer.
45	!
46	! * Redistributions in binary form must reproduce the above copyright
47	! notice, this list of conditions and the following disclaimer listed
48	! in this license in the documentation and/or other materials
49	! provided with the distribution.
50	!
51	! * Neither the name of the copyright holders nor the names of its
52	! contributors may be used to endorse or promote products derived from
53	! this software without specific prior written permission.
54	!
55	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
56	! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
57	! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
58	! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
59	! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
60	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
61	! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
62	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
63	! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
64	! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
65	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
66	!
67
68	!define _TEMPERING_
69
70	!ifdef _DEBUG_
71	!module mt_stream_debug
72	!else
73	module mt_stream
74	!endif
75	use mt_kind_defs
76	implicit none
77	private
78	public :: MTS_SUCCESS
79	public :: MTS_FAIL
80	public :: mt_state
81	public :: new,delete
82	public :: init
83	public :: print
84	public :: save
85	public :: read
86	public :: set_mt19937
87	public :: create_stream
88	public :: genrand_int32
89	public :: genrand_double1
90	public :: genrand_double2
91	public :: genrand_double3
92	public :: genrand_double4
93
94	integer(INT32), parameter :: MTS_SUCCESS = 0
95	integer(INT32), parameter :: MTS_FAIL = -1
96	!#ifdef _DEBUG_
97	! integer(INT32), parameter :: MT_JUMP_DISTANCE_EXP = 16 ! test jump distance (2^16 steps)
98	!#else
99	integer(INT32), parameter :: MT_JUMP_DISTANCE_EXP = 256 ! default jump distance (2^256 steps)
100	!#endif
101
102	!======================================
103	! MT19937 parameters
104	!======================================
105	integer(INT32), parameter :: MT19937_N = 624
106	integer(INT32), parameter :: MT19937_M = 397
107	integer(INT32), parameter :: MT19937_W = 32
108	integer(INT32), parameter :: MT19937_R = 31
109	integer(INT32), parameter :: MT19937_MATA = INT (Z'9908b0df', INT32)
110	integer(INT32), parameter :: MT19937_WMASK = INT (Z'ffffffff', INT32)
111	integer(INT32), parameter :: MT19937_UMASK = INT (Z'80000000', INT32)
112	integer(INT32), parameter :: MT19937_LMASK = INT (Z'7fffffff', INT32)
113	integer(INT32), parameter :: MT19937_SHFT0 = 11
114	integer(INT32), parameter :: MT19937_SHFT1 = 18
115	integer(INT32), parameter :: MT19937_SHFTB = 7
116	integer(INT32), parameter :: MT19937_SHFTC = 15
117	integer(INT32), parameter :: MT19937_MASKB = INT (Z'9d2c5680', INT32)
118	integer(INT32), parameter :: MT19937_MASKC = INT (Z'efc60000', INT32)
119
120	!======================================
121	! MT state
122	! Period = 2^(N*W-R)-1
123	!======================================
124	type mt_state
125	private
126	integer(INT32) :: i = -1 ! state vector index
127	integer(INT32) :: stream_id = -1 ! stream ID
128	integer(INT32) :: istatus = -1 ! initialization status
129	integer(INT32) :: nn = -1 ! MT parameter N
130	integer(INT32) :: mm = -1 ! MT parameter M
131	integer(INT32) :: rr = -1 ! MT parameter R
132	integer(INT32) :: ww = -1 ! MT parameter W (width =32)
133	integer(INT32) :: aaa = 0 ! Companion matrix parameter
134	integer(INT32) :: wmask = 0 ! 32-bit mask
135	integer(INT32) :: umask = 0 ! Twist mask x(1)
136	integer(INT32) :: lmask = 0 ! Twist mask x(0)
137	integer(INT32) :: shift0 = 0 ! Temparing parameters ...
138	integer(INT32) :: shift1 = 0
139	integer(INT32) :: maskB = 0
140	integer(INT32) :: maskC = 0
141	integer(INT32) :: shiftB = 0
142	integer(INT32) :: shiftC = 0
143	integer(INT32) :: mag(0:1) = 0 ! mag(0) = 0, mag(1) = aaa
144	integer(INT32), pointer :: state(:) => NULL() ! state vector
145	end type
146
147	type(mt_state), save :: g_mt_master ! this keeps MT parameters
148	integer(INT32), save :: total_stream = 0
149
150	interface new
151	module procedure mt_new
152	end interface
153
154	interface delete
155	module procedure mt_delete
156	end interface
157
158	interface print
159	module procedure mt_print
160	end interface
161
162	interface read
163	module procedure mt_read
164	end interface
165
166	interface save
167	module procedure mt_save
168	end interface
169
170	interface set_mt19937
171	module procedure mt_set_mt19937
172	end interface
173
174	interface init
175	module procedure mt_init_by_scalar
176	module procedure mt_init_by_array
177	end interface
178
179	interface create_stream
180	module procedure mt_create_stream
181	end interface
182
183	interface genrand_int32
184	! in [0,0xFFFFFFFF]
185	module procedure mt_genrand_int32
186	end interface
187
188	interface genrand_double1
189	! in [0,1] (53-bit resolution)
190	module procedure mt_genrand_double1
191	end interface
192
193	interface genrand_double2
194	! in [0,1) (53-bit resolution)
195	module procedure mt_genrand_double2
196	end interface
197
198	interface genrand_double3
199	! in (0,1) (52-bit resolution)
200	module procedure mt_genrand_double3
201	end interface
202
203	interface genrand_double4
204	! in (0,1] (53-bit resolution)
205	module procedure mt_genrand_double4
206	end interface
207
208	contains
209
210	subroutine mt_new(this,ierr)
211	!
212	!= Allocate mt state
213	!
214	implicit none
215	type(mt_state), intent(inout) :: this
216	integer(INT32), optional, intent(out) :: ierr
217	character(256), parameter :: myname="mt_new"
218	integer(INT32) :: jerr,kerr
219	jerr = MTS_SUCCESS
220	if (0 == g_mt_master%nn) then
221	write(*,'(A,": MT master parameter is not initialized.")')TRIM(myname)
222	write(*,'(A,": Stop!")')TRIM(myname)
223	stop
224	endif
225	call mt_copy(g_mt_master,this)
226	this%stream_id = total_stream
227	total_stream = total_stream + 1
228	if (.not.associated(this%state)) then
229	if (this%nn <= 0) then
230	jerr = MTS_FAIL
231	this%istatus = jerr
232	goto 100
233	endif
234	allocate(this%state(0:this%nn-1),STAT=kerr)
235	if (kerr /= 0) then
236	jerr = MTS_FAIL
237	this%istatus = jerr
238	goto 100
239	endif
240	this%istatus = jerr
241	goto 100
242	endif
243	100 if (present(ierr)) then
244	ierr = jerr
245	return
246	else
247	if (jerr /= MTS_SUCCESS) then
248	write(*,'(A,": Something wrong.")')TRIM(myname)
249	write(*,'(A,": Stop!")')TRIM(myname)
250	stop
251	endif
252	endif
253	end subroutine
254
255	subroutine mt_delete(this)
256	!
257	!= deallocate state vector
258	!
259	implicit none
260	type(mt_state), intent(inout) :: this
261	if (associated(this%state)) then
262	deallocate(this%state)
263	endif
264	total_stream = total_stream - 1
265	return
266	end subroutine
267
268	subroutine mt_set_mt19937
269	!
270	!= set MT19937 parameters on g_mt_master
271	!
272	implicit none
273	g_mt_master%nn = MT19937_N
274	g_mt_master%mm = MT19937_M
275	g_mt_master%ww = MT19937_W
276	g_mt_master%rr = MT19937_R
277	g_mt_master%aaa = MT19937_MATA
278	g_mt_master%wmask = MT19937_WMASK
279	g_mt_master%umask = MT19937_UMASK
280	g_mt_master%lmask = MT19937_LMASK
281	g_mt_master%shift0 = MT19937_SHFT0
282	g_mt_master%shift1 = MT19937_SHFT1
283	g_mt_master%shiftB = MT19937_SHFTB
284	g_mt_master%shiftC = MT19937_SHFTC
285	g_mt_master%maskB = MT19937_MASKB
286	g_mt_master%maskC = MT19937_MASKC
287	g_mt_master%mag(0) = 0
288	g_mt_master%mag(1) = g_mt_master%aaa
289	total_stream = 0
290	return
291	end subroutine
292
293	subroutine mt_copy(this,that)
294	!
295	!= duplicate/copy MT paramters
296	!
297	! that <= this
298	!
299	implicit none
300	type(mt_state), intent(in) :: this
301	type(mt_state), intent(out) :: that
302	that%i = this%i
303	that%stream_id = this%stream_id
304	that%istatus = this%istatus
305	that%nn = this%nn
306	that%nn = this%nn
307	that%mm = this%mm
308	that%ww = this%ww
309	that%aaa = this%aaa
310	that%rr = this%rr
311	that%wmask = this%wmask
312	that%umask = this%umask
313	that%lmask = this%lmask
314	that%shift0 = this%shift0
315	that%shift1 = this%shift1
316	that%shiftB = this%shiftB
317	that%shiftC = this%shiftC
318	that%maskB = this%maskB
319	that%maskC = this%maskC
320	that%mag(0) = 0
321	that%mag(1) = that%aaa
322	return
323	end subroutine
324
325	subroutine mt_save(this,unit)
326	!
327	!= save MT status to file with unit
328	!
329	implicit none
330	type(mt_state), intent(in) :: this
331	integer(INT32), intent(in) :: unit
332	write(unit) this%i
333	write(unit) this%stream_id
334	write(unit) this%istatus
335	write(unit) this%nn
336	write(unit) this%mm
337	write(unit) this%ww
338	write(unit) this%aaa
339	write(unit) this%rr
340	write(unit) this%wmask
341	write(unit) this%umask
342	write(unit) this%lmask
343	write(unit) this%shift0
344	write(unit) this%shift1
345	write(unit) this%shiftB
346	write(unit) this%shiftC
347	write(unit) this%maskB
348	write(unit) this%maskC
349	write(unit) this%mag(:)
350	write(unit) this%state(:)
351	return
352	end subroutine
353
354	subroutine mt_read(this,unit)
355	!
356	!= read MT status from file with unit
357	!
358	implicit none
359	type(mt_state), intent(inout) :: this
360	integer(INT32), intent(in) :: unit
361	read(unit) this%i
362	read(unit) this%stream_id
363	read(unit) this%istatus
364	read(unit) this%nn
365	read(unit) this%mm
366	read(unit) this%ww
367	read(unit) this%aaa
368	read(unit) this%rr
369	read(unit) this%wmask
370	read(unit) this%umask
371	read(unit) this%lmask
372	read(unit) this%shift0
373	read(unit) this%shift1
374	read(unit) this%shiftB
375	read(unit) this%shiftC
376	read(unit) this%maskB
377	read(unit) this%maskC
378	read(unit) this%mag(:)
379	if (associated(this%state)) then
380	deallocate(this%state)
381	allocate(this%state(0:this%nn-1))
382	endif
383	read(unit) this%state(:)
384	return
385	end subroutine
386
387	subroutine mt_print(this)
388	!
389	!= print MT stream parameter and states
390	!
391	implicit none
392	type(mt_state), intent(in) :: this
393	character(256) :: str
394	write(*,'("===============================================================================")')
395	write(*,'(" MT PARAMETERS")')
396	write(str,*)this%stream_id
397	write(*,'(" STREAMID= ",A)')TRIM(ADJUSTL(str))
398	write(str,)this%nnthis%ww-this%rr
399	write(*,'(" PERIOD= 2^",A," - 1")')TRIM(ADJUSTL(str))
400	write(*,'(" NN=",I8," MM=",I8," WW=",I3," RR=",I3)')this%nn,this%mm,this%ww,this%rr
401	write(*,'(" AVEC= 0x",Z8.8)', advance='no' )this%aaa
402	write(*,'(" WMASK= 0x",Z8.8)', advance='no' )this%wmask
403	write(*,'(" UMASK= 0x",Z8.8)', advance='no' )this%umask
404	write(*,'(" LMASK= 0x",Z8.8)') this%lmask
405	write(*,'(" MAG= 0x",Z8.8," 0x",Z8.8)') this%mag
406	write(*,'(" SHFT0=",I3)', advance='no' )this%shift0
407	write(*,'(" SHFT1=",I3)', advance='no' )this%shift1
408	write(*,'(" SHFTB=",I3)', advance='no' )this%shiftB
409	write(*,'(" SHFTC=",I3)', advance='no' )this%shiftC
410	write(*,'(" MASKB= 0x",Z8.8)', advance='no' )this%maskB
411	write(*,'(" MASKC= 0x",Z8.8)') this%maskC
412	write(str,*)this%i
413	write(*,'(" POINTER= ",A)')TRIM(ADJUSTL(str))
414	write(*,'("===============================================================================")')
415	return
416	end subroutine
417
418
419	subroutine mt_init_by_scalar(this,iseed,ierr)
420	!
421	!= initialize MT state by a scalar seed.
422	!
423	implicit none
424	type(mt_state), intent(inout) :: this
425	integer(INT32), intent(in) :: iseed
426	integer(INT32), optional, intent(out) :: ierr
427	character(256), parameter :: myname="mt_init_by_scalar"
428	integer(INT32) :: i,jseed,jerr
429	if (.not.associated(this%state)) then
430	jerr = MTS_FAIL
431	goto 100
432	endif
433	jseed = iseed
434	do i=0,this%nn-1
435	this%state(i) = jseed
436	jseed = 1812433253 * IEOR(jseed, ISHFT(jseed,-30)) + i + 1
437	enddo
438	this%i = this%nn
439	do i=0,this%nn-1
440	this%state(i) = IAND(this%state(i),this%wmask)
441	enddo
442	this%mag(0) = 0
443	this%mag(1) = this%aaa
444	jerr = MTS_SUCCESS
445	100 if (present(ierr)) then
446	ierr = jerr
447	else
448	if (jerr /= MTS_SUCCESS) then
449	write(*,'(A,": State vector allocation fails.")')TRIM(myname)
450	write(*,'(A,": Stop!")')TRIM(myname)
451	stop
452	endif
453	endif
454	return
455	end subroutine
456
457	subroutine mt_init_by_array(this,iseed,ierr)
458	!
459	!= initialize MT state by array seeds.
460	!
461	implicit none
462	type(mt_state), intent(inout) :: this
463	integer(INT32), intent(in) :: iseed(0:)
464	integer(INT32), optional, intent(out) :: ierr
465	integer :: isize,i,j,k,n,jerr
466	character(256), parameter :: myname="mt_init_by_array"
467	integer(INT32), parameter :: ASEED = 19650218
468
469	call mt_init_by_scalar(this,ASEED,jerr)
470	if (jerr /= MTS_SUCCESS) goto 100
471
472	isize = SIZE(iseed)
473	n = this%nn
474	i = 1
475	j = 0
476	do k = MAX(n,isize),1,-1
477	this%state(i) = IEOR(this%state(i),(IEOR(this%state(i-1),ISHFT(this%state(i-1),-30))*1664525)) + iseed(j) + j
478	this%state(i) = IAND(this%state(i),this%wmask) ! for WORDSIZE > 32 machines
479	i = i + 1
480	j = j + 1
481	if (i >= n) then
482	this%state(0) = this%state(n-1)
483	i = 1
484	endif
485	if (j >= isize) j = 0
486	enddo
487	do k = n-1,1,-1
488	this%state(i) = IEOR(this%state(i),(IEOR(this%state(i-1),ISHFT(this%state(i-1),-30))*1566083941)) - i
489	this%state(i) = IAND(this%state(i),this%wmask) ! for WORDSIZE > 32 machines
490	i = i + 1
491	if (i >= n) then
492	this%state(0) = this%state(n-1)
493	i = 1
494	endif
495	enddo
496	this%state(0) = INT (Z'80000000', INT32) ! MSB is 1; assuring non-zero initial array
497	jerr = MTS_SUCCESS
498	100 if (present(ierr)) then
499	ierr = jerr
500	else
501	if (jerr /= MTS_SUCCESS) then
502	write(*,'(A,": mt_init_by_scalar fails.")')TRIM(myname)
503	write(*,'(A,": Stop!")')TRIM(myname)
504	stop
505	endif
506	endif
507	return
508	end subroutine
509
510	function mt_genrand_int32(this) result(ir)
511	!
512	!= return a value in [0,0xFFFFFFFF]
513	!
514	implicit none
515	type(mt_state), intent(inout) :: this
516	integer(INT32) :: ir
517	integer(INT32) :: umask,lmask,n,m,is
518	integer(INT32) :: k,nm,n1
519	if (this%i >= this%nn) then
520	n = this%nn
521	m = this%mm
522	lmask = this%lmask
523	umask = this%umask
524	nm = n - m
525	n1 = n - 1
526	do k=0,nm-1
527	is = IOR(IAND(this%state(k),umask),IAND(this%state(k+1),lmask))
528	this%state(k) = IEOR(IEOR(this%state(k+m),ISHFT(is,-1)),this%mag(IAND(is,1)))
529	enddo
530	do k=nm,n1-1
531	is = IOR(IAND(this%state(k),umask),IAND(this%state(k+1),lmask))
532	this%state(k) = IEOR(IEOR(this%state(k+m-n),ISHFT(is,-1)),this%mag(IAND(is,1)))
533	enddo
534	is = IOR(IAND(this%state(n-1),umask),IAND(this%state(0),lmask))
535	this%state(n-1) = IEOR(IEOR(this%state(m-1),ISHFT(is,-1)),this%mag(IAND(is,1)))
536	this%i = 0
537	endif
538
539	is = this%state(this%i)
540	this%i = this%i + 1
541	!#ifdef _TEMPERING_
542	! is = IEOR(is,ISHFT(is,-this%shift0))
543	! is = IEOR(is,IAND(ISHFT(is, this%shiftB),this%maskB))
544	! is = IEOR(is,IAND(ISHFT(is, this%shiftC),this%maskC))
545	! is = IEOR(is,ISHFT(is,-this%shift1))
546	!#endif
547	ir = is
548	return
549	end function
550
551	subroutine mt_matvec(this,v,w)
552	!
553	!= Multiply transition matrix on a state vector v
554	!
555	! w = B v
556	!
557	! this : MT parameters(transition matrix)
558	! v : input vector
559	! w : output vector
560	!
561	implicit none
562	type(mt_state), intent(in) :: this
563	integer(INT32), intent(in) :: v(0:this%nn-1)
564	integer(INT32), intent(out) :: w(0:this%nn-1)
565	integer(INT32) :: umask,lmask,n,m,is
566	integer(INT32) :: k
567	n = this%nn
568	m = this%mm
569	lmask = this%lmask
570	umask = this%umask
571	w(0) = IAND(v(1),umask)
572	do k=1,n-2
573	w(k) = v(k+1)
574	enddo
575	is = IOR(IAND(v(0),umask),IAND(v(1),lmask))
576	w(n-1) = IEOR(IEOR(v(m),ISHFT(is,-1)),this%mag(IAND(is,1)))
577	return
578	end subroutine
579
580	subroutine mt_create_stream(this,that,id)
581	!
582	!= Create New stream (that) with distance id*2^(jp) from (this)
583	!
584	implicit none
585	type(mt_state), intent(inout) :: this ! input state
586	type(mt_state), intent(inout) :: that ! output state
587	integer(INT32), intent(in) :: id
588	character(256), parameter :: myname="mt_create_stream"
589	integer(INT32), parameter :: jp = MT_JUMP_DISTANCE_EXP
590	if (id < 0) then
591	write(*,'(A,": Positive ID is requried.")')TRIM(myname)
592	write(*,'(A,": Stop!")')TRIM(myname)
593	stop
594	endif
595	call mt_jumpahead(this,that,jp,id)
596	that%stream_id = id
597	total_stream = total_stream + 1
598	this%i = this%nn
599	return
600	end subroutine
601
602	subroutine mt_jumpahead(this,that,jp,id)
603	!
604	!= Jump ahead by (id+1) * 2^jp steps.
605	!
606	! this : input state
607	! that : output state, proceeds by id*2^jp steps.
608	!
609	implicit none
610	type(mt_state), intent(inout) :: this ! input state
611	type(mt_state), intent(inout) :: that ! output state
612	integer(INT32), intent(in) :: jp ! exponent (jump step = id*2^jp)
613	integer(INT32), intent(in) :: id ! id (jump step = id*2^jp)
614	integer(INT32) :: v(0:this%nn-1)
615	integer(INT32) :: w(0:this%nn-1)
616	integer(INT32) :: s(0:this%nn-1)
617	integer(INT32) :: i,iwp,ibp
618	character(256) :: str,str2
619	integer(INT32) :: p(0:this%nn-1)
620	integer(INT32) :: np
621	!
622	! external routine in
623	! jump_ahead_coeff/get_coeff.cxx
624	! written in C++ with NTL and gf2x libraries.
625	!
626	!#ifdef _NTL_
627	! interface get_coeff_interface
628	! subroutine get_coeff(nn,mm,rr,ww,avec,nj,id,pp,np) bind (C, name="get_coeff")
629	! use, intrinsic :: iso_c_binding, only: c_int
630	! integer (c_int), intent(in) :: nn,mm,rr,ww,avec,nj,id
631	! integer (c_int), intent(inout) :: pp(0:nn-1),np
632	! end subroutine get_coeff
633	! end interface get_coeff_interface
634	!#else
635	interface f_get_coeff_interface
636	subroutine f_get_coeff(nn,mm,rr,ww,avec,nj,id,pp,np)
637	integer, intent(in) :: nn,mm,rr,ww,avec,nj,id
638	integer, intent(inout) :: pp(0:nn-1),np
639	end subroutine f_get_coeff
640	end interface f_get_coeff_interface
641	!#endif
642
643	!==================================
644	! state copy: this => that
645	!==================================
646	call mt_copy(this,that)
647
648	write(str,'(I12)')id
649	write(str2,'(I12)')jp
650	write(*,'("# ID ",I12)')id
651	write(,'("# Jump Ahead by (",A,")2^",A)')TRIM(ADJUSTL(str)),TRIM(ADJUSTL(str2))
652	if ( this%i /= this%nn) then
653	write(*,'("Error in jumpahead: input state pointer should point nn.")')
654	write(*,'("this%i = ",I12)')this%i
655	write(*,'("this%nn = ",I12)')this%nn
656	write(str,*)this%i
657	write(str2,*)this%nn-1
658	write(*,'(A,"-",A," random numbers are dropped.")') &
659	& TRIM(ADJUSTL(str)),TRIM(ADJUSTL(str2))
660	write(*,'("forced to point nn")')
661	this%i = this%nn
662	endif
663
664	!==================================
665	! compute jump ahead polynomial
666	! p(x) coefficients
667	! for this MT parameter
668	!==================================
669	!#ifdef _NTL_
670	! call get_coeff(this%nn,this%mm,this%rr,this%ww,this%aaa,jp,id,p,np)
671	!#else
672	call f_get_coeff(this%nn,this%mm,this%rr,this%ww,this%aaa,jp,id,p,np)
673	!#endif
674
675	!==================================
676	! multiply p(B) on a state vector v
677	! p(x) : jump ahead polynomial
678	! B : transition matrix
679	! with simple Horner's method
680	! w = p(B) v = B^(2^jp) v
681	!==================================
682	v(:) = this%state(:)
683	iwp = (np-1)/32
684	ibp = mod(np-1,32)
685	if (BTEST(p(iwp),ibp)) then
686	w(:) = v(:)
687	endif
688	do i=np-2,0,-1
689	iwp = i/32
690	ibp = mod(i,32)
691	call mt_matvec(this,w,s) ! s = B w
692	if (BTEST(p(iwp),ibp)) then
693	w(:) = IEOR(v(:),s(:)) ! w = 1 v + s
694	else
695	w(:) = s(:) ! w = 0 v + s
696	endif
697	enddo
698
699	if (.not. associated(that%state)) then
700	allocate(that%state(0:that%nn-1))
701	endif
702	that%state(:) = w(:)
703	that%i = this%nn
704
705	return
706	end subroutine
707
708	function mt_genrand_double1(this) result(r)
709	!
710	! r in [0,1] (53-bit resolution)
711	!
712	implicit none
713	type(mt_state), intent(inout) :: this
714	real(REAL64) :: r
715	real(REAL64) :: a,b
716	integer(INT32) :: ia,ib
717	ia = mt_genrand_int32(this) ! ia in [0,0xFFFFFFFF]
718	ib = mt_genrand_int32(this) ! ib in [0,0xFFFFFFFF]
719	ia = ISHFT(ia,-5) ! ia in [0,2^27-1]
720	ib = ISHFT(ib,-6) ! ib in [0,2^26-1]
721	a = REAL(ia,kind=KIND(r))
722	b = REAL(ib,kind=KIND(r))
723	!===============================
724	! ( a*2^26 + b ) in [0,2^53-1]
725	! r = ( a*2^26 + b )/(2^53-1)
726	!===============================
727	r = (a67108864.0_REAL64 + b)(1.0_REAL64/9007199254740991.0_REAL64)
728	return
729	end function
730
731	function mt_genrand_double2(this) result(r)
732	!
733	! r in [0,1) (53-bit resolution)
734	!
735	implicit none
736	type(mt_state), intent(inout) :: this
737	real(REAL64) :: r
738	real(REAL64) :: a,b
739	integer(INT32) :: ia,ib
740	ia = mt_genrand_int32(this) ! ia in [0,0xFFFFFFFF]
741	ib = mt_genrand_int32(this) ! ib in [0,0xFFFFFFFF]
742	ia = ISHFT(ia,-5) ! ia in [0,2^27-1]
743	ib = ISHFT(ib,-6) ! ib in [0,2^26-1]
744	a = REAL(ia,kind=KIND(r))
745	b = REAL(ib,kind=KIND(r))
746	!===============================
747	! ( a*2^26 + b ) in [0,2^53-1]
748	! r = ( a*2^26 + b )/(2^53)
749	!===============================
750	r = (a67108864.0_REAL64 + b)(1.0_REAL64/9007199254740992.0_REAL64)
751	return
752	end function
753
754	function mt_genrand_double3(this) result(r)
755	!
756	! r in (0,1) (52-bit resolution)
757	!
758	implicit none
759	type(mt_state), intent(inout) :: this
760	real(REAL64) :: r
761	real(REAL64) :: a,b
762	integer(INT32) :: ia,ib
763	ia = mt_genrand_int32(this) ! ia in [0,0xFFFFFFFF]
764	ib = mt_genrand_int32(this) ! ib in [0,0xFFFFFFFF]
765	ia = ISHFT(ia,-6) ! ia in [0,2^26-1]
766	ib = ISHFT(ib,-6) ! ib in [0,2^26-1]
767	a = REAL(ia,kind=KIND(r))
768	b = REAL(ib,kind=KIND(r))
769	!===============================
770	! ( a*2^26 + b ) in [0,2^52-1]
771	! r = ( a*2^26 + b + 1/2 )/(2^52)
772	!===============================
773	r = (a67108864.0_REAL64 + b + 0.5_REAL64)(1.0_REAL64/4503599627370496.0_REAL64)
774	return
775	end function
776
777	function mt_genrand_double4(this) result(r)
778	!
779	! r in (0,1] (53-bit resolution)
780	!
781	implicit none
782	type(mt_state), intent(inout) :: this
783	real(REAL64) :: r
784	r = 1.0_REAL64 - mt_genrand_double2(this)
785	return
786	end function
787
788	end module

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format