krxrossの雑記帳

2017-10-10

Haskellでラマヌジャン数

見た目がきれいなバージョンがあります。
https://qiita.com/zeal1483/private/69689dd8a68903932170

第4回『Haskellによる関数プログラミングの思考法』読書会に参加してきました。
https://sampou.connpass.com/event/66242/
そこで、ラマヌジャン数をhaskellで求める方法として以下が上がりました。
http://takeichi.ipl-lab.org/~onoue/pub/jssst01.pdf

10 -- Finding Ramanujan numbers
11
12 ramanujan :: [((Int,Int), (Int,Int))]
13 ramanujan = [(x,y)|(x,y)<-zip s (tail s),
14     sumcubes x == sumcubes y]
15   where s = fsort sumcubes 1
16
17 sumcubes :: (Int,Int) -> Int
18 sumcubes (a,b) = a^3 + b^3
20
21 fsort :: ((Int,Int) -> Int) -> Int
22 -> [(Int,Int)]
23 fsort r k
24 = (k,k) : fmerge r [(k,b)|b<-[k+1..]]
25 (fsort r (k+1))
26
27 {-
28 Two argument lists should be infinite,
29 so there is no definition for null list.
30 -}
31
32 fmerge :: (a->Int) -> [a] -> [a] -> [a]
33 fmerge r (x:xs) (y:ys)
34   | r x <= r y = x : fmerge r xs (y:ys)
35   | otherwise = y : fmerge r (x:xs) ys

しかし、下の様なナイーブな方法よりは、早いですが
1000個のラマヌジャン数を求めるなどすると、遅いことがわかります。

rmNum n = [(e,a,b,c,d)|e<-[28..2*n^3]
                      ,a<-[1..n]
                      ,b<-[a+1..n]
                      ,e == a^3 + b^3
                      ,c<-[a+1..n]
                      ,d<-[c..n]
                      ,e == c^3 + d^3]

遅い理由は、下記のフィボナッチ数列の求め方が遅い理由と同じように思えます。

fib n = fibs !! n
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)

Haskell神話
http://d.hatena.ne.jp/kazu-yamamoto/20100624/1277348961
フィボナッチ数列の遅い理由は、上記リンクによると
計算(サンク?、算術木?)が消せないこととあります。

(この遅い理由が、計算が消せないとなっていることについては
自分には、なぜそうなのかが、明瞭には分からないです。
空間計算量や、時間計算量ならば、直感的に分かります。

例えば、メモリの専有領域の問題なら、限界を超えると
スワップが発生し、ディスクIOが保存と取込で2重に発生するから
速いメモリとのやり取りから、遅いディスクとのやり取りが加わり
そのスピード差から、圧倒的に速度が違うことは、明らかです。

また、時間計算量の話は、計算に物理的に時間がかかるので
これも、明らか。
最初のナイーブな書き方で遅い理由だと思われる。

計算が消せないと、なぜ遅くなるのかのメカニズムが見えない。)

これとは別に問題点として
a^2+b^3=c^3+d^3=e^3+f^3=kとなるような
要件を満たす(a,b),(c,d),(e,f)とkが存在することがあります。
例えば、k=87539319,(167,436),(228,423),(255,414)
連続する2つの組で捉える方法では、これを捉えられません。
(これは、groupByやgroupを使用することで、簡単に対応できます。)

読書会の時に、山下さんが、言われた
x+y=kをk=[1..]として進めていく方法の方が
実装してみると明らかに速く
1000個のラマヌジャン数を求めても
苦にはならないくらい十分にスピードがあります。

ここで、実装に入る前に、実装方針の説明をしてみます。

単に、x+y=kをk=[1..]で進めていくと、大小関係がわからない。
(読書会中に、議論になっていたとおもう)

そこで、区間を(1,2]、(2,4]、(4,8]、(8,16]、(16,32]...と区切っていく
その内の(n,2n]と(2n,4n]を考える。
そうすると、下図のように大小関係がはっきりする。

f:id:krxross:20171010133051j:image

図の(n,2n]上のラマヌジャン数を考えると
点線x^3+y^3=n〜x^3+y^3=2nの範囲を考えれば良い

f:id:krxross:20171010133046j:image

A:∫(x,y)dk, x+y=k、ただし積分範囲( n,2n]
B:∫(x,y)dk, x+y=k、ただし積分範囲(2n,4n]
f(x,y)=x^3+y^3
とおく

{(x,y)|(x,y)∈A∪B} <= MAX( {f(x,y)|(x,y)∈A} )
となる(x,y)を取れば、図の△鮠絽造箸垢詢琉茲とれる。

下限のい砲弔い討蓮Aの1つ下の階層O(オー)を考えればよい。

