Method for performing high-speed monte carlo random number generation, the method enabling reproducibility and portability by using only double-precision integer operations to generate random numbers based on multiplicative congruential method whose modulus is a product of two or more pairwise coprime positive odd single-precision integers utilizing the structure of the chinese remainder theorem, and enabling random-number output computation using only real double-precision operations

Abstract

To provide a method for rapidly generating double-precision real random numbers using only double-precision integer operations and double-precision real operations.SOLUTION: The method uses a positive integer d=e1*e2*...*em as a modulus, where m is an integer of 2 or more, and e1, e2, ..., em are pairwise coprime positive single-precision integers. A positive multiplier z<d, coprime with d, and a positive seed n<d, also coprime with d, are used to form a multiplicative congruential (MC) random-number generation mechanism (d, z, n). Starting from the given initial value n, the method produces a sequence of positive integers {x(t) | t=0, 1, 2, ...}, each smaller than d, by the recurrence x(0)=mod(n,d) and x(t+1)=mod(z*x(t),d) for t=0, 1, ... The uniform and independent random-number sequence {v(t):=x(t) / d | t=0, 1, 2, ...} is obtained in this manner.SELECTED DRAWING: None

Description

[Technical Field]

[0001] Positive odd numbers that are relatively prime (have no common prime factors and have a greatest common divisor of 1) m sets e1, e2, ..., e m The product d, d=e1e2…e m , ($) A multiplicative congruential (MC) random number generator (d, z, n) modulo d>0, where d is a prime integer z less than d and called the "multiplier", and d is a prime integer n less than d and called the "initial value or seed", is given by a multiplicative congruential (MC) random number generator. The MC random number generator is given by a recursive congruence relation modulo d for a sequence of integers {z0, z1, z2, ...} in the open interval (0, d). z0=mod(n, d), z k+1 =mod(z×z k ,d),k=0,1,2,… (￥) The well-known definition of the mod function is clearly defined immediately below, and the multiplication symbol × will be omitted from now on. The invention described here points out the important structure that the calculation of modulus d in the form of the congruential arithmetic (￥) always used in MC methods has under Sun Tzu's theorem. The necessity of this usage method is related to the nature of prime numbers in random number problems. Since the random number sequence required for simulation problems is a sample process, statistical uniform independence is a matter of "appearing that way", a matter of hypothesis testing. A full-scale MC random number generation method is based on the famous prime numbers of Fishman and Moore, and the Mersenne prime 2, which is close to the limit of single integer precision. 31 The tests they conducted show that the prime numbers that can give the modulus of an MC random number generator that appears to have good performance are within the integer single precision limit of 2. 32Experience has shown that a size close to d is required. Please see (Reference 1) for details. If the size is small, it is difficult to find random numbers that are deemed "qualified" in testing. Furthermore, simulation experience has shown that sufficient precision can only be achieved using "double-precision real number calculations," and the random number period must also be appropriately large. The way to create random numbers that meet these requirements is to use a modulo that is the product of two or more prime numbers that are close to the limit of single precision. We will begin by explaining why MC random numbers of the form (￥) require a modulus d that is the product of two or more partial modulos that are close to the limit of single precision, due to the requirements of simulation. (Reference 1) Naoya Nakazawa / Hiroshi Nakazawa: “Random Number Generator on Computers”, JENNY STANFORD, proof copy attached. Coming soon.

