フィボナッチ in L

10689番が英文が短かったので解いてみた。
昨日のに比べて簡単過ぎたかもしらん…
解いてる人も多いし。
問題はf(0)=a,f(1)=b,f(n)=f(n-1)+f(n-2)の漸化式で定められる
数列fの何番目かを求めるというもの。
初期値a,bで生成される数列のn番目,下位m桁を出力する。
いわゆるフィボナッチなのだが、喜び勇んで

head $ take n $ fix $ \fib -> a:b:zipWith (+) fib (tail fib)

このようなコードを書いてはこの問題では駄目である。
というのもnの最大が10億なので、O(n)のアルゴリズムでは到底間に合わない。
(というか、Haskell自体が使えないけど…)


f(n-1),f(n) と f(n),f(n+1) の間の関係は
以下のような線形写像になっているので、
\Large\left(\array{f_n\\f_{n+1}}\right)=\left(\array{\0&1\\1&1}\right)\left(\array{f_{n-1}\\f_n}\right)
これを用いて
\Large\left(\array{f_n\\f_{n+1}}\right)=\left(\array{\0&1\\1&1}\right)^n\left(\array{a\\b}\right)
f(n)をこのような式から求めることが出来る。
累乗はO(logn)で計算できるので、結局フィボナッチはO(logn)で計算できる
…というのはそれなりに有名な話である。


ということで、ちゃちゃっと実装してAccept。
非常に気分が良い。
(あと、細かい点でだが、
(a*b) mod n = a mod n * b mod n になることを
適当に使っている)

#include 
using namespace std;

int powi(int b,int n){ return (n==0)?1:b*powi(b,n-1); }

void mulm(int *dest,int *a,int *b,int mod){
  dest[0]=(a[0]*b[0]+a[1]*b[2])%mod;
  dest[1]=(a[0]*b[1]+a[1]*b[3])%mod;
  dest[2]=(a[2]*b[0]+a[3]*b[2])%mod;
  dest[3]=(a[2]*b[1]+a[3]*b[3])%mod;
}

void powm(int *dest,int *base,int n,int mod){
  if (n==0){
    dest[0]=dest[3]=1;
    dest[1]=dest[2]=0;
  }
  else if (n&1){
    int tmp[4];
    powm(tmp,base,n-1,mod);
    mulm(dest,tmp,base,mod);
  }
  else{
    int tmp[4];
    powm(tmp,base,n/2,mod);
    mulm(dest,tmp,tmp,mod);
  }
}

int main(){
  int cases;cin>>cases;
  while(cases--){
    int a,b,n,m;cin>>a>>b>>n>>m;
    int mod=powi(10,m);
    int d[4],s[4]={0,1,1,1};
    powm(d,s,n,mod);
    cout<<(d[0]*a+d[1]*b)%mod<

ついでに(というか、こっちがメインかも)Haskellでも実装。
最近覚えたイディオム、

(f.) . g

をふんだんに使ってみた。

import Data.List

main = getContents >>= mapM_ (print.solve.map read.words).tail.lines

solve :: [Int] -> Int
solve [a,b,n,m] = (head $ powm n [ [0,1],[1,1] ]) `inner` [a,b] where
  powm n m | n == 0 = [ [1,0],[0,1] ]
           | n `mod` 2 == 1 = m `mulm` powm (n-1) m
           | otherwise = let t = powm (n `div` 2) m in t `mulm` t

  mulm a b = [ [la `inner` lb | la <- a] | lb <- transpose b]
  inner = (sum.) . zipWith (((`mod` base).).(*))
  base = powi m 10
  powi = (product.) . replicate

結構短くなってこっちも大満足である。
関数合成演算子(.)は

(.) :: (b -> c) -> (a -> b) -> a -> c

のような型をもつので、2引数関数を合成できない。
a -> b -> c と c -> d を合成して a -> b -> d のような
関数を作りたいときに先のテクニックが使える。

((+1).) . (*)

とかとすると掛け算して、その結果に1足すとかになる。
(.)を使ったテクニックはもっと色々あるらしいが、
あまり頭の中でイメージできない…。
高階プログラミングはまだまだ極められないなぁ。
どこまで勉強しても上が見えてこないような絶望感が。