JPH04199225A - Arithmetic random number generating method in computer - Google Patents

Abstract

PURPOSE:To increase the speed in this method by generating plural pieces of random numbers in a batch and simultaneously performing vector processing for the generation of the plural random numbers. CONSTITUTION:A vector computer is constituted of a scalar processor 400, vector processor 500, input-output processor 600, and main storage device 700 and an extended storage device 800 may be added to the computer. In addition, a function which can generate numerous uniform random numbers by a product combining method is provided and N pieces of powers of multipliers or already generated random numbers are held so that vector processing can be performed by independently performing random number generating calculation on the N pieces of random numbers. Therefore, the processing speed of the vector computer is increased.

Description

DETAILED DESCRIPTION OF THE INVENTION [Industrial Field of Application] The present invention relates to the use of arithmetic random numbers in computers.
r), and particularly relates to a processing method suitable for using a vector computer or an array computer.

[Conventional technology]

Arithmetic random numbers (pseudo-random numbers) are sequentially generated by performing certain numerical operations on given initial values. Arithmetic random numbers that follow various distributions are required, but uniform random numbers are particularly fundamental and important. This is because if a uniform random number is obtained, random numbers following other distributions can be generated based on the uniform random number. The formula for calculating random numbers is
For example, "Monte Carlo Method and Simulation" by Takao Tsuda.

Baifukan (1968), pages 7 to 37 provide detailed information.

The product congruence method, which is a typical method for generating uniform random numbers, will be explained. When the initial value M04 multiplier a, (M, , a is a positive integer) is given, Mr+=aMn-0 (mod L)
-(1) However, n=1. ．． An integer sequence M of 2,3 degrees. 2M2. If you create ..., the interval (0, L
) is obtained, which is a uniformly distributed arithmetic random number. This is the multiplication joint method.

An integer type variable is represented by Ω+1 bits (the first bit represents the sign, and the area is [-2', 2'-1]),
Allow overflow during integer operations (does not cause interrupts)
When using a calculator, if n = 21 and the multiplication a''Mn-1 on the right side of equation (1) is performed and the most significant sign bit is ignored (i.e., treated as 0), the number represented by the lower Q bits is:
This is the remainder when M is divided by L. This method requires processing to take the modulus (generally requires verification, rounding down, multiplication by 2, and subtraction)
It is widely used because it does not require

Uniform random numbers are often used after being standardized so that they are distributed in the interval (0, 1), but for this purpose, they can be divided by Mn (mod L).

The numerical expression range of Mn is [0, L-1], but 1M
If we choose prime numbers other than 2 as 0 and a.

Since Mn=O does not hold and O<M and <L, O<R and <1.

Numerical simulation using random numbers is generally referred to as the Monte Carlo method or Monte Carlo simulation, and this method is applicable to many natural phenomena.

It is used to elucidate social phenomena. Further elucidation of the phenomenon (
This is equivalent to solving an equation of state), so it can also be used as part of a device control device. The application of the Monte Cal method is extremely wide, for example.

・Noise phenomenon, particle collision process, waiting phenomenon, analysis of traffic flow, analysis of electronic circuits, thermal conduction of plasma, macroscopic instability due to noise, and lattice vibration. These are described in detail in Chapters 8 and 9 of Mr. Tsuda's book mentioned above.

[Problem to be solved by the invention]

FIG. 2 is an example of a program using uniform random numbers written in the FORTRAN language. In the figure, the numbers at the right end are card numbers, 100 to 160 are part of the program that calls the subprogram RANDOM that generates uniform random numbers, and 200 to 250 indicate the contents of the subprogram RANDOM.

Here, K and M correspond to M, , and a in equation (1), respectively, and are 97 and 13, respectively. Statement 240 is a process for converting R to a positive number since when the first bit is 1, R becomes a negative number. Note that in one sentence 220, K may be negative, and the values are generally different depending on whether the Q bits below the sign bit are taken or the absolute value, but in general, processing as shown in this example is often sufficient.

As described above, the method of implementing the random number generation unit as a function subprogram or a subroutine subprogram is widely used because it is convenient and easy to provide as a library. However, this method has the disadvantage of requiring a call to an external procedure each time a random number is generated, resulting in a heavy load (and therefore time).

Furthermore, when a vector computer is used, there is a drawback that a D* loop that includes an external procedure call cannot be vectorized. Here, a vector computer is a so-called supercomputer equipped with multiple vector calculation machines and vector registers. The specifications of these computers are shown in Table 1.11 on page 23 of "Supercomputer" by Kenbe Murata et al., Maruzen (published in 1985).

- As an example, if we take Hitachi S-810/20,
Its configuration is as shown in FIG.
0. Vector processor 50o, bathing force processing device 600
．． It consists of a main storage device 700, and an extended storage device 800 may be added.

