Angenommen, Sie haben einen 20-seitigen Würfel. Du fängst an, diesen Würfel zu würfeln und musst ihn ein paar Dutzend Mal würfeln, bevor du schließlich alle 20 Werte würfelst. Sie fragen sich, wie viele Rollen ich brauche, bevor ich eine 50% ige Chance bekomme, alle 20 Werte zu sehen? Und wie viele Rollen eines neinseitigen Würfels muss ich rollen, bevor ich alle nSeiten rolle ?

Nach einigen Recherchen stellen Sie fest, dass es eine Formel zur Berechnung der Wahrscheinlichkeit gibt, alle nWerte nach dem Würfeln zu rwürfeln.

P(r, n) = n! * S(r, n) / n**r

Dabei S(a, b)bezeichnet Stirling-Zahlen der zweiten Art die Anzahl der Möglichkeiten, eine Menge von n Objekten (jede Rolle) in k nicht leere Teilmengen (jede Seite) zu unterteilen.

Sie finden auch die OEIS Sequenz , die wir nennen wollen R(n), dass entspricht dem kleinsten , rwo P(r, n)mindestens 50% beträgt. Die Herausforderung besteht darin, den ndritten Term dieser Sequenz so schnell wie möglich zu berechnen .

Die Herausforderung

• Wenn gegeben n, finden Sie die kleinste, r wo P(r, n)größer oder gleich 0.5oder 50% ist.
• Ihr Code sollte theoretisch jede nicht negative Ganzzahl nals Eingabe behandeln, aber wir werden Ihren Code nur im Bereich von testen 1 <= n <= 1000000.
• Für die Wertung wird die Gesamtzeit benötigt, die erforderlich ist, um die R(n)Eingaben zu 1durchlaufen 10000.
• Wir werden überprüfen, ob Ihre Lösungen korrekt sind, indem wir unsere Version von R(n)auf Ihrer Ausgabe ausführen, um festzustellen, ob P(your_output, n) >= 0.5und P(your_output - 1, n) < 0.5, dh ob Ihre Ausgabe tatsächlich die kleinste rfür eine bestimmte ist n.
• Sie können eine beliebige Definition für S(a, b)in Ihrer Lösung verwenden. Wikipedia hat mehrere Definitionen, die hier hilfreich sein können.
• Sie können in Ihren Lösungen integrierte Funktionen verwenden, einschließlich solcher, die berechnen S(a, b), oder sogar solcher, die P(r, n)direkt berechnen .
• Sie können bis zu 1000 Werte R(n)und eine Million Stirling-Zahlen fest codieren , obwohl keine dieser Grenzwerte hart ist. Sie können geändert werden, wenn Sie ein überzeugendes Argument für das Erhöhen oder Verringern dieser Werte vorbringen können.
• Sie müssen nicht alle möglichen überprüfen müssen rzwischen nund die rwir suchen, aber Sie tun müssen, finden die kleinste rund nicht nur irgendwelche rwo P(r, n) >= 0.5.
• Ihr Programm muss eine Sprache verwenden, die unter Windows 10 frei ausgeführt werden kann.

Die Spezifikationen des Computers, der Ihre Lösungen testet, sind i7 4790k, 8 GB RAM. Vielen Dank an @DJMcMayhem für die Bereitstellung seines Computers für die Tests. Fühlen Sie sich frei, Ihr eigenes inoffizielles Timing als Referenz hinzuzufügen , aber das offizielle Timing wird später bekannt gegeben, sobald DJ es testen kann.

Testfälle

n       R(n)
1       1
2       2
3       5
4       7
5       10
6       13
20      67       # our 20-sided die
52      225      # how many cards from a huge uniformly random pile until we get a full deck
100     497
366     2294     # number of people for to get 366 distinct birthdays
1000    7274
2000    15934
5000    44418
10000   95768
100000  1187943
1000000 14182022

Lassen Sie mich wissen, wenn Sie Fragen oder Anregungen haben. Viel Glück und gute Optimierung!

Python + NumPy, 3,95 Sekunden

from __future__ import division
import numpy as np

def rolls(n):
    if n == 1:
        return 1
    r = n * (np.log(n) - np.log(np.log(2)))
    x = np.log1p(np.arange(n) / -n)
    cx = x.cumsum()
    y = cx[:-1] + cx[-2::-1] - cx[-1]
    while True:
        r0 = np.round(r)
        z = np.exp(y + r0 * x[1:])
        z[::2] *= -1
        r = r0 - (z.sum() + 0.5) / z.dot(x[1:])
        if abs(r - r0) < 0.75:
            return np.ceil(r).astype(int)

for n in [1, 2, 3, 4, 5, 6, 20, 52, 100, 366, 1000, 2000, 5000, 10000, 100000, 1000000]:
    print('R({}) = {}'.format(n, rolls(n)))

import timeit
print('Benchmark: {:.2f}s'.format(timeit.timeit(lambda: sum(map(rolls, range(1, 10001))), number=1)))

Probieren Sie es online aus!

Wie es funktioniert

Dies verwendet die geschlossene Reihe für P ( r , n ) und ihre Ableitung in Bezug auf r , die für numerische Stabilität und Vektorisierung neu angeordnet wurden, um eine Newtonsche Methodensuche nach r durchzuführen, so dass P ( r , n ) = 0,5, Rundung r auf eine ganze Zahl vor jedem Schritt, bis sich der Schritt r um weniger als 3/4 bewegt . Mit einer guten anfänglichen Vermutung dauert dies normalerweise nur ein oder zwei Iterationen.

x _i = log (1 - i / n ) = log (( n - i ) / n )
cx _i = log ( n ! / (( n - i - 1)! ⋅ n ^{i + 1} )
y _i = cx _i + cx _{n - i - 2} - cx _{n - 1} = log binom ( n , i + 1)
z _i = (-1) ^{i + 1} ⋅ binom ( n ,i + 1) ⋅ (( n - i - 1) / n ) ^r
1 + ∑ z _i = n! ⋅ S ( r , n ) / n ^r = P ( r , n )
z _i ⋅ x _{i + 1} = (-1) ^{i + 1} ⋅ Binom ( n , i + 1) ⋅ (( n - i - 1) / n ) ^r log (( n - i - 1) / n)
∑ z _i ⋅ x _{i + 1} = d / d r P ( r , n )