この記事は ひとりアドベントカレンダーRosettaCodeで楽しむプログラミング Advent Calendar 2025の22日めの記事です。
定義
偶数桁の自然数 $n$ で、$n = a \times b$ となる約数を持ち、$a$ と $b$ が以下の条件を満たす数
- $a$ と $b$ はどちらも $n$ の半分の桁数
- 10進数表記で、$a$と$b$ を合わせると $n$ と数字の構成が同じ
- 末尾の桁の0を持つ数はたかだか一つだけ、またその個数は1つまで
(原文: at most one of them has a trailing zero, 強調筆者)
$a$,$b$ を $n$ の fangs と呼ぶ。
タスク
- vampire数の最初の25個を、その fangs とともに示せ。
- 以下の数はvampire数か。もしそうなら、fangsも示せ。
- 16758243290880, 24959017348650, 14593825548650
注意：一つのvampire数に対して、複数のfangsがありうる。
考える
与えられた数がvampire数か判定するには、二つのアプローチが考えられる。
まず、約数を探すときと同じ試し割りを、固定の桁数の全ての数で行い、割り切れたときに数字の構成を確認する方法。
もう一つは、与えられた数に出現する数字を半分個使って作れる数を昇順に生成し、それについて試し割りをする方法。
$n$ 桁の数を任意に選ぶと場合の数は $9 \cdot 10^{n-2}$ 通り、
$2n$ 個の数字から $n$ 個を選んで並べるやり方は、重複を無視して $\prod_{i = n+1}^{2n} i$通り。
$n$ が小さいうちは後者、$n \geq 5$ になると前者が有利だが、個数を管理する配列の操作を伴うので係数が大きい感じ。
スケーラビリティのある後者を選んで
やってみよう。
与えられた整数 x に対して、そのfangsを全て挙げる。ひとつでもある場合、それはvampire数である。
まず、Integer の数字を数える下請け。
import Data.Array.Unboxed
import Data.Char
-- 数字の個数
type Count = UArray Int Int
mkCount :: Integer -> Count
mkCount = accumArray (+) 0 (0,9) . map ((\d -> (d,1)) . digitToInt) . show
Count で使える数字だけを選んで、指定した桁の数を昇順に全て生成する。
-- CountからInt個の数字をちょうど使って作れる数を昇順で
generate :: Count -> Int -> [Integer]
generate cnt0 w0 = recur 0 cnt0 (reverse $ take w0 $ iterate (10 *) 1) []
where
recur :: Integer -- ここまでの結果
-> Count -- 使える数字
-> [Integer] -- 次に使う数を入れる位
-> [Integer] -- これ以降の結果
-> [Integer] -- 答え
recur val _nt [] rest = val : rest
recur val cnt (w:ws) rest = foldr ($) rest
[ recur val1 cnt1 ws
| (d,c) <- assocs cnt, c > 0
, let val1 = val + fromIntegral d * w, val1 > 0 -- 最上位桁に0を置くことを防ぐ
, let cnt1 = cnt // [(d, pred c)] ]
ghci> generate (mkCount 1260) 1
[1,2,6]
ghci> generate (mkCount 1260) 2
[10,12,16,20,21,26,60,61,62]
ghci> generate (mkCount 1260) 3
[102,106,120,126,160,162,201,206,210,216,260,261,601,602,610,612,620,621]
ghci> generate (mkCount 1260) 4
[1026,1062,1206,1260,1602,1620,2016,2061,2106,2160,2601,2610,6012,6021,6102,6120,6201,6210]
x に対して桁半分で generate した結果で割り算して、割り切れた場合の商について、
x と Count の内容が一致する、とすると、generate の内部で作られた cnt が無駄になる。
機能分割の点からは残念だが、generate の生成段階に手を加え、
fangsの対をなすようなものだけを返す、という専用構成にする。
fangs :: Integer -> [(Integer,Integer)]
fangs x
| odd w2 = []
| otherwise = generate1 0 cnt0 (reverse $ take w1 $ iterate (10 *) 1) []
where
cnt0 = mkCount x
w2 = sum $ elems cnt0
w1 = div w2 2
generate1 val cnt [] rest
| cond = (val, q) : rest
| otherwise = rest
where
(q,r) = divMod x val
cond = and
[ r == 0 -- valはxを割りきる
, cnt == mkCount q -- valで使った数字の残りとqの数字が一致
, mod val 10 /= 0 || mod q 10 /= 0 -- 少なくともどちらかは1の位が0でない
, mod val 100 /= 0, mod q 100 /= 0 -- どちらも1,10の位ともに0ではない
]
generate1 val cnt (w:ws) rest = foldr ($) rest
[ generate1 val1 cnt1 ws
| (d,c) <- assocs cnt, c > 0
, let val1 = val + fromIntegral d * w, val1 > 0
, val1 * val1 <= x -- √xを超えない範囲だけ調べる
, let cnt1 = cnt // [(d, pred c)] ]
タスク1は25個見つかるまで数を順に試す。
ただし、奇数桁の数は生成するだけ無駄なので飛ばす。
task1 = take 25
[ (x, fs)
| lb <- iterate (100 *) 10 -- 始 10, 1000, 100000, ...
, x <- [lb .. pred $ 10 * lb] -- 終 99, 9999, 999999, ...
, let fs = fangsof x, not $ null fs
]
task1 = take 25
[ (x, fs)
| lb <- iterate (100 *) 10
, x <- [lb .. pred $ 10 * lb]
, let fs = fangs x, not $ null fs
]
main :: IO ()
main = do
mapM_ print task1
mapM_ (print . \x -> (x, fangs x)) [16758243290880, 24959017348650, 14593825548650]
ghci> main
(1260,[(21,60)])
(1395,[(15,93)])
(1435,[(35,41)])
(1530,[(30,51)])
(1827,[(21,87)])
(2187,[(27,81)])
(6880,[(80,86)])
(102510,[(201,510)])
(104260,[(260,401)])
(105210,[(210,501)])
(105264,[(204,516)])
(105750,[(150,705)])
(108135,[(135,801)])
(110758,[(158,701)])
(115672,[(152,761)])
(116725,[(161,725)])
(117067,[(167,701)])
(118440,[(141,840)])
(123354,[(231,534)])
(124483,[(281,443)])
(125248,[(152,824)])
(125433,[(231,543)])
(125460,[(204,615),(246,510)])
(126027,[(201,627)])
(126846,[(261,486)])
(16758243290880,[(1982736,8452080),(2123856,7890480),(2751840,6089832),(2817360,5948208)])
(24959017348650,[(2947050,8469153),(2949705,8461530),(4125870,6049395),(4129587,6043950),(4230765,5899410)])
(14593825548650,[])
(8.53 secs, 7,093,113,032 bytes)
Rosetta Code に載っている結果と一致しない。他の人の答えだと
120600 = 201 x 600
125500 = 251 x 500
が混入している。Discussionページを見ると
Shouldn't the first 26 vampire numbers be printed?
The trailing zeros condition excludes 126000 = 210 * 600,
which is the first number that breaks only the third condition.
--Choroba (talk) 22:13, 16 August 2013 (UTC)
I agree, that would test one of the exclusion rules.
-- Gerard Schildberger (talk) 22:06, 25 October 2018 (UTC)
というやりとりがあり、どうやら、このページにある答えはこの点で全て間違っているようだ。
（あるいは、trailing zero 条件が後出しで追加されたか。）
ページに掲載されているHaskellプログラム
またしゃらくさいコードでいちいち読みにくい。
ボトムアップに、下請けの、約数を生成している部分から見てみる。
integerFactors :: Int -> [Int]
integerFactors n
| n < 1 = []
| otherwise =
lows ++
(quot n <$>
(if intSquared == n -- A perfect square,
then tail -- and cofactor of square root would be redundant.
else id)
(reverse lows))
where
(intSquared, lows) =
(^ 2) &&& (filter ((0 ==) . rem n) . enumFromTo 1) $
floor (sqrt $ fromIntegral n)
なんか無駄に高度なHaskell使おうとして滑ってない？直してみる。
integerFactors :: Int -> [Int]
integerFactors n
| n < 1 = []
| otherwise = lows ++ map (quot n) ((if rn * rn == n then tail else id) (reverse lows))
where
rn = floor $ sqrt $ fromIntegral n
lows = filter ((0 ==) . rem n) [1 .. rn]
lows の数を作るときに rem n で割り算してからもう一度 quot n しているのが気になる。
こういうのはどうだろう。
factors n
| last ls == head us = ls ++ tail us
| otherwise = ls ++ us
where
(ls, us) = fmap reverse $ unzip $ filter ((n ==) . uncurry (*)) $
takeWhile (uncurry (<=)) $ map (id &&& div n) [1 ..]
剰余をメモリに残してごちゃごちゃするより、掛け算をやり直す方が速かった。メモリ確保はコストが高い。
次はもう本丸。
fangs :: Int -> [(Int, Int)]
fangs n
| odd w = []
| otherwise = ((,) <*> quot n) <$> filter isfang (integerFactors n) -- 約数全てをisfangにかける
where
ndigit :: Int -> Int
ndigit 0 = 0
ndigit n = 1 + ndigit (quot n 10)
w = ndigit n -- nの桁数
xmin = 10 ^ (quot w 2 - 1) -- fangの下限（含…まない？）
xmax = xmin * 10 -- 上限（含まず）
isfang x =
x > xmin &&
x < y && -- わざわざ integerFactors で reverse lows したのは何だったのか
y < xmax && -- same length
(quot x 10 /= 0 || quot y 10 /= 0) && -- not zero-ended
sort (show n) == sort (show x ++ show y) -- 数字の出現が同じ
where
y = quot n x
色々とちぐはぐな感じ。
そもそも、約数を調べるときに桁数の範囲だけ調べれば済むし、後半を捨てるなら初めから返さなくていい。
ndigit で桁数を数えたかと思えば、sort (show n) == sort (show x ++ show y) という極めてズボラな方法を使っている。
どうせ $10^{w/2-1}$ ぴったりな数は規則3で蹴られるだろうとはいえ、x > xmin は x == xmin の場合を捨てている。
そして "not zero-ended" は、それだけでは足らない。
というか quot ではなく rem でないとおかしい。バグっている。
直して、条件3にも対応させる。factors は n 専用で内蔵。
fangs :: Int -> [(Int, Int)]
fangs n =
[ (x,y)
| even w
, x <- factors, let y = quot n x
, rem x 10 /= 0 || rem y 10 /= 0 -- 修正
, rem x 100 /= 0, rem y 100 /= 0 -- 追加
, sort (show n) == sort (show x ++ show y) ]
where
w = length $ show n
xmin = 10 ^ (quot w 2 - 1)
xmax = xmin * 10
factors = takeWhile (xmax >) $ map fst $ filter ((n ==) . uncurry (*)) $
takeWhile (uncurry (<=)) $ map (id &&& div n) [xmin ..]
あとはドライバ部分。
vampires :: [Int]
vampires = filter (not . null . fangs) [1 ..]
main :: IO ()
main =
mapM_
(print . ((,) <*>) fangs)
(take 25 vampires ++ [16758243290880, 24959017348650, 14593825548650])
vampires の中で filter するために fang を使ったのに、もう一度表示のために main で呼びなおすとか、なんかもう！という感じ。奇数桁の数を除けることもしない。
修正版での実行結果：
(2.79 secs, 8,168,476,120 bytes)
それでもこっちの方がずっと速いの、やんなっちゃうな。
桁数が Int で収まらないくらいに増えたら事情は違ってくるだろうけども。
凝った計算
Rosetta Codes のページからリンクがあるvampire search algorithmは消えているが、
webarchiveに残っている。
$a \times b = n$ で、$a,b$ は $n$ のfangsの一つだとする。
同じ数字構成になっているという定義から、当然、数字根を足した結果も同じになる。
数字根を $r(n)$ とすると $r(r(a) + r(b)) = r(n)$ である。
これだけだと情報量が減っているだけにしか見えないが、
数字根が9なら元の数は9の倍数、というよく知られた関係を使うと様子が変わってくる。
つまり $r(x) \equiv x \bmod 9$ である。
最初の式 $a \times b = n$ から普通に $a \times b \equiv n \mod 9$ となるが、
これは $r(a) \times r(b) \equiv r(n) \mod 9$ ということでもある。
つまり $a + b \equiv a * b \mod 9$ が fangs の条件である。
さて、$0 \leq r(a), r(b) < 9$ について、この条件を満たす組み合わせを考えると：
ghci> [(a,b,mod (a+b) 9) | a <- [0 .. 8], b <- [0 .. 8], mod (a * b) 9 == mod (a + b) 9]
[(0,0,0),(2,2,4),(3,6,0),(5,8,4),(6,3,0),(8,5,4)]
$r(n) \in {0,4}$ の場合しか考えなくてよい。
また、$r(n) = 0$ のとき、$r(a) \in {0,3,6}$ の場合だけ、
$r(n) = 4$ のとき、$r(a) \in {2,5,8}$ の場合だけ考えればよい。
これは大きく候補を減らすだろう。
fangs :: Int -> [(Int, Int)]
fangs n
| odd w = []
| elem rn [0,4] = ans -- NEW!
| otherwise = []
where
w = length $ show n
rn = mod n 9 -- NEW!
xs | rn == 0 = [0,3,6] -- NEW!
| rn == 4 = [2,5,8]
xmin = 10 ^ (quot w 2 - 1)
xmax = xmin * 10
factors = takeWhile (xmax >) $ map fst $ filter ((n ==) . uncurry (*)) $
takeWhile (uncurry (<=)) $ map (id &&& div n) [xmin ..]
ans =
[ (x,y)
| x <- factors
, elem (mod x 9) xs -- NEW!
, let y = quot n x
, mod x 10 /= 0 || mod y 10 /= 0
, mod x 100 /= 0, mod y 100 /= 0
, sort (show n) == sort (show x ++ show y) ]
(1.09 secs, 2,373,380,024 bytes)
3倍速くなった。