mt19937ar の一部を Modern Fortran で書いてみた
統計計算を行うのに何かと乱数が必要になります。Fortran には組み込み関数として random_number() がありますが、
コンパイラによって実装が異なります*1*2。環境によって品質が変わると困るので、自分で実装したものを使うのが安全だと思います。私はメルセンヌ ツイスターを用いています。これまでは mt19937 をモジュールとして実装していましたが、この度 mt19937ar をクラスにしてみました。記事の表題に Modern Fortran とあるのはモジュールやクラスを使ったから、くらいの意味です*3。
以下に mt19937ar.c を移植したものを示します*4。都合により、mt19937ar.c にある init_by_array() および genrand_int31() は実装していません。また、関数の名称を変えています。
コメントはできるだけ付けました。しかしながら、私がメルセンヌ ツイスターの数学的な部分を理解していないので、あまり気の利いたコメントが書けませんでした。コメントの書式は Javadoc 風ですが、ドキュメント生成ツールを使うためではなく、変数の説明を書くのにちょうどよい書式だと思ったからです。
実装ミスをはじめ、間違い等ございましたら、ご指摘くださいませ。
! mt19937ar.c を Fortran へ移植したプログラムです。 ! ただし、init_by_array() および genrand_int31() は実装していません。 ! ! [0.0, 1.0]、[0.0, 1.0)、(0.0, 1.0) からの一様乱数を生成します。 ! ! 元のプログラムは以下のリンク先にあります。 ! http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c ! ! ! ----- 移植元ライセンス ここから----- !/* ! A C-program for MT19937, with initialization improved 2002/1/26. ! Coded by Takuji Nishimura and Makoto Matsumoto. ! ! Before using, initialize the state by using init_genrand(seed) ! or init_by_array(init_key, key_length). ! ! Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: ! ! 1. Redistributions of source code must retain the above copyright ! notice, this list of conditions and the following disclaimer. ! ! 2. Redistributions in binary form must reproduce the above copyright ! notice, this list of conditions and the following disclaimer in the ! documentation and/or other materials provided with the distribution. ! ! 3. The names of its contributors may not be used to endorse or promote ! products derived from this software without specific prior written ! permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, ! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, ! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR ! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF ! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING ! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! ! Any feedback is very welcome. ! http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html ! email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) !*/ ! ----- 移植元ライセンス ここまで----- ! module Class_mt19937ar use, intrinsic :: iso_fortran_env implicit none private ! パラメータ integer(int32), parameter :: n = 624 integer(int32), parameter :: m = 397 ! 計算に使うために定義したパラメータです。 integer(int32), parameter :: n_1 = n - 1 integer(int32), parameter :: n_2 = n - 2 integer(int32), parameter :: m_1 = m - 1 integer(int32), parameter :: n_m_1 = n - m - 1 integer(int32), parameter :: n_m = n - m integer(int32), parameter :: m_n = m - n ! 2147483647 は 0x7fffffff の符号付き 10 進表記です。 integer(int32), parameter :: lowerMask = 2147483647 ! 0x80000000 の符号付き 10 進表記です。 integer(int32), parameter :: upperMask = -lowerMask - 1 ! mt19937 による一様乱数生成器です。 type, public :: mt19937ar integer(int32), private :: mt(0:n-1) integer(int32), private :: mti contains ! [0, 1] から 32 bit 精度の一様乱数を生成します。 procedure, pass :: generateOnCC ! [0, 1) から 32 bit 精度の一様乱数を生成します。 procedure, pass :: generateOnCO ! (0, 1) から 32 bit 精度の一様乱数を生成します。 procedure, pass :: generateOnOO ! [0, 1) から 53 bit 精度の一様乱数を生成します。 procedure, pass :: generateOnCO_53bit end type ! コンストラクタです。 interface new_mt19937ar module procedure :: new_mt19937ar end interface ! public にする関数 public new_mt19937ar contains ! コンストラクタです。シード値を指定して乱数生成器を初期化します。 ! ! @param inutSeed 乱数シード ! @return 乱数生成器 function new_mt19937ar(inputSeed) result(newGenerator) ! 入力 integer(int32), optional, intent(in) :: inputSeed ! 変数 integer(int32) seed ! 戻り値 type(mt19937ar) newGenerator ! シード値を決定します。 if(present(inputSeed)) then ! 指定されている場合は、指定の値をシード値とします。 seed = inputSeed else ! 指定されていない場合はデフォルト値 5489 を用います。 seed = 5489 end if ! -1 は 0xffffffffUL の符号あり 10 進表記です。 newGenerator%mt(0) = iand(seed, -1) block ! ループ カウンター integer(int32) i do i = 1, n_1 newGenerator%mt(i) = 1812433253 & & * (ieor(newGenerator%mt(i-1), ishft(newGenerator%mt(i-1), -30))) & & + i ! -1 は 0xffffffffUL の符号あり 10 進表記です。 newGenerator%mt(i) = iand(newGenerator%mt(i), -1); end do newGenerator%mti = n end block end function ! [0x0, 0xffffffff] の範囲から符号あり整数の乱数を返します。 ! ! @param this 乱数生成器 ! @return [0x0, 0xffffffff] の範囲からの符号あり整数の乱数 function generate_int32(this) result(y) ! 引数 class(mt19937ar), intent(inout) :: this ! 戻り値 integer(int32) y if (this%mti >= n) then call generateNWordsAtOneTime(this, y) end if y = this%mt(this%mti) this%mti = this%mti + 1 ! Tempering ! ishft(i, j) で j = n とすると左に n bit シフトします。 ! ishft(i, j) で j = -n とすると右に n bit シフトします。 y = ieor(y, ishft(y, -11)) ! -1658038656 は 9d2c5680 の符号あり 10 進表記です。 y = ieor(y, iand(ishft(y, 7), -1658038656)) ! -272236544 は efc60000 の符号あり 10 進表記です。 y = ieor(y, iand(ishft(y, 15), -272236544)) y = ieor(y, ishft(y, -18)) end function ! mti >= n の時に、内部状態を更新します。 ! ! @param this 乱数生成器 ! @param y 符号あり整数 subroutine generateNWordsAtOneTime(this, y) ! 引数 class(mt19937ar), intent(inout) :: this integer(int32), intent(inout) :: y ! パラメータ ! -1727483681 は 9908b0df の符号なし 10 進表記です。 integer(int32), parameter :: matrixA = -1727483681 integer(int32), parameter :: mag01(0:1) = (/0, matrixA/) ! ループ カウンター integer(int32) kk do kk = 0, n_m_1 y = ior(iand(this%mt(kk) , upperMask), & & iand(this%mt(kk+1), lowerMask)) this%mt(kk) = ieor(ieor(this%mt(kk+m), ishft(y, -1)), & & mag01(iand(y, 1))) end do do kk = n_m, n_2 y = ior(iand(this%mt(kk) , upperMask), & & iand(this%mt(kk+1), lowerMask)) this%mt(kk) = ieor(ieor(this%mt(kk+m_n), ishft(y, -1)), & & mag01(iand(y, 1))) end do y = ior(iand(this%mt(n_1), upperMask), iand(this%mt(0), lowerMask)) this%mt(n_1) = ieor(ieor(this%mt(m_1), ishft(y, -1)), & & mag01(iand(y, 1))) this%mti = 0 end subroutine ! 以下、[0, 1]、[0, 1)、(0, 1) から一様乱数を発生させる関数です。 ! 関数名の on の後につく O は 端点を含まないことを (open の意味) 、 ! C は端点を含むことを (close の意味) 表します。 ! ! 例えば、 ! generateOnCO --- [0, 1) からの一様乱数 ! generateOnOO --- (0, 1) からの一様乱数 ! を表します。 ! [0, 1] から 32 bit 精度の一様乱数を発生させます。 ! ! @param this 乱数生成器 ! @return [0, 1] からの一様乱数 function generateOnCC(this) result(outputRandomNumber) ! 引数 class(mt19937ar), intent(inout) :: this ! 変数 integer(int32) intRandomNumber real(real64) convertedIntRandomNumber ! 戻り値 real(real64) outputRandomNumber ! 符号ありの整数乱数を取得します。 intRandomNumber = generate_int32(this) ! 符号なし整数としての値を、実数として取得します。 convertedIntRandomNumber & & = convertIntegerToUnsignedIntegerAsReal(intRandomNumber) ! [0, 1] の値に変換します。 ! 4294967295 は 2^32-1 です。 outputRandomNumber = convertedIntRandomNumber / 4294967295.0d0 end function ! [0, 1) から 32 bit 精度の一様乱数を発生させます。 ! ! @param this 乱数生成器 ! @return [0, 1) からの一様乱数 function generateOnCO(this) result(outputRandomNumber) ! 引数 class(mt19937ar), intent(inout) :: this ! 変数 integer(int32) intRandomNumber real(real64) convertedIntRandomNumber ! 戻り値 real(real64) outputRandomNumber ! 符号ありの整数乱数を取得します。 intRandomNumber = generate_int32(this) ! 符号なし整数としての値を、実数として取得します。 convertedIntRandomNumber & & = convertIntegerToUnsignedIntegerAsReal(intRandomNumber) ! [0, 1) の値に変換します。 ! 4294967295 は 2^32 です。 outputRandomNumber = convertedIntRandomNumber / 4294967296.0d0 end function ! [0, 1) から 53 bit 精度の一様乱数を発生させます。 ! ! @param this 乱数生成器 ! @return [0, 1) からの一様乱数 function generateOnCO_53bit(this) result(outputRandomNumber) ! 引数 class(mt19937ar), intent(inout) :: this ! 変数 integer(int32) a_int integer(int32) b_int real(real64) a_double real(real64) b_double ! 戻り値 real(real64) outputRandomNumber ! 符号ありの整数乱数を取得します。 a_int = ishft(generate_int32(this), -5) b_int = ishft(generate_int32(this), -6) ! 符号なし整数としての値を、実数として取得します。 a_double = convertIntegerToUnsignedIntegerAsReal(a_int) b_double = convertIntegerToUnsignedIntegerAsReal(b_int) ! [0, 1) の値に変換します。 ! 4294967295 は 2^32 です。 outputRandomNumber = (a_double * 67108864.0 + b_double) & & / 9007199254740992.0d0 end function ! (0, 1) から 32 bit 精度の一様乱数を発生させます。 ! ! @param this 乱数生成器 ! @return (0, 1) からの一様乱数 function generateOnOO(this) result(outputRandomNumber) ! 引数 class(mt19937ar), intent(inout) :: this ! 変数 integer(int32) intRandomNumber real(real64) convertedIntRandomNumber ! 戻り値 real(real64) outputRandomNumber ! 符号ありの整数乱数を取得します。 intRandomNumber = generate_int32(this) ! 符号なし整数としての値を、実数として取得します。 convertedIntRandomNumber = & & convertIntegerToUnsignedIntegerAsReal(intRandomNumber) ! (0, 1) の値に変換します。 ! 4294967295 は 2^32 です。 outputRandomNumber = (convertedIntRandomNumber + 0.5d0) & & / 4294967296.0d0 end function ! 符号あり整数を、符号なし整数を倍精度実数に変換した値にします。 ! Fortran には unsigned long がないため、この関数で unsigned long ! としての値に変換します。 ! ! @param inputInteger 整数 ! @return 符号あり整数を倍精度実数に変換した値 function convertIntegerToUnsignedIntegerAsReal(inputInteger) result(outputReal) ! 引数 integer(int32), intent(in) :: inputInteger ! 戻り値 real(real64) outputReal if (inputInteger < 0) then ! 4294967296 は 2^32 です。 outputReal = dble(inputInteger) + 4294967296.d0 else outputReal = dble(inputInteger) end if end function end module ! メイン関数です。 ! 使用例として、乱数生成器から 4 種類の乱数を生成します。 program MainProgram use, intrinsic :: iso_fortran_env use Class_mt19937ar implicit none ! 変数 ! 乱数生成器 type(mt19937ar) generator integer(int32) file_randomNumber ! ループ カウンター integer(int32) i ! 乱数生成器を生成します。 generator = new_mt19937ar(5489) open(newunit=file_randomNumber, file='fortran_mt19937ar_CC.csv', status='replace') ! 乱数をファイルに出力します。 do i = 1, 10000 write(file_randomNumber, *) generator%generateOnCC() end do close(file_randomNumber) ! 乱数生成器のシードを再設定します。 generator = new_mt19937ar(5489) open(newunit=file_randomNumber, file='fortran_mt19937ar_CO.csv', status='replace') ! 乱数をファイルに出力します。 do i = 1, 10000 write(file_randomNumber, *) generator%generateOnCO() end do close(file_randomNumber) ! 乱数生成器のシードを再設定します。 generator = new_mt19937ar(5489) open(newunit=file_randomNumber, file='fortran_mt19937ar_OO.csv', status='replace') ! 乱数をファイルに出力します。 do i = 1, 10000 write(file_randomNumber, *) generator%generateOnOO() end do close(file_randomNumber) ! 乱数生成器のシードを再設定します。 generator= new_mt19937ar(5489) open(newunit=file_randomNumber, file='fortran_mt19937ar_CO_53bit.csv', status='replace') ! 乱数をファイルに出力します。 do i = 1, 10000 write(file_randomNumber, *) generator%generateOnCO_53bit() end do close(file_randomNumber) end program
*1:例えば、GCC 6.5 では KISS algorithm ですが、GCC 7.4 では xorshift1024* になっています。 gcc.gnu.org gcc.gnu.org
*2:様々な乱数を信頼できるアルゴリズムで実装した標準ライブラリを提供していないのは、 科学技術計算用言語として致命的だと思います。
*3:Modern Fortran の厳密な定義は知りません。
*4:mt19937ar の説明は以下にあります。同リンク先の mt19937ar.c が移植元のプログラムです。 www.math.sci.hiroshima-u.ac.jp