By the way, in order to reduce the load of external procedure calls,
As shown in FIG. 4, the function subprogram can be expanded inline. However, when using the above-mentioned vector calculator, simply performing inline expansion does not speed up the process very much. This is because this calculation involves first-order cyclic operations, and this type of operation is slow compared to addition, subtraction, multiplication, and division.

An object of the present invention is to provide a random number generation method that avoids the above-mentioned problems when using a parallel computer such as a hector computer or an array computer.

[Means to solve the problem]

In order to achieve such an object, the present invention.

Using a computer that has the function of generating a large number of uniform random numbers using the product congruence method, you can generate N random numbers by storing N multipliers raised to the power of - or N pre-created random numbers. It is characterized by making them mutually independent and performing vector processing.

[Effect]

By generating a plurality of random numbers at once, the present invention reduces the load of external procedure calls and enables vector processing of the generation of a plurality of random numbers. Furthermore, a random number sequence M1. , M2. ....

When calculating Mn, . . . , the derivation method is devised so that −th order cyclic operations are not included, thereby increasing the processing speed.

From the definition of equation (1), for any positive numbers m and n, M, +
, =a″+1-M, =-(3)=
a ”-MIl-(3)' ” an-M* ・-(3)
' holds true. Therefore, N successive random numbers M,.

MIL+1. +...If MII+N-□ is found, the relational expression MII+N= a N1Mn
,,, (4) However, m < n < m
+ N, the following N random numbers M tn+N+ M1n
+N+,.

...9M, ya2N-0 can be obtained by simple multiplication. That is, high-speed vector processing is possible.

Here, since aN does not depend on the index r1, it is only necessary to calculate - degrees.

〔Example〕

First, the principle of the present invention will be explained.

FIG. 1 shows a basic method for random number generation according to the present invention.

Step S1 is an initialization process, in which the initial value Mo is multiplied by a to obtain Ml, and further multiplied by a to obtain M2.

Repeat the same procedure to find up to MN. A similar calculation is performed for the initial value 1 to find aN. Step S2 is a state in which initialization is completed, that is, aN and N successive random numbers M. .

Mn and 8. . . . When Mn+N-1 has been found, this is the part to find N random numbers following it. According to formula (4), M n + N - t = a N - Mn
-(4)'where i=
o, 1, 2. It may be set to N-1. However, in order to obtain (M, .) from the second time onwards, it is not necessary to go through the initialization section, so an entry is provided between step S1 and step S2 (see the embodiments from FIG. 4 onwards). In addition.

The standardization process is omitted in FIG.

An embodiment of the present invention will be described below with reference to FIG. The blocks up to block 7 in the figure are FAI) diagrams showing the initialization process for obtaining aN and Myu9M2°...2MN.

The processing contents will be explained below in the order of the boxes in the figure.

1: Entrance to initialization processing.

2: Set the initial value 1 to the parameter AN.

AN is an integer type variable representing N.

3: Perform loop control to repeat procedures 4 to 6.

4: AN (= a'-'(mad L))
Multiply a by a to find a'l+1 and set it as AN again. M.S.
K is a constant whose first bit is O and the following Ω bits are all 1, and IAND is a vector processing function that performs logical product for each bit. As a result, AN
Even if the first bit is set to 1, it is removed.

5: Find Mn according to equation (1). The sign operation is the same as in 4.

6: Multiply and divide Mn to standardize. This division must be performed as an operation between real variables (usually using floating point representation).

FLOAT is a type conversion function that converts an integer type variable to a real number type variable.

To create a random number sequence for the second and subsequent times, it is sufficient to enter from block 8.

9: Performs loop control similar to 3.

10: Calculate aN-Mn according to equation (4), remove the sign bit, and set Mn again. That is, Mn is M n+
It is updated to the value of N.

], standardize Mゎ in the same way as in 6.

In the embodiment shown in FIG. 5, the initialization part is performed in accordance with the definition formula, but it is also possible to perform parallel processing by repeating the same method from the second time onwards.

First, regarding AN, square it in order as a, a'', a', etc. When N is exactly the first power of 2, log
It is sufficient to perform 2N multiplications; if not, find the largest number K of the first power of 2 that does not exceed N, find aK, and then multiply a by N- times. To find K, k, = [
Find k that is log2N ([ ] is the Gagara symbol), and k
= 2'. For example, when N=10, =8, and first a2. a4. Find a8, then 810
=ag・a'a. However, for the part where a is multiplied twice, a2 may be used, and generally when calculating aN-, a2. a4. '...aK2 can be used.

Next, M□2M2. ...Explain how to find MN.

After multiplying Mo by a to create M, proceed as follows.

1) Multiply M□ by a to find M2.