O:∫(x,y)dk, x+y=k、ただし積分範囲(n/2,n] (nが偶数の場合)
O:x+y=0、x^3+y^3=0として、下限値0をとる (nが奇数、すなわち1の場合)

これに対して、下記のように取れば、下限を押さえられる
{(x,y)|(x,y)∈A∪B} > MAX( {f(x,y)|(x,y)∈O} )


以上から、上図の斜線部を1単位として

n = [2^k| k<-[1..]]

を展開して、足し合わせることによって
ラマヌジャン数の探索範囲を敷き詰めることができる。

実装は、以下の通り。
makeSearchAreaで、上図の斜線部を求めている。

import Data.List

-----------------------------------------
-- x + y = n 上の正数値列挙
-- > line 8
-- > [(1,7),(2,6),(3,5),(4,4)]
-----------------------------------------
line n = [(a,n-a)|a<-[1..n-1],a<=n-a]

------------------------------------------------------------------
-- 高さ(a^3 + b^3)付与
-- > estLine 8
-- > [(344,1,7),(224,2,6),(152,3,5),(128,4,4)]
------------------------------------------------------------------
estLine n = map cubeNumNum (line n)
  where cubeNumNum (a,b) = (sumCubes (a,b), a, b)

------------------------------------------------------------------
-- 2n〜4nの区間の等高線を束ねて面にする
-- (これを等高線の層と言う事にする)
-- > sekibun (2, 4)
-- > [(9,1,2),(16,2,2),(28,1,3)]
-- s: 2^n (2の指数乗)
-- e: 2*s
------------------------------------------------------------------
sekibun (s, e) = sort (concatMap estLine [s+1..e])

------------------------------------------------------------------
-- 探索領域をつくる
-- > makeSearchArea 1
-- > [(9,1,2),(16,2,2),(28,1,3)]
-- > makeSearchArea 2
-- > [(35,2,3),(54,3,3),(65,1,4),(72,2,4),(91,3,4),(126,1,5),(128,4,4),(133,2,5),(152,3,5),(189,4,5),(217,1,6),(224,2,6),(243,3,6),(250,5,5),(280,4,6),(341,5,6),(344,1,7)]
------------------------------------------------------------------
makeSearchArea n = filter (\(x,y,z) -> min < x && x <= max) concat2Part
  where -- 上層・中層・下層
        fstPart = sekibun (n,   2*n)
        mdlPart = sekibun (2*n, 4*n)
        lstPart = sekibun (4*n, 8*n)
        -- 上界・下界
        min = fst3 (last fstPart)
        max = fst3 (last mdlPart)
        -- 連続した上位2つの等高線の層を束ね、整列する
        concat2Part = sort (mdlPart ++ lstPart)

------------------------------------------------------------------
-- ラマニュジャン数(等高線の層の指定)
-- > rmNumRange 2
-- > []
-- > rmNumRange 4
-- > [[(1729,1,12),(1729,9,10)]]
-- > rmNumRange 8
-- > [[(4104,2,16),(4104,9,15)],[(13832,2,24),(13832,18,20)],[(20683,10,27),(20683,19,24)]]
------------------------------------------------------------------
rmNumRange n = getRmNums (makeSearchArea n)

------------------------------------------------------------------
-- searchArea上のラマヌジャン数取得
------------------------------------------------------------------
getRmNums searchArea = filter isDuplicate (groupBy eq searchArea)
  where -- 等高線上の点か ?
        eq p q = fst3 p == fst3 q
        -- 複数要素があるか ?
        isDuplicate xs = length xs > 1

------------------------------------------------------------------
-- ラマニュジャン数
-- > take 1 rmNums
-- > [[(1729,1,12),(1729,9,10)]]
-- > take 2 rmNums
-- > [[(1729,1,12),(1729,9,10)],[(4104,2,16),(4104,9,15)]]
-- > take 3 rmNums
-- > [[(1729,1,12),(1729,9,10)],[(4104,2,16),(4104,9,15)],[(13832,2,24),(13832,18,20)]]
------------------------------------------------------------------
rmNums = concatMap rmNumRange ninoshisu



------------------------------------------------------------------
-- 補助関数
-- 2の指数の系列を生成
-- > ninoshisu  8
-- > [2,4,8,16,32,64,128,256...]
------------------------------------------------------------------
ninoshisu = map (\m -> 2^m) [1..]

------------------------------------------------------------------
-- 補助関数
------------------------------------------------------------------
fst3 (x,y,z) = x
sumCubes (a,b) = a^3 + b^3

上層・中層・下層の3つを、それぞれで求めているので
無駄に見えるかもしれない。

  where -- 上層・中層・下層
        fstPart = sekibun (n,   2*n)
        mdlPart = sekibun (2*n, 4*n)
        lstPart = sekibun (4*n, 8*n)

iterateを使えば、上層を、それぞれで求めることになるが
体感的には、速度は上がらないようです。

