Powiedzmy, że masz 20-stronną kostkę. Zaczynasz rzucać tą kością i musisz rzucić ją kilkadziesiąt razy, zanim w końcu rzucisz wszystkie 20 wartości. Zastanawiasz się, ile rzutów potrzebuję, zanim otrzymam 50% szansy na zobaczenie wszystkich 20 wartości? A ile rzutów nkostką jednostronną muszę wykonać, zanim wykonam rzut ze wszystkich nstron?

Po kilku badaniach okazało się, że istnieje wzór do obliczania szansy na wyrzucenie wszystkich nwartości po rrzutach.

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

gdzie S(a, b)oznacza liczby Stirlinga drugiego rodzaju , liczbę sposobów podziału zestawu n obiektów (każda rolka) na k niepustych podzbiorów (z każdej strony).

Znajdziesz również sekwencję OEIS , którą nazwiemy R(n), która odpowiada najmniejszemu, rgdzie P(r, n)wynosi co najmniej 50%. Wyzwanie polega na njak najszybszym obliczeniu th tej sekwencji.

Wyzwanie

• Biorąc pod uwagę an n, znajdź najmniejsze, r gdzie P(r, n)jest większe lub równe 0.5lub 50%.
• Twój kod powinien teoretycznie obsługiwać nieujemną liczbę całkowitą njako dane wejściowe, ale będziemy testować twój kod tylko w zakresie 1 <= n <= 1000000.
• Do punktacji, będziemy brać całkowity czas potrzebny do uruchomienia R(n)na wejściach 1dzięki 10000.
• Sprawdzimy, czy twoje rozwiązania są poprawne, uruchamiając naszą wersję R(n)na twoim wyjściu, aby sprawdzić, czy P(your_output, n) >= 0.5i P(your_output - 1, n) < 0.5, tj. Czy twoje wyjście jest w rzeczywistości najmniejsze rdla danego n.
• Możesz użyć dowolnej definicji S(a, b)w swoim rozwiązaniu. Wikipedia ma kilka definicji, które mogą być tutaj pomocne.
• Możesz używać wbudowanych rozwiązań, w tym tych, które obliczają S(a, b), a nawet tych, które obliczają P(r, n)bezpośrednio.
• Możesz zakodować na stałe do 1000 wartości R(n)i miliona liczb Stirlinga, chociaż żadna z nich nie jest twardym limitem, i można je zmienić, jeśli będziesz w stanie przekonująco argumentować za ich podniesieniem lub obniżeniem.
• Nie musisz sprawdzać wszystkich możliwych rmiędzy ntym r, czego szukamy, ale musisz znaleźć najmniejsze, ra nie tylko rgdziekolwiek P(r, n) >= 0.5.
• Twój program musi używać języka, który można swobodnie uruchamiać w systemie Windows 10.

Specyfikacje komputera, który przetestuje twoje rozwiązania, są i7 4790k, 8 GB RAM. Dzięki @DJMcMayhem za udostępnienie komputera do testowania. Możesz dodać własne nieoficjalne terminy w celach informacyjnych, ale oficjalny harmonogram zostanie podany później, gdy DJ będzie mógł go przetestować.

Przypadki testowe

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

Daj mi znać, jeśli masz jakieś pytania lub sugestie. Powodzenia i dobrej optymalizacji!

Python + NumPy, 3,95 sekundy

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)))

Wypróbuj online!

Jak to działa

Korzysta z serii zamkniętej dla P ( r , n ) i jego pochodnej w odniesieniu do r , przestawionej pod względem stabilności numerycznej i wektoryzacji, aby przeprowadzić metodę Newtona w poszukiwaniu r tak, że P ( r , n ) = 0,5, zaokrąglenie r oznacza liczbę całkowitą przed każdym etapie, do momentu, aż etap r o mniej niż 3/4. Przy dobrym wstępnym przypuszczeniu zwykle zajmuje to tylko jedną lub dwie iteracje.

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 )