2) Multiply Mo, M2 by a2 to find M,,M.

3) Ml, M2. -M, multiplied by a4, M, ・M.

seek.

In the same way, up to Ml, up to M3□, and so on.
Find up to K. Mxya□,..., MN is M n =
a'-M n-K...
(5) However, it can be calculated in parallel using n = K + 1, . . . , N. In the above calculation, sign bit processing like blocks 4 and 5 in FIG. 5 is required every time multiplication is performed.

As shown in formula (3)', Mn and Mn+N (n is any positive integer) have the following formula: Mn+w= a rl-MN (mod L)
-(6) holds true. Therefore, M-engineer 7
M2. ....

a, a2. instead of MN. ..., aN, then N random numbers Mn ya 1. following n are stored. ..., Mn+N can be generated in parallel.

The 6th @ is a PAD diagram showing a random number generation method according to this idea. Blocks 21 to 27 in the same figure are initialization processing, in which a, a'', .
A2. ..., stored in AN. However, Ao=1
shall be. Also, find Mn every time A is found (M, is
Since there is no need to store it, it is represented by a temporary variable t. )
, and then normalize it to obtain a uniform random number R. The closing price of t is
M, and can be carried over to the second and subsequent calculations. In the same figure, blocks 28 to 32 are used for the second and subsequent times, that is, A 1g A 21.
・This is the procedure when the AN has already been created.

The method shown in FIG. 6 is superior to the method shown in FIG. 5 in that it can be changed arbitrarily as long as the number of random numbers generated from the second time onwards is N or less.

Both of the methods shown in FIGS. 5 and 6 employ a method as defined in which when M and l are negative numbers, the sign bit is ignored (if it is 1, it is replaced with O).

Conventionally, other convenient methods have been used, such as a method of taking an absolute value as shown in FIG. 3, and a method of not using 1 as a random number in the case of a negative value. In any case, the sign bit is 1
Special handling is required for these cases.

By the way, if we consider that Mn represents one integer including the sign bit, its range is [-ri.

It is L-1. Therefore, if Mrl is divided by 2L=2"1, -0,5<Mn/2L<0.5. However,
If you choose an appropriate prime number for Mo, a, Mn=-L will never hold, so in reality -0,5<Mn/2L<0.5
It is. (Mll) is replaced with a random number sequence (M
' Make n).

M'p=a'M'n-0n-0(2L)
・= (7) However, M′. =MO L Rl is a uniform random number of (-0, 5, 0, 5) for the above-mentioned reason. Furthermore, R', = o, E5 + R',
-(9), R'n becomes a uniform random number (0, 1).

FIG. 7 shows a uniform random number generation procedure using this numerical expression. The code focus processing that was necessary in blocks 4, 5, and 10 in FIG. 5 is no longer necessary in blocks 44, 45, and 50 in FIG. 7. Instead, an addition of 0.5 is required in block 46.51. Generally, the method shown in FIG. 7 requires less calculation load. The randomness of this method becomes noticeable when uniform random numbers are used to obtain random numbers that follow various distributions.

The above RIo or RITo indicates a random number generation method according to various distributions using a generation method.

(1) A pair of normal random numbers η3'1.η is created from a pair of (0, L) uniform random numbers ξ31' and ξ5''.

Here, Rln is used as a uniform random number. That is, lξ”
= R' n+ to R' n+, +0.5 - (13)
shall be.

cos, 2xξ" = cos(2x R' n+□+
π) = -cos2πR'n+1 Similarly, Therefore, instead of creating a normal random number and taking its cosine and sine, R#', = -M' +1-(16) The uniformity of (-π, π) becomes Create a random number and its cosine.

Just take the sine.

On the other hand, for the terms including η3"1 in equations (10) and (11), equation (12) is used. RJ n, R"i
l and formula (10).

(11) A method for processing the evaluation of the expression within the same iteration loop is shown in Figure 8.9 In the same figure, AN=a' (same as in Figure 1)
ts, tz, +3 are −temporal variables, α=1/2L.

β=π/L, N is an even number.

Block 61 executes the processing of blocks 62 to 68 N/2 times. Represents loop control. n is 1, . 3. ．． 5
．． , N-1.し、はR″.

t2 corresponds to R”n, . Equation (14), (1,
5) Using equations, find the values of equations (10) and (11), and obtain a pair of normal random numbers η. ,η. Ten is found. Since all operations within this loop can be vector processed, a total of N normal random numbers are generated in parallel.

In this way, uniform random numbers and regular random numbers can be created in the same loop,
By doing this, there is a part where there is no need to store uniform random numbers in memory and where there is no need to standardize them to (0, 1). By reducing the load on loop control, faster processing than before is possible.

(2) Triangular random number The random number that follows the triangular distribution f(x) = 1-1 x 1 (-1<x<1) is 1. - a pair of uniform random numbers ξ11. Make it from ξ32゜ as follows.

η=ξ −ξ ・・・(]+
7ξ11', ξ01, R', (Formula (19))
If we decide to use η"R' n R" n+.

"R' n -"R' n +,
・(+8), so (-0, 5, 0 in equation (8)
, 5) Uniform random numbers can be used as is. That is, standardization of equation (9) becomes unnecessary. As in the case of normal random numbers, it is possible to improve the efficiency of making dinner (η.) within the same loop as with uniform random numbers.

(3) Even in calculations according to other distributions, normalization processing can be eliminated or reduced.

a) A random number following the distribution fk(x)=kx” (0<x<1) is η=max(ξ3゛3.ξi21...ξ(kl),
, , j 19). In this case too, R'' is taken as ξ'Jl, and R' n = 0, 5 + R'
. Using the property, 'Q = max(R' fl
l R' n+0. -, R'-ten-t
)+0.5...(],,+9', and the load on normalization processing (addition of 0 and 5) is reduced.

