2012/12/01 〜 2012/12/25TeX & LaTeX Advent Calendar
こちらは(ry
* * *
「変換して TeX プログラミングの例」シリーズの最終回。今回は「ガウス素数の模様」の図を TeX で描いてみる。
ガウス素数の模様って何?
Wikipedia*1を参照。右にある図がその模様である。ということで説明終わり。
……でよいと思うがもう少し説明してみる。素数というのは「自明でない約数を持たない自然数」のことで、列挙すると次のようになるのであった。
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, ……
いま、1 から順に自然数を数えて、それが素数であれば黒、なければ白の四角を並べると次のようになる。
□■■□■□■□□□■□■□□□■□■□□□■□□□□□■□■□□□□□■□□□……
正の実数に対して定義されていた「整数」「素数」を拡張して、複素数の世界について《整数》や《素数》という概念を定義することができ、それぞれ「ガウス整数」「ガウス素数」と呼ばれる。*2その《素数》から上の図に相当するものを作ると、二次元の上下左右に無限に続く「模様」になる。それが「ガウス素数の模様」である。
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
で著すことができる。なので敢えて「関数を手続きに変換する」対象にしなかった。
(続く)