Flashエンジニアのくせに今日もFlashの話題を振らないnaoです。こんにちは。
　クビにならないよう努力はしています。

　今回のお題。何の役に立つのかと問われれば、漢らしく「何の役にもたたん」と答えましょう。しかし、ひょっとしたらいつか、使う日も来るかもしれない。こうして諸葛孔明の様に何事にも用意周到に備える私の部屋はモノで溢れかえっております。みろ、部屋がゴミのようだ。だってさー、PowerBook Duo(270cカラー液晶だぜ)とか、そのうちプレミアつきそうじゃん？いつか本当に必要な時が来たら、その時には声高らかに朗々と「こんな事もあろうかと」と宣言をするのである。

　というわけで今回は、まるで昔使った思い出の品のように、近年すっかり過去の遺物へと変貌し始めているUNIX TIMEです。UNIX TIMEといえば昔Perlのcgiで掲示板スクリプトを作るときなんかにはよく使われましたが、最近ではそこらの無料レンタルサービスでもDBが普通に使えるので、きょう日常用するのはtimestamp型ですよね。人間にも分かりやすいし、なにより遠い未来でも利用できる。それに対してUNIX TIMEは、2038年(2038年問題でググろう)までしか使えません。
　ではなぜ今更UNIX TIMEなのか。それはやっぱり、いつか使うかもしれないじゃん？という個人的なこだわりです。まぁ、使おうと思えばいくらでも使う場面はありますけどね。あえて使う理由も無いんですが。

　UNIX TIMEとは、1970年1月1日0時0分0秒(これをUNIX EPOCHという)より始まるタイムスタンプで、1秒経過すると値が1増えるという非常に扱いやすい性質を持ちます。そのため、例えば1970年1月1日0時1分0秒なら、UNIX TIMEは60というとうに、非常にシンプルです。この性質を利用すれば、例えば現在の時刻が判明すれば、UNIX EPOCHからの秒数を計算して、現在のUNIX TIMEを得ることが出来るのである。今回は、それを計算してみようというわけです。
　勘の良い人はそろそろ気がついたかもしれませんが、実は今回のラボブログ、目的はUNIX TIMEを計算することですが、本質的にはグレゴリオ暦を計算する計算式の話になります。暦が計算できて、UNIX EPOCHからの日数の差分が得られたら、あとはそれを秒数に変換すれば、UNIX TIMEになるのです。

グレゴリオ暦とは
グレゴリオ暦（グレゴリオれき）とは、1582年にローマ教皇グレゴリオ13世がユリウス暦を改良して制定した暦である。現行の太陽暦として世界各国で用いられている。単に新暦（英語：New Style、略称：N.S.、NS）と呼ばれる場合もある。現在使われている西暦はグレゴリオ暦である。
以上、wikipedia グレゴリオ暦より引用。

早い話が、今現在使われている暦の事である。500年近く前に作られたアルゴリズムなんですって。すごいよね。
それでは、さっそくUNIX TIMEの基準となる紀元元年1月1日からUNIX EPOCHまでの経過日数の計算をしてゆきましょう。

閏年の判定の仕方
閏年の判定は、条件式で書くとこんな感じになります。

year % 4 == 0 && year % 100 != 0 || year % 400 == 0

1)西暦を4で割って余りが0(4の倍数)の年は閏年
2)だけど100で割って余りが0(100の倍数)の年は除外して
3)なおかつ400で割った場合に余りが0(400の倍数)の年はやっぱ例外的に閏年

ということですね。

　これだけ理解しておけば、何かしらのプログラミング言語であれば、forとか構造体を使って暦の紀元元年1月1日からの経過日数を計算できますが、そこは公衆の面前に晒されるラボブログの記事。もう一歩突っ込んで、公式を使ってイッパツで計算しましょう。まぁ、誰にも見せないコードの中身だったら、面倒だから何も考えずforとかで回しちゃいますけどね。

