雑記帳

自然数係数の多項式の係数を取り出す

前書き
この記事は、この問題をネットで目にした際に「任意精度演算ができる Haskell を使えばシンプルに実装できそうじゃん」ということで、実際にやってみただけのネタ記事。
問題の内容
自然数を係数に持つ多項式
\[ f(x) = \sum_{i = 0}^n a_i \cdot x^i \]
から、その多項式の係数となる全ての自然数 \(a_i\) を把握することを考えたいとき、その過程で必要になる「その多項式に値を代入する操作の回数」はたった2回だけで済んでしまうというものである。
Haskell で実践的に求めてみる
具体的な解説は、ちゃんとした数学者により行われたものを参考にするのがベストであると思うが、大雑把には
  • まず必要になる1回目の代入操作は、任意の係数 \(a_i\) よりも大きいことが明らかな任意の1つの自然数を得ることを目的として、全ての係数の総和にあたる \(f(1)\) を求めるために行う。
  • \(f(1)\) の10進数表記を行った際の桁数 \(k\) を求める。
  • 続いて2回目の代入操作として \(f(10^k)\) を求める。その数は10進数表記で表示すると、多項式の各々の係数が長さ \(k\) のブロック毎に並べられた「潜在的に全ての係数のリストと見做せる1つの自然数」になる。
  • 数学的にはもう既に十分ではあると思うが、その数から実際に自然数のリストを抽出する。
といった感じで求めることができる。
実際のコードは、ページの終わりに載せておくが、例として
\[ f(x) = 55 x^2 + 2 x^3 + 192 x^4 + 9 x^6 \]
として定まる多項式について処理を行ってみると、以下の結果が得られる。
ghci> main
[0,0,55,2,192,0,9]
無名関数を使っても
ghci> getListOfCoefficients(\x -> 55*x^2 + 2*x^3 + 192*x^4 + 9*x^6)
[0,0,55,2,192,0,9]
というように同様の結果が得られる。
他の例もいくつか挙げておくと、
ghci> getListOfCoefficients(\x -> 0)
[0]
ghci> getListOfCoefficients(\x -> 1)
[1]
ghci> getListOfCoefficients(\x -> 10 + 9*x + 8*x^2 + 5*x^5)
[10,9,8,0,0,5]
ghci> getListOfCoefficients(\x -> x^21)
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1]
因みに、ここではわかりやすいように10進数表記を用いているが、10進数であることが必須であるわけではない。
ソースコード
main :: IO ()
main = do
  print $ getListOfCoefficients(f)


f x = foldr (+) 0 $ do
  i <- [0..(length(listOfCoefficients) - 1)]
  return $ a(i) * (x ^ i)
  where
    a i =
      listOfCoefficients !! i
    listOfCoefficients =
      [ 0
      , 0
      , 55
      , 2
      , 192
      , 0
      , 9 ]

getListOfCoefficients :: (Num a, Integral b, Show b) => (a -> b) -> [Integer]
getListOfCoefficients fnc =
  let
    sumOfAllCoeffients = fnc(1)
    numberWhichWillPotentiallyBeListOfCoefficients = fnc(10 ^ numOfDigitsInDecimal)
    numOfDigitsInDecimal =
      if sumOfAllCoeffients == 0 then
        1
      else
        1 + (floor $ (logBase 10) $ fromIntegral sumOfAllCoeffients)
    xs = reverse . (fmap read) . words . reverse $ do
      (i,d) <- zip [0..] (reverse $ show numberWhichWillPotentiallyBeListOfCoefficients)
      if (i `mod` numOfDigitsInDecimal) == (numOfDigitsInDecimal - 1) then
        [d, ' ']
      else
        return d
  in
    xs