enésimo número de fibonacci em tempo sublinear

Existe algum algoritmo para calcular o enésimo número de fibonacci em tempo sub-linear?

O nésimo número de Fibonacci é dado por

 f(n) = Floor(phi^n / sqrt(5) + 1/2) 

Onde

 phi = (1 + sqrt(5)) / 2 

Assumindo que as operações matemáticas primitivas ( + , - , * e / ) são O(1) você pode usar este resultado para calcular o nésimo número de Fibonacci no tempo O(log n) ( O(log n) por causa da exponenciação em a fórmula).

Em c #:

 static double inverseSqrt5 = 1 / Math.Sqrt(5); static double phi = (1 + Math.Sqrt(5)) / 2; /* should use const double inverseSqrt5 = 0.44721359549995793928183473374626 const double phi = 1.6180339887498948482045868343656 */ static int Fibonacci(int n) { return (int)Math.Floor(Math.Pow(phi, n) * inverseSqrt5 + 0.5); } 

A partir da referência de Pillsy à exponenciação de matriz, tal que para a matriz

 M = [1 1] 
     [1 0] 

então

  fib ( n ) = Mn 1,2 

Levantar matrizes para poderes usando multiplicação repetida não é muito eficiente.

Duas abordagens para a exponenciação da matriz são dividir e conquistar, o que resulta em etapas M n em O ( ln n ), ou decomposição de autovalores, que é um tempo constante, mas pode introduzir erros devido à precisão limitada do ponto flutuante.

Se você quer um valor exato maior que a precisão da sua implementação de ponto flutuante, você tem que usar a abordagem O (ln n) baseada nesta relação:

  M n = ( M n / 2 ) 2 se n igual
    = M · M n -1 se n é ímpar

A decomposição de autovalores em M encontra duas matrizes U e Λ tais que Λ é diagonal e

  M = U Λ U -1 
  M n = ( U Λ U -1 ) n
     = U Λ U -1 U Λ U -1 U Λ U -1 ... n vezes
     = U Λ Λ Λ ... U -1 
     = U Λ n U -1 

Aumentar a matriz diagonal Λ até a enésima potência é uma simples questão de elevar cada elemento em Λ até o enésimo, portanto, isso fornece um método O (1) para elevar M à enésima potência. No entanto, os valores em Λ provavelmente não serão inteiros, portanto, algum erro ocorrerá.

Definindo Λ para nossa matriz 2×2 como

 Λ = [λ 1 0]
   = [0 λ 2 ]

Para encontrar cada λ , nós resolvemos

  |  M - λ I |  = 0 

que dá

  |  M - λ I |  = -λ (1 - λ) - 1

 λ² - λ - 1 = 0

usando a fórmula quadrática

 λ = (-b ± √ (b² - 4ac)) / 2a
      = (1 ± 5) / 2
  {λ 1 , λ 2 } = {Φ, 1-Φ} onde Φ = (1 + √5) / 2

Se você leu a resposta de Jason, pode ver para onde isso está indo.

Resolvendo os autovetores X 1 e X 2 :

 se X 1 = [ X 1,1 , X 1,2 ]

  M.  X 1 1 = λ 1 X 1

  X 1,1 + X 1,2 = 1 x 1,1
  X 1,1 = λ 1 X 1,2

 =>
  X 1 = [Φ, 1]
  X 2 = [1-Φ, 1]

Esses vetores dão a U :

 U = [ X 1,1 , X 2,2 ]
     [ X 1,1 , X 2,2 ]

   = [Φ, 1-Φ]
     [1, 1]

Invertendo U usando

 A = [ab]
       [cd]
 =>
 A -1 = (1 / | A |) [d -b]
                    [-ca]