グレゴリオ暦換算紀元0年1月1日からの経過日数の数え方
ここに「Fairfieldの公式(フェァフィールドのこうしき)」という非常に素晴らしい公式があります。こんなに素晴らしいのに、何故かマイナー。google先生にお伺いしても、あんまり情報がありませんが、今回はコレを使います。

Fairfieldの公式による経過日数の算出
365*（ｙ-1）+[ｙ/4]-[ｙ/100]+[ｙ/400]+31+28+1+[306*(ｍ+1)/10]-122+(d-1)

さぁ、ややこしいモノが出てまいりました。まずは、順を追って見ていきましょう。

365*（ｙ-1）

まずここは簡単。計算したい日の前年までの年数*1年の日数ですね。

+[ｙ/4]-[ｙ/100]+[ｙ/400]

次にこちら。ここは、閏年の処理ですね。[ ]の記号はガウス記号ってヤツで、int()と同じ小数点以下切り捨てですですが、これは高校とか大学で習うのかどうか知らないので、一応説明しておきました。自分中卒なので。 ここまでで、計算したい日の前年までの日数が出ます。

+31+28+1+[306*(ｍ+1)/10]-122+(d-1)

で、問題はこちら。ここが公式の決定的な肝であり、難関でもあります。 ざっくりと説明すると、この公式では1年間を3月〜14月として計算することで、2月の最終日を1年間の最後の日として扱うことが出来ます。 まず「+31+28+1」この部分が1月と2月(公式中の扱いとしては前年の13月と14月)の日数。 次に「[306*(ｍ+1)/10]」ここが、3月〜12月(公式中では1月〜10月)の日数。 最後に「-122+(d-1)」で、13月と14月の補正と、計算したい日の当月の日数になります。

ざっくりとした説明になりましたが、こんなものは数字大好きな人が理解していればいいことなので、私のような凡人は、さらにこの式を発展させて、以下のような式を使います。

365 * y + [y/4] - [y/100] + [y/400] + [306 * (m+1) / 10] + d - 428

「428」って何よ、とか言わないでください。元の式を計算すると、まぁなんとなくこんな感じになります。13月と14月の分と、1年の始まりが3月からになるようにしていた分を全部予め計算しておくと出る数字です。
数学が大の苦手の私的には、こっちがありがたい。

それでは、早速基準となるUNIX EPOCHである1970年1月1日0時0分0秒までの経過日数を計算しましょう。

365 * 1970 + [1970/4)]- [1970/100] + [1970/400] + [306 * (1 + 1) / 10] - 428 + 1 = 719161

　これで、UNIX EPOCHは、グレゴリオ暦換算紀元元年1月1日から719161日経過している、という事がわかりました。 　次は、この719161を起点として、そこから何日経過しているかを計算します。

_EPOC_TIME = 719161;
// Fairfieldの公式で算出した719161日からの経過秒数の算出
diff = ((365 * year + int(year / 4) - int(year / 100) + int(year / 400) + int(306 * (month + 1) / 10) - 428 + day) -  _EPOC_TIME) * 24 *60 * 60;

せっかくなのでFlash Lite 1.1のActionScriptで書いてみました。
これに、時刻を足してあげれば、UNIX TIMEの完成です。

_EPOC_TIME = 719161;
UNIX_TIME = ((365 * year + int(year / 4) - int(year / 100) + int(year / 400) + int(306 * (month + 1) / 10) - 428 + day) -  _EPOC_TIME) * 86400 + (hour * 3600) + (minute * 60) + second;

ほら1行でできた。エンジニアっぽいっ！

　今回はグレゴリオ暦の計算を用いてUNIX TIMEを導き出したわけですが、これを応用すれば、いろいろな事が計算できます。例えば、グレゴリオ暦換算の紀元元年1月1日は月曜日なので、Fairfieldの公式による経過日数を7で割った余りで曜日が判定できます。余り0なら日曜日。6なら土曜日。ちなみに、この割り算の"余り"を求める事をモジュロ演算と言います。演算子の"mod"や"%"のヤツ。サマーウォーズで見たよね？こいつのおかげで、ケンジは得意げになりナツキ先輩はメロメロなわけである。みんな覚えておこう。いつか「こんな事もあろうかと」という日が来るかもしれない。

