ナルシシスト数を出力する（Haskell実装）

こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz です。

はじめに

本稿は『ナルシシスト数（Narcissistic Numbers）』をHaskellで出力するコードを紹介するエントリです。

本稿で紹介する処理は、RossetaCodeのサイトに紹介されているものを使っており、オリジナルではございません。あらかじめご了承ください。

ナルシシスト数

自己愛が強いことをナルシシズムといいますが、それに由来する名称がついている正の整数があります。まずは『ナルシシスト数』の定義を紹介します。

ナルシシスト数の定義

$n$ 桁の正の整数 $N$ について，各桁の $n$ 乗の和がもとの数に等しいとき, $N$ をナルシシスト数と言います。

は、何故この特徴で“自己愛が強い”という意味の名前がつけられているのかは不思議に思いますが、$153 = 1^3 + 5^3 + 3^3$や$370 = 3^3 + 7^3 + 0^3$のような数だそうです。また、１桁の正の整数は全て条件を満たしており、$1 \sim 9 $はナルシシスト数に該当します。
そして、ナルシシスト数が有限個しか存在しない ことは、簡単に証明できる¹らしく、次のサイトに詳しく証明が掲載されております。

この証明によると「ナルシシスト数は60桁以下」ということが示されています。また、Wikipedia では、ナルシシスト数は（0を含めないならば）全部で88個存在し、その最大のものは39桁という紹介もされています。証明では60桁以下は示せていますが、実際はもっと少ない桁数の数（39桁）が最大数ということが、計算で知られているということのようです。
Narcissistic Number | WolframMathWorld のサイトに88個のリストが紹介されておりますので、興味のある方は確認してみてください。

全部で88個しかないのですから、全てのナルシシスト数はプログラムで求められそうです。次にナルシシスト数を求めるプログラムをみていきましょう。

RosettaCode

ナルシシスト数を求めるプログラムは、RosettaCodeにタスクがありました。RosettaCodeはロゼッタストーンに発想を得たサイトです。同じタスクの処理を様々なプログラミング言語で書かれたソースコードが、全て１ページに紹介されているサイトです。
本稿では、Haskell実装を取り上げたいので次のリンクにある処理を用います。

Haskellの実装では２種類のプログラムが紹介されていますが、処理の効率が良いと書かれている "Reduced search (unordered digit combinations)" の方を採用します。次の折りたたみの中でプログラムを紹介します。※ただし、が一部修正しております。

20桁までのナルシシスト数を求めるHaskellのプログラム

RosettaCodeでは7桁までの処理ですが、Integer型に対応して、20桁までを求めるように修正しました。

narcissistic.hs

import Data.Bifunctor (second)

nn :: Integer
nn = 20

narcissiOfLength :: Integer -> [Integer]
narcissiOfLength nDigits = snd <$> go nDigits []
  where
    powers = ((,) <*> (^ nDigits)) <$> [0 .. 9]
    go :: Integer -> [(Integer, Integer)] -> [(Integer, Integer)]
    go n parents
      | 0 < n = go (pred n) (f parents)
      | otherwise = filter (isDaffodil nDigits . snd) parents
      where
        f :: [(Integer, Integer)] -> [(Integer, Integer)]
        f parents
          | null parents = powers
          | otherwise =
            parents >>=
            (\(d, pwrSum) -> second (pwrSum +) <$> take (succ (fromIntegral d)) powers)

isDaffodil :: Integer -> Integer -> Bool
isDaffodil e n =
  (((&&) . (e ==) . (fromIntegral . length)) <*> (n ==) . powerSum e) (digitList n)

powerSum :: Integer -> [Integer] -> Integer
powerSum n = foldr ((+) . (^ n)) 0

digitList :: Integer -> [Integer]
digitList n
 | n < 1 = []
 | otherwise = digitList (n `div` 10) ++ [n `mod` 10]

