MPO 95 : Balayage des rationnels

[gege]

Bonjour,
Ok les maths saoûlent pas mal de gens alors ce MPO évacue la question et on va se concentrer sur la programmation...

Donc, la suite suivante balaye tous les nombres rationnels, on le supposera :
U(1)=1
U(n+1)=1/(2*int(U(n))+1-U(n))

Voici le MPO : calculer le plus vite possible U(500).
Quel est ce rationnel ?
Points additionnels à celui qui trouve U(5000)...

A vos claviers !!
G.E.

Sommaire des MPO

[bernouilli92]

Il manque une parenthèse dans la définition de la suite.

[gege]

Bonjour,
Oups ! C'est corrigé, merci
G.E.

[Hobiecat]

Avec une méthode "lourde", ma 41-CV trouve U(500) en un peu moins de 3mn. Mais le résultat est faux, car elle dévie à cause des erreurs d'arrondis dès U(15) : 3.999999, c'est presque 4, mais quand on prend la partie entière, ça fait seulement 3, donc à partir de U(16), c'est faux....

Il doit falloir des bécanes assez précises pour arriver au bon résultat. J'ai mis un coup d'Excel pour comprendre

Edit : pas mieux sur la 42S...

[C.Ret]

[...]Il doit falloir des bécanes assez précises pour arriver au bon résultat. J'ai mis un coup d'Excel pour comprendre

Ou jouer d'astuces :

Bonjour,
Donc, la suite suivante balaye tous les nombres rationnels, on le supposera :

Comme gégé dit de le supposer, je l'ai pris pour une véritée vraie et toute nue.

Donc, ce que gégé désigne par U, je l'ai tout de suite implémenté comme étant A/B. Et comme il faut utiliser une partie entière, c'est bien pratique, en faisant un peu attention A et B sont toujours des entiers !

La question, faut-il déterminer de temps en temps le PGCD(A,B) ou alors A et B sont des enfants sages qui restent entiers et raisonnables pour nos machines ??


Sinon 3 min c'est pas mal, mon SHARP PC-1211 met encore presque 7 min pour m'afficher U(500) (en fait A(500)/B(500) )

[C.Ret]

Voici le MPO : calculer le plus vite possible U(500).

Ah! Oui, j'ai fais un effort, mon code pour SHARP PC-1211 passe de plus de 9' à moins de 12"6 pour afficher la fraction exacte de U(500).
Il fait 115 bytes avec son label DEF-SPC et le BEEP qui permet de bien chronométrer l'arrivée du résultat.

J'obtiens donc la fraction U(5000) en maintenant 17"9.

L'algorithme utilise astucieusement que des choses triviales pour le BASIC du premier pocket; décodage, fractions, calcul infinitésimal, ... chut!

Je vous donne le code en question, mais n'ayant pas remplacé mon écran j'ai un peu de mal avec l'huile noire:

Code : Tout sélectionner

1:" "AREAD N:A=1,B=0,R=▒,D=▒▒▒▒▒ ▒▒▒▒ N▒▒▒▒ ▒
2:R=▒▒R,F=▒
3:▒▒ ▒▒▒ D▒▒ ▒▒▒ ▒N▒D▒=R▒▒▒ N=N▒RD,D=▒▒D,F=F▒▒:▒▒▒▒ 3
4:W=B▒FA,B=A,A=W:▒▒ D▒=R▒▒▒▒ 2
5:BEEP 1:PRINT A,B

Bon j'aurais pu aussi recopier le code que j'ai noté sur mon bloc-notes, mais je ne veux pas influencer vos propres recherches ou maladroitement algo-spoiler. Je crois avoir trouvé le type d'astuce que gégé cherchait à nous faire trouver.

Avec ce code, je trouve U(1081) = 53/37, U(1990)=37/53, U(14)=3/4, U(9)=4/3 et U(148095)=478/65 en moins de 6"9 sur mon SHARP PC-1360 et 21"5 sur mon PC-1211.

J'ose pas sortir mon HP Prime. Ca va déjà trop vite.

[gege]

Bonjour,
Bravo tes résultats sont exacts, mais bon c'est un peu énorme !
L'idée était d'accélérer sur des machines pas trop puissantes en effet...
G.E.

[Marge]

Ça me rappelle un MPO pour lequel j'avais passé plus de temps à le rédiger que C.Ret à en trouver la solution...

[Hobiecat]