　余談ですが、宇宙戦艦ヤマトの工場長兼技師長である真田志郎氏といえば「こんな事もあろうかと」で有名ですが、wikipediaによると劇中にそんなセリフは無いそうです。

えー、がっかりー。

日時: 20:41 | パーマリンク | コメント (1) | トラックバック (0)

最近、「何故、君はマウスを2つ同時に使っているんだい？」と聞かれることが多くなったbokkoです。VX Revolution RXは右手用ですが、左手で使うならMicrosoftのArcがオススメです。近頃はフットマウスを買うかどうか真剣に悩んでいます。

Software Designには去年にも「ソースを読み，パッチを作成してみよう～GNU GLOBAL，diff，patchの使い方～」という記事を執筆する機会を頂いたので、本誌に執筆するのはこれが二度目になります。

今回の内容は以前当ブログに書いた「diff with C++」の記事をもっと濃くした感じになっていて、編集距離やLCS、SESの解説に始まり、Subversionのdiffエンジンや拙作のdtlで使われているWuのO(NP)差分アルゴリズムについて解説しています。

WuのO(NP)アルゴリズム(以下Wuのアルゴリズム)や、それによく似たMyersのO(ND)アルゴリズムをはじめとするエディットグラフを使って差分を求めるアルゴリズムは普段何気なく使っているdiffコマンドや各種バージョン管理システムで使われていることからもわかるようにとても有用であり、アルゴリズム自体の美しさやスマートさにも目を見張るものがあります。興味がある方はぜひ手にとって読んでみてください。本記事を通してdiffのアルゴリズムのおもしろさが伝われば幸いです。

お詫びと訂正

本記事には訂正箇所が3箇所ほどあります。うち2つは本記事で使用しているdtlのバージョンが執筆後に0.05から0.06に上がったことによるもので、もう一つはWuのアルゴリズムの計算量に対する私の勘違いによるものです。詳しい内容についてはSoftware Design 2009年6月号：サポートページ｜gihyo.jp ... 技術評論社をご覧下さい。

関係者や読者の方にはご迷惑をおかけしました。すみません。以後、気を付けます。

おまけ

本記事にはサンプルコードとしてWuのアルゴリズムを使って2つの文字列間の編集距離を求めるC++のプログラムが載っているのですが、ふと思い立って、同じことをするプログラムおよびLCSとSESも計算するプログラムをLuaで書いてみました。だからどうというわけでもないのですが、何か新しいプログラミング言語を触る際に、こんな風に自分にとってなじみのあるアルゴリズムを記述してみるのはその言語を勉強する上で、とてもいい方法だと再確認した次第です。

実際に書いていてハマったのですが、Luaではほかの大多数の言語と違って配列や文字列のインデックスが0からではなく1から始まるので、別の言語で記述されたアルゴリズムとかのコードを写経する際は注意しましょう。

しかし、この間、プログラマの友達とそんなLuaの独特の挙動について話をしていると、

まあ、(アルゴリズムの)論文に載っているような擬似コードは配列のインデックスが1から始まってるから、論文の擬似コードを写経する分にはちょうどいいんじゃね？

という感じのことを言われて「ああ、言われてみればそうだなあ」と妙に納得した次第です。とりあえず、今後、論文に載ってる疑似コードを実際のコードに落とす際はまずLuaで書いてそれからCとかC++に書き直すということを実践してみようと思います。

話が逸れました。以下にLuaによるWuのアルゴリズムの実装を2種類示します。その1が編集距離のみを計算するバージョン、その2が編集距離に加えてLCS、SESを計算するバージョンになります。なお、dtlと違ってワーストケース(LCSが極端に短くなる場合)への対応は行っていません。

LuaによるWuのアルゴリズムの実装 その1(編集距離のみ)

