Вычисление фибоначчи

Мне была отправлена эта красивая нерекурсивная функция для вычисления последовательности фибоначчи.

Итак, я закодировал немного С# и смог проверить, что все номера до 1474 были правильными.

Проблема возникает при попытке вычислить ее для 1475 и выше. Мои математические навыки в С# просто не справляются с задачей выяснения другого пути. Итак, есть ли у кого-то лучший способ выразить эту конкретную математическую функцию в С#? кроме традиционного способа выполнения рекурсивной функции?

Кстати, я начал использовать BigInteger в качестве возвращаемого типа. Но проблема действительно возникает при попытке повысить (1 + Math.Sqrt(5)/2) до 1475-й мощности. Я просто не вижу, какой тип данных мне нужен (и механизм в этом отношении), чтобы заставить это вернуться с чем-то другим, кроме бесконечности.

Здесь начальная точка.

private Double FibSequence(Int32 input) {
    Double part1 = (1 / Math.Sqrt(5));
    Double part2 = Math.Pow(((1 + Math.Sqrt(5)) / 2), input);
    Double part3 = Math.Pow(((1 - Math.Sqrt(5)) / 2), input);

    return (part1 * part2) - (part1 * part3);
}

И нет, это не домашнее задание. Просто "простая" проблема для медленного дня.

Ответы

Ответ 1

Я не думаю, что у С# есть тип данных с достаточной плавающей точностью и диапазоном, чтобы справиться с этим наивно.

Если вы действительно хотите пойти по этому пути, вы можете заметить, что сопряжение $\Phi=\phi^{-1}=\phi-1=\frac{-1+\sqrt 5}{2}$ меньше единицы, поэтому $-\frac{(-\Phi)^n}{\sqrt 5}$ делает то же самое, что округление до ближайшего целого числа, таким образом вы можете упростить свое решение для поиска $\left\lfloor\frac{\phi^n}{\sqrt 5}+\frac12\right\rfloor$ . Затем используйте биномиальное расширение, чтобы вам приходилось вычислять $\left\lfloor a+b\sqrt 5\right\rfloor$ с соответствующими a и b (которые являются рациональными и могут быть вычислены точно с помощью BigInteger). Если вы все еще вернетесь к Double для этого, вы все равно не получите гораздо больше, чем 1475, но вы должны быть в состоянии выяснить, как сделать эту часть только с точной целочисленной математикой только ☺

$\frac{\phi^n}{\sqrt 5}=\frac{(1+\sqrt 5)^n}{2^n\sqrt 5}=\frac{\sum_{k=0}^n{n\choose k}\sqrt 5^k}{2^n\sqrt 5}$
$=\left(\sum_{k=0}^{\left\lfloor\frac{n-1}2\right\rfloor}\frac{{n\choose 2k+1}5^k}{2^n}\right)+\left(\sum_{k=0}^{\left\lfloor\frac n2\right\rfloor}\frac{{n\choose 2k}5^{k-1}}{2^n}\right)\sqrt 5$

Существует еще один аккуратный метод вычисления чисел Фибоначчи с использованием экспоненциальной матрицы:

$\left(\begin{matrix}1&1\1&0\end{matrix}\right)^n=\left(\begin{matrix}F_n&F_{n-1}\F_{n-1}&F_{n-2}\end{matrix}\right)$

Это можно сделать в O (log n), если вы умны.

Я закончил реализацию этих программ в Haskell. fib1 представляет собой экспонентуцию матрицы, а fib2 - точный целочисленный перевод формулы замкнутой формы, как описано выше. Их соответствующие временные ряды выглядят так, как измерено Criterion при компиляции GHC 7.0.3:

import Control.Arrow
import Data.List
import Data.Ratio

newtype Matrix2 a = Matrix2 (a, a, a, a) deriving (Show, Eq)
instance (Num a) => Num (Matrix2 a) where
    Matrix2 (a, b, c, d) * Matrix2 (e, f, g, h) =
        Matrix2 (a*e+b*g, a*f+b*h, c*e+d*g, c*f+d*h)
    fromInteger x = let y = fromInteger x in Matrix2 (y, 0, 0, y)
fib1 n = let Matrix2 (_, x, _, _) = Matrix2 (1, 1, 1, 0) ^ n in x

binom n =
    scanl (\a (b, c)-> a*b `div` c) 1 $
    takeWhile ((/=) 0 . fst) $ iterate (pred *** succ) (n, 1)
evens (x:_:xs) = x : evens xs
evens xs = xs
odds (_:x:xs) = x : odds xs
odds _ = []
iterate' f x = x : (iterate' f $! f x)
powers b = iterate' (b *) 1
esqrt e n = x where
    (_, x):_ = dropWhile ((<=) e . abs . uncurry (-)) $ zip trials (tail trials)
    trials = iterate (\x -> (x + n / x) / 2) n
fib' n = (a, b) where
    d = 2 ^ n
    a = sum (zipWith (*) (odds $ binom n) (powers 5)) % d
    b = sum (zipWith (*) (evens $ binom n) (powers 5)) % d
fib2 n = numerator r `div` denominator r where
    (a, b) = fib' n
    l = lcm (denominator a) (denominator a)
    r = a + esqrt (1 % max 3 l) (b * b / 5) + 1 % 2

Ответ 2