J'ai refait un programme pour la 41 avec des entiers A/B comme C.Ret. Ayant eu la flemme de le faire sur la pile, j'ai utilisé des mémoires : on obtient U(500)=11/28 en environ 4'30, mais ça doit pouvoir être largement optimisé.

Au moins le résultat est juste.

[C.Ret]

Avec une méthode "lourde", ma 41-CV trouve U(500) en un peu moins de 3mn.

J'ai moi aussi eu envi d'avoir un méthode lourde sur mon HP-41C. Comme je l'avais indiqué, j'ai utilisé deux valeurs A et B pour exprimer sans arrondi chaque rationnel U().

J'obtiens moi aussi U(500)=11/28 en environ 4'12"9 avec cette version n'utilisant que le registre R00: comme compteur. Les rationnels U(n)=a/b sont mémorisés dans les registres X: / Y: (c'est à dire numérateur en X et dénominateur en Y).
C'est assez pénible et source d'erreur car il faut utiliser la touche [x<>y] pour lire le résultat et donc le risque d'inverser X est assez grand.

Je pense que la solution d'Hobicat utilisant des registres est bien plus fiable et n'est pas significativement plus lente. Recommencer tout le calcul parce que l'on sait plus combien de fois on a pressé la touche [X<>Y] fait perdre bien plus que les misèrables secondes qui diffèrent nos deux codes.

Code : Tout sélectionner

01 LBL "WILF
02  STO 00  1  0
05  LBL 00  RCL X  RCL Z  /  INT  ST+ X  1  +  R^  *  -  CHS  R^  DSE 00  GTO 00
20 END

Petite explication sur le calcul réalisé :

La séquence est définie par
avec
Mon code utilise une définition équivalente en utilisant comme terme initial non pas U(1) mais U(0), ce qui simplifie la façon de coder la boucle principale par LBL 00 ... DSE 00 GTO 00 qui est exécutée une fois pour n=1, le code donne un résultat erroné si l'on entre 0 ce qui ne me gène pas car gégé a eut la bonté de définir sa séquence à partir de U(1).


Le code initialise donc U(0) à zéro, c'est à dire Y:1 X:0 c'est ce que font les instructions 03-04.

En posant , la définition devient:

On vérifie bien que U(1)=1 en effectuant une fois le calcul , par contre, comme l'a très justement souligné Hobicat, les fractions au dénominateur vont nous poser rapidement un souci de précision. Il faut donc calculer et de façon à n'avoir que des entiers.

Pour cela il suffit de transformer un peu l'expression de U(n+1) en factorisant par b(n) au dénominateur :

et l'on trouve par identification :

avec