-- editDistance.lua
-- Lua5.1.4で動作確認
ONP = {}
function ONP.new (a, b)           -- ONPクラスのコンストラクタ
   local self = {                           -- メンバ変数
      A = a,
      B = b,
      M = string.len(a),
      N = string.len(b),
   }
   function self.editDistance ()  -- 編集距離を計算する
      offset = self.M + 1
      delta  = self.N - self.M
      size   = self.M + self.N + 3
      fp = {}
      for i = 0, size-1 do
	 fp[i] = -1
      end
      p = -1
      repeat
	 p = p + 1
	 for k=-p, delta-1, 1 do
	    fp[k+offset] = self.snake(k, math.max(fp[k-1+offset]+1, fp[k+1+offset]))
	 end
	 for k=delta+p,delta+1, -1 do
	    fp[k+offset] = self.snake(k, math.max(fp[k-1+offset]+1, fp[k+1+offset]))
	 end
	 fp[delta+offset] = self.snake(delta, math.max(fp[delta-1+offset]+1, fp[delta+1+offset]))
      until fp[delta+offset] >= self.N
      return delta + 2 * p
   end
   function self.snake (k, y)     -- 最遠点のy座標を計算する
      x = y - k
      while (x < self.M and y < self.N and 
	     string.sub(self.A, x+1, x+1) == string.sub(self.B, y+1, y+1)) 
      do
	 x = x + 1
	 y = y + 1
      end
      return y
   end
   if self.M >= self.N then       -- N >= Mになるように調整
      self.A, self.B = self.B, self.A
      self.M, self.N = self.N, self.M
   end
   return self
end
if #arg < 2 then
   error("few argument")
end
a = arg[1]
b = arg[2]
d = ONP.new(a, b)
print("editDistance:" .. d:editDistance())

実行

narazuya@bokkko% lua editDistance.lua abcdef dacfea
6
narazuya@bokkko% lua e.editDistancelua abc abd
2
narazuya@bokkko% 

LuaによるWuのアルゴリズムの実装 その2(編集距離、LCS、SES)

