OverflowError bei Division mit großen Zahlen

[Mr. Q]

Hallo,
ich habe ein Programm geschrieben, dass mit großen Zahlen arbeitet:

Code: Alles auswählen

fac = lambda n:[1,0][n>0] or fac(n-1)*n  #Funktion, die eine Fakulät berechnet
a=366
b=input()                                       #Zahl zwischen 0 und 400
c=fac(a-1)/fac(a-b)/a**(b-1.0)
print c                                           #Ergebnis, sollte zwischen 0 und 1 liegen

Bei Eingabe für b bis 121 gibt das Programm Werte aus, ab 122 wird eine Fehlermeldung ausgegeben:

Traceback (most recent call last):
File "Gebw2.py", line 4, in <module>
c=fac(a-1)/fac(a-b)/a**(b-1.0)
OverflowError: (34, 'Numerical result out of range')

Ich verstehe, dass die Zahlenwerte ein wenig hoch sind aber so sind sie eben. Was kann man denn da machen?

[EyDu]

Das liegt daran, dass das Ergebnis von "a**(b-1.0)" zu groß für einen Float wird.

Umgehen kannst du dass, indem du vor den Divionsionen schon ein wenig kürz. Das ist bei der Fakutlät aber recht einfach.

Das decimal-Modul könnte auch eine Hilfe sein.

[numerix]

Hat eine große Affinität zum sog. "Geburtstagsproblem" ...
b = 400 wird sowieso nix, weil die Fakultät negativer Zahlen nicht definiert ist.

Das liegt daran, dass das Ergebnis von "a**(b-1.0)" zu groß für einen Float wird. Umgehen kannst du dass, indem du vor den Divionsionen schon ein wenig kürz. Das ist bei der Fakutlät aber recht einfach.

Damit allein dürfte es nicht getan sein, weil die Potenz dadurch auch nicht kleiner wird.

Eine Möglichkeit wäre es, auf die Fakultäten und Potenzen ganz zu verzichten, und daraus ein Produkt aus einzelnen Werten x mit 0<x<1 zu machen.

Nimm z.B. mal a = 366 und b = 6 (oder sonst was) und schreib dir für diese Zahlen auf, was genau zu rechnen ist. Dann kürzt du die Fakultäten, schreibst es als Bruch mit Zähler und Nenner als Produkt mit einzelnen Faktoren und das dann als Produkt aus lauter echten Brüchen. Dann kannst du es problemlos multiplizieren.
Wäre dann am Ende nur zu sehen, wie stark die Einbußen bei der Genauigkeit sind.

[Mr. Q]

@ pütone: Goldrichtig, ich versuche das Geburtstagsproblem (Wie groß ist die Wahrscheinlichkeit, dass mindestens 2 von n Personen am gleichen Tag des Jahres Geburtstag haben?)zu lösen

Vielen Dank,
dieses Problem ist gelöst. Dafür habe ich jetzt ein neues.
Ich habe das Programm entsprechend pütones Tipp modifiziert.

Code: Alles auswählen

a=366
g=a
n=input()
Var=1.0
for i in range(1,n):
 g=g-1
 Var=float(str(round(Var*g/a, 9)))
print 1-Var

Nun ergibt '1-Var' aber schon ab n=185 (Personen) 1.
Wenn in Zeile 7 nur ' Var=Var*g/a' stünde sogar schon ab n=135.
Das darf aber erst ab n=366 passieren.

[BlackJack]

Das sind die üblichen Rechenungenauigkeiten bei Fliesskommazahlen. Wenn man das genau, also mit Brüchen, zum Beispiel mit dem `gmpy`-Modul, durchlaufen lässt, gibt's bei ``n = 367`` die 1 als Ergebnis.

[numerix]

Nun ergibt '1-Var' aber schon ab n=185 (Personen) 1.
Wenn in Zeile 7 nur ' Var=Var*g/a' stünde sogar schon ab n=135.
Das darf aber erst ab n=366 passieren.

Du solltest mal reflektieren, WAS du da berechnest. Die Wahrscheinlichkeit nähert sich mit steigender Personenzahl immer mehr der 1 an, ab 58 Personen liegt sie schon bei >0,99.

Die Wahrscheinlichkeit, dass bei Zahlen >100 zwei NICHT am gleichen Tag Geburtstag haben, ist so verschwindend gering, dass man das mit an Sicherheit grenzender Wahrscheinlichkeit ausschließen kann. Dein Ergebnis ist also nicht falsch, sondern im Rahmen der Berechnung mit Fließkommazahlen völlig korrekt.

Wenn du es noch genauer haben willst - sinnvoll oder wirklich brauchbar sind die Ergebnisse allerdings nicht mehr - dann kannst du es z.B. auch so versuchen:

Code: Alles auswählen

from decimal import Decimal, getcontext