[0002] What is the feasibility of technical design for real double precision MC random number generator? (1) For one large odd prime modulo d=p, (2) The product of m distinct large odd prime numbers, d=p1p2…p m in the case of, As an easy case to realize the calculation, in the above (2), in the case of the modulo d of the product of odd prime numbers where m=2, odd prime numbers p1 and p2, which are single precision integers, satisfy the condition "double the product, 2d=2p1p2, is the double precision integer 2 64 The case where the number of products m is 2 or more is an arbitrary case. "The mutually prime submodulo e1, e2, ..., e m There are different single precision There are no conditions other than that they are prime numbers (and therefore mutually prime).” First, we explain why case (1) is not technically suitable. Currently, double-precision real number calculations are often required for computer simulations. To generate double-precision real number random numbers suitable for this purpose using the MC method with a modulus of one prime number d=p, and to make the random number period large enough for normal simulations, we need a size of 2. 64A huge modulus d=p and multiplier z of about double precision integers are required. However, at each stage of the calculation of the MC random number sequence {z0, z1, z2, …}, the congruence recurrence formula of modulus d is z k+1 =mod(z×z k ,d) (*) (The multiplication symbol × will be omitted from now on.) k If is a double-precision integer, this multiplication part must necessarily go through arithmetic using quadruple-precision real numbers, which will result in a significant increase in calculation time (by more than 10 times), and moreover, the maximum period of the modulus of an excellent prime number of the MC random number generator discovered to date (2024) is 2 34 This is too short for practical simulation, and the calculation burden is too great compared to the effect.

[0003] The way to avoid the difficulties pointed out in the previous paragraph is to find a solution where the modulus d is m odd prime numbers "within the single-precision integer limit" (注1) e1, e2, …, e m Product of d=e1e2…e m The goal is to adopt this method, and further to "enable high-speed double-precision integer arithmetic for the single-precision integer partial MC random number generation mechanism provided by each partial modulus, while ensuring both reproducibility and portability." In the previous application, when the double-precision real random number output for simulation satisfies the condition that "m = 2, and twice d = e1e2 remains within the double-precision integer limit," this "high-speed MC real double-precision random number generation mechanism with reproducibility and portability" was obtained by "a method using only double-precision integer arithmetic." Here, we would like to present a method that ensures reproducibility and portability by also using double-precision "real" arithmetic, in a way that "the number of modulus m of the partial single-precision integer has no restrictions."

