マクロツイーター

はてダから移行した記事の表示が崩れてますが、そのうちに直せればいいのに(えっ)

それでも TeX で「ガウス素数の模様」したい人のための何か (1)

【頼りなく宣伝】
2012/12/01 〜 2012/12/25

TeX & LaTeX Advent Calendar

こちらは(ry

*  *  *

「変換して TeX プログラミングの例」シリーズの最終回。今回は「ガウス素数の模様」の図を TeX で描いてみる。

ガウス素数の模様って何?

Wikipedia*1を参照。右にある図がその模様である。ということで説明終わり。

……でよいと思うがもう少し説明してみる。素数というのは「自明でない約数を持たない自然数」のことで、列挙すると次のようになるのであった。

2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, ……

いま、1 から順に自然数を数えて、それが素数であれば黒、なければ白の四角を並べると次のようになる。

□■■□■□■□□□■□■□□□■□■□□□■□□□□□■□■□□□□□■□□□……

正の実数に対して定義されていた「整数」「素数」を拡張して、複素数の世界について《整数》や《素数》という概念を定義することができ、それぞれ「ガウス整数」「ガウス素数」と呼ばれる。*2その《素数》から上の図に相当するものを作ると、二次元の上下左右に無限に続く「模様」になる。それが「ガウス素数の模様」である。

Lua で「ガウス素数の模様」してみた

いつも通り、グローバル変数のみの Lua プログラム。

[pascaltrig0.lua]
-- 整数型: n, rn, x, y, px, py, qx, qy
-- 配列型: a

-- gptexture(整数n)
-- 絶対値がn以下の範囲での模様を出力する.
function gptexture(_1)
  n = _1
  if n < 1 then n = 1 end
  gpt_setrn()
  gpt_init()
  gpt_sieve()
  gpt_show()
end

-- gpt_setrn()
-- nの平方根を整数に切り捨てた値をrnに代入する.
function gpt_setrn()
  rn = math.floor(math.sqrt(n))
end

-- gpt_init()
-- 篩の配列の初期化.
-- ※篩の2次元配列 a の要素は最初は全て 1 にしておき, x+yi が素数で
-- ないと判明したら, a[x][y] を 0 に変える.
function gpt_init()
  a = {}
  for x = 0, n do
    a[x] = {}
    for y = 0, n do
      a[x][y] = 1
    end
  end
  a[0][0] = 0
  a[1][0] = 0; a[0][1] = 0
end

-- gpt_sieve()
-- 篩のアルゴリズムを実行する.
function gpt_sieve()
  for x = 1, rn do
    for y = 0, x do
      if a[x][y] > 0 then
        px, py = x, y
        while px <= n do
          qx, qy = px, py
          while qx >= 0 and qy <= n do
            a[qx][qy] = 0; a[qy][qx] = 0
            qx, qy = qx - y, qy + x
          end
          qx, qy = px, py
          while qx <= n and qy >= 0 do
            a[qx][qy] = 0; a[qy][qx] = 0
            qx, qy = qx + y, qy - x
          end
          px, py = px + x, py + y
        end
        a[x][y] = 1; a[y][x] = 1
      end
    end
  end
end

-- gpt_show()
-- 篩の配列 a にある結果に基づいて模様を出力する.
-- (-n,-n)から(n,n)までの範囲を, (2n+1)行(2n+1)列の文字の
-- 集まりで表す. 絶対値がn以下でかつ素数ならば「#」、
-- それ以外は空白となる.
function gpt_show()
  for x = -n, n do
    for y = -n, n do
      px, py = math.abs(x), math.abs(y)
      if a[px][py] > 0 and x * x + y * y <= n * n then
        io.write("#")
      else
        io.write(" ")
      end
    end
    io.write("\n")
  end
end

-- メイン
gptexture(39)

これを実行すると以下のような「模様」が出力される。

                                                                               
                                #   #     #   #                                
                               #     #   #     #                               
                                  #   # #   #                                  
                       #       #     #   #     #       #                       
                        #   # #   #         #   # #   #                        
                   #           #     #   #     #           #                   
                        # #       # #     # #       # #                        
                   #   #     #   # #   #   # #   #     #   #                   
                #         # #                     # #         #                
                       #     #   # #       # #   #     #                       
              #         # #       #         #       # #         #              
                 # #         #       #   #       #         # #                 
              #   #         # #   #   # #   #   # #         #   #              
           # # # #     # # #     # #       # #     # # #     # # # #           
              #     #             #   # #   #             #     #              
         #       # # #     #   #       #       #   #     # # #       #         
            # # #     # # #       #         #       # # #     # # #            
             #               #     #       #     #               #             
      # #   #   #   #     # #   #   # # # #   #   # #     #   #   #   # #      
               #   #   # #   #   #     #     #   #   # #   #   #               
                #     #         # #         # #         #     #                
                 #   #     # # #     #   #     # # #     #   #                 
    #   # #   #     #         #   #   # #   #   #         #     #   # #   #    
     # #   #     #       #         # #   # #         #       #     #   # #     
              #     #   #   # #       # #       # #   #   #     #              
       # # #     # #       # # #     #   #     # # #       # #     # # #       
              # #     #   #     #             #     #   #     # #              
     #   #   #     #     #       # #   #   # #       #     #     #   #   #     
        # # #     # # #   #   # #   # # # #   # #   #   # # #     # # #        
     #       #         # #   #     #       #     #   # #         #       #     
  # # #         #     #   #     # # #     # # #     #   #     #         # # #  
 #                 # #     # # #     # # #     # # #     # #                 # 
        # #   #     #       #     #   # #   #     #       #     #   # #        
   # # #   # # # #   # #       # # # #   # # # #       # #   # # # #   # # #   
        # #   #   #     #   # #   #   # #   #   # #   #     #   #   # #        
 #     #           #         # #     # # #     # #         #           #     # 
  # # #     #         # # #     # # # # # # # #     # # #         #     # # #  
   #         # #   #   # #   #   # # ## ## # #   #   # #   #   # #         #   
        #       #   #       #   #   #     #   #   #       #   #       #        
   #         # #   #   # #   #   # # ## ## # #   #   # #   #   # #         #   
  # # #     #         # # #     # # # # # # # #     # # #         #     # # #  
 #     #           #         # #     # # #     # #         #           #     # 
        # #   #   #     #   # #   #   # #   #   # #   #     #   #   # #        
   # # #   # # # #   # #       # # # #   # # # #       # #   # # # #   # # #   
        # #   #     #       #     #   # #   #     #       #     #   # #        
 #                 # #     # # #     # # #     # # #     # #                 # 
  # # #         #     #   #     # # #     # # #     #   #     #         # # #  
     #       #         # #   #     #       #     #   # #         #       #     
        # # #     # # #   #   # #   # # # #   # #   #   # # #     # # #        
     #   #   #     #     #       # #   #   # #       #     #     #   #   #     
              # #     #   #     #             #     #   #     # #              
       # # #     # #       # # #     #   #     # # #       # #     # # #       
              #     #   #   # #       # #       # #   #   #     #              
     # #   #     #       #         # #   # #         #       #     #   # #     
    #   # #   #     #         #   #   # #   #   #         #     #   # #   #    
                 #   #     # # #     #   #     # # #     #   #                 
                #     #         # #         # #         #     #                
               #   #   # #   #   #     #     #   #   # #   #   #               
      # #   #   #   #     # #   #   # # # #   #   # #     #   #   #   # #      
             #               #     #       #     #               #             
            # # #     # # #       #         #       # # #     # # #            
         #       # # #     #   #       #       #   #     # # #       #         
              #     #             #   # #   #             #     #              
           # # # #     # # #     # #       # #     # # #     # # # #           
              #   #         # #   #   # #   #   # #         #   #              
                 # #         #       #   #       #         # #                 
              #         # #       #         #       # #         #              
                       #     #   # #       # #   #     #                       
                #         # #                     # #         #                
                   #   #     #   # #   #   # #   #     #   #                   
                        # #       # #     # #       # #                        
                   #           #     #   #     #           #                   
                        #   # #   #         #   # #   #                        
                       #       #     #   #     #       #                       
                                  #   # #   #                                  
                               #     #   #     #                               
                                #   #     #   #                                
                                                                               

このプログラムでは「エラトステネスの篩〈ふるい〉」をガウス整数に拡張したアルゴリズムを使っている。「篩」を表す配列は(複素数なので)2 次元のものが必要になる。これまで配列は全て名前参照として実現していたが、2 次元の配列も名前参照として容易に扱える。例えば次のように変換すればよい。

a[2][3]  →  _G["a/2/3"]
a[x][y]  →  _G["a/"..x.."/"..y]

この後者を TeX に直すと \@nameuse{tcgp@a/\number\tcgp@x/\number\tcgp@y} という非常に長いコードになる。これが何度も出てくるのは嫌なので、TeX ではマクロで表すことにしたい。そこで、「a[x][y]=1 の代入」「a[x][y]=0 の代入」「a[x][y] の参照」を set(x,y)reset(x,y)cur(x,y) という関数で表すことにする。(後掲の gptexture1.lua の (※1) の箇所。)

最終的な Lua 内変換を行った後のプログラムコードを示す。

-- 整数型: n, rn, x, y, px, py, qx, qy, t
-- 配列型: a

-- gptexture(整数n)
function gptexture(_1)
  n = _1
  if n < 1 then n = 1 end
  gpt_setrn()
  gpt_init()
  gpt_sieve()
  gpt_show()
end

-- (※1)
function set(_1,_2)   _G["a/".._1.."/".._2] = 1 end
function reset(_1,_2) _G["a/".._1.."/".._2] = 0 end
function cur(_1,_2)   return _G["a/".._1.."/".._2] end

-- gpt_setrn() (※2)
function gpt_setrn()
  rn = 0; t = 0
  while not (t > n) do
    rn = rn + 1; t = rn; t = t * rn
  end
  rn = rn - 1
end

-- gpt_init()
function gpt_init()
  x = -1
  while x < n do
    x = x + 1; y = -1
    while y < n do
      y = y + 1; set(x, y)
    end
  end
  reset(0, 0); reset(1, 0); reset(0, 1)
end

-- gpt_sieve()
function gpt_sieve()
  x = 0
  while x < rn do
    x = x + 1; y = -1
    while y < x do
      y = y + 1
      if cur(x, y) > 0 then
        px = x; py = y
        while not (px > n) do -- (※3)
          qx = px; qy = py; t = 1
          while t > 0 do  -- (※4)
            reset(qx, qy); reset(qy, qx)
            qx = qx - y; qy = qy + x
            if qx < 0 then t = 0 end
            if qy > n then t = 0 end
          end
          qx = px; qy = py; t = 1
          while t > 0 do
            reset(qx, qy); reset(qy, qx)
            qx = qx + y; qy = qy - x
            if qx > n then t = 0 end
            if qy < 0 then t = 0 end
          end
          px = px + x; py = py + y
        end
        set(x, y); set(y, x)
      end
    end
  end
end

-- gpt_show()
function gpt_show()
  x = -n; x = x - 1
  while x < n do
    x = x + 1; y = -n; y = y - 1
    while y < n do
      y = y + 1
      px = math.abs(x); py = math.abs(y) -- (※5)
      if cur(px, py) > 0 then
        qx = x; qx = qx * x
        qy = y; qy = qy * y; qx = qx + qy
        qy = n; qy = qy * n
        if not (qx > qy) then io.write("#")
        else io.write(" ")
        end
      else io.write(" ")
      end
    end
    io.write("\n")
  end
end

-- メイン
gptexture(39)

特に注意すべきことを挙げる。

  • ※2:gpt_setrn() では math.sqrt() を使っているが、TeX では平方根を自前で求める必要がある。ここでは効率は高くなくてもよい((明らかに gpt_sieve() の処理が所要時間の大部分を占める。))のでニュートン法等ではなく単純に rn をループさせて平方根(の切捨て)を求めている。
  • ※3: TeX には <= の比較が直接できないので、> の否定(not)として表す。ただし not を扱うのはそれほど簡単ではない。(後述)
  • ※4: 元々は while A and B do ... end というループであった。しかし TeX (の \@whilenum)では条件に単純な変数(や定数)同士の比較しか書けないので、何とかしてその形に持ち込む必要がある。ここでは整数変数 t をフラグとして用いて、
    t = 1
    while t > 0 do
      ...
      if not A then t = 0 end
      if not B then t = 0 end
    end
    というパターンに書き換えている。これは元と完全に同値でなくて後判定(repeat ... while A and B に相当する形)に変わっているが、実際は初回は必ずループに入ることが判っているので動作は同じになる。
  • ※5: 絶対値 math.abs(x)
    \ifnum\tcgp@x<0 -\fi\tcgp@x
    で著すことができる。なので敢えて「関数を手続きに変換する」対象にしなかった。

(続く)

*1:ガウス整数」の項目の中の一節。2012年6月8日更新の版。

*2:ガウス整数は「(整数)+(整数)·i」の形の複素数。そして「《自明な約数》を持たないガウス整数」がガウス素数である。