[Défi math] Calculer rapidement une factorielle

[SPH]

mise en applicationde en Purebasic de la formule de STIRLING

Code : Tout sélectionner

e.f=2.7182818
pi.f=3.1415926536

n=4 ; <<<<<<<<< ceci est le chiffredont on calcul la factoriel

debug  Round(sqr(2*pi*n) * Pow((n/e),n),#PB_Round_Nearest)

c'est étonnant ... ça marche

Essaye avec 200! pour voir

[Ollivier]

Tu as compris mon image SPH en fait?

[SPH]

Tu as compris mon image SPH en fait?

Non, je ne t'ai meme pas lu car j'aime des reponses courte et clair et sans images

[Backup]

Essaye avec 200! pour voir

oui je sais , mais je m'étonnais du fonctionnement de cette fomule sur des nombres digèrable par le PureBasic ...

[SPH]

Essaye avec 200! pour voir

oui je sais , mais je m'étonnais du fonctionnement de cette fomule sur des nombres digèrable par le PureBasic ...

Oui moi aussi. C'est fou le nombre de formule utilisant e et pi.

Sinon, pour ma routine, je calcule avec un for i=1 to 200 la factorielle de 200. Et ma routine met.....0 ms !!!

[Backup]

je suppose que tu divise le nombre de depart (200) par x
pour n'avoir que des petits chiffres a gerer , puis ensuite tu reconstitu non ?

[SPH]

Bon ok, je vais vous donner ma solution mais elle va faire bondir plus d'un :

Code : Tout sélectionner

OpenConsole()
;cmb=500000

PrintN("SPH(2009) - Calcul la factoriel d'un nombre")
PrintN("Exemple : Factoriel de 7 = 1*2*3*4*5*6*7")
re:
PrintN(" ")
Print("Factoriel de : ")
a$=Input()
cmb=Val(a$)
If cmb<1
Goto re
EndIf


;Goooooooooooo 
temps = GetTickCount_() 

;========================================
Dim noeud.w(cmb*3)
noeud(0)=1
ici.l=1
reste.l=0
reste2.l=0
For np=1 To cmb
      If np%10000=0
      af.f=ici/np
      PrintN("Factorielle de : "+Str(np)+" -- Produit en "+Str((GetTickCount_()-temps)/1000)+"s  -- "+StrF(af));+" Reste2="+Str(reste2))
      EndIf
premier.l=np
;;;;;;;;;
        !MOV      dword Esi,[a_noeud]
        !MOV      dword Ecx,[v_ici]

                 ! asm:
                 !XOR      Eax, Eax
                 !MOV      word Ax, [Esi]
                 !MUL      dword[v_premier]
                 !ADD      Eax, [v_reste]
                 ;!add      dword [v_reste2],Edx
                                  
                 !MOV      word [Esi],Ax
                 !SHR      Eax,16
                 !add      Eax,edx
                 !MOV      dword [v_reste],Eax
                 
                 !add      Esi,2
                 !loop      asm
                 
                 !CMP      dword [v_reste],0
                 !JE      fasm
                 !INC      Ecx
                 !INC      dword [v_ici]
                 !JMP   asm
                 ! fasm:
                 
               Next
;========================================

temps = GetTickCount_() - temps 
;temps /1000
;MessageRequester("Millisecondes : ",Str(temps))
PrintN("Millisecondes : "+Str(temps))
;STOOOOOOOOOOP 

;;;;
If CreateFile(0, "d:\Factorielle "+Str(cmb)+"_"+Str(temps)+"ms.bank")
Else
MessageRequester("Erreur", "Fichier impossible à créer...")
End
EndIf

For octet=0 To ici
WriteWord(0,noeud(octet)) 
Next

CloseFile(0)

;;;;
PrintN(" ")
PrintN("--FIN--")

For i=0 To 5
Beep_(i*100,60)
Next
Repeat
  Delay(1000)
ForEver

Et ouai, en asm !!! Car c'est rapide et on fait avec de l'asm ce qu'on fait difficilement en PB. J'etais curieux de voir vos codes pour comparer. Il y avait aussi autre chose qui me tracassais : les balaises en ASM peuvent ils faire encore plus rapide ???

[SPH]

Je viens d'essayer entierement en PB mais impossible !! Et si je reussi, ca fera trop de détours pour trop de perte de temps de calcul...