[0004] First, with a positive integer modulus d and a positive integer A in mind, we will clearly prepare the definition of the mod function mod(A, d), which will be used frequently from now on. (Definition 1) mod(A, d) is defined as "the remainder of the integer between 0 and d obtained by subtracting q times the integer quotient of the positive integer d modulo A from A." In terms of a division equation, A=qd+r, where r is an integer greater than or equal to 0 and less than d. Based on this, we define mod(A, d) = r. (End of Definition 1) In practical computer calculations, there is a strict distinction between single precision, double precision, and quadruple precision for all numbers that can be used, but in discussions, (Definition 1) can be used for any of them. Furthermore, in various computational languages, mod functions have already been implemented with definitions that are defined in the same form as (Definition 1) even when the variable A and modulus d are real numbers. Here, considerations of calculation speed necessitate technology that does not require quadruple precision calculations, but since this will be discussed in the construction of the calculation program, all mod functions will be based on (Definition 1) above. What is important is an understanding of Sun Tzu's theorem, which is perfectly suited to the current problem situation. (Theorem 2. Sun Tzu's Theorem) The product of integers d=e1e2…e where the modulus d>0 is positive and coprime (the greatest common factor is 1). m Let's say that m is 2 or more. "Any positive integer A" is A=A1D1D1 -1 +A2D2D2 -1 +…+A m D m D m -1 , A k =mod(A, e k ), A k is 0 or more e k less than, k=1, 2, …, m (*) It is decomposed into the form of "uniquely determined modulo d". However, D k =d / e k , D k -1 is the law k D in k is the positive reciprocal of mod(D k D k -1 , e j )=δ jk (#) where δ jk is the Kronecker delta, which is 1 if j=k and 0 otherwise. (Proof) e1, e2, …, e m are relatively prime, that is, the greatest common factor is 1, so ek Considering e k and D, which represents the product of the remaining partial methods k =e / e k are mutually prime. Therefore, the conclusion of Euclidean algorithm is that there are integers M and N, and the greatest common factor 1 is given by the following "equation": MD k +Ne k =1, MD k =1-Ne k If we look at this as a congruence equation, MD k is the law k 1 and the integer M is 'mod e k In D k The reciprocal of D k -1 In particular, 0 and above e k Can be less than D k ―1 = D k ―1 = x and then change x to 1, 2, ..., e k -1 and D k Multiply by modulo e k Then we can find the value that is 1. From now on, x=D k -1 is taken as such and the integer D k D k -1 For any positive integer A, (*) the right-hand side mod(A, e k ) and (#) (注1) mod(A, e k )=mod(A k D k D k -1 , e k ) =mod(A k , e k )=A k holds true:A k is 0 or more e k Note that it is defined as follows: k =mod(A, e k ) holds for all k = 1, 2, ..., m. Now, a positive integer B can also be expressed in the same form (*). B=A1D1D1 -1 +A2D2D2 -1 +…+A m D m D m -1 , have, i.e. A k =mod(B, e k ), k=1, 2, …, m If this holds, then mod(B, e k )=mod(A, e k ), k=1, 2, …, m So, if we use the division equation to express this, we have quotients q and q', which are integers greater than or equal to 0. B=qe k +A k , and A=q'e k +A k This is BA=(qq')e k , k=1, 2, …, m holds, so the integers BA are e1, e2, ..., e m are multiples of e1, e2, …, e m are relatively prime, that is, the greatest common factor is 1, so BA is the product e1e2…e m = multiple of d, that is, modulo d, BA is 0. Therefore, an integer A greater than or equal to 0 that has the expression (*) is modulo d=e1e2…e m (Theorem 2. Proof of Sun Tzu's Theorem) Although it is the same thing, the following rephrasing and understanding of Theorem 2 above (Sun Tzu's Theorem) is important. (Corollary 3. Sun Tzu's Theorem) The product of integers d=e1e2…e where the modulus d>0 is positive and coprime (the greatest common factor is 1). m Let m be 2 or more. In this case, any positive integer A between 0 and d is modulo e k 0 or more in e k Less than A k =mod(A, e k ), uniquely by a linear combination of k=1, 2, …, m A=mod(A1D1D1 -1 +A2D2D2 -1 +…+Am D m D m -1 , d). This conclusion, Sun Tzu's theorem, may be accepted in theory, but it is difficult to accept in practice. Later, at the end of Figure 1, I would like you to try it with an arbitrarily chosen integer A, using an example where m=3 and e1, e2, and e3 are specifically chosen, to see for yourself how it works. [注1) Note the mod functions for addition and multiplication of integer variables X and Y greater than or equal to 0. (1) d is a positive integer. The mod function mod is the sum of integers X, Y, ... (X+Y, d) is mod(X+Y, d) =mod(mod(X, d)+mod(Y, d), d) is the exact formula. If we think of it as a division equation, let the integer quotients of the division be q and q'. X=qd+mod(X, d), Y=q'd+mod(Y, d), mod(X+Y, d)=mod(mod(X, d)+mod(Y, d), d), is clear; all integer multiples of modulo d may be dropped, but since mod(X, d) + mod(Y, d) may exceed d, an outer "final mod function" is generally required. (2) Similarly, if the variables of the mod function are the product of integers X and Y greater than or equal to 0, mod(XY, d)=mod(mod(X, d)×mod(Y, d), d) Note that the equation X=qd+mod(X, d), Y=q'd+mod(Y, d) holds. q and q' are integer quotients in division. Therefore, XY = integer multiple of d + mod(X, d) × mod(Y, d) and mod(XY, d) =mod(mod(X, d)×mod(Y, d), d) Since mod(X, d) × mod(Y, d) may be larger than d, the outermost mod function is generally necessary.

[0005] Product d=e1e2…e of m mutually prime positive partial modulo m Consider the structure of MC random numbers modulo . In order to express the partial modulus, we generally define the MC random number sequence generated by (d, z, n) as {z (t) =mod(nz t , d)| t=0,1,2,…,} Since j is difficult to read, we generally use t to represent time. This can also be expressed in the form of a recurrence formula: z (0) =mod(n, d), z (t+1) =mod(zz (t) , d), t=0, 1, 2, …(#1 * ) The form d=e1e2…e is written as follows: m When you have Partial multiplier, partial starting value Partial multiplier: {z k = mod(z, e k ), k=1, 2, …, m}, Partial initial value: {n k = mod(n, e k ), k=1, 2, …, m}. These are defined as "part e k MC random number generation mechanism {(e k , z k , n k )|k=1, 2, …, m} and the "MC partial integer random number sequence" generated by them is written as {mod(n k (z k ) t , e k ), t=0,1,2,…, k=1, 2, …, m}(#2) Or in recurrence form: z k (0) =mod(n k , e k ), z k (t+1) =mod(z k zk (t) , e k ), t=0, 1, 2,…, k=1, 2,…, m(#2 * ) The problem is (#1 * ) and (#2 * ) The following holds: (Theorem 4) The modulus d of a positive integer is a product of m positive integer submodulo d=e1e2…e, where d is a positive integer modulus that is relatively prime. m When (e k , z k , n k ) is the MC generating mechanism for the partial MC integer random number sequence (#2 * ) and mod(z (t) , d)=mod(mod(z1 (t) , e1)D1D1 -1 + mod(z2 (t) , e2)D2D2 -1 +…+ mod(z m (t) , e m )D m D m D -1 , d), (#2 * ) The MC random number generator (d, z, n) is z (t) =mod(n(z) t ,d), The actual generation procedure is in the form of a recurrence formula. z (0) =mod(n, d) z (t+1) =mod(zz (t) ,d), t=0,1,2,… The output as a uniform random number is a rational number sequence {x (t) =mod(z (t) , d) / d, t=0, 1, 2, …} is. (Proof) (#1 in recurrence form) * )z (t) In each submethod, each e kcomponent mod(z (t) D k D k -1 , e k ), k=1, 2, …, m These are used in (*) to obtain the original (d, z, n) MC random number output. {z (t) , x (t) =z (t) / d|t=0, 1, 2, …} is a unique construction modulo d by Sun Tzu's theorem. (End of proof)

[0007] The first key to generating (d, z, n) MC random numbers is the modulus d=e1e2…e m odd partial modulo e1, e2, …, e m The key is to select the integers with single precision so that they are relatively prime. This allows us to create a "partial MC random number generator" {(e k , z k , n k )|k=1,…,m} There is no limit to the number m of d'. In this case, partial random numbers can be generated in a reproducible and portable manner using only double-precision integer arithmetic. The synthesis of these into MC random number outputs (d, z, n) using Sun Tzu's theorem can be performed as double-precision real numbers, since there is no need to consider reproducibility, and this reduces some of the difficulties and can be achieved quickly. This is the gist of the current application. However, there is still a limit to the size of the modulus d, and realistically, the possible limit for the number m would be 2 or 3. For example, even if calculations are performed as quadruple-precision real numbers, d is 2 128 It is not possible to exceed this limit, and there is a concern that the size may reach the limit during calculation if used close to this limit. It is difficult to use the full standard even for real numbers. However, the recurrence formula z×z (t) The possibility of performing this in double-precision integers using the small part method is certainly expanded. In the next paragraph, we will show an example of the actual program for the case of m=3, and clearly show which calculation steps are new innovations.

[0008] A random number MC calculation program in FORTRAN format using the modulus d=e1e2e3, which uses three different prime numbers e1, e2, and e3, is shown in Figure 1 on the next page. However, please note that the statistical performance of this MC random number is not excellent. Aiming for the most understandable notation from a current perspective, the program details the size of the variables, such as simple, double, and quadruple. First, we will measure the generation speed. Figure 1. The product of m=3 prime numbers modulo d=e1e2e3 gives 10 million The main program that generates random numbers program main implicit integer * 8(ip),real * 8(ah,rz) common p1,p2,p3,d,iz1,iz2,iz3,d1,d2,d3,invd1,invd2,invd3 p1=134265023! p1 is approximately 2 ∧ (27) prime single-precision integer p2=1000919 ! p2 is approximately 2 ∧ (20) prime single-precision integer p3=257 ! p3 is approximately 2 ∧ (8) Prime single-precision integer num=10000000! Set the number of random numbers to generate !d:=p1 * p2 * p3 ! modulo d is approximately 2 ∧ Integer (54.94) d1=p2 * p3 ! D1:=d / p1(double precision real) d2=p3 * p1 ! D2:=d / p2(double precision real) d3=p1 * p2 ! D3:=d / p3(double precision real) d=p1 * p2 * p3 ! d=p1 * p2 * p3 (double precision real number) iz1=19061252!Primitive root iz1 of p1 iz2=64545 !Primitive root iz2 of p2 iz3=53 !primitive root iz3 of p3 invd1=38789367 !invd1:=mod(D1 ∧ (-1),p1) invd2=388498 !invd2:=mod(D2 ∧ (-1),p2) invd3=83 !invd3:=mod(D3 ∧ (-1), p3) n1=1 ! Partial MC initial value at modulo p1 n2=1 ! Partial MC initial value at modulo p2 n3=1 ! Partial MC initial value at modulo p3 mz1=mod(n1,p1) ! Initial setting of partial MC random number sequence with modulus p1 mz2=mod(n2,p2) ! Initial setting of partial MC random number sequence with modulus p2 mz3=mod(n3,p3) ! Initial setting of partial MC random number sequence with modulus p3 !random subroutine is called num=10 million times to generate and output random numbers. do i=1,num call random(mz1,mz2,mz3,rand)!random subroutine (Figure 2) end do end (End of Figure 1) This main program is used to calculate and fix the necessary values ​​of the variables before generating random numbers (d, z, n). The modulus d=e1e2e3 is 2 55 I will use this opportunity to give you an actual example of Sun Tzu's theorem, so I hope you will understand it. The modulus d is the same as in Figure 1. d=e1e2e3, e1=134265023, e2=1000919, e3=257 A positive integer A=195090965842899 can be decomposed as follows: A=mod(A1D1D1 -1 +A2D2D2 -1 +A3D3D3 -1 ,d), A1=19061252, A2=64545, A3=53, D1=e2e3=257236183, D2=e3e1=34506110911, D3=e1e2=134388412556137, D1 -1 =38789367, D2 -1 =388498, D3 -1 =83, mod(D1D1 -1 , e1)=1, mod(D2D2 -1 , e2)=1, mod(D3D3 -1 , e3)=1. The subroutine in Figure 2 below is related to calculation speed. In searches using two partial MC random numbers, it has been experienced that good combinations cannot be found unless they are "close to the same single-precision integer limit." In future practical calculations, we must aim to discover more prolific combinations. Figure 2. Random number generation subroutine using Sun Tzu's theorem subroutine random(mz1,mz2,mz3,rand) implicit integer * 8(ip),real * 8(ah,qz) common p1,p2,p3,d,iz1,iz2,iz3,d1,d2,d3,invd1,invd2,invd3 !d1:=D1,d2:=D2,d3:=D3 (double precision real number) ! Arguments mz1, mz2, and mz3 are updated with each random number generation and are passed on to the next mz1=mod(mz1 * iz1,p1) ! MC partial integer random number sequence with modulo p1 (double precision integer) mz2=mod(mz2 * iz2,p2) ! MC partial integer random number sequence with modulo p2 (double precision integer) mz3=mod(mz3 * iz3,p3) ! MC partial integer random number sequence with modulo p3 (double precision integer) !(The following process will process the MC random number output with composite number d) rz1=mod(mz1 * invd1,p1)!rz1:=mod(mz1 * D1 ∧ (-1),p1) (double precision real number) rz2=mod(mz2 * invd2,p2)!rz2:=mod(mz2 * D2 ∧ (-1),p2) (double precision real number) rz3=mod(mz3 * invd3,p3)!rz3:=mod(mz3 * D3 ∧ (-1),p3) (double precision real number) az=d1 * rz1+d2 * z2+d3 * rz3 ! Composition of az by Sun Tzu's theorem !az:=D1 * D1 ∧ (-1) * rz1+D2 * D2 ∧ (-1) * rz2+D3 * D3 ∧ (-1) * rz3 rand=mod(az / d,1d0) ! Output: Take the fractional part of az / d (double precision real number) return end (Figure 2. End of subroutine) This random number generation subroutine requires a CPU time of 0.33 to 0.34 seconds per 10 million numbers. While this is longer (approximately 1.2 to 1.3 times) than the CPU time for the relatively prime modulus m=2 described in Patent Application No. 2022-024277, it is not so short that it warrants criticism. However, the CPU time will vary depending on the size of the three prime numbers. This will be precisely understood once a MC generation mechanism that passes testing is obtained. We are optimistic that a modulus with three prime numbers will have more diverse combinations than one with two prime numbers, and that there will be more passing attempts. The number of passing attempts and the CPU time will all be accurately understood once various prime numbers have been tested.

[0009] Figure 3 below shows a program for a primitive method of performing direct multiplication using quadruple precision real arithmetic without decomposing into submodules, using the same modulus d=e1e2e3 as in the previous paragraph

[0008] . Figure 3. Program for the quadruple precision real multiplication method modulo d program main implicit interger * 8(ip),real * 8(ah,rz),real * 16(q) num=10000000! Set the number of random numbers to generate p1=134265023! p1 is approximately 2 ∧ (27) prime single-precision integer p2=1000919 ! p2 is approximately 2 ∧ (20) prime single-precision integer p3=257 ! p3 is approximately 2 ∧ (8) Prime single-precision integer !d:=p1 * p2 * p3 ! modulo d is approximately 2 ∧ Integer (54.94) qd=p1 * p2 * p3 ! qd=p1 * p2 * p3 (quad-precision real number) qz=1950909605842899d0!qz: multiplier (quadruple precision real number) qmz=1 ! MC random number sequence (quadruple precision real number) with modulus d, initial value is 1 ! Call random numbers num=10 million times and generate and output them do i=1,num qmz=mod(qmz * qz,qd)! MC random number sequence modulo d (quadruple precision real number) rand=qmz / qd !output end do end (Figure 3. End) It can be confirmed that this program produces the same output results as those shown in Figures 1 and 2 above. Figure 3 simply repeats multiplication modulo d, so the program structure appears simple at first glance. However, the CPU time is approximately 2.7 seconds per 10 million operations, nearly 10 times longer than the subroutine in Figure 2. This is one important reason why the generation method presented in this invention, which uses the structure of Sun Tzu's Theorem and is based only on double-precision integers and double-precision real numbers, should be put into practical use. Furthermore, the generation method of this invention uses 2 64 It is also possible to output MC random number sequences with a modulus d of a composite number exceeding 1 using only double-precision real number arithmetic, which has great potential and development potential compared to methods that simply use quadruple-precision multiplication.

[0010] Here is an example program with the intention of actually using random numbers in software for simulations, games, etc. Figure 4 below is a basic and simple main program that generates a specified number of random numbers. Figure 4. Main program that generates num random numbers program main implicit integer * 8(ip),real * 8(ah,rz) common p1,p2,p3,d,iz1,iz2,iz3,d1,d2,d3,invd1,invd2,invd3 num=10 ! Number of random numbers to generate p1=134265023! p1 is approximately 2 ∧ (27) prime single-precision integer p2=1000919 ! p2 is approximately 2 ∧ (20) prime single-precision integer p3=257 ! p3 is approximately 2 ∧ (8) Prime single-precision integer d=p1 * p2 * p3 ! modulo d is approximately 2 ∧ Integer (54.94) d1=p2 * p3 ! D1:=d / p1(double precision real) d2=p3 * p1 ! D2:=d / p2(double precision real) d3=p1 * p2 ! D3:=d / p3(double precision real) iz1=19061252!Primitive root of p1 iz2=64545 !Primitive root of p2 iz3=53 !primitive root of p3 invd1=38789367 !invd1:=mod(D1 ∧ (-1),ip1) invd2=388498 !invd2:=mod(D2 ∧ (-1),ip2) invd3=83 !invd3:=mod(D3 ∧ (-1),ip3) n1=10 ! Partial MC initial value at modulo p1 n2=22 ! Partial MC initial value at modulo p2 n3=41 ! Partial MC initial value at modulo p3 mz1=mod(n1,p1)! Initialization of partial MC random number sequence with modulo p1 mz2=mod(n2,p2)! Initialization of partial MC random number sequence with modulo p2 mz3=mod(n3,p3)! Initialization of partial MC random number sequence with modulo p3 !The double-precision integers mz1, mz2, and mz3 are updated every time a random number is generated. Call the !randgenerate subroutine to generate num random numbers. !Arguments mz1, mz2, and mz3 store the num-th partial MC random numbers (double precision integers). !Can be carried over to the next random number generation call randgenerate(mz1,mz2,mz3,num)! randgenerate subroutine (Figure 5) end (Figure 4. End) The basic structure is the same as the main program in Figure 2 in the previous paragraph. However, assuming the development of programs such as simulations and games, it has been modified so that the "required number" of random numbers can be generated as a subroutine at any time. Figure 5 shows a subroutine that generates an arbitrary number of random numbers, num. Figure 5. Subroutine that generates num random numbers subroutine randgenerate(mz1,mz2,mz3,num) implicit integer * 8(ip),real * 8(ah,rz) !Generate random numbers num times Rewrite the partial modulus multipliers mz1, mz2, and mz3 each time you call !random do i=1,num call random(mz1,mz2,mz3,rand)!random subroutine uses Figure 2 end do return end (Figure 5. End) The modulus d=p1 required for the MC random number generator used * p2 * p3 *If the numerical values ​​and initial values ​​for the primitive roots iz1, iz2, and iz3 are given in the main program (Figure 4), then no matter how many times the randgenerate subroutine (Figure 5) is called within the same program, MC random numbers for the total number of steps (elapsed time) of the argument num will be generated without gaps (skips) or overlaps. If you customize the program or subroutine as needed to save the latest information on the partial modulus MC integer random numbers mz1, mz2, and mz3 to a storage medium, and when you interrupt or terminate, set mz1, mz2, and mz3 to their "initial values" in another program, it will be possible to inherit the initial conditions from the "previous program" on the same MC random number generator and continue generating without gaps or overlaps.

Claims

[Claim 1] The mod function mod(A, d) with a positive modulus d is given by mod(A, d) = 'The integer remainder r when A is divided by d, r is greater than or equal to 0 and less than d. and there are two or more positive integers e that do not have m common prime factors, i.e., they are mutually prime. １ , e ２ , ..., e ｍ A positive integer d = e １ e ２ …e ｍ Modulo 'e １ , e ２ , ..., e ｍ is a single-precision integer, and d is a prime positive multiplier z<d, d is a prime positive seed n<d, and In the generation mechanism (d, z, n), a sequence of positive integers less than d obtained sequentially from a given initial value n by a recursive congruence formula {x （ｔ） |t=0, 1, 2, ...}, where x （０） =mod(n、d)、 x （ｔ＋１） =mod(zx （ｔ） 、d)、t=0、1、…、 A uniformly independent random sequence {v （ｔ） :=x （ｔ） / d|t=0, 1, 2, ...}, a form of Sun Tzu's theorem, partial method e １ , e ２ , ..., e ｍ Partial MC (e ｋ , z ｋ , n ｋ ) generation mechanism, i.e., partial recursive congruence x ｋ （０） =mod(n ｋ 、e ｋ )、x ｋ （ｔ＋１） =mod(z ｋ x ｋ （ｔ） 、e ｋ )、 k=1, 2,..., m, t=0, 1,..., z ｋ :=mod(z,e ｋ )、k=1、2、…、m、 n ｋ :=mod(n,e ｋ )、k=1、2、…、m、 Partial MC random number sequence given by {x ｋ （ｔ） |k=1、2、…m、t=0、1、2、…} Using the mutually prime submodulus e １ , e ２ , ..., e ｍ m operator integers defined by {D ｋ D ｋ －１ |D ｋ = d / e ｋ , D ｋ －１ is law e ｋ D in ｋ The reciprocal of That is, mod(D ｋ D ｋ －１ , e ｊ ) is 1 if j = k, If j is different from k, it becomes 0. j, k=1, 2,..., m} Linear combination of m integer operators using x （ｔ） =mod(x １ （ｔ） D １ D １ －１ + +x ２ （ｔ） D ２ D ２ －１ +…+x ｍ （ｔ） D ｍ D ｍ －１ 、d) Let t be a sequence of integers {x （ｔ） |t=0, 1, 2, ...}, MC uniform random number output is real number {v （ｔ） =x （ｔ） / d|t=0、1、2、…} A specific example of a uniform random number sequence generation method for m=3 in FORTRAN format is: Figure 6. Modulus d = e １ e ２ e ３ Random number generation main program program main implicit integer ＊ 8(i-p),real ＊ 8(a-h,r-z) common p1, p2, p3, d, iz1, iz2, iz3, d1, d2, d3, invd1, invd2, invd3 p1 = 134265023! p1 is approximately 2 ∧ Prime single-precision integer of (27) p2 = 1000919! p2 is approximately 2 ∧ Prime single-precision integer of (20) p3 = 257! p3 is approximately 2 ∧ (8) Prime single-precision integer num=10000000! Set the number of random numbers to generate !d:=p1 ＊ p2 ＊ p3 ! Modulo d is approximately 2 ∧ Integer (54.94) d1=p2 ＊ p3! D1:=d / p1 (double precision real number) d2=p3 ＊ p1! D2:=d / p2 (double precision real number) d3=p1 ＊ p2! D3:=d / p3 (double precision real number) d=p1 ＊ p2 ＊ p3 ! d = p1 ＊ p2 ＊ p3 (double precision real number) iz1 = 19061252! Primitive root iz1 of p1 iz2 = 64545! Primitive root iz2 of p2 iz3 = 53! Primitive root iz3 of p3 invd1=38789367.invd1:=aod(D1 ∧ (-1),1) invd2=388498 !invd2:=mod(D2 ∧ (-1) ),p2) P ∧ ()),) n1 = 1 ! Partial MC initial value in modulo p1 n2 = 1 ! Partial MC initial value in modulo p2 n3 = 1 ! Partial MC initial value in modulo p3 mz1 = mod(n1, p1) ! Initialization of partial MC random number sequence with modulo p1 mz2 = mod(n2, p2) ! Initialization of partial MC random number sequence modulo p2 mz3 = mod(n3, p3) ! Initialization of partial MC random number sequence modulo p3 Call the random subroutine num = 10 million times to generate and output random numbers. do i=1,num call random (mz1, mz2, mz3, rand) end do end (End of Figure 6) and Figure 7. Modulus d = e １ e ２ e ３ Random number generation in subroutnine subroutine random (mz1, mz2, mz3, rand) implicit integer ＊ 8(i-p),real ＊ 8(a-h,q-z) common p1, p2, p3, d, iz1, iz2, iz3, d1, d2, d3, invd1, invd2, invd3 ! rd1:=D1, rd2:=D2, rd3:=D3 (double precision real number) mz1=mod(mz1 ＊ iz1, p1) ! MC partial integer random number sequence (double precision integer) modulo p1 mz2=mod(mz2 ＊ iz2, p2) ! MC partial integer random number sequence (double precision integer) modulo p2 mz3=mod(mz3 ＊ iz3, p3) ! MC partial integer random number sequence (double precision integer) modulo p3 rz1=mod(mz1 ＊ invd1, p1)! rz1:=mod(mz1 ＊ D1 ∧ (-1), p1) (double precision real number) rz2=mod(mz2 ＊ invd2, p2)! rz2:=mod(mz2 ＊ D2 ∧ (-1), p2) (double precision real number) rz3=mod(mz3 ＊ invd3, p3)! rz3:=mod(mz3 ＊ D3 ∧ (-1), p3) (double precision real number) a: ＊ ______ ＊ :#) ＊ r: #::#. ＊ .. ∧ ()) ＊ ______________ ＊ .. ∧ ()) ＊ _________________ ＊ ︓ ∧ ()) ＊ r: rand = mod(az / d, 1d0) ! output (double precision real number) return end (Figure 7. End) m=2 or more relatively prime positive single precision integer submodulus e using the following arithmetic: １ , e ２ , ..., e ｍ The product of modulo d = e １ e ２ …e ｍ This method calculates m single-precision multiplication congruential (MC) methods using , and linear combinations (shuffling, mixing) according to Sun Tzu's theorem as double-precision real numbers, thereby quickly generating double-precision real random numbers using only double-precision integer and double-precision real number arithmetic.