b) gk(x)=k(,1-x)'-"(0<k<1
) is a random number that follows the distribution η=min(ξ”' * ξ”Z , ξ”')
It can be found from (20). In this case as well, R
If we take ``0, η=min(R' ny R' n+, t ``'R'
1 for n+)+0.5''(20)', and the load on the standardization process is reduced.

C) There is a method of approximating a -dimensional normal distribution (mean 09 variance 1). In this case, if R′o is taken as ξ(gl), then the normalization process becomes completely unnecessary.

〔Effect of the invention〕

According to the present invention, since a large number of arithmetic random numbers can be generated at once, the load of calling a random number subprogram can be relatively reduced. Furthermore, since the random number calculation formula is not a linear cyclic type, processing that takes advantage of the characteristics of vector computers and parallel computers is possible.

According to the method of the other embodiments, there is no need for calculation load for removing the sign bit, and the value of the modulus can be increased (twice as much as that of the conventional method), so there is an advantage that the period becomes longer. Furthermore, when generating arithmetic random numbers according to various distributions, the uniform random numbers that are the basis thereof can be used without being normalized to (0, 1), making high-speed processing possible.

Furthermore, by configuring uniform random number calculation and random number calculation according to a functional distribution in the same loop, the loop control load is reduced, and the storage of uniform random numbers in memory can be avoided, reducing processing time. Ru.

[Brief explanation of the drawing]

Figure 1 is a PAD diagram showing the basic procedure of the uniform random number generation method of the present invention, Figure 2 is a diagram showing an example of a conventional uniform random number generation program, and Figure 3 is the internal configuration of a vector computer. FIG. 4 is a diagram showing an example of the program of the present invention shown in FIG. 2, FIGS. 5 and 6 are PAD diagrams showing the procedure of an embodiment of the present invention, and FIG. PAD diagram showing the uniform random number generation method when changing the numerical expression, Figure 8 generates uniform random numbers.
95M(), ((T, θ5ndingHEN
/lρZ(I) = ArX(D + Y(1)
t2(ILL S E
t3ρ, Z(In= ρρ
14ρE
NDrF t5ρ
100 C? 1NTINuE
161) R.J.
NC7IγN RANf)'ijM()
200DATA K /q7 Otsu M /13/
'Z/θ=MtK
22θT? = FL7j
Ar(K)/I-L2IAsu (2181)
230RANρiTM=IIBS(R)
240 RETUR
N250E
ND26
θth 3rd stage /i Figure VO/θθ enil, /#
3θθ1, C=to'*M
31θR='rLiAT(K)/r
L11T (21731) 3zOIF (A
BS(R), f(7θ5) THEN
33θK(1)=A! X(J)suY(1)
3i. ELSE35ρ z(1)2ρρ Ground ball θE
Nf) IF 3
7ρ100 C0NTINLAE
3231) Figure 5 v t Figure 7 Figure 3

Claims (1)

[Claims] 1. In a computer that has the function of generating a large number of uniform random numbers using the multiplication product congruence method, by storing N multipliers raised to the power of ■ or N pre-created random numbers. , N random number generation calculations are made mutually independent and vector processing is performed. 2. A method for generating arithmetic random numbers, characterized by generating random numbers using a numerical representation including a sign bit in a computer that calculates a remainder using an overflow tolerance function during calculations for integer type variables. 3. Arithmetic in a computer according to claim 1 or 2, characterized in that the processes of generating uniform random numbers and using the random numbers to generate arithmetic random numbers according to a given distribution are performed within the same iterative loop. How to generate random numbers.