[Ollivier]

Non, je ne t'ai meme pas lu car j'aime des reponses courte et clair et sans images

Spiche, t'es une flèche! Tu me recales la gueule comme une peau de banane. Si Dobro n'en remet pas une petite couche, j'aurais de la chance!!!

Bon ok, je vais vous donner ma solution mais elle va faire bondir plus d'un

C'est l'hystérie là, il y a même des vagues de olas...

les balaises en ASM

Les « gens qui traitent l'ASM » c'est peut-être plus humble comme expression...

Je viens d'essayer entierement en PB mais impossible !!

Chèr SPH,

C'est en PB qu'il faut commencer. Le problème que tu as posé est intéressant. Mais je ne peux pas faire court, ni enlever les images. A la limite, les changer, mais celles-là sont très bien. Et un tel problème est IMPOSSIBLE à expliquer en une phrase ou deux. Alors s'il-te-plaît, essaie de lire ce qui suit pour m'apporter une réponse sur ce que tu comprends...

Programmer c'est comme faire un voyage. On part d'un endroit (par exemple Paris) avec un listing quasi-vide et on aboutit à un autre endroit (donc dans une autre ville), avec un programme à peu près considéré comme fini!!!

Maintenant, si tu as une idée en tête, reste dessus sans trop d'égarer. Essaie d'évaluer ce que ça coute en temps et en effort pour parvenir à la concrétisation de cette idée. Personne, je dis bien PERSONNE n'est à l'abri de se planter sur cette évaluation.

Et c'est justement pour ça que je parle de villes et de voyages. Comme ça les choses sont plus explicables. Si programmer c'est voyager, tu peux quand même comprendre que terminer un programme correctement peut être arriver dans une ville avec toute sa tête et sa santé, pour considérer qu'exécuter un programme c'est découvrir cette ville.

Plus haut, je donne deux exemples de destination : Marseille (avec Yolande) et Amsterdam (avec tout ce qui se fume et qui casse le crâne).

Le but est d'arriver à Amsterdam et éviter d'avoir le crâne en vrac pour découvrir son architecture et son port, et donc permettre d'autres voyages. Par traduction, le but réel c'est une factorielle calculée rapidement pour se permettre d'autres calculs que les mathématiques offrent à l'infini et avec, parfois beaucoup de plaisir.

