Suite de Fibonacci et nombres de Lucas

Même pour $n$ (très) grand ?

e.v.

Je fais joujou. Je vais faire calculer $X^n $ modulo $X^2 -(X+1)$, ce qui fournira $F_nX + F_{n-1}$. On commence petit.

[color=#000000]> Q := RationalField() ;        
> QX<X> := PolynomialRing(Q) ;
> F := X^2 - (X+1) ;
> Modexp(X, 7, F) ;  // X^7 modulo F
13*X + 8
> Fibonacci(7), Fibonacci(7-1) ;
13 8
[/color]

Même chose avec $n = 10^4$. Je ne montre pas le résultat $F_n$ mais seulement son nombre de chiffres en base 10

[color=#000000]> time G := Modexp(X, 10^4, F) ;
Time: 0.000
> time Coefficient(G,1) eq Fibonacci(10^4) ;
true
Time: 0.000
> Ilog(10, Z!Coefficient(G,1)) ;            
2089
[/color]

Play again with $n = 10^7$

[color=#000000]> time G := Modexp(X, 10^7, F) ;            
Time: 0.550
> time Coefficient(G,1) eq Fibonacci(10^7) ;
true
Time: 0.340
> Ilog(10, Z!Coefficient(G,1)) ;            
2089876
[/color]

Voici le $n$ du post précédent et son nombre de chiffres décimaux.

[color=#000000]> n ;
123456789
> Ilog(10, n) ;                             
8
[/color]

Est ce que c'est vraiment raisonnable de vouloir faire calculer $F_n$ ? Qui c'est qui paie ?

[color=#000000]> time Fn := Fibonacci(n) ;
Time: 6.010
> time Ilog(10, Fn) ;
25800942
Time: 10.620
[/color]

Tiens, c'est plus rapide de calculer $F_n$ que de calculer le nombre de chiffres décimaux de $F_n$. Rappel : je suis en mode ``joujou'' donc, en principe, je ne peux pas être exclus pour m'être incrusté dans ce fil. Je dis bien en principe.

Un programme Python (pas commenté, on est en shtam, merdre !) qui crache les décimales (dans l'ordre si vous me posez la question) de $F_{2^{27}}$. Le tout dans un temps raisonnable, sachant que ma babasse commence à accumuler les heures de vol.

a=1
b=1
c=0
f=0
l=27
for k in range(l):
    a,b,c=a*a+b*b,b*(a+c),b*b+c*c
print(b)

e.v.

Pour calculer $F_{123456789}$ un code Python. J'ai eu la flemme d'écrire la conversion en binaire.
Je ne l'ai pas essayé.
Je vais le lancer, si vous voyez de la fumée, c'est sûrement chez moi.

L=[1,1,1,0,1,0,1,1,0,1,1,1,1,0,0,1,1,0,1,0,0,0,1,0,1,0,1]
a=1
b=0
c=1
l=len(L)
for k in range(l):
    a,b,c=a*a+b*b,b*(a+c),b*b+c*c
    if L[k]==1:
        a,b,c=a+b,a,b
print(b)

e.v.

ev
Donc le programme retranscrit en C++ cela devrait prendre moins de 2h... et probablement avec un $F_n$ encore plus grand...

[Inutile de reproduire le message précédent. AD]

@ Wilfrid
Effectivement, on peut rogner un peu de temps de calcul en dissociant les deux cas. Mais comme le calcul dans le cas \( \tt L[k]==1 \) est une bête addition, je ne vois pas trop l'intérêt. Maintenant je vais quand même essayer... avec une valeur raisonnable. Je ne vais pas non plus me servir de ma babasse pour faire un barbecue.

@ Kate.
Woaouh, Sage, ça déchire son langage !
Si je savais l'installer, j'abandonnerais Python et ses farces taquines tout de suite.

amicalement,

e.v.

OK, ça gratte quelques pourcents. D'un autre côté comme les programmes calculent plus vite $F_{1023}$ que $F_{450}$, j'avoue ne pas y comprendre grand chose.

import time

def fibo(n):
    L=[]            
    while n>0:
        L.append(n%2)
        n=n//2
    l=len(L)        
    a=1             
    b=0
    c=1
    for k in range(l):
        if L[l-k-1]==0:
            a,b,c=a*a+b*b,b*(a+c),b*b+c*c   
        else:                 
            a,b,c=a*(a+b)+b*(b+c),a*a+b*b,b*(a+c)
    return(b)

tps1 = time.clock()
print(fibo(1023))
tps2 = time.clock()
print(tps2 - tps1)

e.v.

Pour les plus grandes valeurs de $n$ (pourquoi ?), c'est la première version qui gagne. Si on regarde bien, il y a une addition de moins dans le cas \( \tt L[k]==1 \). Je remballe la deuxième version.

@ WIlfrid.
Je ne comprends rien à ce que tu shtames. Les programmes que j'ai écrits donnent-ils la valeur correcte de $F_n$ oui ou non ?

e.v.

Ah bon. Tu peux donner le code du programme que tu as fait tourner pour obtenir 55 ?

Mon programme à moi donne les mêmes chiffres que mathematica. En tout cas les mêmes premiers et les mêmes derniers.

Je ne sais pas ce que je dois en déduire.

e.v.

C'est sûrement toi qui a raison. (Je ne parle pas le mathematica).

Bon Week end,

e.v.

$\def\N{{\lfloor n/2\rfloor}}$Coucou,
Je n'ai pas lu tous les posts de manière attentive (et peut-être que je n'en ai pas envie) mais je trouve que certains passages sentent un peu la bidouille. Je me demande quelle peut-être la réaction d'un ``matheux pur'' (sic) vis-à-vis d'une telle présentation de l'algorithmique. Du coup, je raconte comment on enseignait une méthode de calcul de $F_n$ (niveau L2, module algorithmique dans un cursus mathématique). Je ne donne pas tous les détails (cf une note de cours attachée pour en savoir plus).

Etape 1. Prendre un peu de recul. Pour un monoïde quelconque $(M, \times, 1_M)$, deux exponentiations dichotomiques fournies sous forme d'un schéma récursif. Je désigne par $n = 2\lfloor n/2\rfloor + r$ la division euclidienne de $n$ par $2$ : $\lfloor n/2\rfloor$ en est le quotient et $r = 0,1$ le reste.
$$
x^n = (x^2)^{\N} \times x^r = \cases {(x^2)^{\N} &si $n$ est pair \cr (x^2)^{\N} \times x &si $n$ est impair }
\qquad\qquad
x^n = (x^\N)^2 \times x^r = \cases {(x^\N)^2&si $n$ est pair \cr (x^\N)^2\times x &si $n$ est impair }
$$
avec bien entendu $x^0 = 1_M$. A gauche, c'est l'exponentation dichotomique la plus répandue. Celle de droite (moins répandue) est parfois nommée ``exponentation dichotomique à la Hörner'' car elle résulte en fait de l'évaluation ``à la Hörner'' de l'exposant $n$ via son écriture en base 2.

Etape 2. Mettre sous forme itérative ces deux schémas récursifs. Interviennent alors les chiffres binaires de $n$ qui sont calculés au fur et à mesure (en commençant par les poids faibles à droite et par les poids forts à gauche). Le nombre d'opérations est grosso-modo $\log_2(n)$ (et pas $n$). Pas de détails ici. Cela vaut le coup une fois dans sa vie de prendre un $n$ ni trop petit ni trop grand et de regarder de manière précise ce que fournissent les schémas : quelles sont exactement les puissances calculées ? ...Etc... Je l'ai fait pour $n = 155$, un extrait d'une page manuscrite (sur deux).

Etape 3. Pour le calcul de $F_n$, considérer l'anneau quotient $\Z[X] / \langle X^2 -(X+1) \rangle$ et noter $\phi$ la classe de $X$ dans le quotient, qui vérifie $\phi^2 = \phi + 1$.
$$
\Z[\phi] =_{\rm def} {\Z[X] \over \langle X^2 - (X+1) \rangle} = \Z \oplus \Z\phi
$$
Le point fondamental est l'égalité $\phi^n = F_{n-1} + F_n\phi$ qui se démontre facilement par récurrence sur $n$. Ici, on joue donc la carte : calculer $F_n$ ou $\phi^n$, c'est kif-kif.

Etape 4. Prendre en charge soi-même l'élévation au carré et la multiplication dans $\Z[\phi]$ en encodant $a + b \phi$ par le couple $(a,b)$.

Etape 5. Implémenter le binz dans son langage préféré. J'ai fait cela un certain nombre de fois dans des langages divers. Je viens de le refaire histoire de jouer le jeu. J'ai choisi l'exponentiation dichotomique à la Hörner car elle présente ici un avantage : la multiplication par $\phi$ est une fausse multiplication car $(a+b\phi) \times \phi = b + (a+b)\phi$.

[color=#000000]FibonacciHorner := function(n)
  // Retourne F(n) en calculant Phi^n via une exponentiation dichotomique à la Horner
  if n eq 0 then return 0 ; end if ;
  p := 2^Ilog(2,n) ;    // p puissance de 2   et   n < 2*p
  // (a,b) <--> a.1 + b.Phi
  a := 1 ; b := 0 ;
  m := n ;
  while true do
    // (a,b) <-- (a + b*Phi)^2 = a^2+b^2 + b*(2*a+b)*Phi
    a_saved := a ;
    a := a^2 + b^2 ;
    b := b*(2*a_saved + b) ;
    if m ge p then
      // (a,b) <-- (a + b*Phi)*Phi = b + (a+b)*Phi 
      a_saved := a ;
      a := b ;
      b := a_saved + b ;
      m := m-p ;
    end if ;
    if p eq 1 then break ; end if ;
    p := ExactQuotient(p,2) ;
  end while ;
  return b ;
end function ;
[/color]

Je n'ai pas cherché à finasser. Pour $n = 123456789$, l'exécution ci-dessus nécessite 8 secondes versus 6 secondes pour celle implémentée en magma.

Note (rien à voir) : je crois voir qu'un noble langage a choisi ``len'' comme nom de fonction prédéfinie pour la longueur. No comment.

@ev
Je viens de comprendre (partiellement). Il me semble (mais je suis peut-être le seul à le penser) que tes posts manquaient d'explications. Et je n'avais pas vu trace, dans tes posts précédents, de cette matrice dont tu parles dans le dernier post : ai je la vue basse ?

Oui les anneaux sont canoniquement isomorphes. De manière précise, si $F \in A[X]$ est un polynôme unitaire ($A$ anneau commutatif) de degré $n \ge 1$ et $C \in M_n(A)$ sa matrice compagnon, alors les $A$-algèbres $A[X]/\langle F\rangle$ et $A[C]$ sont canoniquement isomorphes via $X \bmod F \leftrightarrow C$.

Ci-dessous, la matrice compagnon $C$ (à la française) de $F = X^2 - (X+1)$, qui est conjuguée à la tienne

[color=#000000]> C ;
[0 1]
[1 1]
> S ;
[0 1]
[1 0]
> S*C*S ;
[1 1]
[1 0]
> n := 7 ;
> C^n ;
[ 8 13]
[13 21]
> Fibonacci(n-1), Fibonacci(n), Fibonacci(n+1) ;
8 13 21
[/color]

Par ailleurs, je n'ai pas cherché à comprendre complètement ton code (j'ai passé l'âge). Pas compris non plus les temps de calcul (8 heures ?). Une utilisation brute de l'exponentation de $C$ (pas d'optimisation du fait que $C$ est symétrique) n'est pas si mauvaise.

[color=#000000]> n := 123456789 ;
> time F := C^n ;
Time: 24.270
[/color]

Je ne sais pas quelle exponentiation dichotomique est utilisée (la classique probablement).

$\def\vep{\varepsilon}$@Wilfrid
Je pense pouvoir dire que j'ai lu attentivement ton dernier post et que j'en ai apprécié les explications (quelques légers bémols cependant, cf plus loin). J'ai même été jusqu'à implémenter la méthode de manière à trouver les invariants. De manière à pouvoir en fournir une preuve. Car chez nous autres, il faut faire plus que d'expliquer : un algorithme doit être prouvé.

1. Petits bémols en question. Il est important de préciser dans quel ordre sont consommés les chiffres binaires de $n$ : ils le sont dans le sens poids forts vers poids faibles. Certes, cela se voit sur l'exemple $13 = (1101)_2$ et les explications que tu donnes. Mais après tout, dans les programmes, on ne sait pas trop dans quel sens sont rangés les chiffres binaires. En tout cas, en maple ou en magma, ils ne sont pas rangés dans le même sens que Mathematica.
Il y a aussi le coup d'accés à Li en dehors de la boucle qui définit $i$ : dans certains langages, l'indice d'une boucle n'est défini qu'à l'intérieur de la boucle, pas en dehors. On s'en tire en testant la parité de $n$ (ce qui ne coûte rien). Enfin, quand tu dis que les bits de poids fort et faibles ne sont pas utilisés, c'est à modérer : celui de poids faible est utilisé après la boucle pour décider entre $b$ et $c$. Quant à celui de poids fort, il est égal à 1, alors ...

A noter que l'algorithme ne fonctionne que pour $n \ge 2$ ; ce qui n'est pas un handicap car (petit patch) $F_n = n$ pour $n = 0,1$.

Mais ces remarques ne sont que du détail. Ce que je voudrais essayer de faire passer, c'est le point 2. qui suit.

2. Vers des éléments de preuve. Attention, cela va être du lourd, du très lourd même. Je ne sais pas si tu connais la théorie des invariants en algorithmique. Je vais jouer la carte de l'indexation : certes, c'est lourd (bis) mais sûr. D'abord, pour éviter les conflits avec tes notations, je vais noter la décomposition binaire de $n$ de la manière suivante :
$$
n = \vep_0 + \vep_1 2^1 + \vep_2 2^2 + \cdots + \vep_h 2^h, \qquad \vep_i \in \{0,1\}
$$
Je ne range rien ici dans des listes ou je ne sais trop quoi : pour l'instant, je fais des maths, du moins j'essaie.

J'introduis alors (attention, c'est là le lourd) les sommes partielles :
$$
\left \{
\begin {array} {ccl}
s_h &=& \vep_h (= 1) \\
s_{h-1} &=& \vep_{h-1} + 2 \vep_h \\
& \vdots & \\
s_2 &=& \vep_2 + 2\vep_3 + 2^2 \vep_4 + \cdots \\
s_1 &=& \vep_1 + 2\vep_2 + 2^2 \vep_3 + \cdots \\
s_0 &=& \vep_0 + 2\vep_1 + 2^2 \vep_2 + \cdots \\
\end {array}
\right.
$$
Bien observer les indices.Of course $s_0$ c'est $n$. Ce que je dis, c'est que l'on va calculer les triplets successifs :
$$
T_i = (F_{2s_i-1}, F_{2s_i}, F_{2s_i+1}) \quad \text{pour} \quad i \text { descendant de } h \text { à } 1
$$
Note que chaque triplet est de la forme $(F_{2k-1}, F_{2k}, F_{2k+1})$. Et que l'on descend jusqu'à $s_1$ et pas jusqu'à $s_0 = n$ (qui sera pris en charge comme il faut). Cela commence comment pour $i = h$ ? On a dit que $s_h = 1$ donc :
$$
2s_h -1 = 1, \qquad 2s_h = 2, \qquad 2s_h + 1 = 3 \qquad \qquad \hbox {et} \qquad
(F_1, F_2, F_3) = (1,1,2) \text { que tu reconnais}
$$
Enfin, on passe de $T_i$ à $T_{i-1}$ ``comme tu le penses'' en utilisant $s_{i-1} = 2 s_i + \vep_i$ : regarde par exemple $s_1 = 2s_2 + \vep_1$.

3. Si tu savais tout cela, c'est un coup pour rien. De mon côté, cela m'a permis de verrouiller mon programme. En guise de conclusion, une petite pensée pour quelques anciens (dans le désordre) : Hoare, Floyd, Jacques Arsac ...