Digamos que você tenha um dado de 20 lados. Você começa a rolar esse dado e precisa rolar algumas dezenas de vezes antes de finalmente rolar todos os 20 valores. Você quer saber, quantos rolos eu preciso antes de ter 50% de chance de ver todos os 20 valores? E quantos rolos de nmatriz de um lado eu preciso rolar antes de rolar por todos os nlados?

Após algumas pesquisas, você descobre que existe uma fórmula para calcular a chance de rolar todos os nvalores após as rrolagens.

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

onde S(a, b)denota números Stirling do segundo tipo , o número de maneiras de particionar um conjunto de n objetos (cada rolo) em k subconjuntos não vazios (cada lado).

Você também encontra a sequência OEIS , que chamaremos de R(n), que corresponde à menor rem P(r, n)pelo menos 50%. O desafio é calcular o ntermo th desta sequência o mais rápido possível.

O desafio

• Dado um n, encontre o menor r onde P(r, n)é maior ou igual a 0.5ou 50%.
• Teoricamente, seu código deve manipular qualquer número inteiro não negativo ncomo entrada, mas apenas testaremos seu código no intervalo de 1 <= n <= 1000000.
• Para pontuação, estaremos tomar o tempo total necessário para executar R(n)em entradas 1através 10000.
• Verificaremos se suas soluções estão corretas executando nossa versão de R(n)em sua saída para ver se P(your_output, n) >= 0.5e P(your_output - 1, n) < 0.5, ou seja, se sua saída é realmente a menor rpara um dado n.
• Você pode usar qualquer definição para S(a, b)em sua solução. A Wikipedia possui várias definições que podem ser úteis aqui.
• Você pode usar built-ins em suas soluções, incluindo aquelas que calculam S(a, b), ou mesmo aquelas que calculam P(r, n)diretamente.
• Você pode codificar até 1000 valores R(n)e um milhão de números Stirling, embora nenhum deles seja um limite rígido, e poderá ser alterado se você puder argumentar de maneira convincente para aumentá-los ou diminuí-los.
• Você não precisa verificar todos os possíveis rentre ne o rque estamos procurando, mas precisa encontrar o menor re não apenas ronde P(r, n) >= 0.5.
• Seu programa deve usar um idioma que possa ser executado livremente no Windows 10.

As especificações do computador que testará suas soluções são i7 4790k, 8 GB RAM. Agradecemos a @DJMcMayhem por fornecer seu computador para os testes. Sinta-se à vontade para adicionar seu próprio tempo não oficial para referência, mas o tempo oficial será fornecido mais tarde quando o DJ puder testá-lo.

Casos de teste

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

Entre em contato se tiver alguma dúvida ou sugestão. Boa sorte e boa otimização!

Python + NumPy, 3,95 segundos

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

Experimente online!

Como funciona

Isso usa a série de forma fechada para P ( r , n ) e sua derivada em relação a r , reorganizada para estabilidade numérica e vetorização, para fazer uma pesquisa de método de Newton por r de modo que P ( r , n ) = 0,5, arredondando r para um número inteiro antes de cada etapa, até que a etapa se mova r em menos de 3/4. Com um bom palpite inicial, isso geralmente leva apenas uma ou duas iterações.

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 )