tage = 366
personen = 300
p = Decimal(1)
getcontext().prec = 100
for k in xrange(tage-1,tage-personen,-1):
    p *= Decimal(k)/Decimal(tage)
print "P('Mind. %i Personen am gleichen Tag Geburtstag'):" %personen
print Decimal(1)-p

[sma]

Ist "b" eine Ganzzahl? Dann benutze auch Ganzzahlen in der Berechnung:

Code: Alles auswählen

>>> 366 ** (122-1)
15177887613753647486770963503976599403069898859246477999404393
23912237005154299536838728176049372963106677640549492756096029
35925425499083347624974918807284128222924948002181254069655945
62095414350674285483819146990052099981392178594487189045613846
15398155954378684618053965754616279555858086569634538897119641
6L
>>>

Kann Python auch mit Brüchen rechnen? Das würde bei diesem Problem helfen, "float" als Datentyp komplett zu vermeiden und exakt zu bleiben. Ein gutes Scheme hätte damit ebenso keine Probleme wie Smalltalk. Da muss es doch aber auch etwas für Python geben...

Dieses Smalltalk-Programm liefert "7.83872351432298e-11":

Code: Alles auswählen

| a b |
a := 366. b := 123.
((a - 1) factorial / (a - b) factorial / (a raisedTo: b - 1)) asFloat

Stefan

[BlackJack]

Für Brüche hatte ich ja schon `gmpy`, die Python-Anbindung an die `libgmp` genannt. Damit ist man genau.

Code: Alles auswählen

marc@god:~$ python test.py
366
18598330272217692646750287014083139247278428111290059172714838737443509708825644
09911182717345333073670615793409765697833394537085349808061439017345237170970053
25087179685328965716370924702233728963349855669458306976739703701656268511161756
40507671875284126961554912992040206739337395460707885647158645335403154619433104
66868002839390295680211760730588602200074169511026426054787832951328713449708884
04747710967688582401298256489904367358300531050631661258359136129061743076510129
52078612512191984290930788941954676807279396173544814203001236051247515289984870
00443590720827492715441345205591662818720794486809915778587820110022678106403437
91676902401405655885718585215002706647624914817038404251217784534201297848259868
75674438480231/18598330272217692646750287014083139247278428111290059172714838737
44350970882564409911182717345333073670615793409765697833394537085349808061439017
34523717097015293429396474729535170502122985026783862199794952770407543218822816
36429640690745388000023681348841558545325391500933675732029148591910279141486970
64631014967116151991723638430661601970233411415047942543881937326828704807550997
56227688313258233939049744363350745432733890596700408763807302324336515804151279
99024555723170596422249127811916235307928416228293204056672712377054192253172280
39965151748616481353535708911269260029515488261841753839287049710946404657883581
32577401441954561586418344253331005028819778256636854724711997907609680294591513
28420097451454652926635745856
marc@god:~$ python test.py
367
1

Irgendwo in dem ersten Ergebnis versteckt sich ein '/', der Zähler und Nenner von dem Bruch trennt.

Das Programm:

Code: Alles auswählen

import gmpy

def main():
    a = 366
    n = int(raw_input())
    var = gmpy.mpq(1)
    for g in reversed(xrange(a - n + 1, a)):
        var = var * g / a
    print 1 - var

if __name__ == '__main__':
    main()

[birkenfeld]

Kann Python auch mit Brüchen rechnen?

[mod]fractions[/mod]

Ist allerdings erst neu in 2.6.

[Mr. Q]

Vieln Dank für eure Hilfe. Ich habe das Problem jetzt etwas unsauber, aber (mich) optisch zufriedenstellend gelöst:

Code: Alles auswählen

def Gebwahrsch(Personen): #Die Wahrscheinlichkeit, dass 2 aus n Personen am selben Tag Geburtstag haben
  Auswahl=366
  Prob=1.0
  for f in range(1,Personen):
     Auswahl=Auswahl-1
     Prob=float(str(round(Prob*Auswahl/366, 9)))
 
  if Prob==0 and i < 366: return 0.999999999      #Ergebnisrundung verhindern
  else: return 1-Prob

[BlackJack]

Also ich finde das extrem unsauber. Wenn man Rechenergebnisse so kosmetisch frisiert, kann man sich auch die Rechnung gleich sparen. Wobei diese `float()`/`str()`/`round()`-Kombination schon stinkt. Entweder man akzeptiert die Eigenheiten von Fliesskommazahlen, oder verwendet einen anderen Datentyp.

[EyDu]

Das Verhindern der Ergebnisrundung macht dein Ergebnis übrigens ungenauer:

Für 2 Geburtstage an unterschiedlichen Tagen:

Code: Alles auswählen

366!/((366-n)! * 366^n)

mit: n = 355

366!/(1! * 366^355) < (183/366)^172 = 0.5^172 < 10^(-50)

Und das ist viel dichter an 0 als an 1-0.999999999.