[Défi math] Calculer rapidement une factorielle

[PAPIPP]

Bonjour à tous

Voici à titre d’exemple en PB le début d’une librairie réalisée par Kevin Jasik (CSHW89) (Forum Allemand)

La structure de stockage qu’il propose se situe entre la structure purement binaire de SPH donc difficilement compréhensible pour l’homme mais très efficace pour la Machine et une structure en String très facilement compréhensible pour l’homme mais d’une efficacité faible pour la machine.

La grande difficulté de la structure SPH et la conversion en structure compréhensible pour l’homme
Mais par contre toutes les opérations de bases doivent être les plus rapides puisque en code propre à la machine
Si de plus le langage utilisé est ASM alors il doit être imbattable en vitesse

Dans le PRG ci-dessous comme vous pouvez le voir (CSHW89) a découpé le nombre en 10^9= 1 000 000 000 (1 milliard) et placé les chiffres de poids les plus faibles aux indices faibles (little endian) les chiffres de poids forts aux indices forts (little endian)
Il est à la fois proche de la machine tableau en type .l 2^32 mais limité à 10^9 pour facilité la conversion Homme.-Machine
Ce type de structure est dans l’état actuel des possibilités un bon compromis

On peut continuer à développer la librairie.
Toutefois toutes les valeurs seront passées en pointeurs par adresse dans les procédures f(x,y) au lieu de y=f(x) en valeur
Et pour les opérateurs type addition multiplication etc .. f(x,y,z) ou lieu de z=f(x,y)
Voyez mes remarques ; http://www.purebasic.fr/french/viewtopi ... 45&start=0

Code : Tout sélectionner

; Name: BigDecimal
; Author: Kevin Jasik (CSHW89)
; Date: 14.05.2010
; Description:Calcul avec des nombres infiniment grands et exacts


EnableExplicit


#BDMaxInLong = 1000000000
#BDDigitsInLong = 9

#BDExpNull = 0
#BDExpNaN  = 1
#BDExpInf  = 2