import Data.List

-----------------------------------------
-- x + y = n 上の正数値列挙
-- > line 8
-- > [(1,7),(2,6),(3,5),(4,4)]
-----------------------------------------
line n = [(a,n-a)|a<-[1..n-1],a<=n-a]

------------------------------------------------------------------
-- 高さ(a^3 + b^3)付与
-- > estLine 8
-- > [(344,1,7),(224,2,6),(152,3,5),(128,4,4)]
------------------------------------------------------------------
estLine n = map cubeNumNum (line n)
  where cubeNumNum (a,b) = (sumCubes (a,b), a, b)

------------------------------------------------------------------
-- 等高線を束ねて、面にする(これを等高線の層と言う事にする)
-- > sekibun (2, 4)
-- > [(9,1,2),(16,2,2),(28,1,3)]
------------------------------------------------------------------
sekibun (s, e) = sort (concatMap estLine [s+1..e])

------------------------------------------------------------------
-- 次のラマヌジャン数取得
-- > get6th $ next (1, 0, 0, [], [], [])
-- > []
-- > get6th $ next $ next (1, 0, 0, [], [], [])
-- > []
-- > get6th $ next $ next $ next (1, 0, 0, [], [], [])
-- > [[(1729,1,12),(1729,9,10)]]
-- > get6th $ next $ next $ next $ next (1, 0, 0, [], [], [])
-- > [[(4104,2,16),(4104,9,15)],[(13832,2,24),(13832,18,20)],[(20683,10,27),(20683,19,24)]]
------------------------------------------------------------------
next (n, min, max, fstPnl, mdlPnl, _) = (2*n, max, nextMax, mdlPnl, lstPnl, rmNums)
  where lstPnl = sekibun (4*n, 8*n)
        -- 次の上界
        nextMax = fst3 (last lstPnl)
        -- x^3 + y^3の上下の、上界の等高線を考慮して検索エリア作成
        searchArea = makeSearchArea (min, max) mdlPnl lstPnl
        -- searchArea上のラマヌジャン数
        rmNums = getRmNums searchArea

------------------------------------------------------------------
-- 補助関数
-- 上位2層を結合して、区間で絞り込む
------------------------------------------------------------------
makeSearchArea (min, max) mdlPnl lstPnl = filter p (sort (mdlPnl ++ lstPnl))
  where p (x,y,z) = min < x && x <= max

------------------------------------------------------------------
-- 補助関数
-- searchArea上のラマヌジャン数取得
------------------------------------------------------------------
getRmNums searchArea = filter isDuplicate (groupBy eq searchArea)
  where -- 等高線上の点か ?
        eq p q = fst3 p == fst3 q
        -- 複数要素があるか ?
        isDuplicate xs = length xs > 1

------------------------------------------------------------------
-- ラマニュジャン数
-- > take 1 rmNums
-- > [[(1729,1,12),(1729,9,10)]]
-- > take 2 rmNums
-- > [[(1729,1,12),(1729,9,10)],[(4104,2,16),(4104,9,15)]]
-- > take 3 rmNums
-- > [[(1729,1,12),(1729,9,10)],[(4104,2,16),(4104,9,15)],[(13832,2,24),(13832,18,20)]]
------------------------------------------------------------------
rmNums = concatMap get6th $ drop 3 $ iterate next (1, 0, 0, [], [], [])
get6th (n, min, max, fstPnl, mdlPnl, rmNums) = rmNums

------------------------------------------------------------------
-- 補助関数
------------------------------------------------------------------
fst3 (x,y,z) = x
sumCubes (a,b) = a^3 + b^3

2017-05-20

『ゼロから作るDeep Learning』読書会@高円寺(4)その3

(7shiさんに聞きたいこと)

下記のリンクで、あの時間に見えてきたことは
http://playground.tensorflow.org/


内容的には下リンク内の下記文言のことと
同じと思うけど、あってますか ?

https://rekken.g.hatena.ne.jp/murawaki/20161017/p1

>>引用開始
システムの中身を直感的に説明するのは難しい。LeCun 御大の例えをもじって、ノブを使った説明を試みる。機械翻訳というブラックボックスには、上部と下部に大量の穴があいていて、それぞれ入力と出力に対応する。上部の適切な穴に水を注ぐと、下部の適切な穴から水が出てくる。上部と下部の穴の間には複雑にからまったパイプがあり、途中で分岐 (というより分身) したり合流したりするし、途中に水を貯めている箇所があったりもする。そういう箇所にはノブがあって、水の流れを制御する。実際には、量が増幅されたり、負の値をとったりするので、水で例えるのは微妙だけど。
(中略)
学習とはノブを調整すること。ノブを適切に調整していないと、別の穴からちょろちょろと水が漏れたりする。正しい穴だけから水が出るようにノブを調整する。こういった階層の深いシステムであっても、充分な教師データを与えれば適切にノブを調整できることがわかった。それが深層学習とよばれているもの。とはいえ、途中を流れている水を見ても、何が起きているのか人間には (システム設計者ですら) さっぱりわからない。
そろそろノブの例えが厳しくなってきたのでここでやめにするが、最後に一つ付け加える。システム内部は、水の流れが多いか少ないかという数量で制御されている。確かに入り口と出口は離散的だけど、中は連続値で支配されている。記号の出る幕はない。
引用終了<<