Le passage du rang n au rang n+1 se fait dans la boucle principale LBL 00 en remontant b(n) dans le registre X: (cf R^ de l'instruction 17-) et devient donc le nouveau a(n+1). Ce qui fait glisser le résultat des calculs vers le registre Y: qui devient effectivement b(n+1).


Voilà, voilà, ça mouline, pour U(5000), il faut 5000 itérations soit près de 37 min.
Je posterai la semaine prochaine la version rapide celle qui calcule U(500) en 8 itérations et U(5000) en 12 itérations.

[Hobiecat]

et l'on trouve par identification :


Le passage du rang n au rang n+1 se fait dans la boucle principale LBL 00 en remontant b(n) dans le registre X: (cf R^ de l'instruction 17-) et devient donc le nouveau a(n+1). Ce qui fait glisser le résultat des calculs vers le registre Y: qui devient effectivement b(n+1)

C'est très exactement ce que j'ai fait. En revanche, j'ai conservé le INT de l'énoncé de Gege, car autant A et B sont entiers, autant le rationnel A/B ne l'est pas.

[C.Ret]

OUi, j'ai pas précisé, mais j'utilise aussi INT (cf. instruction 09-) car .

D'ailleurs, le code n'est pas à l'abri d'un souci si a ou b sont très grands, le rapport peut être mal représenté et une erreur d'arrondi peut mettre à mal tout cela
Ce qui n'arrive pas avec les nombres a et b que l'on manipule, il faudrait atteindre les limites de la machine vers 1E10. On en est bien loin.

MPO96 - HP-41C Calkin Wilf full iterative.gif (8.91 Kio) Vu 8372 fois

[Hobiecat]

OUi, j'ai pas précisé, mais j'utilise aussi INT (cf. instruction 09-) car .

Je comprends mieux.

Bizarrement, ma version sur 41 avec registres donne un résultat juste pour U(500), mais pas pour U(1000)... je cherche où est l'erreur...

[C.Ret]

Je trouve U(1000)=11/39 en juste 9" secondes la première fois, puis instantanément les fois suivantes; la méthode utilisée sur mon HP-28S utilise une technique de mnémonisation.
Et pour la première fois sur cette machine, j'utilise non pas un tableau indexé mais bel et bien une mnémonisation avec un dictionnaire.

L'idée initiale était d'économiser du temps de calcul en ne recommençant pas depuis le début de la séquence les calculs des termes U(n).
En mémorisant les valeurs calculées dans un tableau (une liste ou un vecteur sur l'HP-28S) il est possible de gagner du temps.

Par contre, comme il s'agit de rationnels, il va falloir utiliser deux tableaux, l'un pour les numérateurs A(n) et l'autre pour les dénominateurs B(n).

Si je veux pouvoir aller calculer des U(500) ou U(5000), ces deux tableaux seront d'une taille non négligeable (proportionnellement à la capacité mémoire d'une HP-28S). En plus du souci de la mémoire nécessaire, il y a aussi le temps des éventuels pré-calculs.

Je me suis alors rendu compte de deux choses remarquables concernant le problème que pose gégé :
Observons les 32 premiers termes de cette séquence:

Code : Tout sélectionner

U(n) = a_n / b_n
 n    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
a_n   1  1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3  8  5  7  2  7  5  8  3  7  4  5  1...
b_n   1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3  8  5  7  2  7  5  8  3  7  4  5  1  6...

Comme nous l'avons vu précédemment, a(n+1)=b(n), ce qui fait que les deux séries a_n+1 et b_n sont identiques.

Code : Tout sélectionner

U(n) = a_n / b_n
 n    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
a_n   1  1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3  8  5  7  2  7  5  8  3  7  4  5  1...
        /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /
       /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /  /
b_n   1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3  8  5  7  2  7  5  8  3  7  4  5  1  6...

Cela m'arrange bien pour ma mnémonisation, il suffit de n'enregistrer qu'un seul tableau de valeurs. Voilà une bonne économie.

L'idée est donc d'effectuer la détermination de U(n) à partir de deux termes consécutifs de la séquence fusc : U(n) = fusc(n) / fusc(n+1).
(Evidemment, je ne suis pas le premier à avoir eut cette idée et je reprends ici le nom historique de cette séquence).

Code : Tout sélectionner

U(n) = fusc_n / fusc_n+1
 n     0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
fusc_n 0  1  1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3  8  5  7  2  7  5  8  3  7  4  5  1...

pairs     ^  ^  :  ^  !  :  .  ^     !     :     .     ^           !           :           .           ^
impairs        1+1   1+2   2+1   1+3   3+2   2+3 

Il y a beaucoup de répétitions dans cette séquence, je n'ai peut-être pas besoin de mémoriser chacun des termes.
En effet, par exemple, on voit que les valeurs 1 se répètent, on a fusc(1)=fusc(2)=fusc(4)=fusc(8)=fusc(16)=fusc(32)=...=1

Et il en est de même pour toutes les valeurs d'indices pairs :
fusc(2)=fusc(4)=fusc(8)=fusc(16)=fusc(32)=...=1 = fusc(1)
fusc(6)=fusc(12)=fusc(24)=...= 2 = fusc(3)
fusc(10)=fusc(20)=...= 3 = fusc(5)
fusc(14)=fusc(28)= ...=3 = fusc(7)
etc.

Et concernant les valeurs d'indices impairs, elles s'obtiennent par addition des deux valeurs de part et d'autre du demi-indice:
fusc(3)=fusc(1)+fusc(2)=1+1=2
fusc(5)=fusc(2)+fusc(3)=1+2=3
fusc(7)=fusc(3)+fusc(4)=2+1=3
fusc(9)=fusc(4)+fusc(5)=1+3=4
etc.

Bon, j'ai trouvé ça très très bizarre ! Mais en me documentant, j'ai appris que ces propriétés proviennent de la façon dont tout cela est construit. Car l'énumération des rationnels ce fait sous forme d'arbres.
La séquence que gégé nous incite à étudier, ce fait sous la forme d'un arbre binaire. Parcourir systématiquement cet arbre binaire infini produit une séquence très particulière, il s'agit de la suite diatomique de Stern.
Les propriétés de ces séquences et la forme de l'arbre explique aussi pourquoi je peux par ailleurs calculer U(500) et U(500) en respectivement 8 et 12 itération. Mais je reviendrai sur cela car il y a encore beaucoup à découvrir dont le lien de tout cela avec les fractions continues.

J'ai donc programmé dans mon HP-28S le moyen d'accroitre le tableau fusc(n) au fur et à mesure des calculs selon le principe de la mnémonisation en ne mémorisant que les indices impairs, les indices pairs étant systématiquement divisés par 2 afin d'arriver à un nombre impairs (éventuellement 1).

Je comptais ainsi diviser par deux la place mémoire nécessaire pour calculer U(500), il ne faut que 250 registres d'indices impairs.

Mais en testant, je me suis rendu compte que la détermination, de par exemple U(1000) nécessite beaucoup beaucoup moins de valeurs. Et préparer ou préremplir un aussi grand tableau ne sert à rien.

En effet, pour déterminer U(1000) je n'ai besoin que de fusc(1000) et fusc(1001).

La valeur de fusc(1000) d'indice pair dépend donc du premier indice impair issu de la décomposition 1000→500→250→125 :
fusc(1000)=fusc(500)=fusc(250)=fusc(125)
Or, fusc(125) peut être obtenu à l'aide de fusc(62) et fusc(63): fusc(125)=fusc(62)+fusc(63)
L'indice 62 pair s'atomise lui aussi, la valeur de fusc(62)=fusc(31) et fusc(31)=fusc(15)+fusc(16) d'après le principe des indices impairs.
Et ainsi de suite pour déterminer fusc(1000).
Puis faire de même à partir de fusc(1001).

Si bien qu'au total, la détermination n'a besoin que de 10 termes fusc() mémorisés d'indices { 1001 501 251 125 63 31 15 7 3 1 } dont les valeurs sont calculées de façon récurrente et une seule fois car mémorisées au fur et à mesure dans une liste respectivement { 39 28 17 11 6 5 4 3 2 1 }

fusc(1000)=fusc(125)=fusc(31)+fusc(63)=fusc(15)+fusc(1)+fusc(63)=fusc(7)+2*fusc(1)+fusc(63)=fusc(3)+3*fusc(1)+fusc(63)=5*fusc(1)+fusc(63)
fusc(1000)=5*1+fusc(31)+fusc(1)=5+5+1=11

fusc(1001)=fusc(125)+fusc(501)=11+fusc(250)+fusc(251)=11+11+fusc(125)+fusc(63)=11+11+11+6=39

D'où U(1000)=fusc(1000)/fusc(1001)=11/39 déterminé à l'aide des 10 termes mnémonisés d'indice impairs de la suite diatomique de Stern, soit une vintaine (car la réduction des indices pairs se fait de façon systématique au début de la procédure et ne crées aucun niveau d'appels supplémentaires) d'appels récursifs de Fusc, ce qui est bien plus rapide qu'une itération jusqu'à 1000

Le nombre de termes mnémonisés étant très faible, il est rentable d'utiliser un système de mémorisation plus élaboré qu'un simple tableau (qui prend une bien trop grande plage mémoire pour rester à 98% vide). J'utilise donc deux listes, la première FU.d contient les indices impairs mnémonisés et la seconde FU.v la valeur des fusc() correspondants. Evidemment, l'ordre dans les deux listes correspond exactement.
Elles sont toutes les deux initialisées avec { 1 }. Les nouveaux termes sont ajoutés au début (plus facile à programmer dans cet RPL)

Après le calcul de U(1000), on a FU.d={ 1001 501 251 125 63 31 15 7 3 1 } et FU.v= { 39 28 17 11 6 5 4 3 2 1 }.
Les calculs suivants de U(1000) sont instantanés et ne modifient pas ces deux listes : U(1000)=fusc(125)/fusc(1001)=11/39.
Les calculs d'autres U(n) viennent compléter la mnémonisation et seront d'autant plus rapides lorsque les termes sont déjà disponibles.
Par contre, lorsque les listes s'agrandissent au-delà de 500 termes, la recherche dans le dictionnaire, même réalisée avec l'instruction POS, est de plus en plus lente et ralenti sensiblement mêmes les calculs déjà effectués.
La mémoire de l'HP-28S permet d'avoir des listes de plus de 5000 termes précalculés.

Mais bon, la méthode permet de déterminer sur HP-28S que U(479001600)=1314/13345 en moins de 21 secondes

[gege]

Bonjour,
Bravo, plein d'idées que je ne soupconnais pas !

Le dernier qui quitte le forum peut-il éteindre C.ret en partant ?
Sinon il tourne toute la nuit et le lendemain on a un bourrage papier...
[popcorn mode]Sit back and watch the show[/popcorn mode]
G.E.