豊福です。

  算譜の記 2004-04-28 の「n次元の非負整数格子点の列挙(重み
可変版) 」を印刷だけして読まずに忘れていたのが一昨日たまたま
見つかったのでやっと読めました。
  短いプログラムに感心しながらも理解しようとすると頭が爆発し
そうになるのでわかりやすい版(本当?)を書いてみました。
コードは整理されていませんが取りあえず動くということで。

-- R n [i_m, ... ,i2,i1] i f :: 多次元の矩形領域配列
--   (m+n)次元空間中の
--   x_1==i_1, x_2==i_2, ... , x_m==i_m,
--   x_m_1>=i, x_m_2>=0, ... , x_m_n>=0 が表わす n次元矩形
--
--   格子点 (x1,x2, ... ,x_m_n) での値が f(x1,x2, ... ,x_m_n)
--
-- C d1 d2 [i_n, ... ,i2,i1] :: 2つの領域 d1 と d2 を組み合わせた領域
--   [i1,i2, ... , i_n] :: 領域の最小座標
--   i1 == min(d1 の 最小x1座標, d2 の 最小x1座標), ...

data D a = R Int [Int] Int ([Int] -> a) | C (D a) (D a) [Int]

-- head :: 矩形の head は一次元下がった矩形
hd (R n idxs i f) = R (n-1) (i:idxs) 0 f
-- tail :: 矩形の tail は1つずれた矩形
tl (R n idxs i f) = R n idxs (i+1) f

-- 矩形の角(かど)の座標(最小座標)
ind n idxs i = replicate (n-1)  0 ++ (i:idxs)

-- 領域の最小座標
getInd (R n idxs i f) = ind n idxs i
getInd (C _ _ idxs) = idxs

-- 座標の比較
le :: [Int] -> [Int] -> Bool
le x y = and [ix <= iy | (ix,iy) <- zip x y]

-- 1次元直線のソートは簡単
fsort (R 1 idxs i f) = [f (j:idxs) | j <- [i..]]
-- fsort x@(R n idxs i f) = fsort (C (hd x) (tl x) [])
fsort x@(R _ _ _ _) = fsort2 x
-- 最小値が矩形領域側にあるとわかっている場合は特別に処理
fsort (C x@(R _ _ _ _)  y _) =
  if (le (getInd x) (getInd y)) then
    fsort3 x y
  else
    merge (fsort x) (fsort y)
fsort (C x y _) = merge (fsort x) (fsort y)

-- 矩形のソート
fsort2 x = t: fsort r
   where (t, r) = getMinFromR x

fsort3 x y = t: merge (fsort r) (fsort y)
   where (t, r) = getMinFromR x

-- 矩形から最小値だけを取り出して残りの領域と一緒に返す
---  1次元は簡単
getMinFromR (R 1 idxs i f) = ((f (i:idxs)), (R 1 idxs (i+1) f))
getMinFromR x@(R n idxs i f) = (t, C d (tl x) (ind n idxs i))
  where
    (t, d) = getMinFromR (hd x)

-- merge は以前からのと本質的に同じ
merge xs [] = xs
merge [] ys = ys
merge xxs@(x@(_,fx):xs) yys@(y@(_,fy):ys)
  | fx <= fy = x: merge xs yys
  | otherwise = y: merge xxs ys

-- 2次元格子点の列挙
latnds_2dim =
    fsort
  $ R 2 [] 0 (\ idxs -> (idxs, sum idxs)) -- sum で構成される2次元矩形領域

-- ramanujan数
ramanujan =
    filter (not.null.drop 1)
  $ groupBy (\ (_,fx) (_,fy) -> fx == fy)
  $ fsort
  $ tri 0 (\ p@[x,y] -> (p, x^3 + y^3)) -- 三乗和で構成される三角領域
      where -- 1次元直線を繰り返しずらしながら組み合わせたものが三角領域
        tri i f = C (R 1 [i] i f) (tri (i+1) f) [i,i]

---
                        豊福
                        nobu_toyofuku / nifty.com


--
ML: haskell-jp / quickml.com
使い方: http://QuickML.com/