Marseille possède aussi un port. On peut aussi y voyager. L'appart de Yolande à Marseille symbolise la formule de Stirling.
Mais... C'est un piège!!! Et il me faut, je crois autant de ligne d'explications qu'ici, sans image ni métaphore pour t'expliquer pourquoi la formule de Stirling est un piège!!! (Si quelqu'un peut mieux faire, je suis preneur)

Woilà! Donc si tu n'y mets pas du tiens, c'est sûr nous n'aboutirons pas: il faut voir les choses petit à petit. Le défi que tu lances est TRES RAID à expliquer, mais possible. C'est très raid à appliquer aussi, mais possible.

En gros, même si "on arrivait à Amsterdam", pas sûr que l'on atteigne le port, le bon bateau et le bon départ pour un voyage de calculs mathématiques exactes et performant (en ASM MMX): il faut encore éviter toutes les tentations de facilité qui y existent et qui noierait ta concentration dans une ivresse contre-productive.

Alors, question: est-ce que tu me comprends un peu mieux là? Et est-ce que je peux aborder le sujet de Yolande, la femme à barbe de Marseille, ou plus proprement « pourquoi la formule de Stirling ne convient pas à la performance voulue dans ce sujet? »

[SPH]

Alors, question: est-ce que tu me comprends un peu mieux là? Et est-ce que je peux aborder le sujet de Yolande, la femme à barbe de Marseille, ou plus proprement « pourquoi la formule de Stirling ne convient pas à la performance voulue dans ce sujet? »

Seul ta derniere phrase est "utile"
Faire des multiplications en allant de 1 a n (pour faire factorielle de n) reste facile puisqu'on parle de quelque chose d'evident : une banale multiplication.
Maintenant, faire Round(sqr(2*pi*n) * Pow((n/e),n) en traitant chaque terme pour qu'il ne fasse pas un overflow est bien plus compliqué car il y a une racine carré (donc, outre la "complexité", c'est un nombre a virgule), il y a PI (que l'on devra prendre avec une precision consequente selon la factorielle cherchée), il y a une puissance et une super mega division avec "e" qui doit etre (outre la complexité) d'une certaine precision.

Je vais expliquer la facilité de ne faire que des multiplications. Comment ca marche (et pourquoi en PB c'est tres compliqué) :

On procede comme une multiplication sur le papier :
ABCDEFG
* X
=======
les lettres de la ligne du haut seront traité individuellement. Ce sont des nombres en base 256 (comme ca, ca tiendra pile dans un octet). Le X lui est un nombre quelconque.
On multiplie X par G puis on ecrit dans l'octet 0 le resultat. Le depassement, on le garde pour le reporter sur F.
On continue en multipliant X pas F puis on ajoute le depassement obtenu plus haut. Et on ecrit dans l'octet 1.
Bref, on fait comme une multiplication papier. On continuera jusque A et jusqu'a ce qu'il n'y ait plus de reste.

[PAPIPP]

Bonjour a Tous
@ Ollivier
OK
+1
@SPH
Bonjour SPH
Ton programme fonctionne bien j’ai pu le vérifier jusqu’à !20 limite en format .q et précision 100%
Et au de là en format ; d Float double précision jusqu’à !170
!170 =7,2574156153079989673967282111293e+306
Comme la vérification est difficile au dessus de !170 par récurrence on peut admettre que le prg Fonctionne bien
J’ai 2 remarques
1) l’instruction Add Eax,edx (récupération en cas de débordement de la multiplication de
!MUL dword[v_premier] ) est inutile car
Le registre Edx sera toujours à zero en effet :
2^16=65536 capacité max de Ax
Avec 2^16 dans ax =>Ax*Ax= 2^32 capacité max de Eax donc il ne peut y avoir de débordement dans Edx
J’ai fait les essais en conservant les fichiers et en les comparant avec des factorielles >1500, les fichiers sont identiques avec l’instruction ou sans l’instruction
2) La structure des données dans le fichier généré sont difficilement exploitables humainement car elles sont définies en base 2^16=65536 par indice et de surcroît les indices faibles contiennent des valeurs de poids faibles et les indices forts des valeurs de poids forts comme en little-endian
L’exemple de la restitution dans notre écriture pose de sérieux Pb à partir de !170
avec comme exemple res.d=res+Pow (65536.0, mot)*ValD (StrU (noeud (mot), #PB_Word))
J’ai des Prg qui fonctionnent en flottant 80 bits permettant de contrôler jusqu’à !1734
mais pour comparer avec ton Prg il faudrait développer sur ta structure binaire.

Code : Tout sélectionner

; OpenConsole()
;cmb=500000
Procedure.s _mh(*ad,lng)
  For i=0 To lng-1
    num.c=PeekC(*ad+i)
    cart$=cart$+Chr(num)
    hexs.s=hexs.s+RSet(Hex(num),2,"0")+" "
  Next
  sort.s=RSet(Hex(*ad),8,"0")+"  "+hexs+"    "+cart$
  ProcedureReturn sort
EndProcedure
re:
a$=InputRequester("Calcul la factoriel d'un nombre"," Factoriel de 4 = 1*2*3*4","")
cmb=Val(a$)
If cmb<1
  Goto re
EndIf

;Goooooooooooo
temps=GetTickCount_()

;========================================
Dim noeud.w(cmb*3)
noeud(0)=1
ici.l=1
reste.l=0
reste2.l=0
valedx.q=0
For np=1 To cmb
  If np%10000=0
    af.f=ici/np
    ;     PrintN("Factorielle de : "+Str(np)+" -- Produit en "+Str((GetTickCount_()-temps)/1000)+"s  -- "+StrF(af));+" Reste2="+Str(reste2))
  EndIf
  premier.l=np
  ;   6         5         4         3         2         1         0
  ;3210987654321098765432109876543210987654321098765432109876543210
  ;|____________________________Rax_______________________________|  quad word les registres MMx/Fpu et les registres SSE de 0 à 7 sont accessibles en 32 bits
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;|--------------Eax-------------|  Double word
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;|______Ax______|   Word
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;|--Ah---|__Al__|
  ; les registres Rax Rbx Rcx Rdx etc.. ont la même structure
  !MOV dword Esi,[a_noeud] ; chargement du premier élément de la table
  !MOV dword Ecx,[v_ici]   ; valeur max d'incrémentation loop décrémente le registre Cx
  
  ! asm:
  !XOR Eax, Eax
  !MOV word Ax, [Esi]
  !MUL dword[v_premier]
  !ADD Eax, [v_reste]
  ;!add      dword [v_reste2],Edx;
  ;!add Eax,edx  ;  Cette instruction est inutile voir remarque ci-dessous sur les débordement de Eax sur Edx
  
  !MOV word [Esi],Ax
  !SHR Eax,16  ; décalage de 16 bits vers la droite ou division par 65536 valeur max d'un word
  
  ;   !mov edx,dword[v_valedx]
  ;   !pusha
  ;   Debug valedx  ; le registre edx sera toujours à zero en effet WORD(signé)*WORD(signé donne au max 32767^2< 2^32 pas de débordement
  ;   ;;                  et word * word  non signé 65536^2 =2^32 toujours pas de débordement de Eax sur Edx
  ;   !popa
  !MOV dword [v_reste],Eax
  
  !add Esi,2
  !loop asm  ; decremente CX chargé avec ici
  
  !CMP dword [v_reste],0
  !JE fasm
  !INC Ecx
  !INC dword [v_ici]
  !JMP asm
  ! fasm:
  
Next
;========================================

temps=GetTickCount_()-temps
;temps /1000
;MessageRequester("Millisecondes : ",Str(temps))
; PrintN("Millisecondes : "+Str(temps))
;STOOOOOOOOOOP

;;;;
If CreateFile(0,"d:\Factorielle "+Str(cmb)+"_"+Str(temps)+"ms.bank")
Else
  MessageRequester("Erreur","Fichier impossible à créer...")
  End
EndIf

For mot=0 To ici-1
  WriteWord(0,noeud(mot))
  res.d=res+Pow(65536.0,mot)*ValD(StrU(noeud(mot),#PB_Word))
Next
; Debug res

CloseFile(0)
Debug "adresse    0                                               1                                               2                                               3                                               4               "
Debug "adresse    0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F  0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F  0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F  0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F  0  1  2  3  4  5"

Debug _mh(@noeud(0), 64)
Debug _mh(@noeud(0)+64,64)
Debug "!"+a$+"="+StrD(res)

A+

[SPH]

1) l’instruction Add Eax,edx (récupération en cas de débordement de la multiplication de
!MUL dword[v_premier] ) est inutile car
Le registre Edx sera toujours à zero en effet :
2^16=65536 capacité max de Ax
Avec 2^16 dans ax =>Ax*Ax= 2^32 capacité max de Eax donc il ne peut y avoir de débordement dans Edx
J’ai fait les essais en conservant les fichiers et en les comparant avec des factorielles >1500, les fichiers sont identiques avec l’instruction ou sans l’instruction

Non, car mon code permet de calculer des factorielles superieurs a 65536. Il permet de faire une factorielle jusqu'a 4294901760

[PAPIPP]

@SPH
OK c'est vrai
Mais quel est l'intérêt de calculer des factorielles > 65536 avec une précision aussi grande
alors que nous ne pouvons même pas comprendre ou appréhender la factorielle !1734
et avec ta structure >!170
La 2em remarque en fait état
A+

[SPH]

@SPH
OK c'est vrai
Mais quel est l'intérêt de calculer des factorielles > 65536 avec une précision aussi grande
alors que nous ne pouvons même pas comprendre ou appréhender la factorielle !1734
et avec ta structure >!170
La 2em remarque en fait état
A+

Avoir un resultat exact pour le plaisir mathematique uniquement

[Backup]

@SPH
OK c'est vrai
Mais quel est l'intérêt de calculer des factorielles > 65536 avec une précision aussi grande
alors que nous ne pouvons même pas comprendre ou appréhender la factorielle !1734
et avec ta structure >!170
La 2em remarque en fait état
A+

ma reflexion serai legerement differente

puisque SPH arrive a faire afficher un nombre tres grand...

pourquoi Fred n'inclue t'il pas un type de variable qui permettrai
l'affichage de tres grand nombre , comme le Basic 1000D sur Atari ST
http://jean-jacques.labarthe.perso.sfr. ... index.html
je suis sur que cela intéresserai les étudiants ...

apres tout si l'on pouvais faire ça sur un 16 bits , de nos jours ce devrai etre un jeu d'enfant ..