-- onp.lua
-- Lua5.1.4で動作確認
ONP = {}
SES_DELETE = -1
SES_COMMON = 0
SES_ADD    = 1
function ONP.new (a, b)          -- ONPクラスのコンストラクタ
   local self = {                          -- メンバ変数
      A = a,
      B = b,
      M = string.len(a),
      N = string.len(b),
      path       = {},
      pathposi   = {},
      P          = {},
      ses        = {},
      seselem    = {},
      lcs        = "",
      editdis    = 0,
      reverse    = false,
   }
   -- getter
   function self.geteditdistance () 
      return self.editdis
   end
   function self.getlcs ()
      return self.lcs
   end
   function self.getses ()
      return self.ses
   end
   -- constructor
   function self.P.new (x_, y_, k_)
      local self = { x=x_, y=y_, k=k_ }
      return self
   end
   function self.seselem.new (elem_, type_)
      local self = { elem=elem_, type=type_}
      return self
   end
   -- 差分構築
   function self.compose ()
      offset = self.M + 1
      delta  = self.N - self.M
      size   = self.M + self.N + 3
      fp = {}
      for i = 0, size-1 do
	 fp[i]        = -1
	 self.path[i] = -1
      end
      p = -1
      repeat
	 p = p + 1
	 for k=-p, delta-1, 1 do
	    fp[k+offset] = self.snake(k, fp[k-1+offset]+1, fp[k+1+offset])
	 end
	 for k=delta+p,delta+1, -1 do
	    fp[k+offset] = self.snake(k, fp[k-1+offset]+1, fp[k+1+offset])
	 end
	 fp[delta+offset] = self.snake(delta, fp[delta-1+offset]+1, fp[delta+1+offset])
      until fp[delta+offset] >= self.N
      self.editdis = delta + 2 * p
      r    = self.path[delta+offset]
      epc  = {}
      while r ~= -1 do
	 epc[#epc+1] = self.P.new(self.pathposi[r+1].x, self.pathposi[r+1].y, nil)
	 r = self.pathposi[r+1].k
      end
      self.recordseq(epc)
   end
   function self.snake (k, p, pp)     -- 最遠点のy座標を計算する
      r = 0;
      if p > pp then
	 r = self.path[k-1+offset];
      else
	 r = self.path[k+1+offset];
      end
      y = math.max(p, pp);
      x = y - k
      while (x < self.M and y < self.N and 
	     string.sub(self.A, x+1, x+1) == string.sub(self.B, y+1, y+1)) 
      do
	 x = x + 1
	 y = y + 1
      end
      self.path[k+offset] = #self.pathposi
      p = self.P.new(x, y, r)
      self.pathposi[#self.pathposi+1] = p
      return y
   end
   function self.recordseq (epc)          -- LCS、SESを記録する
      x_idx,  y_idx  = 1, 1
      px_idx, py_idx = 0, 0
      for i=#epc, 1, -1 do
	 while (px_idx < epc[i].x or py_idx < epc[i].y) do
	    if (epc[i].y - epc[i].x) > (py_idx - px_idx) then
	       elem = string.sub(self.B, y_idx, y_idx)
	       if self.reverse then 
		  type = SES_DELETE
	       else
		  type = SES_ADD
	       end
	       self.ses[#self.ses+1] = self.seselem.new(elem, type)
	       y_idx  = y_idx  + 1
	       py_idx = py_idx + 1
	    elseif epc[i].y - epc[i].x < py_idx - px_idx then
	       elem = string.sub(self.A, x_idx, x_idx)
	       if self.reverse then 
		  type = SES_ADD
	       else
		  type = SES_DELETE
	       end
	       self.ses[#self.ses+1] = self.seselem.new(elem, type)
	       x_idx  = x_idx  + 1
	       px_idx = px_idx + 1
	    else 
	       elem = string.sub(self.A, x_idx, x_idx)
	       type = SES_COMMON
	       self.lcs = self.lcs .. elem
	       self.ses[#self.ses+1] = self.seselem.new(elem, type)
	       x_idx  = x_idx  + 1
	       y_idx  = y_idx  + 1
	       px_idx = px_idx + 1
	       py_idx = py_idx + 1
	    end
	 end
      end
   end
   if self.M >= self.N then       -- N >= Mになるように調整
      self.A, self.B = self.B, self.A
      self.M, self.N = self.N, self.M
      self.reverse = true
   end
   return self
end
if #arg < 2 then
   error("few argument")
end
a = arg[1]
b = arg[2]
d = ONP.new(a, b)
d:compose()
print("editDistance:" .. d:geteditdistance()) -- 編集距離
print("LCS:"          .. d:getlcs())                    -- Longest Common Subsequence
print("SES")
ses = d:getses()                                          -- Shortest Edit Script
for i=1, #ses do
   if ses[i].type == SES_COMMON then
      print("  " .. ses[i].elem)
   elseif ses[i].type == SES_DELETE then
      print("- " .. ses[i].elem)
   elseif ses[i].type == SES_ADD then
      print("+ " .. ses[i].elem)
   end
end

実行

narazuya@bokkko% lua onp.lua abcdef dacfea
editDistance:6
LCS:acf
SES
+ d
  a
- b
  c
- d
- e
   f
+ e
+ a
narazuya@bokkko% lua a.lua abc abd
editDistance:2
LCS:ab
SES
  a
  b
- c
+ d
narazuya@bokkko%

日時: 10:18 | パーマリンク | コメント (0) | トラックバック (0)

ミートソーススパゲティを作るときは、ミートソースから作るのが信条のbokkoです。それはさておき、今日はdiffのお話です。

diff

diffは指定した2つのファイルの差分を求めるコマンド、もしくはその差分そのものを指します。普段から何気なく使用しているコマンドですが、その中で使われているアルゴリズムは結構難しいです。

差分を計算するということ

差分を計算するというのは以下の3つを求めることに帰結します。

・Levenshtein Distance(Edit Distance)
・LCS(Longest Common Subsequence)
・SES(Shortest Edit Script)

上から順に1つずつ説明していきます。

続きを読む "diff with C++" »

日時: 11:43 | パーマリンク | コメント (3) | トラックバック (0)