LaTeX: メルセンヌ・ツイスタ

mt19937.tex

% https://h20y6m.github.io/posts/2023-12-11-texadvent2023-2/

\documentclass[a4paper,twocolumn]{article}

% https://github.com/h20y6m/latex-biguint
\usepackage{biguint}

\ExplSyntaxOn
% https://en.wikipedia.org/wiki/Mersenne_Twister

% constant n, m
\int_const:Nn \c_mt_n_int { 624 }
\int_const:Nn \c_mt_m_int { 397 }
% constant a
\biguint_const:Nn \c_mt_a_biguint { "9908 B0DF }
% constant u, d
\biguint_const:Nn \c_mt_u_biguint { 11 }
\biguint_const:Nn \c_mt_d_biguint { "FFFF FFFF }
% constant s, b
\biguint_const:Nn \c_mt_s_biguint { 7 }
\biguint_const:Nn \c_mt_b_biguint { "9D2C 5680 }
% constant t, c
\biguint_const:Nn \c_mt_t_biguint { 15 }
\biguint_const:Nn \c_mt_c_biguint { "EFC6 0000 }
% constant l
\biguint_const:Nn \c_mt_l_biguint { 18 }
% constant f
\biguint_const:Nn \c_mt_f_biguint { 1812433253 }

% state MT[0..(n-1)]
\int_step_inline:nnn { 0 } { \c_mt_n_int - 1 }
  { \biguint_new:c { g_mt_MT[ #1 ]_biguint } }
% index
\int_new:N \g_mt_index_int
\int_gset:Nn \g_mt_index_int { \c_mt_n_int + 1}

\biguint_new:N \g_mt_result_biguint

\biguint_const:Nn \c_mt_lower_mask_biguint { "7FFF FFFF }
\biguint_const:Nn \c_mt_upper_mask_biguint { "8000 0000 }
\biguint_const:Nn \c_mt_wbits_mask_biguint { "FFFF FFFF }
\biguint_const:Ne \c_mt_w_minus_two_biguint { \int_eval:n { 32 - 2 } }

\biguint_new:N \l__mt_tempa_biguint
\biguint_new:N \l__mt_tempb_biguint
\biguint_new:N \l__mt_y_biguint
\biguint_new:N \l__mt_x_biguint
\biguint_new:N \l__mt_xA_biguint

\cs_generate_variant:Nn \biguint_shr:NNN { NcN }
\cs_generate_variant:Nn \biguint_bitand:NNN { NcN }
\cs_generate_variant:Nn \biguint_gbitand:NNN { cNN }
\cs_generate_variant:Nn \biguint_bitxor:NNN { NcN }
\cs_generate_variant:Nn \biguint_gbitxor:NNN { ccN }

\msg_new:nnn { mt } { never-seeded }
  { Generator~was~never~seeded. }

\cs_new_protected:Npn \mt_seed_mt:n #1
  {
    \int_gset_eq:NN \g_mt_index_int \c_mt_n_int
    \biguint_gset:cn  { g_mt_MT[ 0 ]_biguint } { #1 }
    \int_step_inline:nnn { 1 } { \c_mt_n_int - 1 }
      {
        % MT[i-1] >> (w - 2)
        \biguint_shr:NcN
          \l__mt_tempa_biguint
          { g_mt_MT[ \int_eval:n { ##1 - 1 } ]_biguint }
          \c_mt_w_minus_two_biguint
        % MT[i-1] xor (...)
        \biguint_bitxor:NcN
          \l__mt_tempa_biguint
          { g_mt_MT[ \int_eval:n { ##1 - 1 } ]_biguint }
          \l__mt_tempa_biguint
        % f * (...)
        \biguint_mul:NNN
          \l__mt_tempa_biguint
          \c_mt_f_biguint
          \l__mt_tempa_biguint
        % ... + i
        \biguint_set:Nn \l__mt_tempb_biguint { ##1 }
        \biguint_add:NNN
          \l__mt_tempa_biguint
          \l__mt_tempa_biguint
          \l__mt_tempb_biguint
        % MT[i] := lowest w bits of ...
        \biguint_gbitand:cNN
          { g_mt_MT[ ##1 ]_biguint }
          \l__mt_tempa_biguint
          \c_mt_wbits_mask_biguint
      }
  }

\cs_new_protected:Npn \mt_extract_number:
  {
    % if index >= n
    \int_compare:nNnF { \g_mt_index_int } < { \c_mt_n_int }
      {
        \int_compare:nNnT { \g_mt_index_int } > { \c_mt_n_int }
          {
            \msg_error:nn { mt } { never-seeded }
          }
        \mt_twist:
      }

    % y := MT[index]
    \biguint_set_eq:Nc
      \l__mt_y_biguint
      { g_mt_MT[ \int_use:N \g_mt_index_int ]_biguint }

    % y := y xor ((y >> u) and d)
    \biguint_shr:NNN
      \l__mt_tempa_biguint \l__mt_y_biguint \c_mt_u_biguint
    \biguint_bitand:NNN
      \l__mt_tempa_biguint \l__mt_tempa_biguint \c_mt_d_biguint
    \biguint_bitxor:NN
      \l__mt_y_biguint \l__mt_tempa_biguint
    % y := y xor ((y << s) and b)
    \biguint_shl:NNN
      \l__mt_tempa_biguint \l__mt_y_biguint \c_mt_s_biguint
    \biguint_bitand:NNN
      \l__mt_tempa_biguint \l__mt_tempa_biguint \c_mt_b_biguint
    \biguint_bitxor:NN
      \l__mt_y_biguint \l__mt_tempa_biguint
    % y := y xor ((y << t) and c)
    \biguint_shl:NNN
      \l__mt_tempa_biguint \l__mt_y_biguint \c_mt_t_biguint
    \biguint_bitand:NNN
      \l__mt_tempa_biguint \l__mt_tempa_biguint \c_mt_c_biguint
    \biguint_bitxor:NN
      \l__mt_y_biguint \l__mt_tempa_biguint
    % y := y xor (y >> l)
    \biguint_shr:NNN
      \l__mt_tempa_biguint \l__mt_y_biguint \c_mt_l_biguint
    \biguint_bitxor:NN
      \l__mt_y_biguint \l__mt_tempa_biguint

    % index := index + 1
    \int_gincr:N \g_mt_index_int

    \biguint_gset_eq:NN
      \g_mt_result_biguint \l__mt_y_biguint
  }

\cs_new_protected:Npn \mt_twist:
  {
    \int_step_inline:nnn { 0 } { \c_mt_n_int - 1 }
      {
        % MT[i] and upper_mask
        \biguint_bitand:NcN
          \l__mt_x_biguint
          { g_mt_MT[ ##1 ]_biguint }
          \c_mt_upper_mask_biguint
        % MT[(i+1) mod n] and lower_mask
        \biguint_bitand:NcN
          \l__mt_tempa_biguint
          { g_mt_MT[ \int_eval:n { \int_mod:nn { ##1 + 1 } { \c_mt_n_int } } ]_biguint }
          \c_mt_lower_mask_biguint
        % x := ... | ...
        \biguint_bitor:NN \l__mt_x_biguint \l__mt_tempa_biguint
        % xA := x >> 1
        \biguint_shr:NNN \l__mt_xA_biguint \l__mt_x_biguint \c_one_biguint
        % if (x mod 2) != 0
        \biguint_if_odd:NT \l__mt_x_biguint
          {
            % xA := xA xor a
            \biguint_bitxor:NN \l__mt_xA_biguint \c_mt_a_biguint
          }
        % MT[i] := MT[(i + m) mod n] xor xA
        \biguint_gbitxor:ccN
          { g_mt_MT[ ##1 ]_biguint }
          { g_mt_MT[ \int_eval:n { \int_mod:nn { ##1 + \c_mt_m_int } { \c_mt_n_int } } ]_biguint }
          \l__mt_xA_biguint
      }
    % index := 0
    \int_gzero:N \g_mt_index_int
  }

\cs_new_protected:Npn \mt_show:
  { \mt__dump:N \tl_show:N }
\cs_new_protected:Npn \mt_log:
  { \mt__dump:N \tl_log:N }
\tl_new:N \l__mt_dump_tl
\cs_new_protected:Npn \mt__dump:N #1
  {
    \tl_clear:N \l__mt_dump_tl
    \int_step_inline:nnn { 0 } { \c_mt_n_int - 1 }
      {
        \tl_put_right:Ne \l__mt_dump_tl
          { \biguint_to_decimal:c { g_mt_MT[ ##1 ]_biguint } ~ }
      }
    \tl_put_right:Ne \l__mt_dump_tl
      { \int_use:N \g_mt_index_int }
    #1 \l__mt_dump_tl
  }

% LaTeX2e interface
\NewDocumentCommand \MTseed { m }
  { \mt_seed_mt:n {#1} }
\NewDocumentCommand \MTrand {}
  { \mt_extract_number: }
\NewExpandableDocumentCommand \MTnumber {}
  { \biguint_to_decimal:N \g_mt_result_biguint }

\ExplSyntaxOff

\usepackage{pgffor}

\begin{document}
\MTseed{5489}
\foreach \x in {1,...,10000}
  {%
    \MTrand
    $ [\x] = \MTnumber $%
    \par
  }%
\end{document}