この例えの文の内容以上のことは
具体に見えるという点以外は
自分には、同じに見えてしまいます。

playground.tensorflow.orgのリンクは
面白いと思います。
ですので、価値云々の話ではなく

自分が話が見えていないのではないかということに
関して、なにが違うかを聞きたいのです。


話は変わります。
上記の機械学習の水での例えの文で
自分が気になる点を抜き出すと
>途中を流れている水を見ても、何が起きているのか人間には
>(システム設計者ですら) さっぱりわからない。
ここの部分が、機械学習の重要な性質というか
限界に思えます。


昔、ninetiesさんの数学講座のあとで
7shiさん、「ハンドル名を覚えていない」さん、alumicanさん、自分
の懇親会で、ninetiesさんが
学習した内容は、誰にもわからない。
引き継いで応用していくことはできないから

人間世界のように、単純に理論を発展展開さ
せることはできない ?

というようなニュアンスの話を聞いた記憶があります。

自分のフィーリング解釈では、ざっくり
機械学習は、人間世界でいうところの
暗黙知しか獲得できない。」
ということなんじゃないかと思いました。

人間なら、暗黙知形式知に変換して
抽象化し、例えば数式などにしてから
書物の形にして、引継ぎ、発展していけます。

ニュートンの「巨人の肩の上に乗る」
といわれるやつです。

更に妄想すると、定理の自動証明なども
暗黙知だけでは、どうしようもできなくて
なにか、ブレークスルーが必要なんじゃないかと?

機械翻訳についても、文脈理解が必要なので
(機械学習で、画像と同じパターン認識
文法の様なものも、見えてくるようですが)
文法をとばして、文脈まで見えるように
ならないんじゃないかと思います。

googleは、少なくとも、日本語文法は無視する
宣言していて、アルゴリズムの力と火力で
単語列のパターンを対応付けて
google翻訳している模様。

ブレークスルーを見つければ
彼らを、出し抜けるかも !?

2017-05-17

『ゼロから作るDeep Learning』読書会@高円寺(4)その2

⇒発言しなかったが、思ったこと
ニューロンというと、言葉の意味が文脈でかなり
変わるので、掲題のような説明に使われるのは
なんとなく嫌な感じがした。

神経学者や、認知科学者のような
実用観点ではなく、人間の現実的モデルに関心を
持っている層は、ニューラルネットワークと(脳細胞という
意味での)ニューロンとは、マッタク違うというのが
共通認識(常識的)になっていると、何かで見たことがあります。
ニューラルネットワークでは、どうも人間の神経系のモデル
としては、役にたたないとのこと。
昔は、同じと思って研究していたことは確かなようですが。

直感的に、人の思考と機会学習では動きが違うし
機械学習は、学習効率がかなり悪いと思います。
機械翻訳などについて言えば、最低100万コーパス位の
教師データが無いと、まともな結果が得られないとありますが
中学・高校の英語で習うテキストで100万コーパスは無いです。
学生用の英語辞書でさえ、100万コーパス無いのではないでしょうか?
これで、一昔前の学生は、受験英語ラッセルなどの難解な文章
を読まされ、テストに回答させられたんです。
人間の方が、相当凄くないですか?

また、囲碁の世界チャンピオンに勝った話ありますが
(裏付けのある確かな情報ではないですが
たぶん電気代だけで)2億円(ドルだったかな?)を溶かして
半年とか1年とかをかけてgoogleの火力のすごさで
やっと到達したと聞きました。

ニューロンニューラルネットワークは別ものという話は
直感的に十分納得できる話だと思います。
実用的な人工知能を作りたいという機械学習系の話位でしか
ニューラルネットワークが登場しないので
文脈依存にならなく、一意に意味がとれるので
ニューロンという言葉を重要な説明の箇所に
使ってほしくないと個人的に思いました。

『ゼロから作るDeep Learning』読書会@高円寺(4)

『ゼロから作るDeep Learning読書会高円寺(4)
https://koenjidl.connpass.com/event/56139/
に参加してきた。

参加者で、気になった発言を列記してみる。

以下、概ね時系列で列挙する
(1)なぜ、線形領域しか扱えない線形分離器を重ね合わせると
  非線形領域を扱えるようになるのか ?
 ⇒自分の不規則発言で、答えのつもりで(2)を発言
