2 進表現

\begin{array}{rcl} E_{m+n} &=& (F_{3m+3n},F_{3m+3n+1})\\ &=& (F_{3m+1}F_{3n}+F_{3m}F_{3n-1},F_{3m+1}F{3n+1}+F{3m}F{3n})\\ &=& (F_{3m+1}F_{3n}+F_{3m}(F_{3n+1}-F_{3n}),F_{3m+1}F_{3n+1}+F_{3m}F_{3n}) \end{array}

これは,E_{m+n}F_{3m}F_{3m+1}F_{3n}F_{3n+1}で構成できる.ということは,E_{m+n}E_{m}E_{n}で構成できる.つまり,
E_{m} \oplus E_{n} = E_{m+n}
となるような演算子\oplusを定義できる.そうするとE_{n}は,nの 2 進表現を逆順にならべたものと,
E_{2^0},E_{2^1},E_{2^2},\cdots,E_{2^i},\cdots
を対応させて構成できる.このときE_{n}O(\log n)で構成できる.これをHaskellのコードにすると以下のようになる.これで4\cdot10^{10^6} でもなんとかなりそう.

p0002 u = (a+b-1) `div` 2
  where
    (a,b) = foldr1 f . takeWhile ((u>=).fst) $  bfib3s
    x `f` y | fst z > u = y
            | otherwise = z
                    where z = x `oplus` y

oplus :: (Integer,Integer) -> (Integer,Integer) -> (Integer,Integer)
(a,b) `oplus` (c,d) = (b*c + a*d - ac, b*d + ac) where ac = a*c

bfib3s :: [(Integer,Integer)]
bfib3s = iterate dbl (2,3)

dbl :: (Integer,Integer) -> (Integer,Integer)
dbl x = x `oplus` x