G'day all.

Quoting [EMAIL PROTECTED]:

There is one solution missing there (unless I skipped it) fib
n=((1+s)/2)^n-((1-s)/2)^n)/s where s=sqrt 5 If some of you complain
that this is real, not integer, please remember that
Leonardo of Pisa thought of applying this to rabbits. Well, rabbits are
not integers, they eat carrots and have long ears. They are real thing.

As noted, floating-point arithmetic diverges from integer arithmetic
fairly quickly in this case.

Of course, we can avoid this by doing computations in the
field extension Q[sqrt 5]:

    data QS5 = QS5 Rational Rational

    infixl 7 <*>,</>
    infixl 6 <->,<+>

    conjugate :: QS5 -> QS5
    conjugate (QS5 a1 a2) = QS5 a1 (negate a2)

    (<+>),(<->),(<*>),(</>) :: QS5 -> QS5 -> QS5
    (QS5 a1 a2) <+> (QS5 b1 b2) = QS5 (a1+b1) (a2+b2)
    (QS5 a1 a2) <-> (QS5 b1 b2) = QS5 (a1-b1) (a2-b2)
    (QS5 a1 a2) <*> (QS5 b1 b2) = QS5 (a1*b1 + 5*a2*b2) (a1*b2 + a2*b1)
    a@(QS5 a1 a2) </> b@(QS5 b1 b2)
        = let QS5 c1 c2 = a <*> conjugate b
              s = (b1*b1 - 5*b2*b2)
          in QS5 (c1 / s) (c2 / s)

    qpow :: QS5 -> Integer -> QS5
    qpow q n
        | n < 3 = case n of
                      0 -> QS5 1 0
                      1 -> q
                      2 -> q <*> q
        | even n = let q' = qpow q (n `div` 2) in q' <*> q'
        | otherwise = let q' = qpow q (n `div` 2) in q' <*> q' <*> q

    fib ::Integer -> Integer
    fib n
        = let (QS5 fn _) = (qpow phi n <-> qpow phi' n) </> s5 in numerator fn
        where
            phi = QS5 (1%2) (1%2)
            phi' = QS5 (1%2) (negate 1%2)
            s5 = QS5 0 1

However, this is still an O(log n) algorithm, because that's the
complexity of raising-to-the-power-of.  And it's slower than the
simpler integer-only algorithms.  It might be amusing to see if this
could be transformed into one of the simpler algorithms, though.

Cheers,
Andrew Bromage
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Reply via email to