知財求人 - 知財ポータルサイト「IP Force」
2023-111796２個の正で互いに素な単精度整数の積を法とする乗算合同法乱数の整数倍精度と実数倍精度のみでの計算方法
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023111796
(43)【公開日】2023-08-10
(54)【発明の名称】２個の正で互いに素な単精度整数の積を法とする乗算合同法乱数の整数倍精度と実数倍精度のみでの計算方法
(51)【国際特許分類】
G06F 7/58 20060101AFI20230803BHJP
【ＦＩ】
G06F7/58 660
【審査請求】未請求
【請求項の数】1
【出願形態】書面
(21)【出願番号】P 2022024277
(22)【出願日】2022-01-31
(71)【出願人】
【識別番号】508260898
【氏名又は名称】中澤 直也
(72)【発明者】
【氏名】中澤 直也
(72)【発明者】
【氏名】中澤 宏
(57)【要約】 （修正有）
【要約】
【課題】倍精度整数と倍精度実数の演算だけで乱数を生成する計算方法を提供する。
【解決手段】方法は、２つの共通素因数を持たない、互いに素な、正の整数ｅとｆの積であるｄ＝ｅｆを法とし、ｄとは素な乗数ｚ、ｄとは素なｓｅｅｄ種ｎで構成される乗算合同法乱数（ｄ、ｚ、ｎ）において、整数の任意に与えられた初期値ｎから順次得られる漸化合同式による一様独立乱数｛ｖｋ：＝ｚｋ／ｄ｜ｋ＝０、１、２、…｝を得るための正の整数列｛ｚｋ｜ｋ＝０、１、２、…｝、ｚ０＝ｍｏｄ（ｎ、ｄ）、ｚｋ＋１＝ｍｏｄ（ｚｚｋ、ｄ）、ｋ＝０、１、２、…、を得る。
【選択図】なし
【課題】倍精度整数と倍精度実数の演算だけで乱数を生成する計算方法を提供する。
【解決手段】方法は、２つの共通素因数を持たない、互いに素な、正の整数ｅとｆの積であるｄ＝ｅｆを法とし、ｄとは素な乗数ｚ、ｄとは素なｓｅｅｄ種ｎで構成される乗算合同法乱数（ｄ、ｚ、ｎ）において、整数の任意に与えられた初期値ｎから順次得られる漸化合同式による一様独立乱数｛ｖｋ：＝ｚｋ／ｄ｜ｋ＝０、１、２、…｝を得るための正の整数列｛ｚｋ｜ｋ＝０、１、２、…｝、ｚ０＝ｍｏｄ（ｎ、ｄ）、ｚｋ＋１＝ｍｏｄ（ｚｚｋ、ｄ）、ｋ＝０、１、２、…、を得る。
【選択図】なし
【特許請求の範囲】
【請求項1】
２つの共通素因数を持たない、即ち互いに素な、正の整数ｅとｆの積であるｄ＝ｅｆを法とし、ｄとは素な乗数ｚ、ｄとは素なｓｅｅｄ種ｎで構成される乗算合同法乱数（ｄ、ｚ、ｎ）において、整数の任意に与えられた初期値ｎから順次得られる漸化合同式による一様独立乱数｛ｖｋ：＝ｚｋ／ｄ｜ｋ＝０、１、２、…｝を得るための正の整数列｛ｚｋ｜ｋ＝０、１、２、…｝、
ｚ０＝ｍｏｄ（ｎ、ｄ）、
ｚｋ＋１＝ｍｏｄ（ｚｚｋ、ｄ）、 ｋ＝０、１、２、…、
を得る計算方法であって、孫子の定理の一つの形，部分法ｅとｆだけに対する成分
ｘ：＝ｍｏｄ（ｚ，ｅ）、
ｘｋ：＝ｍｏｄ（ｚｋ，ｅ）、ｋ＝０、１、２、…、
ｙ：＝ｍｏｄ（ｚ，ｆ）、
ｙｋ：＝ｍｏｄ（ｚｋ，ｆ）、ｋ＝０、１、２、…、
ａ：＝ｍｏｄ（ｎ，ｅ）
ｂ：＝ｍｏｄ（ｎ，ｆ）
を用いて定義される２組のＭＣ部分乱数列、ＭＣ生成機構（ｅ、ｘ、ａ）からの法ｅの合同漸化式の解
ｘ０＝ｍｏｄ（ａ、ｅ）、
ｘｋ＋１＝ｍｏｄ（ｘｘｋ、ｅ）、 ｋ＝０、１、２、…、
及びＭＣ生成機構（ｆ、ｙ、ｂ）からの法ｆの合同漸化式の解
ｙ０＝ｍｏｄ（ｂ、ｆ）、
ｙｋ＋１＝ｍｏｄ（ｙｙｋ、ｆ）、 ｋ＝０、１、２、…、
から得られる一次結合｛ｚｋ｜ｋ＝０、１、２、…｝：
ｚｋ：＝ｍｏｄ｛ｆ ｍｏｄ（ｘｋＦ－１、ｅ）＋
ｅ ｍｏｄ（ｙｋＥ－１、ｆ）、ｅｆ｝、
Ｆ－１：は法ｅでのｆの逆数、 ０＜Ｆ－１＜ｅ、
Ｅ－１：は法ｆでのｅの逆数、 ０＜Ｅ－１＜ｆ、
で与えてよい、という形の孫子の定理に基づき、さらに任意の正の整数Ａ、Ｂに対する関係
ｍｏｄ（ｅＡ，ｅｆ）＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）、
ｍｏｄ（ｆＢ，ｅｆ）＝ｍｏｄ（ｆ ｍｏｄ（Ｂ、ｅ）、ｅｆ）、
によってすべての演算を倍精度整数、倍精度実数だけの範囲内で正確高速に得る計算方法。
ｚ０＝ｍｏｄ（ｎ、ｄ）、
ｚｋ＋１＝ｍｏｄ（ｚｚｋ、ｄ）、 ｋ＝０、１、２、…、
を得る計算方法であって、孫子の定理の一つの形，部分法ｅとｆだけに対する成分
ｘ：＝ｍｏｄ（ｚ，ｅ）、
ｘｋ：＝ｍｏｄ（ｚｋ，ｅ）、ｋ＝０、１、２、…、
ｙ：＝ｍｏｄ（ｚ，ｆ）、
ｙｋ：＝ｍｏｄ（ｚｋ，ｆ）、ｋ＝０、１、２、…、
ａ：＝ｍｏｄ（ｎ，ｅ）
ｂ：＝ｍｏｄ（ｎ，ｆ）
を用いて定義される２組のＭＣ部分乱数列、ＭＣ生成機構（ｅ、ｘ、ａ）からの法ｅの合同漸化式の解
ｘ０＝ｍｏｄ（ａ、ｅ）、
ｘｋ＋１＝ｍｏｄ（ｘｘｋ、ｅ）、 ｋ＝０、１、２、…、
及びＭＣ生成機構（ｆ、ｙ、ｂ）からの法ｆの合同漸化式の解
ｙ０＝ｍｏｄ（ｂ、ｆ）、
ｙｋ＋１＝ｍｏｄ（ｙｙｋ、ｆ）、 ｋ＝０、１、２、…、
から得られる一次結合｛ｚｋ｜ｋ＝０、１、２、…｝：
ｚｋ：＝ｍｏｄ｛ｆ ｍｏｄ（ｘｋＦ－１、ｅ）＋
ｅ ｍｏｄ（ｙｋＥ－１、ｆ）、ｅｆ｝、
Ｆ－１：は法ｅでのｆの逆数、 ０＜Ｆ－１＜ｅ、
Ｅ－１：は法ｆでのｅの逆数、 ０＜Ｅ－１＜ｆ、
で与えてよい、という形の孫子の定理に基づき、さらに任意の正の整数Ａ、Ｂに対する関係
ｍｏｄ（ｅＡ，ｅｆ）＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）、
ｍｏｄ（ｆＢ，ｅｆ）＝ｍｏｄ（ｆ ｍｏｄ（Ｂ、ｅ）、ｅｆ）、
によってすべての演算を倍精度整数、倍精度実数だけの範囲内で正確高速に得る計算方法。
【発明の詳細な説明】
【技術分野】
【0001】
整数の３つ組（ｄ、ｚ、ｎ）、ｄ＞０、ｄとは素で乗数と呼ばれるｚ＞０、ｄとは素なｎで初期値或いはｓｅｅｄと呼ばれるもの、が与える乗算合同ｍｕｌｔｉｐｌｉｃａｔｉｖｅ ｃｏｎｇｒｕｅｎｔｉａｌ（ＭＣ）法と呼ばれる乱数生成機構は、開区間（０、ｄ）の整数の列｛ｚ０，ｚ１，ｚ２，…｝をｄを法とする再帰的合同関係
ｚ０＝ｍｏｄ（ｎ，ｄ）， ｚｋ＋１＝ｍｏｄ（ｚｚｋ，ｄ）
ｋ＝１，２，…
によって生成する。即ちＭＣ（ｄ、ｚ、ｎ）乱数列は、区間（０、ｄ）内の整数列｛ａｋ＝ｍｏｄ（ｎｚｋ，ｄ）、ｋ＝０、１、２、…｝を生成する。ここでｍｏｄ（ａ、ｄ）は整数ｄとは素な整数ａを法ｄで割った余り、０とｄの間の剰余、を与える関数である。開区間（０、１）内に分布したＭＣ一様有理乱数列は
｛ｘｋ：＝ｚｋ／ｄ｜ ｋ＝０，１，２，…、０＜ｘｋ＜１｝
で得られ、実際に乱数として用いられる。
ｚ０＝ｍｏｄ（ｎ，ｄ）， ｚｋ＋１＝ｍｏｄ（ｚｚｋ，ｄ）
ｋ＝１，２，…
によって生成する。即ちＭＣ（ｄ、ｚ、ｎ）乱数列は、区間（０、ｄ）内の整数列｛ａｋ＝ｍｏｄ（ｎｚｋ，ｄ）、ｋ＝０、１、２、…｝を生成する。ここでｍｏｄ（ａ、ｄ）は整数ｄとは素な整数ａを法ｄで割った余り、０とｄの間の剰余、を与える関数である。開区間（０、１）内に分布したＭＣ一様有理乱数列は
｛ｘｋ：＝ｚｋ／ｄ｜ ｋ＝０，１，２，…、０＜ｘｋ＜１｝
で得られ、実際に乱数として用いられる。
【0002】
優れた乱数列を得るためには、条件を満たす任意の整数の組（ｄ、ｚ、ｎ）を取ってよいわけではなく、様々な法ｄ、乗数ｚを試し（検定し）それらが優れた一様独立性と見える出力を与えるものを採用する事が唯一の技術的手段である。この最も簡単で、ある意味では原始的な方法は、しかし乱数生成方法の中で最も広汎強力なもので、驚くべきことにどの様な方式に基づく有限乱数列であってもＭＣ乱数で乗数ｚの逆数の精度で近似する。これは２１世紀に漸く発見された事実であり、（文献１）または（文献２）を参照して頂きたい。そしてＭＣ乱数生成方法は最も信頼性の高い検定方法を持つ事が既知である。これも参考文献に委ねて、再度の説明は省く。この様な優れたＭＣ乱数は既に実在する。ここではそれらの実用の新しい技法を報告する。
（文献１）中澤 直也／中澤 宏：“コンピュータで生成する優れた乱数とは何か”、東京図書、出版予定。近似定理は定理１２に記されている。
（文献２）Ｎａｏｙａ Ｎａｋａｚａｗａ／Ｈｉｒｏｓｈｉ Ｎａｋａｚａｗａ：
“Ｒａｎｄｏｍ Ｎｕｍｂｅｒ Ｇｅｎｅｒａｔｏｒｓ ｏｎ Ｃｏｍｐｕｔｅｒｓ”、ＪＥＮＮＹ ＳＴＡＮＦＯＲＤ、２０２２年出版予定。
（文献１）中澤 直也／中澤 宏：“コンピュータで生成する優れた乱数とは何か”、東京図書、出版予定。近似定理は定理１２に記されている。
（文献２）Ｎａｏｙａ Ｎａｋａｚａｗａ／Ｈｉｒｏｓｈｉ Ｎａｋａｚａｗａ：
“Ｒａｎｄｏｍ Ｎｕｍｂｅｒ Ｇｅｎｅｒａｔｏｒｓ ｏｎ Ｃｏｍｐｕｔｅｒｓ”、ＪＥＮＮＹ ＳＴＡＮＦＯＲＤ、２０２２年出版予定。
【0003】
ＭＣ生成機構（ｄ、ｚ、ｎ）を設計するとき、第一の問題は生成乱数列の周期Ｔである。これには種々の場合があるが、基本は法ｄが奇素数でｚがその法で最長の周期Ｔ＝（ｄ－１）／２を与える原始根の場合であり、他のすべての場合でＴ／ｄはそれより小さくなる。シミュレーションは単一のＭＣ乱数生成機構で行なわれなければならず、２つ以上のＭＣ乱数生成機構は、たとえそれら個々の統計性がどんなに優れていても、相次いで或いは並列に、合成された全体乱数列の検定なしに用いる事はできない。その意味で検定を経て優良性が確認されて現存する２つの素数とそれらの原始根達の孫子の定理による組み合わせの合成ＭＣ乱数生成機構は、一様独立性の生成方式として最も優れている。整数ｐが奇素数、乗数ｚがその原始根であれば、その実際に用いる事のできる周期Ｔ＝（ｄ－１）／２は現在のシミュレーションでは２４０以上である事を求められるから、対応して素数ｐと原始根ｚとに基づくＭＣ乱数生成機構はｄが２４１以上である事を求め、奇素数ｐと原子根ｚのＭＣ乱数生成機構では法ｐは倍精度整数でなければならない。この必要は素数の法ｐとその原始根ｚで構成されるＭＣ乱数生成機構（ｐ、ｚ）の発見にも乱数としての演算にも大きな困難をもたらす。なぜならＭＣ乱数の基本は、前段階の整数出力ｚｋ＝ｍｏｄ（ｚｋｚ，ｐ）にさらにｚを掛けて法ｐを計算するから、必ず
ｍｏｄ（ｚ・ｚｋ、ｐ）、
を求めなければならない。この途中でのｚ・ｚｋの計算は倍精度を一時的に超えるので、計算時間を非常に多く要する４倍精度実数として計算し、整数部分を取らなければならない。しかし２４０程度の小さな周期ではごく短時間で使い切ってしまい、この様な実装は効果に対して余りに過大な負担となる。この事は後の段落で具体例で示される。
ｍｏｄ（ｚ・ｚｋ、ｐ）、
を求めなければならない。この途中でのｚ・ｚｋの計算は倍精度を一時的に超えるので、計算時間を非常に多く要する４倍精度実数として計算し、整数部分を取らなければならない。しかし２４０程度の小さな周期ではごく短時間で使い切ってしまい、この様な実装は効果に対して余りに過大な負担となる。この事は後の段落で具体例で示される。
【0004】
前段落、素数の法と原始根の乗数での計算や周期の短さの困難には、ＭＣ乱数が深く考えられるようになって初めて、２つの素数の法を採用すると孫子の定理が魔術的な解決を与えていることが発見された。現在の発明はこの発見に基づくものである。基礎の理解から再考する様、孫子の定理から述べる。法ｄが２つの互いに素な２つの正整数の積ｄ＝ｅｆであるとする。孫子の定理の実現は様々な形で与えられるが、現在の議論に最も便利なのは下の（Ｂ）の形である。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の式（＊＊）で与える事ができる：
ｘ＝ｍｏｄ（ａｆＦ－１、ｅ）＋ｍｏｄ（ｂｅＥ－１、ｆ） （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。
（証明）先に（Ｂ）の形の解（＊＊）がある事を示し、後にそれが法ｄ＝ｅｆでは一意である事を示す。（＊＊）の式を法ｅ、法ｆを取って見ると直ちに：
ｍｏｄ（ｘ、ｅ）＝ａ， ｍｏｄ（ｘ、ｆ）＝ｂ
の成立がわかる。（Ａ）の（＊）の解がもう一つの形ｙ：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
も持つとする。このとき差ｘ－ｙについて
０＝ｍｏｄ（ｘ、ｅ）－ｍｏｄ（ｙ、ｅ）＝ｍｏｄ（ｘ―ｙ、ｅ）
０＝ｍｏｄ（ｘ、ｆ）－ｍｏｄ（ｙ、ｆ）＝ｍｏｄ（ｘ―ｙ、ｆ）
が成り立つ。故にｘ－ｙはｅの倍数、かつｆの倍数である。正のｅとｆは互いに素な（共通素因数のない）整数だから、それはｘ－ｙがｅｆ＝ｄの倍数である事を意味し、法ｄ＝ｅｆでｘとｙとは同一である。（証明終り）
（注１）逆数Ｆ－１とＥ－１について注釈する。法ｅの既約剰余類群の言葉を使わず説明する。整数ｅ＞０とは素な（共通因数を持たない）ｅ未満の整数がオイラーの関数０＜φ（ｅ）＜ｅ個だけある事は既知である。大切なのはそれらの整数の積もすべて法ｅと素で、積に関して法ｅと素な整数の中（法ｅの既約剰余類と呼ばれる）の集合を出ない事である。このφ（ｅ）個の整数を｛ｉ、ｊ、…、ｍ｝と１列に並べる。これらにｅとは素な（ｅとは共通素因数を持たない）整数ｆを法ｅで掛けた列を法ｅで１からｅ－１に引き戻した集合
｛ｍｏｄ（ｉｆ、ｅ）、ｍｏｄ（ｊｆ、ｅ）、…、ｍｏｄ（ｍｆ、ｅ）｝（％）を作ると、これらは再び１以上ｅ未満のφ（ｅ）個の整数で、その中には必ず１が存在する。それをｍｏｄ（ｋｆ，ｅ）＝１とすると、このｋが法ｅでのｆの逆数ｋ＝Ｆ－１である。Ｆ－１を求めるのには（ｅは単精度整数だから）列（％）を順次生成して探しても困難ではない。乱数を使う前に行なえばよい。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の式（＊＊）で与える事ができる：
ｘ＝ｍｏｄ（ａｆＦ－１、ｅ）＋ｍｏｄ（ｂｅＥ－１、ｆ） （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。
（証明）先に（Ｂ）の形の解（＊＊）がある事を示し、後にそれが法ｄ＝ｅｆでは一意である事を示す。（＊＊）の式を法ｅ、法ｆを取って見ると直ちに：
ｍｏｄ（ｘ、ｅ）＝ａ， ｍｏｄ（ｘ、ｆ）＝ｂ
の成立がわかる。（Ａ）の（＊）の解がもう一つの形ｙ：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
も持つとする。このとき差ｘ－ｙについて
０＝ｍｏｄ（ｘ、ｅ）－ｍｏｄ（ｙ、ｅ）＝ｍｏｄ（ｘ―ｙ、ｅ）
０＝ｍｏｄ（ｘ、ｆ）－ｍｏｄ（ｙ、ｆ）＝ｍｏｄ（ｘ―ｙ、ｆ）
が成り立つ。故にｘ－ｙはｅの倍数、かつｆの倍数である。正のｅとｆは互いに素な（共通素因数のない）整数だから、それはｘ－ｙがｅｆ＝ｄの倍数である事を意味し、法ｄ＝ｅｆでｘとｙとは同一である。（証明終り）
（注１）逆数Ｆ－１とＥ－１について注釈する。法ｅの既約剰余類群の言葉を使わず説明する。整数ｅ＞０とは素な（共通因数を持たない）ｅ未満の整数がオイラーの関数０＜φ（ｅ）＜ｅ個だけある事は既知である。大切なのはそれらの整数の積もすべて法ｅと素で、積に関して法ｅと素な整数の中（法ｅの既約剰余類と呼ばれる）の集合を出ない事である。このφ（ｅ）個の整数を｛ｉ、ｊ、…、ｍ｝と１列に並べる。これらにｅとは素な（ｅとは共通素因数を持たない）整数ｆを法ｅで掛けた列を法ｅで１からｅ－１に引き戻した集合
｛ｍｏｄ（ｉｆ、ｅ）、ｍｏｄ（ｊｆ、ｅ）、…、ｍｏｄ（ｍｆ、ｅ）｝（％）を作ると、これらは再び１以上ｅ未満のφ（ｅ）個の整数で、その中には必ず１が存在する。それをｍｏｄ（ｋｆ，ｅ）＝１とすると、このｋが法ｅでのｆの逆数ｋ＝Ｆ－１である。Ｆ－１を求めるのには（ｅは単精度整数だから）列（％）を順次生成して探しても困難ではない。乱数を使う前に行なえばよい。
【0005】
この明快な孫子の定理は、ＭＣ乱数を考えることになって始めて、我々の認識に達した次の魔術的な補題を意味する。
（補題：２つの互いに異なる素数の部分法が与える孫子縮小）法ｄ＝ｅｆは互いに素な２つの部分法ｅ＞０とｆ＞０の積、とする。任意の整数Ａ＞０対して、
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）
＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）
＝ｅ ｍｏｄ（Ａ、ｆ）
であって、ｍｏｄ（ｅＡ、ｆ）は０以上ｄ＝ｅｆ未満である。
（証明）整数Ａ＞０を整数ｆ＞０で割った等式は
Ａ＝ｑｆ＋ｒ、 ｑは商、ｒは余り、 ｒは０以上ｆ未満
である。これから
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）＝
ｍｏｄ（ｅｆｑ＋ｅｒ、ｅｆ）＝ｅｒ＝ｅ ｍｏｄ（Ａ、ｆ）
だから、補題が示された。（証明終り）
法ｄ＝ｅｆが互いに素な２つの部分法ｅ＞０とｆ＞０の積であり、ｅとｆは単精度整数の限界未満、すなわち２３１未満だとする。２ｄ＝２ｅｆは倍精度整数の上限値２６４より小さい。法ｄでのＭＣ乱数列｛ｚ０、ｚ１、…、ｚｋ、ｚｋ＋１、…｝の漸化式の解は、整数単精度内の部分ＭＣ漸化式、途中に倍精度整数の積ａ・ａｋ、ｂ・ｂｋ、を要するが、の解達
ａｋ＋１＝ｍｏｄ（ａ・ａｋ、ｅ）
ｂｋ＋１＝ｍｏｄ（ｂ・ｂｋ、ｆ）
によって
ｚｋ＋１＝ｍｏｄ（ａｋｆＦ－１、ｅ）＋ｍｏｄ（ｂｋｅＥ－１、ｆ） （＊＊）と解かれる。（孫子縮小の補題）より、
ｚｋ＋１＝ｍｏｄ（ｆ ｍｏｄ（ａｋＦ－１、ｅ）、ｅｆ）＋
ｍｏｄ（ｅ ｍｏｄ（ｂｋＥ－１、ｆ）、ｅｆ）。
この右辺の内側の２つの法ｅと法ｆでのｍｏｄ関数はそれぞれ、ｅとｆと未満の正の整数を与え、外側のｍｏｄ関数の内部は２ｅｆ＝２６４を越えない。だから上のすべての計算は倍精度整数内で行なわれ、それを実数の有理数乱数とする計算も倍精度実数内で行なう事ができる。４倍精度計算を必要としないこのＭＣ乱数生成方法は、実用に高速の大きな利便を与える。
（補題：２つの互いに異なる素数の部分法が与える孫子縮小）法ｄ＝ｅｆは互いに素な２つの部分法ｅ＞０とｆ＞０の積、とする。任意の整数Ａ＞０対して、
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）
＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）
＝ｅ ｍｏｄ（Ａ、ｆ）
であって、ｍｏｄ（ｅＡ、ｆ）は０以上ｄ＝ｅｆ未満である。
（証明）整数Ａ＞０を整数ｆ＞０で割った等式は
Ａ＝ｑｆ＋ｒ、 ｑは商、ｒは余り、 ｒは０以上ｆ未満
である。これから
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）＝
ｍｏｄ（ｅｆｑ＋ｅｒ、ｅｆ）＝ｅｒ＝ｅ ｍｏｄ（Ａ、ｆ）
だから、補題が示された。（証明終り）
法ｄ＝ｅｆが互いに素な２つの部分法ｅ＞０とｆ＞０の積であり、ｅとｆは単精度整数の限界未満、すなわち２３１未満だとする。２ｄ＝２ｅｆは倍精度整数の上限値２６４より小さい。法ｄでのＭＣ乱数列｛ｚ０、ｚ１、…、ｚｋ、ｚｋ＋１、…｝の漸化式の解は、整数単精度内の部分ＭＣ漸化式、途中に倍精度整数の積ａ・ａｋ、ｂ・ｂｋ、を要するが、の解達
ａｋ＋１＝ｍｏｄ（ａ・ａｋ、ｅ）
ｂｋ＋１＝ｍｏｄ（ｂ・ｂｋ、ｆ）
によって
ｚｋ＋１＝ｍｏｄ（ａｋｆＦ－１、ｅ）＋ｍｏｄ（ｂｋｅＥ－１、ｆ） （＊＊）と解かれる。（孫子縮小の補題）より、
ｚｋ＋１＝ｍｏｄ（ｆ ｍｏｄ（ａｋＦ－１、ｅ）、ｅｆ）＋
ｍｏｄ（ｅ ｍｏｄ（ｂｋＥ－１、ｆ）、ｅｆ）。
この右辺の内側の２つの法ｅと法ｆでのｍｏｄ関数はそれぞれ、ｅとｆと未満の正の整数を与え、外側のｍｏｄ関数の内部は２ｅｆ＝２６４を越えない。だから上のすべての計算は倍精度整数内で行なわれ、それを実数の有理数乱数とする計算も倍精度実数内で行なう事ができる。４倍精度計算を必要としないこのＭＣ乱数生成方法は、実用に高速の大きな利便を与える。
【0006】
用いられた（孫子の縮小）の内容は殆ど自明だが、与える便宜は魔術的で大きい。ここではこれを用いた異なる２つの素数の積を法とするＭＣ乱数の大きな利点の実証として、良い機会なので、中澤直也が発見し（文献１）（文献２）で公開される優れたＭＣ乱数生成機構＃００１について具体的に述べる。この乱数は２つの素数ｅとの積ｄ＝ｅｆ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５２．７８乗
で生成される。これらの詳しい組成は（文献１）（文献２）に開示されていて、下の図１の乱数計算ｍａｉｎプログラムと図２のｓｕｂｒｏｕｔｉｎｅ、にも開示されているので見て頂きたい。
図１ ＃００１生成機構の乱数発生ｍａｉｎ ｐｒｏｇｒａｍ
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｉｐ１＝１３４２６５２０３
ｉｐ２＝１３４４７５８２７
ｉｄ＝ｉｐ１＊ｉｐ２ ！ｉｄはおよそ２５４です
ａｄ＝ｉｄ
ｉｚ１＝１９０６１２５２
ｉｚ２＝７７６００５２５
ｎ１＝１０ ！初期値１は１０です
ｎ２＝１３ ！初期値２は１３です
ｉｐ２ｍｐ１＝５２５７７００７
ｉｐ１ｍｐ２＝８１８１６２７１
ｄｏ ｉ＝１，１０００００００
ｃａｌｌ ｒａｎｄｏｍ（ｒａｎｄ）
ｅｎｄ ｄｏ
ｅｎｄ
図２ ＃００１生成機構の乱数発生ｓｕｂｒｏｕｔｎｉｎｅ
ｓｕｂｒｏｕｔｉｎｅ ｒａｎｄｏｍ（ｒａｎｄ）
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｍｚ１＝ｍｏｄ（ｍｚ１＊ｉｚ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｍｚ２＊ｉｚ２，ｉｐ２）
ｍｚ１ａ＝ｍｏｄ（ｍｚ１＊ｉｐ２ｍｐ１，ｉｐ１）
ｍｚ２ａ＝ｍｏｄ（ｍｚ２＊ｉｐ１ｍｐ２，ｉｐ２）
ａｚ＝ｍｏｄ（ｉｐ２＊ｍｚ１ａ＋ｉｐ１＊ｍｚ２ａ，ｉｄ）
ｒａｎｄ＝ａｚ／ａｄ
ｒｅｔｕｒｎ
ｅｎｄ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５２．７８乗
で生成される。これらの詳しい組成は（文献１）（文献２）に開示されていて、下の図１の乱数計算ｍａｉｎプログラムと図２のｓｕｂｒｏｕｔｉｎｅ、にも開示されているので見て頂きたい。
図１ ＃００１生成機構の乱数発生ｍａｉｎ ｐｒｏｇｒａｍ
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｉｐ１＝１３４２６５２０３
ｉｐ２＝１３４４７５８２７
ｉｄ＝ｉｐ１＊ｉｐ２ ！ｉｄはおよそ２５４です
ａｄ＝ｉｄ
ｉｚ１＝１９０６１２５２
ｉｚ２＝７７６００５２５
ｎ１＝１０ ！初期値１は１０です
ｎ２＝１３ ！初期値２は１３です
ｉｐ２ｍｐ１＝５２５７７００７
ｉｐ１ｍｐ２＝８１８１６２７１
ｄｏ ｉ＝１，１０００００００
ｃａｌｌ ｒａｎｄｏｍ（ｒａｎｄ）
ｅｎｄ ｄｏ
ｅｎｄ
図２ ＃００１生成機構の乱数発生ｓｕｂｒｏｕｔｎｉｎｅ
ｓｕｂｒｏｕｔｉｎｅ ｒａｎｄｏｍ（ｒａｎｄ）
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｍｚ１＝ｍｏｄ（ｍｚ１＊ｉｚ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｍｚ２＊ｉｚ２，ｉｐ２）
ｍｚ１ａ＝ｍｏｄ（ｍｚ１＊ｉｐ２ｍｐ１，ｉｐ１）
ｍｚ２ａ＝ｍｏｄ（ｍｚ２＊ｉｐ１ｍｐ２，ｉｐ２）
ａｚ＝ｍｏｄ（ｉｐ２＊ｍｚ１ａ＋ｉｐ１＊ｍｚ２ａ，ｉｄ）
ｒａｎｄ＝ａｚ／ａｄ
ｒｅｔｕｒｎ
ｅｎｄ
【0007】
前ページ図１と図２のプログラムが必要とするＣＰＵ時間は０．２８１秒／１０００万個である。出力は１／ｄ、約２－５４の精度（注）だから倍精度実数として十分な桁数を持つ。この計算速度と倍精度とに比べる事は、＃Ｍ００１生成機構の持つ絶対的に不利な性能の位置を与える。下の図３は＃Ｍ００１で同じ１０００万個の乱数を生成するプログラムである。構造は一見大変簡単である。
図３ ＃Ｍ００１で１０００万個を出力するプログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｇ，ｏ－ｐ，ｒ－ｚ），ｒｅａｌ＊１６（ｑ）
ｉｐ＝１７１７９８６９９８９ｄ０ ！ｉｐはおよそ２３４．００
ｒｐ＝ｉｐ
ｉｚ＝７９２８４１００７２ｄ０
ｉｎ＝１０ ！初期値は１０
ｑｍｚ＝ｉｎ，
ｄｏ ｉ＝１，１０００００００
ｒａｎｄ＝ｍｏｄ（ｑｍｚ，ｒｐ）／ｒｐ
ｑｍｚ＝ｍｏｄ（ｑｍｚ＊ｉｚ，ｒｐ）
ｅｎｄ ｄｏ
ｅｎｄ
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数の計算時間のため大幅に増大する。上の実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は（注）１／ｄ程度でほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないものの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。
（注）（文献１）の近似定理（定理１２）によればＭＣ乱数は１／ｚの精度で任意のコンピュータ上の有限長さ一様乱数を近似している。その意味ではＭＣ乱数の数値精度として１／ｚを取るべきだが、ＭＣ乱数として法ｄで整数１からｄ－１までの区間内の整数を割った有理数列としては１／ｄの精度と考える方が（用いる周期Ｔによって区間内の有理数を１／Ｔおきに取るとしても）自然だろう。
図３ ＃Ｍ００１で１０００万個を出力するプログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｇ，ｏ－ｐ，ｒ－ｚ），ｒｅａｌ＊１６（ｑ）
ｉｐ＝１７１７９８６９９８９ｄ０ ！ｉｐはおよそ２３４．００
ｒｐ＝ｉｐ
ｉｚ＝７９２８４１００７２ｄ０
ｉｎ＝１０ ！初期値は１０
ｑｍｚ＝ｉｎ，
ｄｏ ｉ＝１，１０００００００
ｒａｎｄ＝ｍｏｄ（ｑｍｚ，ｒｐ）／ｒｐ
ｑｍｚ＝ｍｏｄ（ｑｍｚ＊ｉｚ，ｒｐ）
ｅｎｄ ｄｏ
ｅｎｄ
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数の計算時間のため大幅に増大する。上の実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は（注）１／ｄ程度でほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないものの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。
（注）（文献１）の近似定理（定理１２）によればＭＣ乱数は１／ｚの精度で任意のコンピュータ上の有限長さ一様乱数を近似している。その意味ではＭＣ乱数の数値精度として１／ｚを取るべきだが、ＭＣ乱数として法ｄで整数１からｄ－１までの区間内の整数を割った有理数列としては１／ｄの精度と考える方が（用いる周期Ｔによって区間内の有理数を１／Ｔおきに取るとしても）自然だろう。
【技術的展望】
【0008】
孫子の定理は一般にｋ個、ｋ＞２、の互いに素な正の数の積の法についても成り立つので、２個の積ｄ＝ｅｆに限るものではない。だから、３個以上の部分法を持つものでも考える事ができる。ただその様な場合に数値的に倍精度計算だけへの設計をどの様に選ぶべきか、撰ぶ事ができるか、は大きな問題である。それは必然的に２成分のｄ＝ｅｆの場合よりは複雑な考察を要するだろう。２成分の場合がきれいな構成を得て倍精度整数演算、倍精度実数演算だけで実現している事実、多分この構成が部分法での成功を最も得やすい形ではないかという予感、からは発明者達はこのより進んだ形は将来への課題ではないか、より大きな乱数周期への需要、より多い桁数精度の必要の生起、を見たとき、自然に得られるのでははないか、と考える。
【図面の簡単な説明】
【0009】
図１はＭＣ乱数生成機構＃００１での１０００万個の乱数発生Ｆｏｒｔｒａｎ主プログラムの一例である。この図はまた＃００１の詳しい情報、部分法や部分乗数の開示、も与えている。図２はそれに対応するＭＣ乱数生成機構＃００１の乱数発生サブプログラム例である。そして図３はＭＣ乱数生成機構＃Ｍ００１で１０００万個の乱数を出力するＦｏｒｔｒａｎプログラムを示す図であり、＃Ｍ００１の構成の詳しい情報も開示している。
【手続補正書】
【提出日】2022-04-04
【手続補正1】
【補正対象書類名】明細書
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【発明の詳細な説明】
【技術分野】
【0001】
整数の３つ組（ｄ、ｚ、ｎ）、ｄ＞０、ｄとは素で乗数と呼ばれるｚ＞０、ｄとは素なｎで初期値或いはｓｅｅｄと呼ばれるもの、が与える乗算合同ｍｕｌｔｉｐｌｉｃａｔｉｖｅ ｃｏｎｇｒｕｅｎｔｉａｌ（ＭＣ）法と呼ばれる乱数生成機構は、開区間（０、ｄ）の整数の列｛ｚ０，ｚ１，ｚ２，…｝をｄを法とする再帰的合同関係
ｚ０＝ｍｏｄ（ｎ，ｄ）， ｚｋ＋１＝ｍｏｄ（ｚｚｋ，ｄ）
ｋ＝１，２，…
によって生成する。即ちＭＣ（ｄ、ｚ、ｎ）乱数列は、区間（０、ｄ）内の整数列｛ａｋ＝ｍｏｄ（ｎｚｋ，ｄ）、ｋ＝０、１、２、…｝を生成する。ここでｍｏｄ（ａ、ｄ）は整数ｄとは素な整数ａを法ｄで割った余り、０とｄの間の剰余、を与える関数である。開区間（０、１）内に分布したＭＣ一様有理乱数列は
｛ｘｋ：＝ｚｋ／ｄ｜ ｋ＝０，１，２，…、０＜ｘｋ＜１｝
で得られ、実際に乱数として用いられる。
ｚ０＝ｍｏｄ（ｎ，ｄ）， ｚｋ＋１＝ｍｏｄ（ｚｚｋ，ｄ）
ｋ＝１，２，…
によって生成する。即ちＭＣ（ｄ、ｚ、ｎ）乱数列は、区間（０、ｄ）内の整数列｛ａｋ＝ｍｏｄ（ｎｚｋ，ｄ）、ｋ＝０、１、２、…｝を生成する。ここでｍｏｄ（ａ、ｄ）は整数ｄとは素な整数ａを法ｄで割った余り、０とｄの間の剰余、を与える関数である。開区間（０、１）内に分布したＭＣ一様有理乱数列は
｛ｘｋ：＝ｚｋ／ｄ｜ ｋ＝０，１，２，…、０＜ｘｋ＜１｝
で得られ、実際に乱数として用いられる。
【0002】
優れた乱数列を得るためには、条件を満たす任意の整数の組（ｄ、ｚ、ｎ）を取ってよいわけではなく、様々な法ｄ、乗数ｚを試し（検定し）それらが優れた一様独立性と見える出力を与えるものを採用する事が唯一の技術的手段である。この最も簡単で、ある意味では原始的な方法は、しかし乱数生成方法の中で最も広汎強力なもので、驚くべきことにどの様な方式に基づく有限乱数列であってもＭＣ乱数で乗数ｚの逆数の精度で近似する。これは２１世紀に漸く発見された事実であり、［参考文献］の（文献１）または（文献２）を参照して頂きたい。そしてＭＣ乱数生成方法はすべての生成方法の中で最も信頼性の高い検定方法を持つ事が既知である。これも参考文献に委ねて、再度の説明は省く。この様な優れたＭＣ乱数は既に実在する。ここではそれらの実用の新しい技法発明を報告する。
参考文献
（文献１）中澤 直也／中澤 宏：コンピュータで生成する優れた乱数とはなに か、東京図書、出版予定。
（文献２）Ｎａｏｙａ Ｎａｋａｚｗａ／Ｈ．Ｎａｋａｚａｗａ：Ｒａｎｄｏｍ Ｎｕｍｂｅｒ Ｇｅｎｅｒａｔｏｒｓ ｏｎ Ｃｏｍｐｕｔｅｒｓ，Ｊｅｎｎｙ Ｓｔａｎｆｏｒｄ，２０２２年出版予定。
参考文献
（文献１）中澤 直也／中澤 宏：コンピュータで生成する優れた乱数とはなに か、東京図書、出版予定。
（文献２）Ｎａｏｙａ Ｎａｋａｚｗａ／Ｈ．Ｎａｋａｚａｗａ：Ｒａｎｄｏｍ Ｎｕｍｂｅｒ Ｇｅｎｅｒａｔｏｒｓ ｏｎ Ｃｏｍｐｕｔｅｒｓ，Ｊｅｎｎｙ Ｓｔａｎｆｏｒｄ，２０２２年出版予定。
【0003】
ＭＣ生成機構（ｄ、ｚ、ｎ）を設計するとき、第一の問題は生成乱数列の周期Ｔである。これには種々の場合があるが、基本は法ｄが奇素数でｚがその法で最長の周期Ｔ＝（ｄ－１）／２を与える原始根の場合であり、他のすべての場合でＴ／ｄはそれより小さくなる。シミュレーションは単一のＭＣ乱数生成機構で行なわれなければならず、２つ以上のＭＣ乱数生成機構は、たとえそれら個々の統計性がどんなに優れていても、相次いで或いは並列に、合成された全体乱数列の検定なしに用いる事はできない。その意味で検定を経て優良性が確認されて現存する２つの素数とそれらの原始根達の孫子の定理による組み合わせの合成ＭＣ乱数生成機構は、一様独立性の生成方式として最も優れている。整数ｄ＝ｐが奇素数、乗数ｚがその原始根であれば、その実際に用いる事のできる周期Ｔ＝（ｄ－１）／２は現在のシミュレーションでは２４０以上である事を求められるから、対応して素数ｐと原始根ｚとに基づくＭＣ乱数生成機構はｄが２４１以上である事を求め、奇素数ｐと原子根ｚのＭＣ乱数生成機構では法ｐは倍精度整数でなければならない。この必要は素数の法ｐとその原始根ｚで構成されるＭＣ乱数生成機構（ｐ、ｚ）の発見にも乱数としての演算にも大きな困難をもたらす。なぜならＭＣ乱数の基本は、前段階の整数出力ｚｋ＝ｍｏｄ（ｎｚ
ｋ
，ｐ）にさらにｚを掛けて法ｐを計算するから、必ず
ｍｏｄ（ｚｚ ｋ、ｐ）、
を求めなければならない。この途中でのｚ・ｚｋの計算は倍精度を一時的に超えるので、計算時間を非常に多く要する４倍精度実数として計算し、整数部分を取らなければならない。しかし２４０程度の小さな周期ではごく短時間で使い切ってしまい、この様な実装は効果に対して余りに過大な負担となる。この事は後の段落で具体例で示される。
ｍｏｄ（ｚｚ ｋ、ｐ）、
を求めなければならない。この途中でのｚ・ｚｋの計算は倍精度を一時的に超えるので、計算時間を非常に多く要する４倍精度実数として計算し、整数部分を取らなければならない。しかし２４０程度の小さな周期ではごく短時間で使い切ってしまい、この様な実装は効果に対して余りに過大な負担となる。この事は後の段落で具体例で示される。
【0004】
素数の法と原始根の乗数での計算や周期の短さの困難には、ＭＣ乱数が深く考えられるようになって初めて、２つの素数の法を採用すると孫子の定理が魔術的な解決を与えていることが発見された。現在の発明はこの発見に基づくものである。基礎の理解から再考する様、孫子の定理から述べる。法ｄが２つの互いに素な２つの正整数の積ｄ＝ｅｆであるとする。孫子の定理の実現は様々な形で与えられるが、現在の議論に最も便利なのは下の（Ｂ）の形である。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の式（＊＊）で与える事ができる：
ｘ＝ａｆＦ －１ ＋ｂｅＥ －１ （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。
（証明）先に（Ｂ）の形の解（＊＊）がある事を示し、後にそれが法ｄ＝ｅｆでは一意である事を示す。（＊＊）の式を法ｅ、法ｆを取って見ると直ちに：
ｍｏｄ（ｘ、ｅ）＝ａ， ｍｏｄ（ｘ、ｆ）＝ｂ
の成立がわかる。（Ａ）の（＊）の解がもう一つの形ｙ：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
も持つとする。このとき差ｘ－ｙについて
０＝ｍｏｄ（ｘ、ｅ）－ｍｏｄ（ｙ、ｅ）＝ｍｏｄ（ｘ―ｙ、ｅ）
０＝ｍｏｄ（ｘ、ｆ）－ｍｏｄ（ｙ、ｆ）＝ｍｏｄ（ｘ―ｙ、ｆ）
が成り立つ。故にｘ－ｙはｅの倍数、かつｆの倍数である。正のｅとｆは互いに素な（共通素因数のない）整数だから、それはｘ－ｙがｅｆ＝ｄの倍数である事を意味し、法ｄ＝ｅｆでｘとｙとは同一である。（証明終）
（注１）逆数Ｆ－１とＥ－１について注釈する。法ｅの既約剰余類群の言葉を使わず説明する。整数ｅ＞０とは素な（共通因数を持たない）ｅ未満の整数がオイラーの関数０＜φ（ｅ）＜ｅ個だけある事は既知である。大切なのはそれらの整数の積もすべて法ｅと素で、積に関して法ｅと素な整数の中（法ｅの既約剰余類と呼ばれる）の集合を出ない事である。このφ（ｅ）個の整数を｛ｉ、ｊ、…、ｍ｝と１列に並べる。これらにｅとは素な（ｅとは共通素因数を持たない）整数ｆを法ｅで掛けた列を法ｅで１からｅ－１に引き戻した集合
｛ｍｏｄ（ｉｆ、ｅ）、ｍｏｄ（ｊｆ、ｅ）、…、ｍｏｄ（ｍｆ、ｅ）｝（％）
を作ると、これらは再び１以上ｅ未満のφ（ｅ）個の整数で、その中には必ず１が存在する。それをｍｏｄ（ｋｆ，ｅ）＝１とすると、このｋが法ｅでのｆの逆数ｋ＝Ｆ－１である。Ｆ－１を求めるのには（ｅは単精度整数だから）列（％）を順次生成して探しても困難ではない。乱数生成機構を作成するに先立って行なえばよい。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の式（＊＊）で与える事ができる：
ｘ＝ａｆＦ －１ ＋ｂｅＥ －１ （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。
（証明）先に（Ｂ）の形の解（＊＊）がある事を示し、後にそれが法ｄ＝ｅｆでは一意である事を示す。（＊＊）の式を法ｅ、法ｆを取って見ると直ちに：
ｍｏｄ（ｘ、ｅ）＝ａ， ｍｏｄ（ｘ、ｆ）＝ｂ
の成立がわかる。（Ａ）の（＊）の解がもう一つの形ｙ：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
も持つとする。このとき差ｘ－ｙについて
０＝ｍｏｄ（ｘ、ｅ）－ｍｏｄ（ｙ、ｅ）＝ｍｏｄ（ｘ―ｙ、ｅ）
０＝ｍｏｄ（ｘ、ｆ）－ｍｏｄ（ｙ、ｆ）＝ｍｏｄ（ｘ―ｙ、ｆ）
が成り立つ。故にｘ－ｙはｅの倍数、かつｆの倍数である。正のｅとｆは互いに素な（共通素因数のない）整数だから、それはｘ－ｙがｅｆ＝ｄの倍数である事を意味し、法ｄ＝ｅｆでｘとｙとは同一である。（証明終）
（注１）逆数Ｆ－１とＥ－１について注釈する。法ｅの既約剰余類群の言葉を使わず説明する。整数ｅ＞０とは素な（共通因数を持たない）ｅ未満の整数がオイラーの関数０＜φ（ｅ）＜ｅ個だけある事は既知である。大切なのはそれらの整数の積もすべて法ｅと素で、積に関して法ｅと素な整数の中（法ｅの既約剰余類と呼ばれる）の集合を出ない事である。このφ（ｅ）個の整数を｛ｉ、ｊ、…、ｍ｝と１列に並べる。これらにｅとは素な（ｅとは共通素因数を持たない）整数ｆを法ｅで掛けた列を法ｅで１からｅ－１に引き戻した集合
｛ｍｏｄ（ｉｆ、ｅ）、ｍｏｄ（ｊｆ、ｅ）、…、ｍｏｄ（ｍｆ、ｅ）｝（％）
を作ると、これらは再び１以上ｅ未満のφ（ｅ）個の整数で、その中には必ず１が存在する。それをｍｏｄ（ｋｆ，ｅ）＝１とすると、このｋが法ｅでのｆの逆数ｋ＝Ｆ－１である。Ｆ－１を求めるのには（ｅは単精度整数だから）列（％）を順次生成して探しても困難ではない。乱数生成機構を作成するに先立って行なえばよい。
【0005】
明快な孫子の定理は、ＭＣ乱数を考えることになって始めて、我々の認識に達した次の魔術的な補題を意味する。
（補題：２つの互いに異なる素なの部分法が与える孫子縮小）法ｄ＝ｅｆは互いに素な２つの部分法ｅ＞０とｆ＞０の積、とする。任意の整数Ａ＞０対して、
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）
＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）
＝ｅ ｍｏｄ（Ａ、ｆ）
であって、ｍｏｄ（ｅＡ、ｄ）は０以上ｄ＝ｅｆ未満である。
（証明）整数Ａ＞０を整数ｆ＞０で割った等式は
Ａ＝ｑｆ＋ｒ、 ｑは商、ｒは余り、 ｒは０以上ｆ未満
である。これから
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）＝
ｍｏｄ（ｅｆｑ＋ｅｒ、ｅｆ）＝ｅｒ＝ｅ ｍｏｄ（Ａ、ｆ）
だから、補題が示された。（証明終り）
法ｄ＝ｅｆが互いに素な２つの部分法ｅ＞０とｆ＞０の積であり、ｅとｆは単精度整数の限界未満、すなわち２３１未満だとする。２ｄ＝２ｅｆは倍精度整数の上限値２６４より小さい。法ｄでのＭＣ乱数列｛ｚ０、ｚ１、…、ｚｋ、ｚｋ＋１、…｝の漸化式の解は、整数単精度内の部分ＭＣ漸化式、途中に倍精度整数の積、ａとａ ｋ の積、ｂとｂ ｋ の積、を要するが、その解達
ａｋ＋１＝ｍｏｄ（ａａ ｋ、ｅ）
ｂｋ＋１＝ｍｏｄ（ｂｂ ｋ、ｆ）
によって、法ｄで一意に
ｚｋ＋１＝ａ ｋ＋１ ｆＦ －１ ＋ｂ ｋ＋１ ｅＥ －１ （＊＊）
と解かれる。（孫子縮小の補題）より、
ｚｋ＋１＝ｍｏｄ（ｆ ｍｏｄ（ａ ｋ＋１Ｆ－１、ｅ）、ｅｆ）＋
ｍｏｄ（ｅ ｍｏｄ（ｂ ｋ＋１Ｅ－１、ｆ）、ｅｆ）。
この右辺の内側の２つの法ｅと法ｆでのｍｏｄ関数はそれぞれ、ｅとｆと未満の正の整数を与え、外側のｍｏｄ関数の内部は２ｅｆ＝２６４を越えない。だから上のすべての計算は倍精度整数内で行なわれ、それを実数の有理数乱数とする計算も倍精度実数内で行なう事ができる。４倍精度計算を必要としないこのＭＣ乱数生成方法は、実用計算に高速の大きな利便を与える。
（補題：２つの互いに異なる素なの部分法が与える孫子縮小）法ｄ＝ｅｆは互いに素な２つの部分法ｅ＞０とｆ＞０の積、とする。任意の整数Ａ＞０対して、
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）
＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）
＝ｅ ｍｏｄ（Ａ、ｆ）
であって、ｍｏｄ（ｅＡ、ｄ）は０以上ｄ＝ｅｆ未満である。
（証明）整数Ａ＞０を整数ｆ＞０で割った等式は
Ａ＝ｑｆ＋ｒ、 ｑは商、ｒは余り、 ｒは０以上ｆ未満
である。これから
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）＝
ｍｏｄ（ｅｆｑ＋ｅｒ、ｅｆ）＝ｅｒ＝ｅ ｍｏｄ（Ａ、ｆ）
だから、補題が示された。（証明終り）
法ｄ＝ｅｆが互いに素な２つの部分法ｅ＞０とｆ＞０の積であり、ｅとｆは単精度整数の限界未満、すなわち２３１未満だとする。２ｄ＝２ｅｆは倍精度整数の上限値２６４より小さい。法ｄでのＭＣ乱数列｛ｚ０、ｚ１、…、ｚｋ、ｚｋ＋１、…｝の漸化式の解は、整数単精度内の部分ＭＣ漸化式、途中に倍精度整数の積、ａとａ ｋ の積、ｂとｂ ｋ の積、を要するが、その解達
ａｋ＋１＝ｍｏｄ（ａａ ｋ、ｅ）
ｂｋ＋１＝ｍｏｄ（ｂｂ ｋ、ｆ）
によって、法ｄで一意に
ｚｋ＋１＝ａ ｋ＋１ ｆＦ －１ ＋ｂ ｋ＋１ ｅＥ －１ （＊＊）
と解かれる。（孫子縮小の補題）より、
ｚｋ＋１＝ｍｏｄ（ｆ ｍｏｄ（ａ ｋ＋１Ｆ－１、ｅ）、ｅｆ）＋
ｍｏｄ（ｅ ｍｏｄ（ｂ ｋ＋１Ｅ－１、ｆ）、ｅｆ）。
この右辺の内側の２つの法ｅと法ｆでのｍｏｄ関数はそれぞれ、ｅとｆと未満の正の整数を与え、外側のｍｏｄ関数の内部は２ｅｆ＝２６４を越えない。だから上のすべての計算は倍精度整数内で行なわれ、それを実数の有理数乱数とする計算も倍精度実数内で行なう事ができる。４倍精度計算を必要としないこのＭＣ乱数生成方法は、実用計算に高速の大きな利便を与える。
【0006】
（孫子縮小）の内容は殆ど自明だが、与える便宜は魔術的で大きい。
ここではこれを用いた異なる２つの素数の積を法とするＭＣ乱数の大きな利点の実証として、中澤直也が発見し（文献１）と（文献２）で公開される優れたＭＣ乱数生成機構＃００１について具体的に述べる。この乱数は２つの素数ｅとｆとの積ｄ＝ｅｆ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５３乗
で生成される。これらの詳しい組成は（文献１）と（文献２）に開示されている。下の乱数計算主プログラムと副ルーチンにも開示されるので見て頂きたい。
！＃００１生成機構の乱数１０００万個発生主プログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｉｐ１＝１３４２６５０２３
ｉｐ２＝１３４４７５８２７
ｉｄ＝ｉｐ１＊ｉｐ２ ！ｉｄはおよそ２５４です
ａｄ＝ｉｄ
ｉｚ１＝１９０６１２５２
ｉｚ２＝７７６００５２５
ｎ１＝１０ ！初期値１は１０です
ｎ２＝１３ ！初期値２は１３です
ｉｐ２ｍｐ１＝５２５７７００７
ｉｐ１ｍｐ２＝８１８１６２７１
ｄｏ ｉ＝１，１０００００００
ｃａｌｌ ｒａｎｄｏｍ（ｒａｎｄ）
ｅｎｄ ｄｏ
ｅｎｄ
！＃００１生成機構の乱数発生副プログラム
ｓｕｂｒｏｕｔｉｎｅ ｒａｎｄｏｍ（ｒａｎｄ）
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｍｚ１＝ｍｏｄ（ｍｚ１＊ｉｚ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｍｚ２＊ｉｚ２，ｉｐ２）
ｍｚ１ａ＝ｍｏｄ（ｍｚ１＊ｉｐ２ｍｐ１，ｉｐ１）
ｍｚ２ａ＝ｍｏｄ（ｍｚ２＊ｉｐ１ｍｐ２，ｉｐ２）
ａｚ＝ｍｏｄ（ｉｐ２＊ｍｚ１ａ＋ｉｐ１＊ｍｚ２ａ，ｉｄ）
ｒａｎｄ＝ａｚ／ａｄ
ｒｅｔｕｒｎ
ｅｎｄ
ここではこれを用いた異なる２つの素数の積を法とするＭＣ乱数の大きな利点の実証として、中澤直也が発見し（文献１）と（文献２）で公開される優れたＭＣ乱数生成機構＃００１について具体的に述べる。この乱数は２つの素数ｅとｆとの積ｄ＝ｅｆ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５３乗
で生成される。これらの詳しい組成は（文献１）と（文献２）に開示されている。下の乱数計算主プログラムと副ルーチンにも開示されるので見て頂きたい。
！＃００１生成機構の乱数１０００万個発生主プログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｉｐ１＝１３４２６５０２３
ｉｐ２＝１３４４７５８２７
ｉｄ＝ｉｐ１＊ｉｐ２ ！ｉｄはおよそ２５４です
ａｄ＝ｉｄ
ｉｚ１＝１９０６１２５２
ｉｚ２＝７７６００５２５
ｎ１＝１０ ！初期値１は１０です
ｎ２＝１３ ！初期値２は１３です
ｉｐ２ｍｐ１＝５２５７７００７
ｉｐ１ｍｐ２＝８１８１６２７１
ｄｏ ｉ＝１，１０００００００
ｃａｌｌ ｒａｎｄｏｍ（ｒａｎｄ）
ｅｎｄ ｄｏ
ｅｎｄ
！＃００１生成機構の乱数発生副プログラム
ｓｕｂｒｏｕｔｉｎｅ ｒａｎｄｏｍ（ｒａｎｄ）
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ．ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｍｚ１＝ｍｏｄ（ｍｚ１＊ｉｚ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｍｚ２＊ｉｚ２，ｉｐ２）
ｍｚ１ａ＝ｍｏｄ（ｍｚ１＊ｉｐ２ｍｐ１，ｉｐ１）
ｍｚ２ａ＝ｍｏｄ（ｍｚ２＊ｉｐ１ｍｐ２，ｉｐ２）
ａｚ＝ｍｏｄ（ｉｐ２＊ｍｚ１ａ＋ｉｐ１＊ｍｚ２ａ，ｉｄ）
ｒａｎｄ＝ａｚ／ａｄ
ｒｅｔｕｒｎ
ｅｎｄ
【0007】
＃００１ＭＣ乱数の前出実行プログラムが必要とするＣＰＵ時間は０．２８１秒／１０００万個である。出力は１／ｄ、約２－５４の精度（注）だから倍精度として十分な桁数を持つ。この計算速度と倍精度とに比べる事は、＃Ｍ００１生成機構の持つ絶対的に不利な性能の位置を明らかにする。下は＃Ｍ００１で同じ１０００万個の乱数を生成するプログラムである。構造は一見大変簡単である。
！＃Ｍ００１で１０００万個を出力するプログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｇ，ｏ－ｐ，ｒ－ｚ），ｒｅａｌ＊１６（ｑ）
ｉｐ＝１７１７９８６９９８９ｄ０ ！ｉｐはおよそ２３４．００
ｒｐ＝ｉｐ
ｉｚ＝７９２８４１００７２ｄ０
ｉｎ＝１０ ！初期値は１０
ｑｍｚ＝ｉｎ，
ｄｏ ｉ＝１，１０００００００
ｒａｎｄ＝ｍｏｄ（ｑｍｚ，ｒｐ）／ｒｐ
ｑｍｚ＝ｍｏｄ（ｑｍｚ＊ｉｚ，ｒｐ）
ｅｎｄ ｄｏ
ｅｎｄ
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数のため大幅にその計算時間を増大する。上の実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は（注）１／ｄ程度でほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないものの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。
のコンピュータ上の有限長さ一様乱数を近似している。その意味ではＭＣ乱数の数値精度として１／ｚを取るべきだが、ＭＣ乱数として法ｄで整数１からｄ－１までの区間内の整数を割った有理数列としては１／ｄの精度と考える方が（用いる周期Ｔによって区間内の有理数を１／Ｔおきに取るとしても）自然だろう。
！＃Ｍ００１で１０００万個を出力するプログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｇ，ｏ－ｐ，ｒ－ｚ），ｒｅａｌ＊１６（ｑ）
ｉｐ＝１７１７９８６９９８９ｄ０ ！ｉｐはおよそ２３４．００
ｒｐ＝ｉｐ
ｉｚ＝７９２８４１００７２ｄ０
ｉｎ＝１０ ！初期値は１０
ｑｍｚ＝ｉｎ，
ｄｏ ｉ＝１，１０００００００
ｒａｎｄ＝ｍｏｄ（ｑｍｚ，ｒｐ）／ｒｐ
ｑｍｚ＝ｍｏｄ（ｑｍｚ＊ｉｚ，ｒｐ）
ｅｎｄ ｄｏ
ｅｎｄ
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数のため大幅にその計算時間を増大する。上の実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は（注）１／ｄ程度でほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないものの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。
【手続補正書】
【提出日】2022-05-31
【手続補正1】
【補正対象書類名】明細書
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【発明の詳細な説明】
【技術分野】
【0001】
整数の３つ組（ｄ、ｚ、ｎ）、ｄ＞０、ｄとは素で乗数と呼ばれるｚ＞０、ｄとは素なｎで初期値或いはｓｅｅｄと呼ばれるもの、が与える乗算合同ｍｕｌｔｉｐｌｉｃａｔｉｖｅ ｃｏｎｇｒｕｅｎｔｉａｌ（ＭＣ）法と呼ばれる乱数生成機構は、開区間（０、ｄ）の整数の列｛ｚ０，ｚ１，ｚ２，…｝をｄを法とする再帰的合同関係
ｚ０＝ｍｏｄ（ｎ，ｄ）， ｚｋ＋１＝ｍｏｄ（ｚｚｋ，ｄ）
ｋ＝１，２，…
によって生成する。即ちＭＣ（ｄ、ｚ、ｎ）乱数列は、区間（０、ｄ）内の整数列｛ａｋ＝ｍｏｄ（ｎｚｋ，ｄ）、ｋ＝０、１、２、…｝を生成する。ここでｍｏｄ（ａ、ｄ）は整数ｄとは素な整数ａを法ｄで割った余り、０とｄの間の剰余、を与える関数である。開区間（０、１）内に分布したＭＣ一様有理乱数列は
｛ｘｋ：＝ｚｋ／ｄ｜ ｋ＝０，１，２，…、０＜ｘｋ＜１｝
で得られ、実際に乱数として用いられる。
ｚ０＝ｍｏｄ（ｎ，ｄ）， ｚｋ＋１＝ｍｏｄ（ｚｚｋ，ｄ）
ｋ＝１，２，…
によって生成する。即ちＭＣ（ｄ、ｚ、ｎ）乱数列は、区間（０、ｄ）内の整数列｛ａｋ＝ｍｏｄ（ｎｚｋ，ｄ）、ｋ＝０、１、２、…｝を生成する。ここでｍｏｄ（ａ、ｄ）は整数ｄとは素な整数ａを法ｄで割った余り、０とｄの間の剰余、を与える関数である。開区間（０、１）内に分布したＭＣ一様有理乱数列は
｛ｘｋ：＝ｚｋ／ｄ｜ ｋ＝０，１，２，…、０＜ｘｋ＜１｝
で得られ、実際に乱数として用いられる。
【0002】
優れた乱数列を得るためには、条件を満たす任意の整数の組（ｄ、ｚ、ｎ）を取ってよいわけではなく、様々な法ｄ、乗数ｚを試し（検定し）それらが優れた一様独立性と見える出力を与えるものを採用する事が唯一の技術的手段である。この最も簡単で、ある意味では原始的な方法は、しかし乱数生成方法の中で最も広汎強力なもので、驚くべきことにどの様な方式に基づく有限乱数列であってもＭＣ乱数で乗数ｚの逆数の精度で近似する。これは２１世紀に漸く発見された事実であり、（文献１）または（文献２）を参照して頂きたい。そしてＭＣ乱数生成方法は最も信頼性の高い検定方法を持つ事が既知である。これも参考文献に委ねて、再度の説明は省く。この様な優れたＭＣ乱数は既に実在する。ここではそれらの実用の新しい技術を報告する。
（文献１）中澤 直也／中澤 宏：“コンピュ－タで生成する優れた乱数とは何か”、東京図書、出版予定。
（文献２）Ｎａｏｙａ Ｎａｋａｚａｗａ／Ｈｉｒｏｓｈｉ Ｎａｋａｚａｗａ：“Ｒａｎｄｏｍ Ｎｕｍｂｅｒ Ｇｅｎｅｒａｔｏｒｓ ｏｎ Ｃｏｍｐｕｔｅｒｓ”、ＪＥＮＮＹ ＳＴＡＮＦＯＲＤ、２０２２年出版予定。
（文献１）中澤 直也／中澤 宏：“コンピュ－タで生成する優れた乱数とは何か”、東京図書、出版予定。
（文献２）Ｎａｏｙａ Ｎａｋａｚａｗａ／Ｈｉｒｏｓｈｉ Ｎａｋａｚａｗａ：“Ｒａｎｄｏｍ Ｎｕｍｂｅｒ Ｇｅｎｅｒａｔｏｒｓ ｏｎ Ｃｏｍｐｕｔｅｒｓ”、ＪＥＮＮＹ ＳＴＡＮＦＯＲＤ、２０２２年出版予定。
【0003】
ＭＣ生成機構（ｄ、ｚ、ｎ）を設計するとき、第一の問題は生成乱数列の周期Ｔである。これには種々の場合があるが、基本は法ｄが奇素数でｚがその法で最長の周期Ｔ＝（ｄ－１）／２を与える原始根の場合であり、他のすべての場合でＴ／ｄはそれより小さくなる。シミュレーションは単一のＭＣ乱数生成機構で行なわれなければならず、２つ以上のＭＣ乱数生成機構は、たとえそれら個々の統計性がどんなに優れていても、相次いで或いは並列に、合成された全体乱数列の検定なしに用いる事はできない。その意味で検定を経て優良性が確認されて存在する２つの素数とそれらの原始根達の孫子の定理による組み合わせの合成ＭＣ乱数生成機構は、一様独立性の生成方式として最も有力である。整数ｐが奇素数、乗数ｚがその原始根であれば、その実際に用いる事のできる周期Ｔ＝（ｄ－１）／２は現在のシミュレーションでは２４０以上である事を求められるから、対応して素数ｐと原始根ｚとに基づくＭＣ乱数生成機構はｄが２４１以上である事を求める。即ち奇素数ｐと原子根ｚのＭＣ乱数生成機構では法ｐは倍精度整数でなければならない。この事は素数の法ｐとその原始根ｚで構成されるＭＣ乱数生成機構（ｐ、ｚ）の算出演算に大きな困難をもたらす。ＭＣ乱数の基本は、前段階の整数出力ｚｋ＝ｍｏｄ（ｎｚｋ，ｐ）にさらにｚを掛けて法ｐを計算する、のであるから、必ず
ｍｏｄ（ｚ・ｚｋ、ｐ）、
を求めなければならない（積を表す・は今後省略する）。この途中でのｚｚｋの計算は倍精度を一時的に超えるので、用いる事が困難で計算時間を非常に多く要する４倍精度実数として計算し、整数部分を取らなければならない。しかし２４０程度の小さな周期ではごく短時間で使い切ってしまい、実装は効果に対して余りにも過大な負担となる。この事は後の段落で具体例で示される。
ｍｏｄ（ｚ・ｚｋ、ｐ）、
を求めなければならない（積を表す・は今後省略する）。この途中でのｚｚｋの計算は倍精度を一時的に超えるので、用いる事が困難で計算時間を非常に多く要する４倍精度実数として計算し、整数部分を取らなければならない。しかし２４０程度の小さな周期ではごく短時間で使い切ってしまい、実装は効果に対して余りにも過大な負担となる。この事は後の段落で具体例で示される。
【0004】
前段落、素数の法と原始根の乗数の計算や周期の短さの困難は、ＭＣ乱数が深く考えられるようになって初めて、孫子の定理が魔術的な解決を与えていることが発見された。現在の発明はこの発見に基づくものである。理解把握から再考する様、孫子の定理の再認識から述べる。法ｄが２つの互いに素な２つの正整数の積ｄ＝ｅｆであるとする。孫子の定理の実現は様々な形を与えられるが、現在の議論に最も便利なのは下の（Ｂ）の形である。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の等式（＊＊）で与えられる：
ｘ＝ｍｏｄ（ａｆＦ－１＋ｂｅＥ－１、ｅｆ） （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。
（証明）先に（Ｂ）の形の解（＊＊）は（＊）を満たす事を示す。最も簡単で紛れのない初等的証明として（＊＊）を等式で表す。（＊＊）の右辺はｄ＝ｅｆで割った余りなのだから、ある整数ｋを取って
ｘ＝ａｆＦ－１＋ｂｅＥ－１－ｋｅｆ
と表わす事ができる。（注２）これを法ｅで見ればｆＦ－１＝１となって
ｍｏｄ（ｘ、ｅ）＝ａ
が得られる。同様に法ｆで見ればｅＥ－１＝１となって、
ｍｏｄ（ｘ、ｅ）＝ｂ
であるとわかる。
（Ａ）２つの合同式（＊）がもう一つの解ｙも持つとする：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
このとき差ｘ－ｙについて
０＝ｍｏｄ（ｘ、ｅ）－ｍｏｄ（ｙ、ｅ）＝ｍｏｄ（ｘ―ｙ、ｅ）
０＝ｍｏｄ（ｘ、ｆ）－ｍｏｄ（ｙ、ｆ）＝ｍｏｄ（ｘ―ｙ、ｆ）
が成り立つ。故にｘ－ｙはｅの倍数、かつｆの倍数である。正のｅとｆは互いに素な（共通素因数のない）整数だから、それはｘ－ｙがｅｆ＝ｄの倍数である事を意味し、法ｄ＝ｅｆでｘとｙとは同一である。（証明終り）
（注１）逆数Ｆ－１とＥ－１について詳しく述べる。整数ｅ＞０とは素な（共通因数を持たない）ｅ未満の整数がオイラーの関数０＜φ（ｅ）＜ｅ個だけある事は既知である。大切なのはそれらの整数の積もすべて法ｅと素で積に関して法ｅの（既約剰余類）群を作る事である。このφ（ｅ）個の整数を｛ｉ、ｊ、…、ｍ｝と１列に並べる。これらにｅとは素な（ｅとは共通素因数を持たない）整数ｆを法ｅで掛けた列
｛ｍｏｄ（ｆｉ、ｅ）、ｍｏｄ（ｆｊ、ｅ）、…、ｍｏｄ（ｆｍ、ｅ）｝（％）
は再び１以上ｅ未満の同じ物を含まないφ（ｅ）個の整数で、その中には必ず１が存在する。それをｍｏｄ（ｆｋ，ｅ）＝１とすると、このｋが法ｅでのｆの逆数ｋ＝Ｆ－１である。Ｆ－１を求めるには、ｅは単精度整数だから、（％）を順次生成して探しても困難ではない。乱数を使う前に十分時間を掛けて行なえばよい。
（注２）整数ｋは法ｄ＝ｅｆで割った商に当たるが、この詳細はどうでもよい。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の等式（＊＊）で与えられる：
ｘ＝ｍｏｄ（ａｆＦ－１＋ｂｅＥ－１、ｅｆ） （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。
（証明）先に（Ｂ）の形の解（＊＊）は（＊）を満たす事を示す。最も簡単で紛れのない初等的証明として（＊＊）を等式で表す。（＊＊）の右辺はｄ＝ｅｆで割った余りなのだから、ある整数ｋを取って
ｘ＝ａｆＦ－１＋ｂｅＥ－１－ｋｅｆ
と表わす事ができる。（注２）これを法ｅで見ればｆＦ－１＝１となって
ｍｏｄ（ｘ、ｅ）＝ａ
が得られる。同様に法ｆで見ればｅＥ－１＝１となって、
ｍｏｄ（ｘ、ｅ）＝ｂ
であるとわかる。
（Ａ）２つの合同式（＊）がもう一つの解ｙも持つとする：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
このとき差ｘ－ｙについて
０＝ｍｏｄ（ｘ、ｅ）－ｍｏｄ（ｙ、ｅ）＝ｍｏｄ（ｘ―ｙ、ｅ）
０＝ｍｏｄ（ｘ、ｆ）－ｍｏｄ（ｙ、ｆ）＝ｍｏｄ（ｘ―ｙ、ｆ）
が成り立つ。故にｘ－ｙはｅの倍数、かつｆの倍数である。正のｅとｆは互いに素な（共通素因数のない）整数だから、それはｘ－ｙがｅｆ＝ｄの倍数である事を意味し、法ｄ＝ｅｆでｘとｙとは同一である。（証明終り）
（注１）逆数Ｆ－１とＥ－１について詳しく述べる。整数ｅ＞０とは素な（共通因数を持たない）ｅ未満の整数がオイラーの関数０＜φ（ｅ）＜ｅ個だけある事は既知である。大切なのはそれらの整数の積もすべて法ｅと素で積に関して法ｅの（既約剰余類）群を作る事である。このφ（ｅ）個の整数を｛ｉ、ｊ、…、ｍ｝と１列に並べる。これらにｅとは素な（ｅとは共通素因数を持たない）整数ｆを法ｅで掛けた列
｛ｍｏｄ（ｆｉ、ｅ）、ｍｏｄ（ｆｊ、ｅ）、…、ｍｏｄ（ｆｍ、ｅ）｝（％）
は再び１以上ｅ未満の同じ物を含まないφ（ｅ）個の整数で、その中には必ず１が存在する。それをｍｏｄ（ｆｋ，ｅ）＝１とすると、このｋが法ｅでのｆの逆数ｋ＝Ｆ－１である。Ｆ－１を求めるには、ｅは単精度整数だから、（％）を順次生成して探しても困難ではない。乱数を使う前に十分時間を掛けて行なえばよい。
（注２）整数ｋは法ｄ＝ｅｆで割った商に当たるが、この詳細はどうでもよい。
【0005】
この明快な孫子の定理は、ＭＣ乱数を考えることになって始めて、我々の認識に達した次の魔術的な補題を意味する。
（補題：２つの互いに素な部分法が与える孫子縮小）法ｄ＝ｅｆは互いに素な２つの部分法ｅ＞０とｆ＞０の積、とする。任意の整数Ａ＞０対して、
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）
＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）
＝ｅ ｍｏｄ（Ａ、ｆ）
０＜ｍｏｄ（ｅＡ、ｆ）＜ｅｆ＝ｄ
および０＜ｍｏｄ（ｆＢ，ｄ）＜ｅｆが成り立つ。
（証明）整数Ａ＞０を整数ｆ＞０で割った等式は
Ａ＝Ｑｆ＋Ｒ、 Ｑは商、Ｒは余り、 ０＜Ｒ＜ｆ
である。これから直ちに補題が従う。同様の証明は明らかに
０＜ｍｏｄ（ｆＢ，ｄ）＜ｅｆ
についても成り立つ。（証明終り）
法ｄ＝ｅｆが互いに素な２つの部分法ｅ＞０とｆ＞０の積であり、ｅとｆは単精度整数の限界未満、すなわち２３２未満だとし、２ｄ＝２ｅｆは倍精度整数の上限値２６４より小さいと仮定する。法ｄでのＭＣ乱数列｛ｚ０、ｚ１、…、ｚｋ、ｚｋ＋１、…｝の漸化式の解は、整数単精度内の部分ＭＣ漸化式、途中に倍精度整数の積ａａｋ、ｂｂｋ、を要するが、の解達
ａｋ＋１＝ｍｏｄ（ａａｋ、ｅ）
ｂｋ＋１＝ｍｏｄ（ｂｂｋ、ｆ）
によって
ｚｋ＋１＝ｍｏｄ（ｍｏｄ（ａｋ＋１ｆＦ－１、ｅ）＋
ｍｏｄ（ｂｋ＋１ｅＥ－１、ｆ））、ｅｆ） （＊＊）
と解かれる。（孫子縮小の補題）より、
ｚｋ＋１＝ｍｏｄ（ｆ ｍｏｄ（ａｋ＋１Ｆ－１、ｅ）＋
ｅ ｍｏｄ（ｂｋ＋１Ｅ－１、ｆ）、ｅｆ）
である。この右辺の内側の２つのｍｏｄ関数、法ｅと法ｆでのｍｏｄ関数、はそれぞれｆとｅと未満の正の整数を与え、外側のｍｏｄ関数の内部は２ｅｆ＝２６４を越えない。だから上のすべての計算は倍精度整数内で行なわれ、それを実数の有理数乱数とする計算も倍精度実数内で行なう事ができる。４倍精度計算を必要としないこのＭＣ乱数生成方法は、実用に高速の大きな利便を与える。
（補題：２つの互いに素な部分法が与える孫子縮小）法ｄ＝ｅｆは互いに素な２つの部分法ｅ＞０とｆ＞０の積、とする。任意の整数Ａ＞０対して、
ｍｏｄ（ｅＡ，ｄ）＝ｍｏｄ（ｅＡ、ｅｆ）
＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）
＝ｅ ｍｏｄ（Ａ、ｆ）
０＜ｍｏｄ（ｅＡ、ｆ）＜ｅｆ＝ｄ
および０＜ｍｏｄ（ｆＢ，ｄ）＜ｅｆが成り立つ。
（証明）整数Ａ＞０を整数ｆ＞０で割った等式は
Ａ＝Ｑｆ＋Ｒ、 Ｑは商、Ｒは余り、 ０＜Ｒ＜ｆ
である。これから直ちに補題が従う。同様の証明は明らかに
０＜ｍｏｄ（ｆＢ，ｄ）＜ｅｆ
についても成り立つ。（証明終り）
法ｄ＝ｅｆが互いに素な２つの部分法ｅ＞０とｆ＞０の積であり、ｅとｆは単精度整数の限界未満、すなわち２３２未満だとし、２ｄ＝２ｅｆは倍精度整数の上限値２６４より小さいと仮定する。法ｄでのＭＣ乱数列｛ｚ０、ｚ１、…、ｚｋ、ｚｋ＋１、…｝の漸化式の解は、整数単精度内の部分ＭＣ漸化式、途中に倍精度整数の積ａａｋ、ｂｂｋ、を要するが、の解達
ａｋ＋１＝ｍｏｄ（ａａｋ、ｅ）
ｂｋ＋１＝ｍｏｄ（ｂｂｋ、ｆ）
によって
ｚｋ＋１＝ｍｏｄ（ｍｏｄ（ａｋ＋１ｆＦ－１、ｅ）＋
ｍｏｄ（ｂｋ＋１ｅＥ－１、ｆ））、ｅｆ） （＊＊）
と解かれる。（孫子縮小の補題）より、
ｚｋ＋１＝ｍｏｄ（ｆ ｍｏｄ（ａｋ＋１Ｆ－１、ｅ）＋
ｅ ｍｏｄ（ｂｋ＋１Ｅ－１、ｆ）、ｅｆ）
である。この右辺の内側の２つのｍｏｄ関数、法ｅと法ｆでのｍｏｄ関数、はそれぞれｆとｅと未満の正の整数を与え、外側のｍｏｄ関数の内部は２ｅｆ＝２６４を越えない。だから上のすべての計算は倍精度整数内で行なわれ、それを実数の有理数乱数とする計算も倍精度実数内で行なう事ができる。４倍精度計算を必要としないこのＭＣ乱数生成方法は、実用に高速の大きな利便を与える。
【0006】
用いられた（孫子の縮小）の内容は殆ど自明だが、与える技術便宜は魔術的で大きい。ここではこれを用いた異なる２つの素数の積を法とするＭＣ乱数の大きな利点の実証として、良い機会なので、中澤直也が発見し（文献１）（文献２）で公開された優れたＭＣ乱数生成機構＃００１について具体的に述べる。この乱数は２つの素数ｅとｆの積ｄ＝ｅｆ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５２．７８乗
で生成される。これらの詳しい組成は（文献１）（文献２）に開示されていて、下の乱数計算ｍａｉｎプログラムとｓｕｂｒｏｕｔｉｎｅも見て頂きたい。
＃００１生成機構の乱数発生ｍａｉｎ ｐｒｏｇｒａｍ
＃００１生成機構の乱数発生ｓｕｂｒｏｕｔｎｉｎｅ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５２．７８乗
で生成される。これらの詳しい組成は（文献１）（文献２）に開示されていて、下の乱数計算ｍａｉｎプログラムとｓｕｂｒｏｕｔｉｎｅも見て頂きたい。
＃００１生成機構の乱数発生ｍａｉｎ ｐｒｏｇｒａｍ
【0007】
前掲のプログラムが必要とするＣＰＵ時間は０．２８１秒／１０００万個である。出力は１／ｄ、約２－５４の精度だから倍精度実数として十分な桁数を持つ。この計算速度と倍精度とに比べる事は、次の＃Ｍ００１生成機構の持つ絶対的な性能の位置を与える。下は＃Ｍ００１で同じ１０００万個の乱数を生成するプログラムである。構造は一見大変簡単である。
＃Ｍ００１で１０００万個を出力するプログラム
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数の計算時間を大幅に増大させる。上のプログラムの実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は１／ｄ程度ほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないもの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。
＃Ｍ００１で１０００万個を出力するプログラム
【技術的展望】
【0008】
孫子の定理は一般にｋ個、ｋ＞２、の互いに素な正の数の積の法についても成り立つので、２個の積ｄ＝ｅｆに限るものではない。だから、３個以上の部分法を持つものでも考える事ができる。ただその様な場合に数値的に倍精度計算だけへの設計をどの様に選ぶべきか、撰ぶ事ができるか、は大きな問題である。それは必然的に２成分のｄ＝ｅｆの場合よりは複雑な考察を要するだろう。２成分の場合がきれいな構成を得て倍精度整数演算、倍精度実数演算だけで実現している事実、多分この構成が部分法での成功を最も得やすい形ではないかという予感、からは発明者達はこのより進んだ形は将来への課題ではないか、より大きな乱数周期への需要、より多い桁数精度の必要の生起、を見れば自然に得られるのでははないか、と考える。
【手続補正2】
【補正対象書類名】特許請求の範囲
【補正対象項目名】全文
【補正方法】変更
【補正の内容】
【特許請求の範囲】
【請求項1】
２つの共通素因数を持たない、即ち互いに素な、正の整数ｅとｆの積であるｄ＝ｅｆを法とし、ｄとは素な乗数ｚ、ｄとは素なｓｅｅｄ種ｎで構成される乗算合同法乱数生成機構（ｄ、ｚ、ｎ）において、整数の任意に与えられた初期値ｎから順次得られる漸化合同式による一様独立乱数｛ｖｋ：＝ｚｋ／ｄ｜ｋ＝０、１、２、…｝を得るための正の整数列｛ｚｋ｜ｋ＝０、１、２、…｝、
ｚ０＝ｍｏｄ（ｎ、ｄ）、
ｚｋ＋１＝ｍｏｄ（ｚｚｋ、ｄ）、 ｋ＝０、１、２、…、
を得る計算方法であって、孫子の定理の一つの形，部分法ｅとｆだけに対する成分
ｘ：＝ｍｏｄ（ｚ，ｅ）、
ｘｋ：＝ｍｏｄ（ｚｋ，ｅ）、ｋ＝０、１、２、…、
ｙ：＝ｍｏｄ（ｚ，ｆ）、
ｙｋ：＝ｍｏｄ（ｚｋ，ｆ）、ｋ＝０、１、２、…、
ａ：＝ｍｏｄ（ｎ，ｅ）
ｂ：＝ｍｏｄ（ｎ，ｆ）
を用いて定義される２組のＭＣ部分乱数列、ＭＣ生成機構（ｅ、ｘ、ａ）からの法ｅの合同漸化式の解
ｘ０＝ｍｏｄ（ａ、ｅ）、
ｘｋ＋１＝ｍｏｄ（ｘｘｋ、ｅ）、 ｋ＝０、１、２、…、
及びＭＣ生成機構（ｆ、ｙ、ｂ）からの法ｆの合同漸化式の解
ｙ０＝ｍｏｄ（ｂ、ｆ）、
ｙｋ＋１＝ｍｏｄ（ｙｙｋ、ｆ）、 ｋ＝０、１、２、…、
から得られる一次結合｛ｚｋ｜ｋ＝０、１、２、…｝：
ｚｋ：＝ｍｏｄ｛ｆ ｍｏｄ（ｘｋＦ－１、ｅ）＋
ｅ ｍｏｄ（ｙｋＥ－１、ｆ）、ｅｆ｝、
Ｆ－１：は法ｅでのｆの逆数、 ０＜Ｆ－１＜ｅ、
Ｅ－１：は法ｆでのｅの逆数、 ０＜Ｅ－１＜ｆ、
で与えてよい、という形の孫子の定理に基づき、さらに任意の正の整数Ａ、Ｂに対する関係
ｍｏｄ（ｅＡ，ｅｆ）＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）、
ｍｏｄ（ｆＢ，ｅｆ）＝ｍｏｄ（ｆ ｍｏｄ（Ｂ、ｅ）、ｅｆ）、
によって、２ｄ＝２ｅｆが倍精度整数限界２６４未満の場合に、すべての演算を倍精度整数、倍精度実数だけの範囲内で正確高速に実現する計算方法。
ｚ０＝ｍｏｄ（ｎ、ｄ）、
ｚｋ＋１＝ｍｏｄ（ｚｚｋ、ｄ）、 ｋ＝０、１、２、…、
を得る計算方法であって、孫子の定理の一つの形，部分法ｅとｆだけに対する成分
ｘ：＝ｍｏｄ（ｚ，ｅ）、
ｘｋ：＝ｍｏｄ（ｚｋ，ｅ）、ｋ＝０、１、２、…、
ｙ：＝ｍｏｄ（ｚ，ｆ）、
ｙｋ：＝ｍｏｄ（ｚｋ，ｆ）、ｋ＝０、１、２、…、
ａ：＝ｍｏｄ（ｎ，ｅ）
ｂ：＝ｍｏｄ（ｎ，ｆ）
を用いて定義される２組のＭＣ部分乱数列、ＭＣ生成機構（ｅ、ｘ、ａ）からの法ｅの合同漸化式の解
ｘ０＝ｍｏｄ（ａ、ｅ）、
ｘｋ＋１＝ｍｏｄ（ｘｘｋ、ｅ）、 ｋ＝０、１、２、…、
及びＭＣ生成機構（ｆ、ｙ、ｂ）からの法ｆの合同漸化式の解
ｙ０＝ｍｏｄ（ｂ、ｆ）、
ｙｋ＋１＝ｍｏｄ（ｙｙｋ、ｆ）、 ｋ＝０、１、２、…、
から得られる一次結合｛ｚｋ｜ｋ＝０、１、２、…｝：
ｚｋ：＝ｍｏｄ｛ｆ ｍｏｄ（ｘｋＦ－１、ｅ）＋
ｅ ｍｏｄ（ｙｋＥ－１、ｆ）、ｅｆ｝、
Ｆ－１：は法ｅでのｆの逆数、 ０＜Ｆ－１＜ｅ、
Ｅ－１：は法ｆでのｅの逆数、 ０＜Ｅ－１＜ｆ、
で与えてよい、という形の孫子の定理に基づき、さらに任意の正の整数Ａ、Ｂに対する関係
ｍｏｄ（ｅＡ，ｅｆ）＝ｍｏｄ（ｅ ｍｏｄ（Ａ、ｆ）、ｅｆ）、
ｍｏｄ（ｆＢ，ｅｆ）＝ｍｏｄ（ｆ ｍｏｄ（Ｂ、ｅ）、ｅｆ）、
によって、２ｄ＝２ｅｆが倍精度整数限界２６４未満の場合に、すべての演算を倍精度整数、倍精度実数だけの範囲内で正確高速に実現する計算方法。
【手続補正書】
【提出日】2022-06-03
【手続補正1】
【補正対象書類名】明細書
【補正対象項目名】０００４
【補正方法】変更
【補正の内容】
【0004】
前段落、素数の法と原始根の乗数の計算や周期の短さの困難は、ＭＣ乱数が深く考えられるようになって初めて、孫子の定理が魔術的な解決を与えていることが発見された。現在の発明はこの発見に基づくものである。理解、把握から再考する様、孫子の定理のｍｏｄ関数を用いた形から述べる。法ｄが２つの互いに素な２つの正整数の積ｄ＝ｅｆであるとする。孫子の定理の実現には様々な形を与えられるが、現在の議論に最も便利なのは下の（Ｂ）の形である。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の等式（＊＊）で与えられる：
ｘ＝ｍｏｄ（ａｆＦ－１＋ｂｅＥ－１、ｅｆ） （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。（注１）（証明）先に（Ｂ）の形の解（＊＊）は（＊）を満たす事を示す。最も紛れのない初等的証明として（＊＊）を等式で表す。（＊＊）の右辺はｄ＝ｅｆで割った余りなのだから、ある整数ｋを取って
ｘ＝ａｆＦ－１＋ｂｅＥ－１－ｋｅｆ
と表わす事ができる。（注２）これを法ｅで見ればｍｏｄ（ｆＦ－１、ｅ）＝１となって ｍｏｄ（ｘ、ｅ）＝ａ
が得られる。同様に法ｆで見ればｍｏｄ（ｅＥ－１、ｆ）＝１となって、
ｍｏｄ（ｘ、ｆ）＝ｂ
であるとわかる。
（Ａ）２つの合同式（＊）がもう一つの解ｙも持つとする：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
（＊）とこの式とは次の等式をある整数ｊ、ｋ、ｐ、ｑに対して意味する。
ｘ＝ａ＋ｊｅ、ｙ＝ａ＋ｋｅ
ｘ＝ｂ＋ｐｆ、ｙ＝ｂ＋ｑｆ
これらはまず等式ｘ－ｙ＝（ｊ－ｋ）ｅを与える。即ちｘ－ｙはｅの倍数である。同様にｘ－ｙ＝（ｐ－ｑ）ｆ、ｘ－ｙはｅとは素なｆの倍数でもあって、ｘとｙは法ｄで同一である。
（証明終り）
（注１）逆数Ｅ－１とＦ－１について説明する。整数ｅ＞０とは素な任意の正の整数ｘをｅで割るとき、その余りは１からｅ－１までの中にある。その総個数はオイラーの関数φ（ｅ）と呼ばれ、１からｅまでの整数からｅと共通な素因数を持つものを除いた個数、であり１＜φ（ｅ）＜ｅだが、これらの詳細はどうでもよい。これらを１列に、例えば昇順に
｛１、ａ、ｂ、…、ｆ、…、ｅ－１｝ （＄）
と並べる。これらの中の任意の整数ｘをｆに法ｅで掛けたものの列
｛ｘ、ｍｏｄ（ａｘ、ｅ）、ｍｏｄ（ｂｘ、ｅ）、…、
ｍｏｄ（ｆｘ、ｅ）、…，ｍｏｄ（（ｅ－１）ｘ）｝ （＃）
を考えると、これらは法ｅですべて異なる整数である。実際帰謬法で仮に
Ａ＝ｍｏｄ（ｐｘ、ｅ）＝ｍｏｄ（ｑｘ、ｅ）
となると仮定すると、ｐｘ＝Ａ＋ｊｅ，ｑｘ＝Ａ＋ｋｅとなる整数ｊ、ｋがある事から
ｐｘ－ｑｘ＝（ｐ－ｑ）ｘ＝（ｋ－ｊ）ｅ
でなければならず、ｘはｅとは共通素因数を持たないからｐ－ｑがｅの倍数でなければならない、即ちｐとｑが法ｅで等しくなければならない矛盾に至るからである。故に（＃）は元の（＄）の並べ替え（順列）であり、（＄）の中には必ずｍｏｄ（ｘｆ、ｅ）＝１となるものがある、このｘがｆの法ｅでの逆数Ｆ－１である。ｆの法ｅでの逆数Ｆ－１を求めるには、ｅ、ｆは単精度整数だから、
ｍｏｄ（ｆ、ｅ）、ｍｏｄ（ｆａ、ｅ）、
ｍｏｄ（ｆｂ、ｅ）、 ．．．、ｍｏｄ（ｆ（ｅ－１）、ｅ）
を順次生成して探しても困難ではない。乱数を使う前に十分時間を掛けて求めればよい。（注２）整数ｋは法ｄ＝ｅｆで割った商に当たるが、この詳細はどうでもよい。
（孫子の定理）（Ａ）正の整数ｅ、ｆが互いに素であれば、０以上ｄ＝ｅｆ未満の整数ｘに対する連立合同式
ａ＝ｍｏｄ（ｘ、ｅ） ｂ＝ｍｏｄ（ｘ、ｆ） （＊）
は整数ｘを法ｄ＝ｅｆで（０以上ｄ未満で）一意に定める。
（Ｂ）上の連立合同式（＊）に対して、その法ｄでの一意の解である整数ｘは次の等式（＊＊）で与えられる：
ｘ＝ｍｏｄ（ａｆＦ－１＋ｂｅＥ－１、ｅｆ） （＊＊）
ここでＦ－１は法ｅでｆの逆数になるもの、即ちｍｏｄ（Ｆ－１ｆ、ｅ）＝１となるものであり、Ｅ－１は法ｆでのｅの逆数、即ちｍｏｄ（Ｅ－１ｅ、ｆ）＝１となるものである。（注１）（証明）先に（Ｂ）の形の解（＊＊）は（＊）を満たす事を示す。最も紛れのない初等的証明として（＊＊）を等式で表す。（＊＊）の右辺はｄ＝ｅｆで割った余りなのだから、ある整数ｋを取って
ｘ＝ａｆＦ－１＋ｂｅＥ－１－ｋｅｆ
と表わす事ができる。（注２）これを法ｅで見ればｍｏｄ（ｆＦ－１、ｅ）＝１となって ｍｏｄ（ｘ、ｅ）＝ａ
が得られる。同様に法ｆで見ればｍｏｄ（ｅＥ－１、ｆ）＝１となって、
ｍｏｄ（ｘ、ｆ）＝ｂ
であるとわかる。
（Ａ）２つの合同式（＊）がもう一つの解ｙも持つとする：
ａ＝ｍｏｄ（ｙ、ｅ） ｂ＝ｍｏｄ（ｙ、ｆ）
（＊）とこの式とは次の等式をある整数ｊ、ｋ、ｐ、ｑに対して意味する。
ｘ＝ａ＋ｊｅ、ｙ＝ａ＋ｋｅ
ｘ＝ｂ＋ｐｆ、ｙ＝ｂ＋ｑｆ
これらはまず等式ｘ－ｙ＝（ｊ－ｋ）ｅを与える。即ちｘ－ｙはｅの倍数である。同様にｘ－ｙ＝（ｐ－ｑ）ｆ、ｘ－ｙはｅとは素なｆの倍数でもあって、ｘとｙは法ｄで同一である。
（証明終り）
（注１）逆数Ｅ－１とＦ－１について説明する。整数ｅ＞０とは素な任意の正の整数ｘをｅで割るとき、その余りは１からｅ－１までの中にある。その総個数はオイラーの関数φ（ｅ）と呼ばれ、１からｅまでの整数からｅと共通な素因数を持つものを除いた個数、であり１＜φ（ｅ）＜ｅだが、これらの詳細はどうでもよい。これらを１列に、例えば昇順に
｛１、ａ、ｂ、…、ｆ、…、ｅ－１｝ （＄）
と並べる。これらの中の任意の整数ｘをｆに法ｅで掛けたものの列
｛ｘ、ｍｏｄ（ａｘ、ｅ）、ｍｏｄ（ｂｘ、ｅ）、…、
ｍｏｄ（ｆｘ、ｅ）、…，ｍｏｄ（（ｅ－１）ｘ）｝ （＃）
を考えると、これらは法ｅですべて異なる整数である。実際帰謬法で仮に
Ａ＝ｍｏｄ（ｐｘ、ｅ）＝ｍｏｄ（ｑｘ、ｅ）
となると仮定すると、ｐｘ＝Ａ＋ｊｅ，ｑｘ＝Ａ＋ｋｅとなる整数ｊ、ｋがある事から
ｐｘ－ｑｘ＝（ｐ－ｑ）ｘ＝（ｋ－ｊ）ｅ
でなければならず、ｘはｅとは共通素因数を持たないからｐ－ｑがｅの倍数でなければならない、即ちｐとｑが法ｅで等しくなければならない矛盾に至るからである。故に（＃）は元の（＄）の並べ替え（順列）であり、（＄）の中には必ずｍｏｄ（ｘｆ、ｅ）＝１となるものがある、このｘがｆの法ｅでの逆数Ｆ－１である。ｆの法ｅでの逆数Ｆ－１を求めるには、ｅ、ｆは単精度整数だから、
ｍｏｄ（ｆ、ｅ）、ｍｏｄ（ｆａ、ｅ）、
ｍｏｄ（ｆｂ、ｅ）、 ．．．、ｍｏｄ（ｆ（ｅ－１）、ｅ）
を順次生成して探しても困難ではない。乱数を使う前に十分時間を掛けて求めればよい。（注２）整数ｋは法ｄ＝ｅｆで割った商に当たるが、この詳細はどうでもよい。
【手続補正書】
【提出日】2022-06-08
【手続補正1】
【補正対象書類名】明細書
【補正対象項目名】０００６
【補正方法】変更
【補正の内容】
【0006】
用いられた（孫子の縮小）の内容は殆ど自明だが、与える技術便宜は魔術的で大きい。ここではこれを用いた異なる２つの素数の積を法とするＭＣ乱数の大きな利点の実証として、良い機会なので、中澤直也が発見し（文献１）（文献２）で公開された優れたＭＣ乱数生成機構＃００１について具体的に述べる。この乱数は２つの素数ｅとｆの積ｄ＝ｅｆ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５２．７８乗
で生成される。これらの詳しい組成は（文献１）（文献２）に開示されていて、下の乱数計算ｍａｉｎプログラムとｓｕｂｒｏｕｔｉｎｅも見て頂きたい。
＃００１生成機構の乱数発生ｍａｉｎ ｐｒｏｇｒａｍ
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ，ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｉｐ１＝１３４２６５０２３
ｉｐ２＝１３４４７５８２７
ｉｄ＝ｉｐ１＊ｉｐ２ ！ｉｄはおよそ２５４です
ａｄ＝ｉｄ
ｉｚ１＝１９０６１２５２
ｉｚ２＝７７６００５２５
ｎ１＝１０ ！初期値１は１０です
ｎ２＝１３ ！初期値２は１３です
ｉｐ２ｍｐ１＝５２５７７００７
ｉｐ１ｍｐ２＝８１８１６２７１
ｍｚ１＝ｍｏｄ（ｎ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｎ２，ｉｐ２）
ｄｏ ｉ＝１，１０００００００
ｃａｌｌ ｒａｎｄｏｍ（ｒａｎｄ）
ｅｎｄ ｄｏ
ｅｎｄ
＃００１生成機構の乱数発生ｓｕｂｒｏｕｔｎｉｎｅ
ｓｕｂｒｏｕｔｉｎｅ ｒａｎｄｏｍ（ｒａｎｄ）
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ，ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｍｚ１＝ｍｏｄ（ｍｚ１＊ｉｚ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｍｚ２＊ｉｚ２，ｉｐ２）
ｍｚ１ａ＝ｍｏｄ（ｍｚ１＊ｉｐ２ｍｐ１，ｉｐ１）
ｍｚ２ａ＝ｍｏｄ（ｍｚ２＊ｉｐ１ｍｐ２，ｉｐ２）
ａｚ＝ｍｏｄ（ｉｐ２＊ｍｚ１ａ＋ｉｐ１＊ｍｚ２ａ，ｉｄ）
ｒａｎｄ＝ａｚ／ａｄ
ｒｅｔｕｒｎ
ｅｎｄ
ｄ＝１８０５５４００００５０９９０２１、約２５４乗
を法とし、乗数
ｚ＝７７５９０９７９５８７８２９３５、約２５２．７８乗
で生成される。これらの詳しい組成は（文献１）（文献２）に開示されていて、下の乱数計算ｍａｉｎプログラムとｓｕｂｒｏｕｔｉｎｅも見て頂きたい。
＃００１生成機構の乱数発生ｍａｉｎ ｐｒｏｇｒａｍ
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ，ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｉｐ１＝１３４２６５０２３
ｉｐ２＝１３４４７５８２７
ｉｄ＝ｉｐ１＊ｉｐ２ ！ｉｄはおよそ２５４です
ａｄ＝ｉｄ
ｉｚ１＝１９０６１２５２
ｉｚ２＝７７６００５２５
ｎ１＝１０ ！初期値１は１０です
ｎ２＝１３ ！初期値２は１３です
ｉｐ２ｍｐ１＝５２５７７００７
ｉｐ１ｍｐ２＝８１８１６２７１
ｍｚ１＝ｍｏｄ（ｎ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｎ２，ｉｐ２）
ｄｏ ｉ＝１，１０００００００
ｃａｌｌ ｒａｎｄｏｍ（ｒａｎｄ）
ｅｎｄ ｄｏ
ｅｎｄ
＃００１生成機構の乱数発生ｓｕｂｒｏｕｔｎｉｎｅ
ｓｕｂｒｏｕｔｉｎｅ ｒａｎｄｏｍ（ｒａｎｄ）
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｈ，ｏ－ｚ）
ｃｏｍｍｏｎ ｉｐ１，ｉｐ２，ｉｄ，ａｄ，ｉｚ１，ｉｚ２，ｍｚ１，ｍｚ２，ｉｐ２ｍｐ１，ｉｐ１ｍｐ２
ｍｚ１＝ｍｏｄ（ｍｚ１＊ｉｚ１，ｉｐ１）
ｍｚ２＝ｍｏｄ（ｍｚ２＊ｉｚ２，ｉｐ２）
ｍｚ１ａ＝ｍｏｄ（ｍｚ１＊ｉｐ２ｍｐ１，ｉｐ１）
ｍｚ２ａ＝ｍｏｄ（ｍｚ２＊ｉｐ１ｍｐ２，ｉｐ２）
ａｚ＝ｍｏｄ（ｉｐ２＊ｍｚ１ａ＋ｉｐ１＊ｍｚ２ａ，ｉｄ）
ｒａｎｄ＝ａｚ／ａｄ
ｒｅｔｕｒｎ
ｅｎｄ
【手続補正2】
【補正対象書類名】明細書
【補正対象項目名】０００７
【補正方法】変更
【補正の内容】
【0007】
前掲のプログラムが必要とするＣＰＵ時間は０．２８１秒／１０００万個である。出力は１／ｄ、約２－５４の精度だから倍精度実数として十分な桁数を持つ。この計算速度と倍精度とに比べる事は、次の＃Ｍ００１生成機構の持つ絶対的な性能の位置を与える。下は＃Ｍ００１で同じ１０００万個の乱数を生成するプログラムである。構造は一見大変簡単である。
＃Ｍ００１で１０００万個を出力するプログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｇ，ｏ－ｐ，ｒ－ｚ），ｒｅａｌ＊１６（ｑ）
ｉｐ＝１７１７９８６９９８９ｄ０ ！ｉｐはおよそ２３４．００
ｒｐ＝ｉｐ
ｉｚ＝７９２８４１００７２ｄ０
ｉｎ＝１０ ！初期値は１０
ｑｍｚ＝ｉｎ
ｄｏ ｉ＝１，１０００００００
ｒａｎｄ＝ｍｏｄ（ｑｍｚ，ｒｐ）／ｒｐ
ｑｍｚ＝ｍｏｄ（ｑｍｚ＊ｉｚ，ｒｐ）
ｅｎｄ ｄｏ
ｅｎｄ
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数の計算時間を大幅に増大させる。上のプログラムの実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は１／ｄ程度ほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないもの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。
＃Ｍ００１で１０００万個を出力するプログラム
ｍａｉｎ ｐｒｏｇｒａｍ
ｉｍｐｌｉｃｉｔ ｉｎｔｅｇｅｒ＊８（ｉ－ｎ），ｒｅａｌ＊８（ａ－ｇ，ｏ－ｐ，ｒ－ｚ），ｒｅａｌ＊１６（ｑ）
ｉｐ＝１７１７９８６９９８９ｄ０ ！ｉｐはおよそ２３４．００
ｒｐ＝ｉｐ
ｉｚ＝７９２８４１００７２ｄ０
ｉｎ＝１０ ！初期値は１０
ｑｍｚ＝ｉｎ
ｄｏ ｉ＝１，１０００００００
ｒａｎｄ＝ｍｏｄ（ｑｍｚ，ｒｐ）／ｒｐ
ｑｍｚ＝ｍｏｄ（ｑｍｚ＊ｉｚ，ｒｐ）
ｅｎｄ ｄｏ
ｅｎｄ
この簡潔なプログラムは、しかし、ｒｅａｌ＊１６（４倍精度実数）を含むｍｏｄ関数の計算時間を大幅に増大させる。上のプログラムの実行ＣＰＵ時間は＃００１のほぼ１０倍、２．６２５秒となる。出力される乱数精度は１／ｄ程度ほぼ２－３４、周期はほぼｄ／２で２３３、このＭＣ乱数専用にコンピュータを動かすと６分ほどで使い切られる。現在のシミュレーションに用いるには性能も精度も不足である。孫子の定理に基づく＃００１生成機構や、２つの異なる素数、或いは２つの素数のマイナス、で２３１を超えないもの積を法とするＭＣ乱数生成機構がこの問題の現在の最適解である。