então U- 1 é dado por

 U -1 = (1 / (Φ - (1 - Φ)) [1 Φ-1]
                                [-1 Φ]
 U -1 = (√5) -1 [1 Φ-1]
                [-1 Φ]

Verificação de sanidade:

 UΛU -1 = (√5) -1 [Φ 1-Φ].  [Φ 0]  [1 Φ-1] 
                      [1 1] [0 1-Φ] [-1 Φ]

 seja Ψ = 1-Φ, o outro autovalor

 como Φ é uma raiz de λ²-λ-1 = 0 
 então -ΨΦ = Φ²-Φ = 1
 e Ψ + Φ = 1

 UΛU -1 = (√5) -1 [Φ Ψ].  [Φ 0]  [1 -Ψ] 
                  [1 1] [0 Ψ] [-1 Φ]

        = (√5) -1 [Φ Ψ].  [Φ -ΨΦ] 
                  [1 1] [-Ψ ΨΦ]

        = (√5) -1 [Φ Ψ].  [Φ 1] 
                  [1 1] [-Ψ -1]

        = (√5) -1 [²-² Φ Ψ] 
                   [Φ-0]

        = [Φ + Ψ 1]    
          [1 0]

        = [1 1] 
          [1 0]

        = M 

Então a verificação de sanidade é válida.

Agora temos tudo o que precisamos para calcular M n 1,2 :

 M n = U n U -1
    = (√5) -1 [Φ Ψ].  [ N 0]  [1 -Ψ] 
               [1 1] [0 Ψ n ] [-1 Φ]

    = (√5) -1 [Φ Ψ].  [Φ n- n ] 
               [1 1] [-Ψ n Ψ n Φ]

    = (√5) -1 [Φ Ψ].  [Φ n Φ n -1 ] 
               [1 1] [-Ψ n- Ψ n -1 ] como ΨΦ = -1

    = (√5) -1 [ +1 n +1n +1 Φ n- Ψ n ]
               [Φ n- Ψ n Φ n -1 Ψ n -1 ]

assim

  fib ( n ) = Mn 1,2
         = (Φ n - (1-Φ) n ) / √5

O que concorda com a fórmula dada em outro lugar.

Você pode derivá-lo de uma relação de recorrência, mas na computação de engenharia e simulação calcular os autovalores e autovetores de matrizes grandes é uma atividade importante, pois dá estabilidade e harmônicos de sistemas de equações, bem como permite levantar matrizes para altas potências de forma eficiente.

Se você quer o número exato (que é um “bignum”, ao invés de um int / float), então eu estou com medo de que

É impossível!

Como dito acima, a fórmula para os números de Fibonacci é:

fib n = andar (phi n / √5 + 1/2 )

fib n ~ = phi n / √5

Quantos dígitos é fib n ?

numDigits (fib n) = log (fib n) = log (phi n / √5) = log phi n – log √5 = n * log phi-log √5

numDigits (fib n) = n * const + const

é O ( n )

Como o resultado solicitado é O ( n ), ele não pode ser calculado em menos de O ( n ) tempo.

Se você quer apenas os dígitos mais baixos da resposta, então é possível calcular em tempo sub-linear usando o método de exponenciação de matriz.

Um dos exercícios do SICP é sobre isso, que tem a resposta descrita aqui.

No estilo imperativo, o programa seria algo como

 Função Fib ( contagem )
     um ← 1
     b ← 0
     p ← 0
     q ← 1

     Enquanto contagem > 0 Do
         Se Even ( contar ) Então
              p + 
              q ← 2 pq + 
              contarcontar ÷ 2
         Outro
              abq + aq + ap
              bbp + aq
              contagemcontagem - 1
         Fim se
     End While

     Devolução b
 Função final

Você pode fazer isso exponenciando uma matriz de inteiros também. Se você tem a matriz

  / 1 1 \ M = | | \ 1 0 / 

então (M^n)[1, 2] vai ser igual ao nésimo número de Fibonacci, se [] é um subscrito de matriz e ^ é exponenciação de matriz. Para uma matriz de tamanho fixo, a exponenciação para uma potência integral positiva pode ser feita no tempo O (log n) da mesma maneira que com números reais.

EDIT: Claro, dependendo do tipo de resposta que você quer, você pode ser capaz de fugir com um algoritmo de tempo constante. Como as outras fórmulas mostram, o número n Fibonacci cresce exponencialmente com n . Mesmo com números inteiros não assinados de 64 bits, você precisará apenas de uma tabela de consulta de 94 inputs para cobrir todo o intervalo.

SEGUNDO EDIT: Fazer a matriz exponencial com uma eigendecomposição primeiro é exatamente equivalente à solução de JDunkerly abaixo. Os autovalores dessa matriz são (1 + sqrt(5))/2 e (1 - sqrt(5))/2 .

A Wikipedia tem uma solução de formulário fechado http://en.wikipedia.org/wiki/Fibonacci_number

Ou em c #:

  public static int Fibonacci(int N) { double sqrt5 = Math.Sqrt(5); double phi = (1 + sqrt5) / 2.0; double fn = (Math.Pow(phi, N) - Math.Pow(1 - phi, N)) / sqrt5; return (int)fn; } 

Para realmente grandes, esta function recursiva funciona. Ele usa as seguintes equações:

 F(2n-1) = F(n-1)^2 + F(n)^2 F(2n) = (2*F(n-1) + F(n)) * F(n) 

Você precisa de uma biblioteca que permita trabalhar com inteiros grandes. Eu uso a biblioteca BigInteger de https://mattmccutchen.net/bigint/ .

Comece com uma série de números de fibonacci. Use fibs [0] = 0, fibs [1] = 1, fibs [2] = 1, fibs [3] = 2, fibs [4] = 3, etc. Neste exemplo, eu uso uma matriz dos primeiros 501 (contando 0). Você pode encontrar os primeiros 500 números não-zero de Fibonacci aqui: http://home.hiwaay.net/~jalison/Fib500.html . É preciso um pouco de edição para colocá-lo no formato certo, mas isso não é muito difícil.

Então você pode encontrar qualquer número de Fibonacci usando esta function (em C):

 BigUnsigned GetFib(int numfib) { int n; BigUnsigned x, y, fib; if (numfib < 501) // Just get the Fibonacci number from the fibs array { fib=(stringToBigUnsigned(fibs[numfib])); } else if (numfib%2) // numfib is odd { n=(numfib+1)/2; x=GetFib(n-1); y=GetFib(n); fib=((x*x)+(y*y)); } else // numfib is even { n=numfib/2; x=GetFib(n-1); y=GetFib(n); fib=(((big2*x)+y)*y); } return(fib); } 

Eu testei isso para o número 25.000 Fibonacci e similares.

Aqui está minha versão recursiva que recursa log (n) vezes. Eu acho que é mais fácil de ler na forma recursiva:

 def my_fib(x): if x < 2: return x else: return my_fib_helper(x)[0] def my_fib_helper(x): if x == 1: return (1, 0) if x % 2 == 1: (p,q) = my_fib_helper(x-1) return (p+q,p) else: (p,q) = my_fib_helper(x/2) return (p*p+2*p*q,p*p+q*q) 

Ele funciona porque você pode calcular fib(n),fib(n-1) usando fib(n-1),fib(n-2) se n é ímpar e se n for par, você pode calcular fib(n),fib(n-1) usando fib(n/2),fib(n/2-1) .

O case base e o case ímpar são simples. Para derivar o caso par, comece com a, b, c como valores consecutivos de fibonacci (por exemplo, 8,5,3) e escreva-os em uma matriz, com a = b + c. Aviso prévio:

 [1 1] * [ab] = [a+ba] [1 0] [bc] [ab] 

A partir disso, vemos que uma matriz dos três primeiros números de fibonacci, vezes uma matriz de três números de fibonacci consecutivos, é igual à próxima. Então nós sabemos que:

  n [1 1] = [fib(n+1) fib(n) ] [1 0] [fib(n) fib(n-1)] 

Assim:

  2n 2 [1 1] = [fib(n+1) fib(n) ] [1 0] [fib(n) fib(n-1)] 

Simplificar o lado direito leva ao caso par.

usando R

 l1 < - (1+sqrt(5))/2 l2 <- (1-sqrt(5))/2 P <- matrix(c(0,1,1,0),nrow=2) #permutation matrix S <- matrix(c(l1,1,l2,1),nrow=2) L <- matrix(c(l1,0,0,l2),nrow=2) C <- c(-1/(l2-l1),1/(l2-l1)) k<-20 ; (S %*% L^k %*% C)[2] [1] 6765 

veja dividir e conquistar o algoritmo aqui

O link tem pseudocódigo para a exponenciação da matriz mencionada em algumas das outras respostas para essa questão.

A aritmética de ponto fixo é imprecisa. Código c # de Jason dá resposta incorreta para n = 71 (308061521170130 em vez de 308061521170129) e além.

Para resposta correta, use um sistema de álgebra computacional. O Sympy é uma biblioteca desse tipo para o Python. Há um console interativo no http://live.sympy.org/ . Copie e cole esta function

 phi = (1 + sqrt(5)) / 2 def f(n): return floor(phi**n / sqrt(5) + 1/2) 

Então calcule

 >>> f(10) 55 >>> f(71) 308061521170129 

Você pode querer tentar inspecionar phi .

Você pode usar a equação esquisita quadrada para obter uma resposta exata. A razão é que o $ \ sqrt (5) $ cai no final, você só tem que manter o controle dos coeficientes com o seu próprio formato de multiplicação.

 def rootiply(a1,b1,a2,b2,c): ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b''' return a1*a2 + b1*b2*c, a1*b2 + a2*b1 def rootipower(a,b,c,n): ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format''' ar,br = 1,0 while n != 0: if n%2: ar,br = rootiply(ar,br,a,b,c) a,b = rootiply(a,b,a,b,c) n /= 2 return ar,br def fib(k): ''' the kth fibonacci number''' a1,b1 = rootipower(1,1,5,k) a2,b2 = rootipower(1,-1,5,k) a = a1-a2 b = b1-b2 a,b = rootiply(0,1,a,b,5) # b should be 0! assert b == 0 return a/2**k/5 if __name__ == "__main__": assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3) assert fib(10)==55 

Aqui está um one-liner que calcula F (n), usando inteiros de tamanho O (n), em operações aritméticas de O (log n):

 for i in range(1, 50): print(i, pow(2<  

Usando números inteiros de tamanho O (n) é razoável, pois isso é comparável ao tamanho da resposta.

Para entender isso, seja phi a razão áurea (a maior solução para x ^ 2 = x + 1) e F (n) o n'ésimo número de Fibonacci, onde F (0) = 0, F (1) = F (2) = 1

Agora, phi n = F (n-1) + F (n) phi.

Prova por indução: phi ^ 1 = 0 + 1 * phi = F (0) + F (1) phi. E se phi n = F (n-1) + F (n) phi, então phi ^ (n + 1) = F (n-1) phi + F (n) phi ^ 2 = F (n-1) phi + F (n) (phi + 1) = F (n) + (F (n) + F (n1)) phi = F (n) + F (n + 1) phi. O único passo complicado nesse cálculo é aquele que substitui phi ^ 2 por (1 + phi), o que se segue porque phi é a proporção áurea.

Também números da forma (a + b * phi), onde a, b são inteiros são fechados em multiplicação.

Prova: (p0 + p1 * phi) (q0 + q1 * phi) = p0q0 + (p0q1 + q1p0) phi + p1q1 * phi ^ 2 = p0q0 + (p0q1 + q1p0) phi + p1q1 * (phi + 1) = ( p0q0 + p1q1) + (p0q1 + q1p0 + p1q1) * phi.

Usando esta representação, pode-se calcular phi ^ n em operações inteiras O (log n) usando exponenciação por quadratura. O resultado será F (n-1) + F (n) phi, a partir do qual se pode ler o n-ésimo número de Fibonacci.

 def mul(p, q): return p[0]*q[0]+p[1]*q[1], p[0]*q[1]+p[1]*q[0]+p[1]*q[1] def pow(p, n): r=1,0 while n: if n&1: r=mul(r, p) p=mul(p, p) n=n>>1 return r for i in range(1, 50): print(i, pow((0, 1), i)[1]) 

Observe que a maioria desse código é uma function padrão de exponenciação por esquadramento.

Para chegar ao one-liner que inicia esta resposta, pode-se notar que representar phi por um inteiro X grande o suficiente, pode-se executar (a+b*phi)(c+d*phi) como a operação inteira (a+bX)(c+dX) modulo (X^2-X-1) . Então a function pow pode ser substituída pela function Python pow padrão (que convenientemente inclui um terceiro argumento z que calcula o resultado do módulo z . O X escolhido é 2< .