--------------------------- TEST ---------------------------
main :: IO ()
main =
  putStrLn $
  fTable
    ("Narcissistic decimal numbers of length 1-" ++ (show nn) ++ ":\n")
    show
    show
    narcissiOfLength
    [1 .. nn]

fTable :: String -> (a -> String) -> (b -> String) -> (a -> b) -> [a] -> String
fTable s xShow fxShow f xs =
  let rjust n c = drop . length <*> (replicate n c ++)
      w = maximum (length . xShow <$> xs)
  in unlines $
     s : fmap (((++) . rjust w ' ' . xShow) <*> ((" -> " ++) . fxShow . f)) xs

以下のように、プログラムをコンパイルします。

% ghc narcissistic.hs
[1 of 2] Compiling Main             ( narcissistic.hs, narcissistic.o ) [Source file changed]
[2 of 2] Linking narcissistic [Objects changed]

実行結果

% time ./narcissistic
Narcissistic decimal numbers of length 1-20:

 1 -> [1,2,3,4,5,6,7,8,9]
 2 -> []
 3 -> [153,370,371,407]
 4 -> [1634,8208,9474]
 5 -> [54748,92727,93084]
 6 -> [548834]
 7 -> [1741725,4210818,9800817,9926315]
 8 -> [24678050,24678051,88593477]
 9 -> [146511208,472335975,534494836,912985153]
10 -> [4679307774]
11 -> [32164049650,32164049651,40028394225,42678290603,44708635679,49388550606,82693916578,94204591914]
12 -> []
13 -> []
14 -> [28116440335967]
15 -> []
16 -> [4338281769391370,4338281769391371]
17 -> [21897142587612075,35641594208964132,35875699062250035]
18 -> []
19 -> [1517841543307505039,3289582984443187032,4498128791164624869,4929273885928088826]
20 -> [63105425988599693916]

./narcissistic  235.38s user 0.21s system 99% cpu 3:57.03 total

処理時間も計測しましたが、4分ほど待つと計算結果が出力されました。結果があっているかを確認してみましょう。WolframMathのリストと比較すると合致していることが確認できました。

RosettaCodeのままでは、19桁までは合致していますが、20桁のナルシシスト数が出力されない結果 になりました。HaskellではInteger型は多倍長整数を扱えるため桁数を気にしない（処理時間は膨大になる可能性があるため注意）で計算されますが、RosettaCodeのサンプルプログラムではInt型で計算されておりmaxBound :: Intの値9223372036854775807（19桁）より大きな数が扱えなくなるためです。

おわりに

ご一読いただきまして有り難うございます。

本来ならば、39桁までリストを作りたいところです。プログラム中のnn = 20をnn = 39に書き換えてコンパイル＆実行すれば出力は可能だと思います。が、かなりの時間を要すると思われるため未検証です。
ちなみに、21桁は試してみました。

21 -> [128468643043731391252,449177399146038697307]

が出力され、処理時間は385.76s(6分半)ほどかかりました。

公開後にその先も試しましたが、

22 -> []
23 -> [21887696841122916288858,27879694893054074471405,27907865009977052567814,28361281321319229463398,35452590104031691935943]
24 -> [174088005938065293023722,188451485447897896036875,239313664430041569350093]
25 -> [1550475334214501539088894,1553242162893771850669378,3706907995955475988644380,3706907995955475988644381,4422095118095899619457938]

1時間ほどで25桁までは計算できました。
39桁までを試すには、何か処理を工夫したいところですが、本稿ではここまでとさせて頂きます。

最後に、本稿を記載するために検証したHaskell環境を記しておきます。お手元の環境で検証する際に動作が異なる場合などは、参考になるかもしれません。

本稿の環境

本稿のために使用した環境は以下となります。
macOS: Sonoma 14.4 (chip: Apple M1)
GHCup: 0.1.22.0
GHC: 9.6.4

（●）（●） Happy Hacking!
/""　＿＿""＼

1. 中学数学・高校数学の範囲でもかなりのところまで理解できるそうです（https://www.eishinkan.net/creation/ikkan_class/blog/math/column_4.php） ↩