using System;
using Nat = System.Numerics.BigInteger; // needs a reference to System.Numerics

class Program
{
    static void Main()
    {
        Console.WriteLine(Fibonacci(1000));
    }

    static Nat Fibonacci(Nat n)
    {
        if (n == 0) return 0;
        Nat _, fibonacci = MatrixPower(1, 1, 1, 0, Nat.Abs(n) - 1, out _, out _, out _);
        return n < 0 && n.IsEven ? -fibonacci : fibonacci;
    }

    /// <summary>Calculates matrix power B = A^n of a 2x2 matrix.</summary>
    /// <returns>b11</returns>
    static Nat MatrixPower(
        Nat a11, Nat a12, Nat a21, Nat a22, Nat n,
        out Nat b12, out Nat b21, out Nat b22)
    {
        if (n == 0)
        {
            b12 = b21 = 0; return b22 = 1;
        }

        Nat c12, c21, c22, c11 = MatrixPower(
            a11, a12, a21, a22,
            n.IsEven ? n / 2 : n - 1,
            out c12, out c21, out c22);

        if (n.IsEven)
        {
            a11 = c11; a12 = c12; a21 = c21; a22 = c22;
        }

        b12 = c11 * a12 + c12 * a22;
        b21 = c21 * a11 + c22 * a21;
        b22 = c21 * a12 + c22 * a22;
        return c11 * a11 + c12 * a21;
    }
}

Ответ 3

Тип данных Double имеет верхний предел значения 1.7 x 10 ^ 308

Расчет для 1474 включает в себя значение ~ 1.1 x 10 ^ 308 на одном шаге.

Итак, к 1475 году вы определенно превысите то, что может представлять Double. К сожалению, С# только более простой примитив, Decimal (128-битное число) рассчитан на очень высокую прецессию, но с относительно небольшим диапазоном (всего около 10 ^ 28).

Без создания настраиваемого типа данных, который может обрабатывать числа, превышающие 10 ^ 308 с некоторой степенью десятичной точности, я не вижу способа сделать это. Тем не менее, кто-то, возможно, уже сделал такой класс, как я могу себе представить сценарии, где это может пригодиться.

См. double: http://msdn.microsoft.com/en-us/library/678hzkk9(v=VS.80).aspx

и decimal: http://msdn.microsoft.com/en-us/library/364x0z75(v=VS.80).aspx

Ответ 4

"Библиотека Solver Foundation", как представляется, включает некоторые "большие" типы номеров. Возможно, что его тип Rational может предоставить вам точность и диапазон, который вам нужен. Он представляет собой рациональное соотношение двух значений BigInteger. (Он приносит свой собственный BigInteger - я думаю, он был написан до отправки .NET 4.)

В теории это позволяет ему представлять очень большие числа, но также и представлять большую точность. (Очевидно, что ваша формула не имеет отношения к рациональности, но тогда с плавающей точкой также приближается.)

Он предлагает метод для повышения Rational в силе чего-то еще: http://msdn.microsoft.com/en-us/library/microsoft.solverfoundation.common.rational.power(v=VS.93).aspx

Ответ 5

Как вы правильно указали, для BigInteger не существует метода Sqrt. Вы можете реализовать его самостоятельно, хотя:

Вычислить квадратный корень BigInteger (System.Numerics.BigInteger)

Тем не менее точность должна быть проблемой для вашего кода.

Ответ 6

Многие ответы здесь показывают, что сложность может быть минимизирована до O (log (n)). Почему бы не попробовать целочисленную реализацию подхода log (n)?

Сначала рассмотрим, что вам даны два термина из последовательности Фибоначчи: F(n) и F(n+1). Логично, что более крупные члены F(n+k) могут быть записаны как линейная функция от F(n) и F(n+1) как

 F(n+k) = Ck1*F(n) + Ck2*F(n+1)

Вы можете просто вычислить эти коэффициенты (которые будут зависеть только от k) (интересно, они тоже являются последовательностью Фибоначчи!) и использовать их для ускорения, а затем снова вычислять их для больших значений k иметь возможность продвигаться еще быстрее и т.д.

Ответ 7

Самый быстрый (и самый грязный)?: D

private Double dirty_math_function(Int32 input){
       Double part1 = (1 / Math.Sqrt(5));
       Double part2 = Math.Pow(((1 + Math.Sqrt(5)) / 2), input);
       Double part3 = Math.Pow(((1 - Math.Sqrt(5)) / 2), input);
       return (part1 * part2) - (part1 * part3);
 }

private Double FibSequence(Int32 input) {
  if(input < 1475)
       return dirty_math_function(input);
  else{
       return (FibSequence(input -1) + FibSequence(intput -2));
  }
}

Ответ 8

проблема состоит в том, что (5 ^ (1/2) ^ 1475) легко переполняет int. Что вам нужно сделать, так это написать библиотеку "большой математики", чтобы обрабатывать математику из памяти (бит для бит), а не использовать типы жестких типизированных данных. это боль в но, я знаю. Посмотрите на квадрат и метод умножения.

Ответ 9

Это не точная формула, она даст вам только оценочное значение. И, поскольку арифметика с плавающей запятой ограничена 6-8 байтами на число, отклонение будет увеличиваться с большими числами.
Почему бы не использовать большое целочисленное добавление в цикле, оно должно работать нормально. Гораздо лучше, чем с плавающей точкой.