(2)行列で線形分離をする決定面(y = ax + b)を変換かけるのだから
  非線形にできるのは不思議ではない。

  (「例はたくさんあるから不思議におもわなくても
   良いんじゃなかな?」という話
   他のジャンルでも変換で複雑な式に変換する例は
   あると思う。

   例えば
    線形ではないが、地動説流の惑星の楕円起動から
    天動説流の時間tを含むtの6次式だったか
    8次式だったか定かでないが
    複雑な方程式に座標変換できたと思う。
    その変換がガリレイ変換だったかローレンツ
    変換だったか覚えてないです。

    天動説が間違っているかについては
    相対論とのからみで下リンクに説明がある
    http://fj.sci.physics.narkive.com/tR5DC096/4-009
    要するに、どの慣性系をえらんでもよい。
    だけど、太陽を中心に計算する座標系を選んだ方が
    「計算が楽だから、みんな当然そうするよね。」
    という話。
    天動説 vs 地動説の話は、相対論のある現代に
    おいては死んだ話題。
   
   細かく見ると、上記変換には
   時間tがパラメータに含まれているので
   ただ線形和をとってる訳ではない
   
   線形和だけでは、非線形には
   ならないと思われます。
   例えば、ABCD...を行列として
   行列の積の結果Eを
   E = ABCD...
   の様に考える
   Eとベクトル(x,y)の積は
   Eの要素がE11,E12,E21,E22で
   あると考えると
   (E11*x + E12*y,E21*x + E22*y)になり
   変換後も、線形であることには、
   変わりがないことになる。

   線形和のみではなく、活性化関数hが含まれるため
   この状況が変わると思われます。
   h関数は、x、yをとるので、h(x,y)と考えても
   良いと思います。
   また、この話の流れとは直接関係ないですが
   重みの行列をWとすると
   hとW(x y)の積でhとWは交換可能でない
   つまり、hW-Wh≠0(hW≠Wh)です。

   話は変わりますが、もっと俯瞰した話として
   ニューラルネットワークには、任意の関数を表現できる話もあります。
   http://yuki-sato.com/wordpress/2014/09/03/performingfunction/

   もっと、正確な議論は下リンク内にある論文など
   見ることになるんですかね?
   http://qiita.com/HirofumiYashima/items/774f8b41489e2622e1db

(3)パーセプトロンニューラルネットワーク
  言っているが、なぜパーセプトロン
  対してニューラルネットワークでなく
  ニューロンを使わないのはなぜ ?

  ⇒回答として結論はでていて
    この本では、ニューロンは、グラフのノードと同じ意味
    で使用しているので、その流儀では、ニューラルネットワーク
    が正しい。
  ⇒発言しなかったが、思ったこと
    ニューロンという言葉は、意味が文脈でかなり
    変わるので、掲題のような説明に使われるのは
    なんとなく嫌な感じがした。
    ⇒いなや理由
    http://d.hatena.ne.jp/krxross/20170517/1495032611

(4)ローカルミニマムに落ちない保証は ?
  ⇒初っ端に思ったこと、上手く説明できないので
   発言しなかったが、端的に言うと
   「あるワケ無い」とフィーリングで思った。
   任意性の高い関数超曲面の最小値をどうやって
   求められるだろうか ?

   話が違いますが、数理最適化問題
   巡回セールス問題を遺伝的アルゴリズム
   探す手法があるが、最適解でないばかりか
   最悪どこまで、最適解とかけ離れているかも
   保証できなかったとおもいます。
   上記の問題よりも、離散値でなく
   連続値であることや、線形でなく
   任意のグラフが書けてしまうことを
   考えると、問題はさらに難しいと
   フィーリングでは思ってしまいます。
   
   (質問者は、数式で論理的に分かりたいのだと
   思うので、質問者にとって上記はただの
   ポエムでしかなく、なんの価値もないですが。)

   ところで、問題が解けると言う時に
   嬉しい順は、次の順のように思う。
   (a)代数的に解ける。
   (b)近似で解け、最適解に限りなく
     近づくことができる。
   (c)近似で解け、最適解から
     かけ離れている率が、高々決まる。
   (d)近似で解け、最適解から
     どれだけ、かけ離れているのか
     推定すらできない。

   アルゴリズムの性能上(d)を選択するしか
   ない問題が、最適化問題には
   多くあると思います。
   それと、ニューラルネットワーク
   ローカルミニマム問題は
   変わらないと思いますがどうでしょうか?
   (所詮近似でしかないですし)






■自分が7shiさんに聞きたいこと。
http://d.hatena.ne.jp/krxross/20170520/1495250231

2017-02-05

SQLで組合せ最適化(メモ2)

あまり、整理できていないので
あとで整理する予定だが、列記する

以下SQLServerで確認



素因数分解
nullやゴミを取り除いていない


WITH Input(iStart, iEnd) AS (
  SELECT 1, 10
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber INNER JOIN Input
    ON NaturalNumber.num + 1 <= Input.iEnd
),
prime_factorization (num, value, prime, success) AS (
   SELECT CAST(num AS INTEGER)
         ,CAST(num AS INTEGER)
         ,CAST(2 AS INTEGER)
         ,NULL
   FROM NaturalNumber
   UNION ALL
   SELECT
     num,
     CASE WHEN (value % prime) = 0
       THEN value/prime ELSE value
     END,
     CASE WHEN (value % prime) = 0
       THEN prime ELSE prime + 1
     END,
     CASE WHEN (value % prime) = 0
       THEN 1 ELSE 0
     END
   FROM prime_factorization
   WHERE value > 1 OR (value = prime)
)
SELECT
  num,
   (SELECT ',' + CAST(prime AS VARCHAR(MAX)) 
    FROM prime_factorization PF
    WHERE success = 1
         AND PF.num = NN.num
    FOR XML PATH('')
  ) AS primes
FROM NaturalNumber NN


ごみを取り除くと


WITH Input(iStart, iEnd) AS (
  SELECT 1, 50
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber INNER JOIN Input
    ON NaturalNumber.num + 1 <= Input.iEnd
),
primes (num, value, prime, success) AS (
   SELECT CAST(num AS INTEGER)
         ,CAST(num AS INTEGER)
         ,CAST(2 AS INTEGER)
         ,NULL
   FROM NaturalNumber
   UNION ALL
   SELECT
     num,
     CASE WHEN (value % prime) = 0
       THEN value/prime ELSE value
     END,
     CASE WHEN (value % prime) = 0
       THEN prime ELSE prime + 1
     END,
     CASE WHEN (value % prime) = 0
       THEN 1 ELSE 0
     END
   FROM primes
   WHERE value > 1 OR (value = prime)
),
primes2 AS (
SELECT
  num,
   (SELECT ',' + CAST(prime AS VARCHAR(MAX)) 
    FROM primes PF
    WHERE success = 1
         AND PF.num = NN.num
    FOR XML PATH('')
  ) AS primes,
   (SELECT count(prime) 
    FROM primes PF
    WHERE success = 1
          AND PF.num = NN.num
  ) AS cnt
FROM NaturalNumber NN
)
SELECT num, SUBSTRING(primes, 2, LEN(primes))
FROM primes2
WHERE cnt > 0



完全数


WITH Input(iStart, iEnd) AS (
  SELECT 1, 100
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber
  INNER JOIN Input
    ON NaturalNumber.num + 1 <= Input.iEnd
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND (A.num % B.num) = 0
)
SELECT NN.num "完全数",
   (SELECT '+' + CAST(PF.prime AS VARCHAR(MAX)) 
    FROM Primes PF
    WHERE PF.num = NN.num
    FOR XML PATH('')
  ) AS "計算式"
FROM Primes NN
GROUP BY num
HAVING SUM(prime) = num



不足数
完全数と条件の符号が異なるだけである


WITH Input(iStart, iEnd) AS (
  SELECT 1, 100
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber
  INNER JOIN Input
    ON NaturalNumber.num + 1 <= Input.iEnd
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND (A.num % B.num) = 0
)
SELECT NN.num "不足数",
   (SELECT '+' + CAST(PF.prime AS VARCHAR(MAX)) 
    FROM Primes PF
    WHERE PF.num = NN.num
    FOR XML PATH('')
  ) AS "計算式"
FROM Primes NN
GROUP BY num
HAVING SUM(prime) < num



■過剰数
完全数と条件の符号が異なるだけである


WITH Input(iStart, iEnd) AS (
  SELECT 1, 100
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber
  INNER JOIN Input
    ON NaturalNumber.num + 1 <= Input.iEnd
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND (A.num % B.num) = 0
)
SELECT NN.num "過剰数",
   (SELECT '+' + CAST(PF.prime AS VARCHAR(MAX)) 
    FROM Primes PF
    WHERE PF.num = NN.num
    FOR XML PATH('')
  ) AS "計算式"
FROM Primes NN
GROUP BY num
HAVING SUM(prime) > num




完全数でない擬似完全数
擬似完全数の前に、完全数をのぞいたもの


WITH Input(iStart, iEnd) AS (
  SELECT 1, 100
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber
  INNER JOIN Input
    ON NaturalNumber.num + 1 <= Input.iEnd
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND (A.num % B.num) = 0
),
Abundants AS ( 
  SELECT NN.num,
    SUM(prime) - num AS abundant
  FROM Primes NN
  GROUP BY num
  HAVING SUM(prime) > num
),
Goods AS (
SELECT P.num, prime, A.abundant
FROM Abundants A
JOIN Primes P
    ON P.num = A.num
)
,Knapsack(num, prime, pLst, sumVal) AS (
  SELECT num, prime
       , CAST(prime AS varchar(max))
       , prime
  FROM Goods
  UNION ALL
  SELECT a.num
       , b.prime
       , CAST(a.pLst + ','
            + CAST(b.prime AS varchar(max))
         AS varchar(max))
       , a.sumVal+b.prime
  FROM Knapsack a JOIN goods b
    ON a.num = b.num
   AND a.prime < b.prime
   AND a.prime <= b.abundant
) 
SELECT num "擬似完全数", MAX(pLst) "式"
FROM Knapsack
WHERE num = sumVal
GROUP BY num
ORDER BY num


最後のグループ化は無駄に
重複計算してしまうので
なんとかしたい


■擬似完全数


WITH Input(iStart, iEnd) AS (
  SELECT 1, 100
),
NaturalNumber(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM NaturalNumber
  INNER JOIN Input
       ON NaturalNumber.num + 1 <= Input.iEnd
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND (A.num % B.num) = 0
),
Parfects AS (
  SELECT NN.num
  FROM Primes NN
  GROUP BY num
  HAVING SUM(prime) = num
),
Abundants AS ( 
  SELECT num, SUM(prime) - num AS abundant
  FROM Primes
  GROUP BY num
  HAVING SUM(prime) > num
),
Goods AS (
  SELECT P.num, prime, A.abundant
  FROM Abundants A
  JOIN Primes P
    ON P.num = A.num
),
Knapsack(num, prime, pLst, sumVal) AS (
  SELECT num, prime
       , CAST(prime AS varchar(max))
       , prime
  FROM Goods
  UNION ALL
  SELECT a.num
       , b.prime
       , CAST(a.pLst + '+'
            + CAST(b.prime AS varchar(max))
         AS varchar(max))
       , a.sumVal+b.prime
  FROM Knapsack a JOIN goods b
    ON a.num = b.num
   AND a.prime < b.prime
   AND a.prime <= b.abundant
) 
SELECT num "擬似完全数", MAX(pLst) "式"
FROM Knapsack
WHERE num = sumVal
GROUP BY num
UNION
SELECT PF.num,
   (SELECT '+' + CAST(PR.prime AS VARCHAR(MAX)) 
    FROM Primes PR
    WHERE PF.num = PR.num
    FOR XML PATH('')
  ) AS "計算式"
FROM Parfects PF
ORDER BY num



ORACLEでは、上がそのまま動かない
テーブルを作成しながら作成を試みる


create table natural_numbers as 
WITH Input(iStart, iEnd) AS (
  SELECT 0, 9 FROM DUAL
),
numb(num) AS (
  SELECT iStart FROM Input
  UNION ALL
  SELECT num + 1 FROM numb
  INNER JOIN Input
    ON numb.num + 1 <= Input.iEnd
)
SELECT TO_NUMBER(k.num || l.num || m.num || n.num) num
FROM numb k, numb l, numb m, numb n
ORDER BY k.num, l.num, m.num, n.num



完全数でない擬似完全数(ORACLE)


WITH
NaturalNumber AS (
  SELECT num FROM Natural_Numbers WHERE num <= 100
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND MOD(A.num, B.num) = 0
),
Parfects AS (
  SELECT NN.num
  FROM Primes NN
  GROUP BY num
  HAVING SUM(prime) = num
),
Abundants AS ( 
  SELECT num, SUM(prime) - num AS abundant
  FROM Primes
  GROUP BY num
  HAVING SUM(prime) > num
),
Goods AS (
  SELECT P.num, prime, A.abundant
  FROM Abundants A
  JOIN Primes P
    ON P.num = A.num
),
Knapsack(num, prime, pLst, sumVal) AS (
  SELECT num, prime
       , TO_CHAR(prime)
       , prime
  FROM Goods
  UNION ALL
  SELECT a.num
       , b.prime
       , a.pLst || '+' || TO_CHAR(b.prime)
       , a.sumVal+b.prime
  FROM Knapsack a JOIN goods b
    ON a.num = b.num
   AND a.prime < b.prime
   AND a.prime <= b.abundant
)
SELECT num "擬似完全数", MAX(pLst) "式"
FROM Knapsack
WHERE num = sumVal
GROUP BY num
ORDER BY num



■擬似完全数(ORACLE)


WITH
NaturalNumber AS (
  SELECT num FROM Natural_Numbers WHERE num BETWEEN 1 AND 200
),
Primes AS (
  SELECT A.num, B.num prime
  FROM NaturalNumber A
  JOIN NaturalNumber B
    ON B.num <= A.num / 2
   AND MOD(A.num, B.num) = 0
),
Parfects AS (
  SELECT NN.num
  FROM Primes NN
  GROUP BY num
  HAVING SUM(prime) = num
),
Abundants AS ( 
  SELECT num, SUM(prime) - num AS abundant
  FROM Primes
  GROUP BY num
  HAVING SUM(prime) > num
),
Goods AS (
  SELECT P.num, prime, A.abundant
  FROM Abundants A
  JOIN Primes P
    ON P.num = A.num
),
Knapsack(num, prime, pLst, sumVal) AS (
  SELECT num, prime
       , TO_CHAR(prime)
       , prime
  FROM Goods
  UNION ALL
  SELECT a.num
       , b.prime
       , a.pLst || '+' || TO_CHAR(b.prime)
       , a.sumVal+b.prime
  FROM Knapsack a JOIN goods b
    ON a.num = b.num
   AND a.prime < b.prime
   AND a.prime <= b.abundant
)
SELECT num "擬似完全数", MAX(pLst) "式"
FROM Knapsack
WHERE num = sumVal
GROUP BY num
UNION
SELECT PF.num,
(SELECT LISTAGG(PR.prime, '+') WITHIN GROUP (order by PR.prime) FROM Primes PR  WHERE PF.num = PR.num)
FROM Parfects PF
ORDER BY "擬似完全数"

2017-01-11

Haskellで組合せ最適化(メモ)

SQLで作成したコードに対応するHaskellコードを作成した
あまり、カッコ良く書けていないですが。。。

巡回セールスマン全探索


import Data.List (maximumBy,minimumBy,permutations)
import Data.Function (on)

-- エッジのデータ型
data Edge = Edge { path :: [Int], name :: String, cost :: Integer } deriving (Show)  

-- 入力データ
nodes = zip [1..] ["Zurich", "London",  "Berlin", "Roma", "Madrid"]
edges = [(Edge [1, 2] "Zurich,London"  476)
        ,(Edge [2, 3] "London,Berlin"  569)
        ,(Edge [3, 1] "Berlin,Zurich"  408)
        ,(Edge [4, 3] "Roma,Berlin"    736)
        ,(Edge [4, 1] "Roma,Zurich"    434)
        ,(Edge [4, 2] "Roma,London"    894)
        ,(Edge [5, 4] "Madrid,Roma"    852)
        ,(Edge [5, 1] "Madrid,Zurich"  774)
        ,(Edge [5, 2] "Madrid,London"  786)
        ,(Edge [5, 3] "Madrid,Berlin"  1154)]

-- 最大・最小
-- minRoot len or maxRoot len
-- minRoot (length nodes) => ([1,3,2,5,4,1],3049)
-- maxRoot (length nodes) => ([1,5,3,4,2,1],4034)
minRoot n = minimumBy (compare `on` snd) (costRoots (enumRoots [1..n])) 
maxRoot n = maximumBy (compare `on` snd) (costRoots (enumRoots [1..n])) 

-- ルート列挙
-- enumRoots [1,2,3] => [[1,2,3,1],[1,3,2,1]]
-- enumRoots [1,2,3,4] => [[1,2,3,4,1],[1,3,2,4,1],[1,3,4,2,1],[1,2,4,3,1],[1,4,2,3,1],[1,4,3,2,1]]
enumRoots (x:xs) = map (\ys -> x:ys ++ [x]) (permutations xs)

-- ラベル付きコスト算出を展開
-- costRoots [[1,2,3,4,1],[1,2,4,3,1]] => [([1,2,3,4,1],2215),([1,2,4,3,1],2514)]
costRoots roots = map costRoot roots

-- ルート(例:[1,2,3,4,1])のラベル付きコスト算出
-- costRoot [1,2,3,4,1] => ([1,2,3,4,1],2215)
costRoot rt = (rt, sumRoot rt)

-- アトム化したルート(例:[[1,2],[2,3][3,4][4,1]])のコスト算出
-- sumRoot [1,2,3,4,1] => 2215
sumRoot rt = foldl (\acc pt -> acc + (atomCost pt)) 0 (atomRoot rt)

--パス(例:[2,3])のコスト算出
-- atomCost [2,3] => 569
atomCost pt = cost (findEdge pt)

-- パスからエッジを取得
-- formalPath [3, 2] => [2, 3]
-- formalPath [2, 3] => [2, 3]
-- findEdge [2,3] => Edge {path = [2,3], name = "London,Berlin", cost = 569}
formalPath [a, b] = if a < b then [a, b] else [b, a]
findEdge pt = head (filter (\e -> formalPath (path e) == formalPath pt) edges)

-- ルート(例:[1,2,3,1])をパス(例:[1,2])に分解
-- atomRoot [1,2,3,1] => [[1,2],[2,3],[3,1]]
atomRoot (a:b:as) = [a,b] : atomRoot (b:as)
atomRoot as = []

Connection: close