Global Dim BDMultiExp.l(#BDDigitsInLong-1)


Enumeration 0
  #BDRoundDown     ; zur 0 runden
  #BDRoundUp       ; weg von der 0 runden
  #BDRoundCeiling  ; zu Infinity runden
  #BDRoundFloor    ; zu -Infinity runden
  #BDRoundHalfDown ; zur nächsten ganzen Zahl runden (bei .5 zur 0 runden)
  #BDRoundHalfUp   ; zur nächsten ganzen Zahl runden (bei .5 weg von der 0 runden)
EndEnumeration

Structure BigDecimal
  sgn.b  ; Signum
  exp.i  ; Zahl multiplizieren mit (10^9)^exp
  size.i ; wieviele Mantissenteile gibt es
  Array man.l(1) ; Mantisse
EndStructure

DeclareDLL BDAdd(*bda.BigDecimal, *bdb.BigDecimal, *result.BigDecimal, digits=#PB_Default)
DeclareDLL BDSub(*bda.BigDecimal, *bdb.BigDecimal, *result.BigDecimal, digits=#PB_Default)
DeclareDLL BDFromQuad(quad.q, *result.BigDecimal)


Global BDRoundMode
Global BDSpecValNegOne.BigDecimal
Global Dim BDSpecValue.BigDecimal(10)


; Initializations-Methode (für Tailbite)
ProcedureDLL BigDecimal_Init()
  Protected i
  BDRoundMode = #BDRoundDown
  BDMultiExp(0) = 1
  For i = 1 To #BDDigitsInLong-1
    BDMultiExp(i) = BDMultiExp(i-1)*10
  Next
  BDFromQuad(-1, BDSpecValNegOne)
  For i = 0 To 10
    BDFromQuad(i, BDSpecValue(i))
  Next
EndProcedure
BigDecimal_Init()


; Verkleinert eine BigDecimal-Zahl auf 'digits'-Nachkomma-
; stellen (rundet dabei) und normalisiert sie
; (d.h. keine Nullen am Anfang und am Ende)
Procedure _BDNorm(*bd.BigDecimal, digits=#PB_Default)
  Protected count, eman, sman, i, stmax, exp, long, rup
  If (*bd\size = 0)
    ; ist 0, Infinity oder NaN
    ProcedureReturn 
  EndIf
  i = 0
  If (digits >= 0)
    ; Zahl wird verkleinert
    ; Alle überflüssigen Teile der Mantisse auf 0 setzen
    exp = *bd\exp - *bd\size
    While ((exp+1)*#BDDigitsInLong < -digits)
      *bd\man(i) = 0
      i + 1
      exp + 1
    Wend
    If (exp*#BDDigitsInLong < -digits)
      ; Alle überflüssigen Stellen ignorieren, außer die die
      ; zum Runden wichtig ist (long%10)
      long = *bd\man(i) / BDMultiExp(#BDDigitsInLong-digits%#BDDigitsInLong-1)
      ; Wenn rup = 1, dann wird aufgerundet
      Select BDRoundMode
      Case #BDRoundDown
      Case #BDRoundUp
        rup = 1
      Case #BDRoundCeiling
        If (*bd\sgn > 0)
          rup = 1
        EndIf
      Case #BDRoundFloor
        If (*bd\sgn < 0)
          rup = 1
        EndIf
      Case #BDRoundHalfDown, #BDRoundHalfUp
        If (long%10 > 5)
          rup = 1
        ElseIf (long%10 = 5)
          If (BDRoundMode = #BDRoundHalfUp)
            rup = 1
          EndIf
        EndIf
      EndSelect
      If (digits%#BDDigitsInLong = 0)
        ; Spezialfall, wenn die Stelle, die zum Runden wichtig ist,
        ; im nächsten Mantissenteil steht
        *bd\man(i) = 0
        i + 1
        If (i = *bd\size)
          If (rup = 0)
            ClearStructure(*bd, BigDecimal)
            *bd\exp = #BDExpNull
            ProcedureReturn 
          EndIf
          ; Für das Aufrunden ist kein Platz, Mantisse muss vergrößert werden
          ReDim *bd\man(i)
          *bd\size + 1
          *bd\exp + 1
        EndIf
        *bd\man(i) + rup
      Else
        ; Die Zahl ohne die überflüssigen Stellen wird wieder
        ; in die Mantisse geschrieben (ggf. wird aufgerundet)
        *bd\man(i) = (long/10 + rup) * BDMultiExp(#BDDigitsInLong-digits%#BDDigitsInLong)
      EndIf
      If (rup = 1)
        ; Durch das Aufrunden ist die Zahl im Mantissenteil zu groß
        While (*bd\man(i) = #BDMaxInLong)
          *bd\man(i) = 0
          i + 1
          If (i = *bd\size)
            ; Für das Aufrunden ist kein Platz, Mantisse muss vergrößert werden
            ReDim *bd\man(i)
            *bd\size + 1
            *bd\exp + 1
          EndIf
          *bd\man(i) + 1
        Wend
      EndIf
    EndIf
  EndIf
  ; Ist der Mantissenteil am Ende gleich 0
  eman = i
  While (*bd\man(eman) = 0)
    eman + 1
    If (eman = *bd\size)
      ; Die gesammte Mantisse ist gleich 0 -> Zahl ist 0
      ClearStructure(*bd, BigDecimal)
      *bd\exp = #BDExpNull
      ProcedureReturn 
    EndIf
  Wend
  ; Ist der Mantissenteil am Anfang gleich 0
  sman = 0
  While (*bd\man(*bd\size-sman-1) = 0)
    sman + 1
  Wend
  If (eman+sman > 0)
    ; Mantisse muss verkleinert werden
    Protected Dim help.l(*bd\size-eman-sman-1)
    For i = eman To *bd\size-sman-1
      help(i-eman) = *bd\man(i)
    Next
    CopyArray(help(), *bd\man())
    *bd\size = *bd\size-eman-sman
    *bd\exp - sman
  EndIf
EndProcedure


; Kopiert eine BigDecimal-Zahl und speichert das Ergebnis
; in '*result', ggf. wird '*result' ungenauer durch Angabe von 'digits'
ProcedureDLL BDCopy(*bd.BigDecimal, *result.BigDecimal, digits=#PB_Default)
  If (*bd <> *result)
    If (*bd\man() <> #Null)
      Dim *result\man(0)
    EndIf
    *result\exp  = *bd\exp
    *result\sgn  = *bd\sgn
    *result\size = *bd\size
    CopyArray(*bd\man(), *result\man())
    ;   CopyStructure(*bd, *result, BigDecimal)
  EndIf
  If (digits >= 0)
    _BDNorm(*result, digits)
  EndIf
  ProcedureReturn *result
EndProcedure

; Negiert eine BigDecimal-Zahl und speichert das Ergebnis
; in '*result', ggf. wird '*result' ungenauer durch Angabe von 'digits'
ProcedureDLL BDNegative(*bd.BigDecimal, *result.BigDecimal, digits=#PB_Default)
  BDCopy(*bd, *result, #PB_Default)
  *result\sgn * (-1)
  If (digits >= 0)
    _BDNorm(*result, digits)
  EndIf
  ProcedureReturn *result
EndProcedure


; Setzt den Modus, wie gerundet wird
ProcedureDLL BDRoundMode(mode)
  If (mode >= 0) And (mode <= #BDRoundHalfUp)
    BDRoundMode = mode
  EndIf
EndProcedure


; Konvertiert ein String in eine BigDecimal-Zahl und
; speichert das Ergebnis in '*result'
ProcedureDLL BDFromString(STR.s, *result.BigDecimal)
  Protected *string.Character, pos, i, cntvman, cntnman, count
  ClearStructure(*result, BigDecimal)
  *result\sgn = 1
  *string = @str
  If (PeekC(*string) = '-')
    ; Zahl ist negativ
    *string + SizeOf(Character)
    *result\sgn = -1
  ElseIf (PeekC(*string) = '+')
    *string + SizeOf(Character)
  EndIf
  If (CompareMemoryString(*string, @"infinity", #PB_String_NoCase, 8) = #PB_String_Equal)
    ; Zahl ist Infinity
    *result\exp = #BDExpInf
    ProcedureReturn *result
  ElseIf (CompareMemoryString(*string, @"nan", #PB_String_NoCase, 8) = #PB_String_Equal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  EndIf
  STR = PeekS(*string)
  *string = @str
  pos = FindString(STR, ".", 0)
  If (pos = 0)
    pos = Len(STR)+1
  EndIf
  ; Anzahl der Mantissenteile hinter/vor dem Punkt
  cntnman = Round((Len(STR)-pos)/#BDDigitsInLong,1)
  cntvman = Round((pos-1)/#BDDigitsInLong,1)
  If (cntnman+cntvman <= 0)
    *result\sgn = 0
    *result\exp = #BDExpNull
    ProcedureReturn *result
  EndIf
  Dim *result\man(cntvman+cntnman-1)
  count = #BDDigitsInLong - (cntnman*#BDDigitsInLong-(Len(STR)-pos))
  For i = 0 To cntnman-1
    ; Mantisse erstellen für Ziffern hinterm Punkt
    *result\man(i) = Val(PeekS(*string+(pos+(cntnman-1-i)*#BDDigitsInLong)*SizeOf(Character), count)) * BDMultiExp(#BDDigitsInLong-count)
    count = #BDDigitsInLong
  Next
  pos - 1
  count = #BDDigitsInLong
  For i = cntnman To cntnman+cntvman-1
    ; Mantisse erstellen für Ziffern vorm Punkt
    pos - #BDDigitsInLong
    If (pos < 0)
      count = #BDDigitsInLong+pos
      pos = 0
    EndIf
    *result\man(i) = Val(PeekS(*string+pos*SizeOf(Character), count))
  Next
  *result\size = cntnman+cntvman
  *result\exp = cntvman
  _BDNorm(*result, #PB_Default)
  ProcedureReturn *result
EndProcedure


; Konvertiert ein Double in eine BigDecimal-Zahl und
; speichert das Ergebnis in '*result'
ProcedureDLL BDFromDouble(double.d, *result.BigDecimal, digits=#PB_Default)
  ClearStructure(*result, BigDecimal)
  If IsNAN(double)
    *result\exp = #BDExpNaN
  ElseIf (double = 0)
    *result\sgn = 0
    *result\exp = #BDExpNull
  Else
    If (digits < 0)
      digits = 15*#BDDigitsInLong-1
    EndIf
    BDFromString(StrD(double, digits+1), *result)
  EndIf
  _BDNorm(*result, digits)
  ProcedureReturn *result
EndProcedure


; Konvertiert ein Long in eine BigDecimal-Zahl und
; speichert das Ergebnis in '*result'
ProcedureDLL BDFromQuad(quad.q, *result.BigDecimal)
  ClearStructure(*result, BigDecimal)
  If (quad < 0)
    quad * (-1)
    *result\sgn = -1
  Else
    *result\sgn = 1
  EndIf
  If (quad = 0)
    *result\sgn = 0
    *result\exp = #BDExpNull
    ProcedureReturn *result
  ElseIf (quad < #BDMaxInLong)
    Dim *result\man(0)
    *result\man(0) = quad
    *result\exp = 1
  ElseIf (quad < #BDMaxInLong*#BDMaxInLong)
    Dim *result\man(1)
    *result\man(0) = quad % #BDMaxInLong
    *result\man(1) = quad / #BDMaxInLong
    *result\exp = 2
  Else
    Dim *result\man(2)
    *result\man(0) =  quad % #BDMaxInLong
    *result\man(1) = (quad / #BDMaxInLong) % #BDMaxInLong
    *result\man(2) =  quad / #BDMaxInLong  / #BDMaxInLong
    *result\exp = 3
  EndIf
  *result\size = *result\exp
  _BDNorm(*result, #PB_Default)
  ProcedureReturn *result
EndProcedure


; Konvertiert die BigDecimal-Zahl in ein String
ProcedureDLL.s BDStr(*bd.BigDecimal, digits=#PB_Default)
  Protected STR.s
  Protected i, cntdig
  If (digits >= 0)
    BDCopy(*bd, *bd, digits)
  EndIf
  If (*bd\sgn = -1)
    STR = "-"
;   ElseIf (*bd\sgn = 1)
;     str = "+"
  EndIf
  If (*bd\size = 0)
    ; Spezielle Werte
    Select *bd\exp
    Case #BDExpNaN
      ProcedureReturn "NaN"
    Case #BDExpInf
      ProcedureReturn STR+"Infinity"
    EndSelect
  EndIf
  cntdig = -1
  If (*bd\exp <= 0)
    ; Zahl ist kleiner 0
    If (digits = 0)
      ProcedureReturn STR+"0"
    EndIf
    cntdig = -*bd\exp*#BDDigitsInLong
    STR + "0."
    If (cntdig <> 0)
      ; Schreibt Nullen, die nicht in der Mantisse gespeichert sind
      STR + RSet("", cntdig, "0")
    EndIf
  EndIf
  i = *bd\size
  While (#True)
    i - 1
    If (digits >= 0) And (cntdig >= digits)
      ; Zu viele Ziffern sind im String -> String wird gekürzt
      STR = PeekS(@str, Len(STR)-(cntdig-digits))
      Break
    EndIf
    If (cntdig = -1)
      If (*bd\size-i-1 = *bd\exp)
        STR + "."
        cntdig = #BDDigitsInLong
      EndIf
    Else
      cntdig + #BDDigitsInLong
    EndIf
    If (i >= 0)
      ; Schreibt die Zahl aus der Mantisse (ggf. mit Nullen vorweg)
      If (i = *bd\size-1) And (*bd\exp > 0)
        STR + Str(*bd\man(i))
      Else
        STR + RSet(Str(*bd\man(i)), #BDDigitsInLong, "0")
      EndIf
    Else
      ; Die Mantisse ist durchlaufen
      If (digits < 0) And (cntdig <> -1)
        ; Wenn 'digits' nicht angegeben wurde, sind wir fertig
        Break
      EndIf
      STR + RSet("", #BDDigitsInLong, "0")
    EndIf
  Wend
  ProcedureReturn STR
EndProcedure


; Vergleicht zwei BigDecimal-Zahlen
; -1: *bda < *bdb, 0: *bda = *bdb, 1: *bda > *bdb
ProcedureDLL BDCompare(*bda.BigDecimal, *bdb.BigDecimal)
  If ((*bda\size = 0) And (*bda\exp = #BDExpNaN)) Or ((*bdb\size = 0) And (*bdb\exp = #BDExpNaN))
    ; Wenn eine Zahl 'NaN' ist, ist das Ergebnis 0
    ProcedureReturn 0
  EndIf
  ; Signum vergleichen
  If (*bda\sgn < *bdb\sgn)
    ProcedureReturn -1
  ElseIf (*bda\sgn > *bdb\sgn)
    ProcedureReturn 1
  ElseIf (*bda\sgn = 0)
    ; Beide Zahlen sind gleich 0
    ProcedureReturn 0
  EndIf
  ; Vergleich mit Infinity
  If (*bda\size = 0) And (*bda\exp = #BDExpInf)
    If (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
      ProcedureReturn 0
    EndIf
    ProcedureReturn 1 * *bda\sgn
  ElseIf (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
    ProcedureReturn (-1) * *bda\sgn
  EndIf
  ; Vergleich zweier normaler Zahlen
  Protected ia, ib, CMP
  CMP = *bda\exp - *bdb\exp
  ; Vergleich der Exponenten
  If (CMP < 0)
    ProcedureReturn (-1) * *bda\sgn
  ElseIf (CMP > 0)
    ProcedureReturn 1 * *bda\sgn
  EndIf
  ia = *bda\size
  ib = *bdb\size
  Repeat
    ia - 1
    ib - 1
    If (ia < 0) And (ib < 0)
      ; Beide Zahlen enden gleichzeitig
      Break
    EndIf
    ; Eine Zahl endet, die andere Zahl ist größer (bei +) bzw. kleiner (bei -)
    If (ia < 0)
      ProcedureReturn (-1) * *bda\sgn
    ElseIf (ib < 0)
      ProcedureReturn 1 * *bda\sgn
    EndIf
    ; Vergleich der Mantissenteile
    CMP = *bda\man(ia) - *bdb\man(ib)
    If (CMP < 0)
      ProcedureReturn (-1) * *bda\sgn
    ElseIf (CMP > 0)
      ProcedureReturn 1 * *bda\sgn
    EndIf
  ForEver
  ProcedureReturn 0
EndProcedure


; Addiert zwei BigDecimal-Zahlen und speichert das Ergebnis
; in '*result'
ProcedureDLL BDAdd(*bda.BigDecimal, *bdb.BigDecimal, *result.BigDecimal, digits=#PB_Default)
  If ((*bda\size = 0) And (*bda\exp = #BDExpNaN)) Or ((*bdb\size = 0) And (*bdb\exp = #BDExpNaN))
    ; Eine Zahl ist 'NaN'
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  ElseIf (*bda\size = 0) And (*bda\exp = #BDExpNull)
    ; '*bda' ist 0, das Ergebnis ist '*bdb'
    ProcedureReturn BDCopy(*bdb, *result, digits)
  ElseIf (*bdb\size = 0) And (*bdb\exp = #BDExpNull)
    ; '*bdb' ist 0, das Ergebnis ist '*bda'
    ProcedureReturn BDCopy(*bda, *result, digits)
  EndIf
  Protected SUB.BigDecimal
  ; Zwei verschiedene Signums -> es wird subrahiert
  If (*bdb\sgn = -1) And (*bda\sgn = 1)
    BDNegative(*bdb, SUB, #PB_Default)
    ProcedureReturn BDSub(*bda, SUB, *result, digits)
  ElseIf (*bda\sgn = -1) And (*bdb\sgn = 1)
    BDNegative(*bda, SUB, #PB_Default)
    ProcedureReturn BDSub(*bdb, SUB, *result, digits)
  EndIf
  ; Addition mit Infinity
  If (*bda\size = 0) And (*bda\exp = #BDExpInf)
    BDCopy(*bda, *result, digits)
    ProcedureReturn *result
  ElseIf (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
    BDCopy(*bdb, *result, digits)
    ProcedureReturn *result
  EndIf
  ; Addition zweier normaler Zahlen
  Protected sman, eman, i, ia, ib, longa, longb, ub, sgn
  sgn = *bda\sgn
  ; Suche den größeren Exponenten
  sman = *bda\exp
  If (sman < *bdb\exp)
    sman = *bdb\exp
  EndIf
  ; Suche den kleineren Exponenten (bzw. Grenze ist 'digits')
  If (digits < 0)
    eman = *bda\exp - *bda\size
    If (eman > *bdb\exp - *bdb\size)
      eman = *bdb\exp - *bdb\size
    EndIf
  Else
    eman = -Round((digits+1)/#BDDigitsInLong, 1)
  EndIf
  Protected Dim help.l(sman-eman)
  i = eman
  ub = 0
  While (i < sman)
    ; Addiere beide Mantissenteile
    longa = 0
    longb = 0
    ia = i + *bda\size - *bda\exp
    ib = i + *bdb\size - *bdb\exp
    If (ia < 0)
    ElseIf (ia < *bda\size)
      longa = *bda\man(ia)
    EndIf
    If (ib < 0)
    ElseIf (ib < *bdb\size)
      longb = *bdb\man(ib)
    EndIf
    help(i-eman) = (longa+longb+ub) % #BDMaxInLong
    ; Übertrag der Addition
    ub = (longa+longb+ub) / #BDMaxInLong
    i + 1
  Wend
  ; Speichere letzten Übertrag
  help(sman-eman) = ub
  ClearStructure(*result, BigDecimal)
  Dim *result\man(sman-eman)
  CopyArray(help(), *result\man())
  *result\exp = sman+1
  *result\sgn = sgn
  *result\size = sman-eman+1
  _BDNorm(*result, digits)
  ProcedureReturn *result
EndProcedure


; Subtrahiert zwei BigDecimal-Zahlen und speichert das Ergebnis
; in '*result'
ProcedureDLL BDSub(*bda.BigDecimal, *bdb.BigDecimal, *result.BigDecimal, digits=#PB_Default)
  If ((*bda\size = 0) And (*bda\exp = #BDExpNaN)) Or ((*bdb\size = 0) And (*bdb\exp = #BDExpNaN))
    ; Eine Zahl ist 'NaN'
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  EndIf
  Protected ADD.BigDecimal
  If (*bda\sgn * *bdb\sgn < 1)
    ; Zahlen haben verscheidene Signums, oder eine Zahl ist 0
    BDNegative(*bdb, ADD, #PB_Default)
    ProcedureReturn BDAdd(*bda, ADD, *result, digits)
  EndIf
  ; Subtraktion mit Infinity
  If (*bda\size = 0) And (*bda\exp = #BDExpInf)
    If (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
      ; Infinity-Infinity = NaN
      ClearStructure(*result, BigDecimal)
      *result\exp = #BDExpNaN
      ProcedureReturn *result
    EndIf
    BDCopy(*bda, *result, digits)
    ProcedureReturn *result
  ElseIf (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
    BDNegative(*bdb, *result, digits)
    ProcedureReturn *result
  EndIf
  ; Subtraktion mit normalen Zahlen
  Protected CMP, sman, eman, i, ia, ib, longa, longb, long, ub
  ; 'cmp' ist Signum der neuen Zahl
  CMP = BDCompare(*bda, *bdb)
  If (CMP = 0)
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNull
    ProcedureReturn *result
  ElseIf (CMP * *bda\sgn < 0)
  ; Vertausche die Zahlen, falls |*bda|<|*bdb|
    Swap *bda, *bdb
  EndIf
  sman = *bda\exp
  ; Suche den kleineren Exponenten (bzw. Grenze ist 'digits')
  If (digits < 0)
    eman = *bda\exp - *bda\size
    If (eman > *bdb\exp - *bdb\size)
      eman = *bdb\exp - *bdb\size
    EndIf
  Else
    eman = -Round((digits+1)/#BDDigitsInLong, 1)
  EndIf
  Protected Dim help.l(sman-eman-1)
  i = eman
  ub = 0
  While (i < sman)
    ; Subtrahiere beide Mantissenteile
    longa = 0
    longb = 0
    ia = i + *bda\size - *bda\exp
    ib = i + *bdb\size - *bdb\exp
    If (ia < 0)
    ElseIf (ia < *bda\size)
      longa = *bda\man(ia)
    EndIf
    If (ib < 0)
    ElseIf (ib < *bdb\size)
      longb = *bdb\man(ib)
    EndIf
    long = longa-longb-ub
    If (long < 0)
      ; Erste Zahl ist kleiner als zweite -> Übertrag
      long + #BDMaxInLong
      ub = 1
    EndIf
    help(i-eman) = long
    i + 1
  Wend
  ClearStructure(*result, BigDecimal)
  Dim *result\man(sman-eman-1)
  CopyArray(help(), *result\man())
  *result\exp = sman
  *result\sgn = CMP
  *result\size = sman-eman
  _BDNorm(*result, digits)
  ProcedureReturn *result
EndProcedure


; Multipliziert zwei BigDecimal-Zahlen und speichert das Ergebnis
; in '*result'
ProcedureDLL BDMul(*bda.BigDecimal, *bdb.BigDecimal, *result.BigDecimal, digits=#PB_Default)
  If ((*bda\size = 0) And (*bda\exp = #BDExpNaN)) Or ((*bdb\size = 0) And (*bdb\exp = #BDExpNaN))
    ; Eine Zahl ist 'NaN'
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  EndIf
  Protected sgn
  sgn = *bda\sgn * *bdb\sgn
  If ((*bda\size = 0) And (*bda\exp = #BDExpInf)) Or ((*bdb\size = 0) And (*bdb\exp = #BDExpInf))
    ; Multiplikation mit Infinity
    ClearStructure(*result, BigDecimal)
    *result\sgn = sgn
    If (*result\sgn = 0)
      ; 0*Infinity = NaN
      *result\exp = #BDExpNaN
    Else
      *result\exp = #BDExpInf
    EndIf
    ProcedureReturn *result
  EndIf
  ; Multiplikation mit normalen Zahlen
  Protected size, exp, stman, i, k, rsi, ub, quad.q
  If (sgn = 0)
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNull
    ProcedureReturn *result
  EndIf
  stman = 0
  ; Neuer Exponent und Länge der Mantisse
  exp  = *bda\exp  + *bdb\exp
  size = *bda\size + *bdb\size
  If (digits >= 0)
    ; Wenn 'digits' angegeben, werden die hinteren Teile nicht berechnet
    i = Round((digits+2)/#BDDigitsInLong, 1)
    If (size > i+exp)
      stman = size-(i+exp)
      size = i+exp
    EndIf
  EndIf
  If (size <= 0)
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNull
    ProcedureReturn *result
  EndIf
  Protected Dim help.l(size-1)
  For i = 0 To *bda\size-1
    ; Gehe erste Zahl durch
    ub = 0
    For k = 0 To *bdb\size-1
      ; Gehe zweite Zahl durch
      rsi = i+k-stman
      If (rsi >= 0)
        ; Multipliziere beide Mantissenteile, addiere dazu
        ; den vorherigen Übertrag und ggf. den alten Wert
        ; im Ergebnis (schriftliches Multiplizieren)
        If (i > 0)
          quad = *bda\man(i) * *bdb\man(k) + ub + help(rsi)
        Else
          quad = *bda\man(i) * *bdb\man(k) + ub
        EndIf
        help(rsi) = quad % #BDMaxInLong
        ub = quad / #BDMaxInLong
      EndIf
    Next
    rsi = i+*bdb\size-stman
    ; Speichere letzten Übertrag
    If (rsi >= 0)
      help(rsi) + ub
    EndIf
  Next
  ClearStructure(*result, BigDecimal)
  Dim *result\man(size-1)
  CopyArray(help(), *result\man())
  *result\exp = exp
  *result\sgn = sgn
  *result\size = size
  _BDNorm(*result, digits)
  ProcedureReturn *result
EndProcedure


; Dividiert zwei BigDecimal-Zahlen und speichert das Ergebnis
; in '*result'
ProcedureDLL BDDiv(*bda.BigDecimal, *bdb.BigDecimal, *result.BigDecimal, digits=#PB_Default)
  If ((*bda\size = 0) And (*bda\exp = #BDExpNaN)) Or ((*bdb\size = 0) And (*bdb\exp = #BDExpNaN))
    ; Eine Zahl ist 'NaN'
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  EndIf
  Protected sgn
  sgn = *bda\sgn * *bdb\sgn
  If (*bda\size = 0) And (*bda\exp = #BDExpInf)
    If (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
      ; Infinity/Infinity = NaN
      ClearStructure(*result, BigDecimal)
      *result\exp = #BDExpNaN
      ProcedureReturn *result
    EndIf
    If (sgn = 0)
      sgn = *bda\sgn
    EndIf
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpInf
    *result\sgn = sgn
    ProcedureReturn *result
  ElseIf (*bdb\size = 0) And (*bdb\exp = #BDExpInf)
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNull
    ProcedureReturn *result
  EndIf
  If (*bdb\size = 0) And (*bdb\exp = #BDExpNull)
    ; Division durch 0
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  ElseIf (*bda\size = 0) And (*bda\exp = #BDExpNull)
    ; Divident ist 0
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNull
    ProcedureReturn *result
  EndIf
  Protected size, exp, i, ia, k, ka, asize
  Protected diva.q, quad.q, quadb.q, ub, qhat, qrem, double.d
  ; Neuer Exponent und Länge der Mantisse
  If (digits < 0)
    digits = ((*bda\size-*bda\exp) + (*bdb\size-*bdb\exp) + 1) * #BDDigitsInLong
  EndIf
  exp  = *bda\exp  - *bdb\exp + 1
  size = exp+Round((digits+2)/#BDDigitsInLong, 1)
  If (size <= 0)
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNull
    ProcedureReturn *result
  EndIf
  ; Divident braucht neue Länge
  asize = size + *bdb\size
  If (asize <= *bda\size)
    asize = *bda\size+1
  EndIf
  Dim help(size-1)
  Protected Dim diva.l(asize-1)
  Protected Dim divb.l(*bdb\size-1)
  ; Divisor normieren (erste Mantisse darf nicht mit einer 0 starten)
  k = 0
  qhat = *bdb\man(*bdb\size-1)
  Repeat
    qhat / 10
    k + 1
  Until (qhat = 0)
  ; Divident speichern (dabei normieren)
  ub = 0
  ia = asize-*bda\size-1
  For i = 0 To *bda\size-1
    quad = *bda\man(i) * BDMultiExp(#BDDigitsInLong-k)
    diva(ia) = quad % #BDMaxInLong + ub
    ub = quad / #BDMaxInLong
    ia + 1
  Next
  diva(ia) = ub
  ; Divisor speichern (dabei normieren)
  ub = 0
  For i = 0 To *bdb\size-1
    quad = *bdb\man(i) * BDMultiExp(#BDDigitsInLong-k)
    divb(i) = quad % #BDMaxInLong + ub
    ub = quad / #BDMaxInLong
  Next
  ; Berechnung starten
  ia = asize-1
  For i = size-1 To 0 Step -1
    If (ia > 0)
      diva = diva(ia) * #BDMaxInLong + diva(ia-1)
    Else
      diva = diva(ia) * #BDMaxInLong
    EndIf
    qhat = diva / divb(*bdb\size-1)
    qrem = diva % divb(*bdb\size-1)
    If (*bdb\size >= 2)
      ; Divisor hat mehr als einen Mantissenteil
      ; Ergebnis der Divisorn korrigieren
      If (ia >= 2)
        quad = qrem * #BDMaxInLong + diva(ia-2)
      Else
        quad = qrem * #BDMaxInLong
      EndIf
      quadb = divb(*bdb\size-2) * qhat
      If (quad < quadb)
        qhat - 1
      EndIf
    EndIf
    If (qhat > 0)
      ; Vom Divident 'qhat*Divisor' subtrahieren
      ka = ia-*bdb\size
      ub = 0
      For k = 0 To *bdb\size-1
        quad = diva(ka) - divb(k) * qhat - ub
        If (quad < 0)
          double = Round((-quad) / #BDMaxInLong, 1)
          ub = double
          quad + ub * #BDMaxInLong
        Else
          ub = 0
        EndIf
        diva(ka) = quad
        ka + 1
      Next
      If (ub > diva(ka))
        ; Ergebnis der Divisorn nochmal korrigieren, da bei
        ; der Subraktion eine negative Zahl rauskam
        diva(ka) - ub
        Repeat
          qhat - 1
          ; Addiere zum Divident einmal den Divisor
          ka = ia-*bdb\size
          ub = 0
          For k = 0 To *bdb\size-1
            diva(ka) + divb(k) + ub
            If (diva(ka) >= #BDMaxInLong)
              diva(ka) - #BDMaxInLong
              ub = 1
            Else: ub = 0
            EndIf
            ka + 1
          Next
          diva(ka) + ub
        Until (diva(ka) >= 0)
      EndIf
      ; Speichere Ergebnis der Division
      diva(ka) = 0
      help(i) = qhat
    EndIf
    ia - 1
  Next
  ClearStructure(*result, BigDecimal)
  Dim *result\man(size-1)
  CopyArray(help(), *result\man())
  *result\exp = exp
  *result\sgn = sgn
  *result\size = size
  _BDNorm(*result, digits)
  ProcedureReturn *result
EndProcedure


; Berechnet das Ergebnis der Exponentialfunktion
; und speichert es in '*result'
ProcedureDLL BDExp(*bd.BigDecimal, *result.BigDecimal, digits)
  Protected.BigDecimal num, den, sum, bdi, DIV
  Protected i
  If (*bd\size = 0)
    Select *bd\exp
    Case #BDExpNaN
      BDCopy(*bd, *result)
    Case #BDExpInf
      If (*bd\sgn = 1)
        ; e^inf = inf
        BDCopy(*bd, *result)
      Else
        ; e^(-inf) = 0
        ClearStructure(*result, BigDecimal)
        *result\exp = #BDExpNull
      EndIf
    Case #BDExpNull
      ; e^0 = 1
      BDCopy(BDSpecValue(1), *result)
    EndSelect
    ProcedureReturn *result
  EndIf
  If (*bd\sgn < 0)
    ; e^(-a) = 1/e^a
    BDNegative(*bd, num)
    BDExp(num, den, digits)
    BDDiv(BDSpecValue(1), den, *result, digits)
    ProcedureReturn *result
  EndIf
  ; e^z = sum(z^k/k!, k=0..infinity)
  BDCopy(*bd, num)
  BDCopy(BDSpecValue(1), den)
  BDCopy(BDSpecValue(1), sum)
  i = 1
  Repeat
    BDDiv(num, den, DIV, digits)
    If (DIV\size = 0) And (DIV\exp = #BDExpNull)
      ; Die Summanten sind zu klein geworden
      Break
    EndIf
    BDAdd(sum, DIV, sum)
    i + 1
    BDMul(num, *bd, num)
    BDFromQuad(i, bdi)
    BDMul(den, bdi, den)
  ForEver
  ProcedureReturn BDCopy(sum, *result, digits)
EndProcedure


; Berechnet das Ergebnis der Logarithmusfunktion
; und speichert es in '*result'
ProcedureDLL BDLog(*bd.BigDecimal, *result.BigDecimal, digits)
  Protected.BigDecimal num, den, qnum, qden, DIV, MUL, sum, sumb, bdi
  Protected long, multi, i, ub, bdsize, CMP
  If ((*bd\size = 0) And (*bd\exp = #BDExpNaN)) Or (*bd\sgn < 0)
    ; Zahl ist NaN oder negativ
    ClearStructure(*result, BigDecimal)
    *result\exp = #BDExpNaN
    ProcedureReturn *result
  EndIf
  If (*bd\size = 0)
    Select *bd\exp
    Case #BDExpInf
      ; ln(inf) = inf
      BDCopy(*bd, *result)
    Case #BDExpNull
      ; log(0) = -inf
      ClearStructure(*result, BigDecimal)
      *result\exp = #BDExpInf
      *result\sgn = -1
    EndSelect
    ProcedureReturn *result
  EndIf
  If (BDCompare(*bd, BDSpecValue(10)) > 0)
    ; Wenn die Zahl größer 10 ist
    ; log(a) = log(a/10^multi) + log(10)*multi
    long = *bd\man(*bd\size-1) / 10
    While (long <> 0)
      multi + 1
      long / 10
    Wend
    bdsize = *bd\size
    Dim num\man(bdsize)
    num\size = bdsize+1
    num\exp = 1
    num\sgn = 1
    For i = bdsize-1 To 0 Step -1
      If (multi = 0)
        num\man(i+1) = *bd\man(i)
      Else
        num\man(i+1) = *bd\man(i) / BDMultiExp(multi) + ub
        ub = (*bd\man(i) % BDMultiExp(multi)) * BDMultiExp(#BDDigitsInLong-multi)
      EndIf
    Next
    num\man(0) = ub
    _BDNorm(num)
    BDLog(num, sum, digits)
    multi + (*bd\exp-1) * #BDDigitsInLong
    BDFromQuad(multi, MUL)
    BDLog(BDSpecValue(10), sumb, digits+multi/2)
    BDMul(MUL, sumb, sumb)
    BDAdd(sum, sumb, *result, digits)
    ProcedureReturn *result
  ElseIf (BDCompare(*bd, BDSpecValue(2)) > 0)
    ; Wenn die Zahl größer 2 ist
    ; log(a) = log(a/2^multi) + log(2)*multi
    long = *bd\man(*bd\size-1)
    If (long >= 8)
      BDDiv(*bd, BDSpecValue(8), num, digits)
      BDLog(num, sum, digits)
      BDLog(BDSpecValue(2), sumb, digits)
      BDMul(BDSpecValue(3), sumb, sumb)
    ElseIf (long >= 4)
      BDDiv(*bd, BDSpecValue(4), num, digits)
      BDLog(num, sum, digits)
      BDLog(BDSpecValue(2), sumb, digits)
      BDMul(BDSpecValue(2), sumb, sumb)
    Else
      BDDiv(*bd, BDSpecValue(2), num, digits)
      BDLog(num, sum, digits)
      BDLog(BDSpecValue(2), sumb, digits)
    EndIf
    BDAdd(sum, sumb, *result, digits)
    ProcedureReturn *result
  EndIf
  CMP = BDCompare(*bd, BDSpecValue(1))
  If (CMP > 0)
    ; Wenn die Zahl größer 1 ist
    ; log(a) = sum((2/(2*k-1))*((a-1)/(a+1))^(2*i-1), i=1..infinity)
    BDCopy(BDSpecValue(0), sum)
    BDSub(*bd, BDSpecValue(1), num)
    BDAdd(*bd, BDSpecValue(1), den)
    BDMul(num, num, qnum)
    BDMul(den, den, qden)
    BDMul(num, BDSpecValue(2), num)
    i = 1
    Repeat
      BDFromQuad(i*2-1,bdi)
      BDDiv(num,den,DIV,digits)
      BDDiv(DIV,bdi,DIV,digits)
      If (DIV\size=0) And (DIV\exp=#BDExpNull)
        ; Die Summanten sind zu klein geworden
        Break
      EndIf
      BDAdd(sum,DIV,sum)
      i+1
      BDMul(num,qnum,num)
      BDMul(den,qden,den)
    ForEver
    ProcedureReturn BDCopy(sum,*result,digits)
  ElseIf (CMP=0)
    ; log(1) = 0
    BDCopy(BDSpecValue(0),*result)
    ProcedureReturn *result
  Else
    ; Wenn die Zahl kleiner 1 ist
    ; log(a) = -log(1/a)
    BDDiv(BDSpecValue(1),*bd,DIV,digits)
    BDLog(DIV,num,digits)
    BDNegative(num,*result,digits)
    ProcedureReturn *result
  EndIf
EndProcedure
Procedure BDPow(*a.BigDecimal,*b.BigDecimal,*c.BigDecimal,d)
  BDLog(*a,*c,d)
  BDMul(*b,*c,*c,d)
  BDExp(*c,*c,d)
EndProcedure
Define.BigDecimal bda,bdb,result

BDFromQuad(2947,bda)
BDFromString("0.000000000000000001",bdb)
BDAdd(bda,bdb,result)
Debug BDStr(result)
Debug BDStr(result,10)

BDFromString("10",bda)
BDFromString("3",bdb)
BDDiv(bda,bdb,result,200)
Debug BDStr(result)

BDFromString("1",bda)
BDExp(bda,result,200)
Debug BDStr(result)

BDFromString("150",bda)
BDLog(bda,result,200)
Debug BDStr(result)

BDFromDouble(-Infinity(),bda)
BDExp(bda,result,20)
Debug BDStr(result)

BDFromString("12.74",bda)
BDFromString("12.75",bdb)
Debug BDStr(bda,1)
Debug BDStr(bdb,1)

BDRoundMode(#BDRoundHalfUp)

BDFromString("12.74",bda)
BDFromString("12.75",bdb)
Debug BDStr(bda,1)
Debug BDStr(bdb,1)
Debug _n(#PB_Default)

; IDE Options = PureBasic 4.50 Beta 4 (Windows - x86)

A+

[Backup]

interressant

[Starwolf20]

Voici mon code pour les factorielles en multiprécision.
NB :Il est possible d'optimiser le code (au détriment de la lisibilité) en réduisant le nombre de multiplications au niveau de la boucle
for i= 2 to n : au lieu de multiplier par 2 puis par 3 puis par 4 ..... on multiplie par 2*3*4.... d'un coup enfin c'est l'idée de départ.....

Il existe pourtant des algorithmes beaucoup plus efficaces qu'une simple boucle de multiplication :
pour ceux que ca intéresses, je conseille le site FastFactorialFunctions
A+

Code : Tout sélectionner

;PB : Quad .q 8 octets -9223372036854775808 à +9223372036854775807 
; --> Max  = 9223372036854775807
; --> MaxB = 1000000000000000000
; --> MaxB = 1.000000000.000000000 

;site avec des grandes factorielles précalculées pour verifier les calculs.
;http://factorielle.free.fr/index_fr.html

Global Base = 1000000000

#BNum_Imax = 10000
Global Dim bn1.q(#BNum_Imax)
Global Dim bn2.q(#BNum_Imax)
Global Dim bn3.q(#BNum_Imax)

Procedure BigMult_MP(Array b1.q(1), Array b2.q(1), Array b3.q(1))
  Protected i, j
  Protected t.q, r.q

;=== Raz du tableau b3 : Resultat de la multiplication.
For i=0 To b1(0) + b2(0) +2
    b3(i)=0
Next

;===Boucle de multiplication en multiprecision.
For j=1 To b2(0)
    r=0
    For i = 1 To b1(0)
      t = b3(i+j-1) + b1(i) * b2(j) + r
      
        b3(i+j-1) = t%Base ; modulo
        r = Int(t/Base)
    Next
    ;=== Traitement de la retenue intermediaire
    If r <> 0 
        b3(b1(0)+j) = r
    EndIf
Next

;=== Calcul de la taille du resultat
b3(0) = b1(0)+b2(0)-1
If r <> 0 
    b3(0) = b3(0)+1
EndIf

EndProcedure

;=============================== Calcul de Factorielle en MultiPrecision
; Starwolf : 2/10/2010

bn1(0)=1
bn1(1)=1

For i =2 To 100
  ; ok pour des factorielle < a la base
    bn2(0)=1
    bn2(1)=i

    BigMult_MP(bn1(), bn2(), bn3())
    CopyArray(bn3(), bn1())
Next

;=== Affichage du resultat
For i=bn3(0) To 1 Step -1
    Debug bn3(i)
Next

[SPH]

Voici mon code pour les factorielles en multiprécision.
NB :Il est possible d'optimiser le code (au détriment de la lisibilité) en réduisant le nombre de multiplications au niveau de la boucle
for i= 2 to n : au lieu de multiplier par 2 puis par 3 puis par 4 ..... on multiplie par 2*3*4.... d'un coup enfin c'est l'idée de départ.....

Il existe pourtant des algorithmes beaucoup plus efficaces qu'une simple boucle de multiplication :
pour ceux que ca intéresses, je conseille le site FastFactorialFunctions
A+

Code : Tout sélectionner

;PB : Quad .q 8 octets -9223372036854775808 à +9223372036854775807 
; --> Max  = 9223372036854775807
; --> MaxB = 1000000000000000000
; --> MaxB = 1.000000000.000000000 

;site avec des grandes factorielles précalculées pour verifier les calculs.
;http://factorielle.free.fr/index_fr.html

Global Base = 1000000000

#BNum_Imax = 10000
Global Dim bn1.q(#BNum_Imax)
Global Dim bn2.q(#BNum_Imax)
Global Dim bn3.q(#BNum_Imax)

Procedure BigMult_MP(Array b1.q(1), Array b2.q(1), Array b3.q(1))
  Protected i, j
  Protected t.q, r.q

;=== Raz du tableau b3 : Resultat de la multiplication.
For i=0 To b1(0) + b2(0) +2
    b3(i)=0
Next

;===Boucle de multiplication en multiprecision.
For j=1 To b2(0)
    r=0
    For i = 1 To b1(0)
      t = b3(i+j-1) + b1(i) * b2(j) + r
      
        b3(i+j-1) = t%Base ; modulo
        r = Int(t/Base)
    Next
    ;=== Traitement de la retenue intermediaire
    If r <> 0 
        b3(b1(0)+j) = r
    EndIf
Next

;=== Calcul de la taille du resultat
b3(0) = b1(0)+b2(0)-1
If r <> 0 
    b3(0) = b3(0)+1
EndIf

EndProcedure

;=============================== Calcul de Factorielle en MultiPrecision
; Starwolf : 2/10/2010

bn1(0)=1
bn1(1)=1

For i =2 To 100
  ; ok pour des factorielle < a la base
    bn2(0)=1
    bn2(1)=i

    BigMult_MP(bn1(), bn2(), bn3())
    CopyArray(bn3(), bn1())
Next

;=== Affichage du resultat
For i=bn3(0) To 1 Step -1
    Debug bn3(i)
Next

Approche interessante et propre. Ca marche bien mais ca met un certain temps a calculer

[Backup]

...

[Frenchy Pilou]

ou plus simplement
922337203622337203692233720369223372036^2
Il sert à tout ce Google
Bon là, pour le cas Dobro il va plier les genoux à partir de ^8 (puissance 8 )

[Backup]

oui c'est pratique la fonction calculette

[Starwolf20]

Finalement,ce "defi" était très intéressant. Comme je l'indiquais dans mon post précédent, il existe des algo plus efficaces qu'une simple boucle de multiplication.
En utilisant la library MPIR qui permet d'effectuer des opérations mathématiques en multiprecision, j'ai pu comparer l'efficacité des algo dont voici les résultats assez renversants :

il s'agit de calculer la factorielle de 1 million qui représente un nombre de 5 565 709 chiffres.
* Programme de SPH (Boucle de multiplication ) t=1220 Sec. (le plus rapide proposé dans ce post. NB : Resultat non verifié)
* Même algo en utilisant la library MPIR t=540 sec
* Algo de factorisation et Split Binaire avec multiplication par MPIR (Explication sur le site Fast Factorial Fonctions) t=3,75 sec.
* J'étais assez content de moi jusqu'au moment ou j'ai essayé la fonction Factorielle de MPIR : t=1,61 sec.

La conclusion est assez évidente mais j'ajoute que je suis (encore!!) resté sur le Q en voyant les performances du PureBasic en facilité d'écriture comme en performance. Définitivement un grand langage de programmation.

[SPH]

il s'agit de calculer la factorielle de 1 million qui représente un nombre de 5 565 709 chiffres.
* Programme de SPH (Boucle de multiplication ) t=1220 Sec. (le plus rapide proposé dans ce post. NB : Resultat non verifié)
* Même algo en utilisant la library MPIR t=540 sec

Allez, je me met au defit de battre la library MPIR. Il y a en effet un truc qui peux accelerer le calcul...


Cette version ne donne pas le bon resultat et je ne sais pas pourqoui mais regardez le principe qui boost le calcul :

Code : Tout sélectionner

OpenConsole()
cmb=1

PrintN("SPH(2009) - Calcul la factorielle d'un nombre (cette version buggue)")
PrintN("Exemple : Factorielle de 7 = 1*2*3*4*5*6*7")
re:
PrintN(" ")
Print("Factorielle 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/2
;       If np%10000=0
;       af.f=ici/np
;       PrintN("Factoriel2le de : "+Str(np)+" -- Produit en "+Str((GetTickCount_()-temps)/1000)+"s  -- "+StrF(af));+" Reste2="+Str(reste2))
;       EndIf

; Debug np
; Debug (cmb+1-np)
; Debug("===")
premier.l=np*(cmb+1-np)
; If premier<0
;   Debug ("negatif")
; EndIf

;;;;;;;;;
        !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]
                                  
                 !MOV      word [Esi],Ax
                 !SHR      Eax,16
                 
                 !add      Eax,edx
                 ;!XOr       eax,eax
                 !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(v2_bug) "+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

[PAPIPP]

Bonjour à tous

Ici 2 tableaux de comparaisons entre GMP et MPIR

http://www.luschny.de/math/factorial/Benchmark.html

et ce n'est pas les algos qui manquent.
http://www.luschny.de/math/factorial/conclusions.html
A+

[Starwolf20]

2 pb, SPH

1) l'astuce ne marche pas pour les factorielle impaires
2) a la fin de la boucle, la multiplication vaut environ (n/2) au carre : dans le cas de 1000000!, ca fait #500 000^2 = 25 000 000 000 000
d'ou overflow

[SPH]

2 pb, SPH

1) l'astuce ne marche pas pour les factorielle impaires
2) a la fin de la boucle, la multiplication vaut environ (n/2) au carre : dans le cas de 1000000!, ca fait #500 000^2 = 25 000 000 000 000
d'ou overflow

1> oui
2> meme de tout petits nombres comme factorielle 1000 (voir meme moins) donne un mauvais resultat et je ne sais pas pourquoi...

[graph100]

Ce topic super interressant sur les plans des maths et programmation m'a fait retourner sur un code que j'avais fait pour calculer Pi et les racines carré et d'autres truc.
Enfin, pour tester des algo de calcul et les comparer entre eux. (après je suis pas doué en optimisation !)

Et comme Dobro a sorti le sujet des nombres en string je poste mon code (900 ligne hélas =)

j'ai abandonné l'idée de calculer 200 000! parceque pour 2000! ca met déjà 200sec xD

Code : Tout sélectionner

If OpenWindow(0, 0, 0, 700, 280, "Calculatrice Décimale", #PB_Window_ScreenCentered | #PB_Window_MinimizeGadget)
    If CreateStatusBar(0, WindowID(0))
        AddStatusBarField(700)
    EndIf
    
    TextGadget(0, 10, 10, 20, 20, "A=")
    TextGadget(1, 10, 40, 20, 20, "B=")
    StringGadget(2, 30, 10, 660, 20, "")
    StringGadget(3, 30, 40, 660, 20, "")
    StringGadget(5, 30, 70, 660, 20, "")
    TextGadget(4, 10, 70, 20, 20, "R=")
    ButtonGadget(6, 30, 100, 150, 30, "R = A + B")
    ButtonGadget(9, 190, 100, 160, 30, "R = A - B")
    ButtonGadget(8, 360, 100, 160, 30, "R = A * B")
    ButtonGadget(7, 530, 100, 160, 30, "R = A / B")
    ButtonGadget(10, 30, 140, 150, 30, "R = A ^ B")
    ButtonGadget(11, 190, 140, 160, 30, "R = Racine de A")
    ButtonGadget(12, 360, 140, 160, 30, "A en facteur 1er")
    ButtonGadget(13, 530, 140, 160, 30, "Comparaison")
    ButtonGadget(14, 30, 180, 150, 30, "Ordre de grandeur")
    ButtonGadget(15, 190, 180, 160, 30, "Pi - Formule de Wallis")
    ButtonGadget(16, 360, 180, 160, 30, "Pi - Formule de Leibnis")
    ButtonGadget(17, 530, 180, 160, 30, "Pi - Formule BBP")
    ButtonGadget(18, 30, 220, 150, 30, "Nombres Premiers")
    ButtonGadget(19, 190, 220, 160, 30, "Factorielle")
    ButtonGadget(20, 360, 220, 160, 30, "")
    ButtonGadget(21, 530, 220, 160, 30, "")
Else
    End
EndIf

Declare.s PasZeroDebFin(str1.s)
Declare.s CompareNumString(str1.s, str2.s)
Declare.s OrdreGrandeurString(str1.s, str2.s)
Declare.s FracString(str1.s)

Declare.s AddString(str1.s, str2.s)
Declare.s SubstracString(str1.s, str2.s)
Declare.s MultipliString(str1.s, str2.s)
Declare.s DivideString(str1.s, str2.s, decimale.l)

Declare.s PowString(str1.s, str2.s)
Declare.s RacineCarreString(str1.s, decimale.l)
Declare.s DecomposeString(str1.s, gad.l)

; fonctions finies

; Comparaison avec décimales et signes
; Addition avec décimales et signes
; Soustraction avec décimales et signes
; Multiplication avec décimales et signes
; Division avec décimales et signes
; Racine carre avec décimales et signes
; Decomposition en nombres premiers de nombre entier

Procedure.s PasZeroDebFin(str1.s)
    ; debut
    signe.s = ""
    If Left(str1, 1) = "-" : signe.s = "-" : str1 = Right(str1, Len(str1) - 1) : EndIf
    
    For a = 1 To Len(str1)
        If Left(str1, 1) = "0" And Left(str1, 2) <> "0."
            str1 = Right(str1, Len(str1) - 1)
            a - 1
        Else
            Break
        EndIf
    Next
    
    If FindString(str1, ".", 0)
        For a = Len(str1) To 1 Step -1
            If Right(str1, 1) = "."
                str1 = Left(str1, a - 1)
                Break
            ElseIf Right(str1, 1) = "0"
                str1 = Left(str1, a - 1)
            Else
                Break
            EndIf
        Next
    EndIf
    
    ProcedureReturn signe + str1
EndProcedure

Procedure.s CompareNumString(str1.s, str2.s)
    str1 = PasZeroDebFin(str1)
    str2 = PasZeroDebFin(str2)
    
    ; signes
    s1.s = Left(str1, 1)
    s2.s = Left(str2, 1)
    
    If s1 = "-" : str1 = Right(str1, Len(str1) - 1) : EndIf
    If s2 = "-" : str2 = Right(str2, Len(str2) - 1) : EndIf
    
    If s1 = "-" And s2 <> "-"
        ProcedureReturn "<"
    ElseIf s1 <> "-" And s2 = "-"
        ProcedureReturn ">"
    ElseIf s1 = "-" And s2 = "-"
        str.s = str1
        str1 = str2
        str2 = str
    EndIf
    
    ; virgule
    If str1 = str2
        final.s = "="
    Else
        vir1 = FindString(str1, ".", 0)
        vir2 = FindString(str2, ".", 0)
        
        If vir1 Or vir2
            If vir1 = 0 : vir1 = Len(str1) : EndIf
            If vir2 = 0 : vir2 = Len(str2) : EndIf
            
            vir1 = Len(str1) - vir1
            vir2 = Len(str2) - vir2
            
            vir = vir1
            If vir1 < vir2 : vir = vir2 : EndIf
            
            str1 = RemoveString(str1, ".") + ReplaceString(Space(vir - vir1), " ", "0")
            str2 = RemoveString(str2, ".") + ReplaceString(Space(vir - vir2), " ", "0")
        EndIf
        
        If Len(str1) < Len(str2)
            final.s = "<"
        ElseIf Len(str1) > Len(str2)
            final.s = ">"
        Else
            For a = 1 To Len(str1)
                char1 = Val(Mid(str1, a, 1))
                char2 = Val(Mid(str2, a, 1))
                If char1 > char2
                    final.s = ">"
                    Break
                ElseIf char1 < char2
                    final.s = "<"
                    Break
                EndIf
            Next
        EndIf
    EndIf
    
    ProcedureReturn final
EndProcedure

Procedure.s OrdreGrandeurString(str1.s, str2.s); donne un ordre de grandeur de sorte que : str1 * ordre < str2
    final.s = ""
    
    str1 = PasZeroDebFin(str1)
    str2 = PasZeroDebFin(str2)
    
    ; signes
    s1.s = Left(str1, 1)
    s2.s = Left(str2, 1)
    
    If s1 = "-" : str1 = Right(str1, Len(str1) - 1) : EndIf
    If s2 = "-" : str2 = Right(str2, Len(str2) - 1) : EndIf
    
    signe.s = ""
    
    If (s1 = "-" Or s2 = "-") And Not (s1 ="-" And s2 = "-")
        signe = "-"
    EndIf
    
    ; virgule
    If str1 = str2
        final.s = "1"
    Else
        vir1 = FindString(str1, ".", 0)
        vir2 = FindString(str2, ".", 0)
        
        If vir1 Or vir2
            If vir1 = 0 : vir1 = Len(str1) : EndIf
            If vir2 = 0 : vir2 = Len(str2) : EndIf
            
            vir1 = Len(str1) - vir1
            vir2 = Len(str2) - vir2
            
            vir = vir1
            If vir1 < vir2 : vir = vir2 : EndIf
            
            str1 = RemoveString(str1, ".") + ReplaceString(Space(vir - vir1), " ", "0")
            str2 = RemoveString(str2, ".") + ReplaceString(Space(vir - vir2), " ", "0")
            
            str1 = PasZeroDebFin(str1)
            str2 = PasZeroDebFin(str2)
        EndIf
        
        dif = Len(str2) - Len(str1)
        
        nb.s = ReplaceString(Space(Abs(dif)), " ", "0")
        
        If dif < 0 : nb = "0." + Left(nb, -dif - 1) + "1"
        Else : nb = "1" + nb
        EndIf
        
        resultd.s = PasZeroDebFin(MultipliString(nb, str1))
        
        If ComparenumString(resultd, str2) = ">"
            If dif < 0
                dif + 1
                nb.s = ReplaceString(Space(Abs(dif)), " ", "0")
                
                If dif < 0 : nb = "0." + Left(nb, -dif - 1) + "1"
                Else : nb = "1" + nb
                EndIf
            Else
                dif - 1
                nb.s = ReplaceString(Space(Abs(dif)), " ", "0")
                
                If dif < 0 : nb = "0." + Left(nb, -dif - 1) + "1"
                Else : nb = "1" + nb
                EndIf
            EndIf
        EndIf
        final = nb
    EndIf
    
    final = signe + final
    
    ProcedureReturn final.s
EndProcedure

Procedure.s FracString(str1.s)
    str1 = PasZeroDebFin(str1)
    pos = FindString(str1, ".", 0)
    If pos = 0 : ProcedureReturn "0" : EndIf
    
    ProcedureReturn "0" + Mid(str1, pos, Len(str1))
EndProcedure

Procedure.s IntString(str1.s)
    str1 = PasZeroDebFin(str1)
    pos = FindString(str1, ".", 0)
    If pos = 0 : ProcedureReturn str1 : EndIf
    
    ProcedureReturn Left(str1, pos - 1)
EndProcedure

Procedure.s AddString(str1.s, str2.s)
    If Left(Str1, 1) <> "-" And Left(str2, 1) <> "-"
        
        ; traitement de la décimale on ramène le problème sur un terrain sûr !
        vir1 = FindString(str1, ".", 0)
        vir2 = FindString(str2, ".", 0)
        
        If vir1 = 0 And vir2 = 0
            vir = 0
        Else
            If vir1 = 0 : vir1 = Len(str1) : EndIf
            If vir2 = 0 : vir2 = Len(str2) : EndIf
            
            vir1 = Len(str1) - vir1
            vir2 = Len(str2) - vir2
            
            vir = vir1
            If vir1 < vir2 : vir = vir2 : EndIf
            
            str1 = RemoveString(str1, ".") + ReplaceString(Space(vir - vir1), " ", "0")
            str2 = RemoveString(str2, ".") + ReplaceString(Space(vir - vir2), " ", "0")
        EndIf
        
        ; addition
        If Len(str1) < Len(str2) : str.s = str1 : str1 = str2 : str2 = str : EndIf
        
        len1 = Len(str1)
        len2 = Len(str2)
        ret = 0
        final.s = ""
        
        For pos = 0 To len1 - 1
            If pos < len2
                resultat.s = Str(Val(Mid(str1, len1 - pos, 1)) + Val(Mid(str2, len2 - pos, 1)))
            Else
                resultat.s = Mid(str1, len1 - pos, 1)
            EndIf
            
            resultat = Str(Val(resultat) + ret)
            
            If Len(resultat) = 2
                final = Right(resultat, 1) + final
                ret = Val(Left(resultat, 1))
            Else
                ret = 0
                final = resultat + final
            EndIf
        Next
        
        If ret
            final = Str(ret) + final
        EndIf
        
        If vir
            final = Left(final, Len(final) - vir) + "." + Right(final, vir)
        EndIf
        
    Else
        If Left(str1, 1) = "-"
            str = Right(str1, Len(str1) - 1)
            str1 = str2
            str2 = str
            If Left(str1, 1) = "-"
                final = "-" + AddString(Right(str1, Len(str1) - 1), str2)
            Else
                final = SubstracString(str1, str2)
            EndIf
        Else
            str2 = Right(str2, Len(str2) - 1)
            final = SubstracString(str1, str2)
        EndIf
    EndIf
    
    ProcedureReturn final
EndProcedure

Procedure.s SubstracString(str1.s, str2.s)
    If Left(Str1, 1) <> "-" And Left(str2, 1) <> "-"
        
        ; traitement de la décimale on ramène le problème sur un terrain sûr !
        vir1 = FindString(str1, ".", 0)
        vir2 = FindString(str2, ".", 0)
        
        If vir1 = 0 And vir2 = 0
            vir = 0
        Else
            If vir1 = 0 : vir1 = Len(str1) : EndIf
            If vir2 = 0 : vir2 = Len(str2) : EndIf
            
            vir1 = Len(str1) - vir1
            vir2 = Len(str2) - vir2
            
            vir = vir1
            If vir1 < vir2 : vir = vir2 : EndIf
            
            str1 = RemoveString(str1, ".") + ReplaceString(Space(vir - vir1), " ", "0")
            str2 = RemoveString(str2, ".") + ReplaceString(Space(vir - vir2), " ", "0")
        EndIf
        
        ; signe
        signe = 1
        If Len(str2) > Len(str1) Or (Len(Str2) = Len(str1) And Val(Left(Str2, 1)) > Val(Left(str1, 1)))
            signe = -1
            str.s = str1
            str1 = str2
            str2 = str
        EndIf
        
        ; soustraction
        If Len(str1) < Len(str2) : str.s = str1 : str1 = str2 : str2 = str : EndIf
        
        len1 = Len(str1)
        len2 = Len(str2)
        ret = 0
        final.s = ""
        
        sign = 0
        
        For pos = 0 To len1 - 1
            If pos < len2
                resultat.s = Str(Val(Mid(str1, len1 - pos, 1)) - Val(Mid(str2, len2 - pos, 1)) + sign)
            Else
                resultat.s = Str(Val(Mid(str1, len1 - pos, 1)) + sign)
            EndIf
            
            If Left(resultat, 1) = "-"
                sign = -1
                resultat = Str(10 + Val(resultat))
            Else
                sign = 0
            EndIf
            
            final = resultat + final
        Next
        
        If signe = -1
            final = "-" + final
        EndIf
        
        If vir
            final = Left(final, Len(final) - vir) + "." + Right(final, vir)
        EndIf
    Else
        If Left(str1, 1) = "-"
            If Left(str2, 1) = "-"
                final = AddString(str1, Right(str2, Len(str2) - 1))
            Else
                final = "-" + AddString(Right(str1, Len(str1) - 1), str2)
            EndIf
        Else
            str2 = Right(str2, Len(str2) - 1)
            final = AddString(str1, str2)
        EndIf
    EndIf
    
    ProcedureReturn final
EndProcedure

Procedure.s MultipliString(str1.s, str2.s)
    str1 = PasZeroDebFin(str1)
    str2 = PasZeroDebFin(str2)
    
    If str1 = "" : str1 = "0" : EndIf
    If str2 = "" : str2 = "0" : EndIf
    
    ; traitement de la décimale on ramène le problème sur un terrain sûr !
    vir1 = FindString(str1, ".", 0)
    vir2 = FindString(str2, ".", 0)
    
    If vir1 : vir1 = Len(str1) - vir1 : EndIf
    If vir2 : vir2 = Len(str2) - vir2 : EndIf
    
    str1 = RemoveString(str1, ".")
    str2 = RemoveString(str2, ".")
    
    ; traitement du signe
    signe$ = "+"
    If Left(str1, 1) = "-" : signe$ = "-" : EndIf
    If Left(str2, 1) = "-" : signe$ = "-" : EndIf
    If Left(str1, 1) = "-" And Left(str2, 1) = "-" : signe$ = "+" : EndIf
    
    ; muliplication
    
    If Len(str2) > Len(str1)
        str.s = str1
        str1 = str2
        str2 = str
    EndIf
    
    len1 = Len(str1)
    len2 = Len(str2)
    
    Dim Mul.s(len2 - 1)
    
    For pos2 = 0 To len2 - 1
        char = Val(Mid(str2, len2 - pos2, 1))
        ret = 0
        For pos = 0 To len1 - 1
            resultat.s = Str(Val(Mid(str1, len1 - pos, 1)) * char + ret)
            
            If Len(resultat) = 2 And pos <> len1 - 1
                Mul(pos2) = Right(resultat, 1) + Mul(pos2)
                ret = Val(Left(resultat, 1))
            Else
                ret = 0
                Mul(pos2) = resultat + Mul(pos2)
            EndIf
        Next
    Next
    
    
    final.s = ""
    
    For pos = 0 To len2 - 1
        final = AddString(final, Mul(pos) + ReplaceString(Space(pos), " ", "0"))
    Next
    
    ; Finalisations
    If vir1 + vir2 > 0
        final = Left(final, Len(final) - (vir1 + vir2)) + "." + Right(final, vir1 + vir2)
    EndIf
    
    If signe$ = "-" : final = "-" + final : EndIf
    
    ProcedureReturn final
EndProcedure

Procedure.s DivideString(str1.s, str2.s, decimale.l)
    str1 = PasZeroDebFin(str1)
    str2 = PasZeroDebFin(str2)
    
    ; traitement du signe
    signe$ = "+"
    If Left(str1, 1) = "-" : signe$ = "-" : EndIf
    If Left(str2, 1) = "-" : signe$ = "-" : EndIf
    If Left(str1, 1) = "-" And Left(str2, 1) = "-" : signe$ = "+" : EndIf
    
    ; virgule
    vir1 = FindString(str1, ".", 0)
    vir2 = FindString(str2, ".", 0)
    
    If vir1 Or vir2
        If vir1 = 0 : vir1 = Len(str1) : EndIf
        If vir2 = 0 : vir2 = Len(str2) : EndIf
        
        vir1 = Len(str1) - vir1
        vir2 = Len(str2) - vir2
        
        vir = vir1
        If vir1 < vir2 : vir = vir2 : EndIf
        
        str1 = RemoveString(str1, ".") + ReplaceString(Space(vir - vir1), " ", "0")
        str2 = RemoveString(str2, ".") + ReplaceString(Space(vir - vir2), " ", "0")
    EndIf
    
    ;Division
    final.s
    
    des.s = OrdreGrandeurString(str2, str1)
    
    res.s = des
    
    nb_decimal = FindString(res, ".", 0)
    If nb_decimal = 0 : nb_decimal = Len(res) : EndIf
    nb_decimal = Len(res) - nb_decimal
    
    Repeat
        mul.s = MultipliString(res, str2)
        comp.s = CompareNumString(mul, str1)
        
        If comp = "="
            Break
        ElseIf comp = ">"
            res = SubstracString(res, des)
            des = PasZeroDebFin(MultipliString(des, "0.1"))
            If Mid(des, 2, 1) = "." : nb_decimal + 1 : EndIf
        EndIf
        res = AddString(res, des)
        If nb_decimal > decimale : res = Left(res, Len(res) - 1) : EndIf
    Until nb_decimal > decimale
    
    final = res
    
    ; Finalisations
    If signe$ = "-" : final = "-" + final : EndIf
    
    ProcedureReturn final
EndProcedure

Procedure.s PowString(str1.s, str2.s)
    final.s = "1"
    
    If Val(str2) = 0 : ProcedureReturn final : EndIf
    
    For a = 1 To Val(str2)
        final = MultipliString(final, str1)
    Next
    
    ProcedureReturn final
EndProcedure

Procedure.s RacineCarreString(str1.s, decimale.l)
    str1 = PasZeroDebFin(str1)
    
    ; traitement du signe
    If Left(str1, 1) = "-"
        ProcedureReturn "Erreur: nombre négatif sous la racine"
    EndIf
    
    ;Division
    final.s
    
    nbdec = Len(str1) / 2 - 1
    If nbdec < 0 Or Left(str1, 2) = "0." : nbdec = 0 : EndIf
    
    des.s = "1" + ReplaceString(Space(nbdec), " ", "0")
    
    res.s = des
    
    nb_decimal = FindString(res, ".", 0)
    If nb_decimal = 0 : nb_decimal = Len(res) : EndIf
    nb_decimal = Len(res) - nb_decimal
    
    Repeat
        mul.s = MultipliString(res, res)
        
        comp.s = CompareNumString(mul, str1)
        
        If comp = "="
            Break
        ElseIf comp = ">"
            res = SubstracString(res, des)
            des = PasZeroDebFin(MultipliString(des, "0.1"))
            If Mid(des, 2, 1) = "." : nb_decimal + 1 : EndIf
        EndIf
        res = AddString(res, des)
        If nb_decimal > decimale : res = Left(res, Len(res) - 1) : EndIf
    Until nb_decimal > decimale
    
    final = res
    
    ProcedureReturn final
EndProcedure

Procedure.s DecomposeString(str1.s, gad.l); 61932510657186624000
    If FindString(str1, ".", 0) : ProcedureReturn "Erreur: un nombre à virgule ne peut être décomposé en nombre premier" : EndIf
    
    nb_dec = 10
    
    final.s = ""
    
    A.s = str1
    B.s = ""
    C.s = ""
    
    While FracString(PasZeroDebFin(MultipliString(A, "0.5"))) = "0"
        ; ###### point de mesure 2 compose A debut
        final + "*2"
        
        A = PasZeroDebFin(MultipliString(A, "0.5"))
        If A = "1"
            final = Right(final, Len(final) - 1)
            If Left(str1, 1) = "-" : final = "-" + final : EndIf
            ProcedureReturn final
        EndIf
    Wend
    
    B = "3"
    
    Repeat
        C = AddString(RacineCarreString(A, nb_dec), "1")
        
        D:
        SetGadgetText(gad, MultipliString(SubstracString(IntString(C), B), "0.5"))
        WindowEvent()
        
        If CompareNumString(B, C) = ">" Or CompareNumString(B, C) = "="
            ; ###### point de mesure A compose A debut
            final + "*" + A
            
            final = Right(final, Len(final) - 1)
            If Left(str1, 1) = "-" : final = "-" + final : EndIf
            ProcedureReturn final
        EndIf
        
        If FracString(DivideString(A, B, nb_dec)) <> "0"
            B = AddString(B, "2")
            Goto D
        EndIf
        
        ; ###### point de mesure B compose A debut
        final + "*" + B
        
        A = PasZeroDebFin(DivideString(A, B, nb_dec))
    Until 0 = 1
EndProcedure

Procedure.s Cal_Pi_1(iteration, decimale) ; méthode du produit de Wallis
    pi.s = "2"
    
    For a = 0 To iteration * 2 Step 2
        b.s = DivideString(Str(a + 2), Str(a + 1), decimale)
        c.s = DivideString(Str(a + 2), Str(a + 3), decimale)
        d.s = MultipliString(b, c)
        pi = Left(MultipliString(pi, d), 2 + decimale)
        
        SetGadgetText(5, Str(a) + " - " + pi)
        If WindowEvent() = #PB_Event_CloseWindow : End : EndIf
        
    Next
    
    ProcedureReturn pi
EndProcedure

Procedure.s Cal_Pi_2(iteration, decimale) ; méthode de la formule de Leibniz
    pi.s = "0"
    
    For a = 0 To iteration
        c.f = a / 2
        b.s = DivideString("4", Str(a * 2 + 1), decimale)
        If c - Int(c) = 0
            pi = AddString(pi, b)
        Else
            pi = SubstracString(pi, b)
        EndIf
        
        SetGadgetText(5, Str(a) + " - " + pi)
        If WindowEvent() = #PB_Event_CloseWindow : End : EndIf
        
    Next
    ProcedureReturn pi
EndProcedure

Procedure.s Cal_Pi_3(iteration, decimale) ; méthode de la formule BBP
    pi.s = "0"
    
    decimale + 2
    
    For a = 0 To iteration + 2
        b.s = DivideString("1", PowString("16", Str(a)), decimale)
        c.s = DivideString("4", Str(a * 8 + 1), decimale)
        d.s = DivideString("2", Str(a * 8 + 4), decimale)
        e.s = DivideString("1", Str(a * 8 + 5), decimale)
        f.s = DivideString("1", Str(a * 8 + 6), decimale)
        
        d = SubstracString(c, d)
        e = SubstracString(d, e)
        f = SubstracString(e, f)
        
        old_pi.s = pi
        
        pi = AddString(Left(MultipliString(b, f), 2 + decimale), pi)
        
        ; Debug Str(a) + " - " + pi
        
        SetGadgetText(5, Str(a) + " - " + pi)
        If WindowEvent() = #PB_Event_CloseWindow : End : EndIf
        Delay(10)
        
        If pi = old_pi : Break : EndIf
        
    Next
    ProcedureReturn Left(pi, Len(pi) - 2)
EndProcedure

Procedure.s Nombre_Premier(iteration, decimale)
    B.s = "1"
    
    Dim_prem = 100
    Dim Prem.s(Dim_prem)
    
    Prem(0) = "2"
    
    
    deb = ElapsedMilliseconds()
    
    ; MeasureHiResIntervalStart()
    
    Repeat
        B = AddString(B, "2")
        C = 0
        
        Repeat
            If FracString(DivideString(B, Prem(C), decimale)) = "0"
                B = AddString(B, "1")
                C = -1
            EndIf
            
            C + 1
        Until prem(C) = ""
        
        Prem(C) = B
        
        ;   Debug Str(C) + " - " + B + "     temp = " + StrD(1000 * MeasureHiResIntervalStop()) + "ms"
        
        If C = Dim_prem And C <> iteration
            Dim_prem + 100
            ReDim Prem.s(Dim_prem)
        EndIf
        
        ;   MeasureHiResIntervalStart()
        
    Until C = iteration
    
    fin = ElapsedMilliseconds()
    
    Debug ""
    Debug "Total pour " + Str(iteration) + " nombre premier"
    Debug Str(fin - deb) + "ms"
    
    Debug ""
    
EndProcedure

Procedure.s Factorielle(str1.s)
    str1 = PasZeroDebFin(str1)
    
    res.s ="1"
    
    For k = 1 To Val(str1)
        res = MultipliString(res, Str(k))
    Next
    
    ProcedureReturn res
EndProcedure


; resultat :
; 1: tres lent (3 decimales heures)
; 2: tres lent (3 decimales heures)
; 3: efficace, mais les 2 derniere décimales sont fausses (2 decimales a chaque boucle)

nb_dec = 100


Repeat
    event = WaitWindowEvent()
    
    ;{ Gadget
    
    If event = #PB_Event_Gadget
        Select EventGadget()
            Case 6
                temp = ElapsedMilliseconds()
                SetGadgetText(5, PasZeroDebFin(AddString(GetGadgetText(2), GetGadgetText(3))))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "R = A + B          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 7
                nb = Val(InputRequester("Division", "Nombres de décimales ?", Str(nb_dec)))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, PasZeroDebFin(DivideString(GetGadgetText(2), GetGadgetText(3), nb)))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "R = A / B          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 8
                temp = ElapsedMilliseconds()
                SetGadgetText(5, PasZeroDebFin(MultipliString(GetGadgetText(2), GetGadgetText(3))))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "R = A * B          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 9
                temp = ElapsedMilliseconds()
                SetGadgetText(5, PasZeroDebFin(SubstracString(GetGadgetText(2), GetGadgetText(3))))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "R = A - B          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 10
                temp = ElapsedMilliseconds()
                SetGadgetText(5, PasZeroDebFin(PowString(GetGadgetText(2), GetGadgetText(3))))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "R = A ^ B          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 11
                nb = Val(InputRequester("Calcul de Racine", "Nombres de décimales ?", Str(nb_dec)))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, PasZeroDebFin(RacineCarreString(GetGadgetText(2), nb)))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "R = Racine Bième de A          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 12
                temp = ElapsedMilliseconds()
                SetGadgetText(5, DecomposeString(GetGadgetText(2), 5))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "A en facteur 1er          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 13
                temp = ElapsedMilliseconds()
                SetGadgetText(5, CompareNumString(GetGadgetText(2), GetGadgetText(3)))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Comparaison          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 14
                temp = ElapsedMilliseconds()
                SetGadgetText(5, OrdreGrandeurString(GetGadgetText(2), GetGadgetText(3)))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 15
                iteration = Val(InputRequester("Calcul de PI", "Nombres d'iteration ?", ""))
                nb = Val(InputRequester("Calcul de PI", "Nombres de décimales ?", ""))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, Cal_Pi_1(iteration, nb))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 16
                iteration = Val(InputRequester("Calcul de PI", "Nombres d'iteration ?", ""))
                nb = Val(InputRequester("Calcul de PI", "Nombres de décimales ?", ""))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, Cal_Pi_2(iteration, nb))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 17
                iteration = Val(InputRequester("Calcul de PI", "Nombres d'iteration ?", ""))
                nb = Val(InputRequester("Calcul de PI", "Nombres de décimales ?", ""))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, Cal_Pi_3(iteration, nb))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 18
                iteration = Val(InputRequester("Nombres premiers", "NB de premier ?", ""))
                nb = Val(InputRequester("", "Nombres de décimales ?", ""))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, Nombre_Premier(iteration, nb))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 19
                temp = ElapsedMilliseconds()
                SetGadgetText(5, Factorielle(GetGadgetText(2)))
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 20
                iteration = Val(InputRequester("", "Nombres d'iteration ?", ""))
                nb = Val(InputRequester("", "Nombres de décimales ?", ""))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, "")
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
            Case 21
                iteration = Val(InputRequester("", "Nombres d'iteration ?", ""))
                nb = Val(InputRequester("", "Nombres de décimales ?", ""))
                
                temp = ElapsedMilliseconds()
                SetGadgetText(5, "")
                temp = ElapsedMilliseconds() - temp
                StatusBarText(0, 0, "Ordre de grandeur          card(A) = " + Str(Len(GetGadgetText(2))) + "    card(B) = " + Str(Len(GetGadgetText(3))) + "    card(R) = " + Str(Len(GetGadgetText(5))) + "      Delay de reponse = " + Str(temp) + " ms")
                
                
        EndSelect
    EndIf
    
    ;}
    
Until event = #PB_Event_CloseWindow

End