¡Gira para ver todos los lados!

10

Digamos que tienes un dado de 20 lados. Empiezas a tirar ese dado y tienes que tirarlo unas docenas de veces antes de tirar finalmente los 20 valores. Te preguntas, ¿cuántos rollos necesito antes de tener un 50% de posibilidades de ver los 20 valores? ¿Y cuántas tiradas de ndado muero necesito lanzar antes de tirar todos los nlados?

Después de investigar un poco, descubre que existe una fórmula para calcular la posibilidad de obtener todos los nvalores después de los resultados r.

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

donde S(a, b)denota los números de Stirling del segundo tipo , el número de formas de dividir un conjunto de n objetos (cada rollo) en k subconjuntos no vacíos (cada lado).

También encontrará la secuencia OEIS , que llamaremos R(n), que corresponde a la más pequeña rdonde P(r, n)es al menos el 50%. El desafío es calcular el ntérmino th de esta secuencia lo más rápido posible.

El reto

  • Dado un n, encuentre el más pequeño r donde P(r, n)sea ​​mayor o igual 0.5o 50%.
  • Teóricamente, su código debería manejar cualquier número entero no negativo ncomo entrada, pero solo probaremos su código en el rango de 1 <= n <= 1000000.
  • Para la puntuación, estaremos tomar el tiempo total necesario para funcionar R(n)en las entradas 1a través 10000.
  • Verificaremos si sus soluciones son correctas ejecutando nuestra versión de R(n)en su salida para ver si P(your_output, n) >= 0.5y P(your_output - 1, n) < 0.5, es decir, que su salida es realmente la más pequeña rpara un determinado n.
  • Puede usar cualquier definición para S(a, b)en su solución. Wikipedia tiene varias definiciones que pueden ser útiles aquí.
  • Puede utilizar las funciones integradas en sus soluciones, incluidas las que calculan S(a, b)o incluso las que calculan P(r, n)directamente.
  • Puede codificar hasta 1000 valores R(n)y un millón de números de Stirling, aunque ninguno de estos son límites estrictos, y puede cambiarse si puede presentar un argumento convincente para aumentarlos o disminuirlos.
  • No es necesario que compruebes cada posible rentre ny lo rque estamos buscando, pero sí necesitas encontrar el más pequeño ry no cualquier rlugar P(r, n) >= 0.5.
  • Su programa debe usar un lenguaje que se pueda ejecutar libremente en Windows 10.

Las especificaciones de la computadora que probará sus soluciones son i7 4790k, 8 GB RAM. Gracias a @DJMcMayhem por proporcionar su computadora para la prueba. Siéntase libre de agregar su propio tiempo no oficial para referencia, pero el tiempo oficial se proporcionará más tarde una vez que DJ pueda probarlo.

Casos de prueba

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

Avíseme si tiene alguna pregunta o sugerencia. ¡Buena suerte y buena optimización!

Sherlock9
fuente
1
@ JonathanAllan sabía que debería haber elegido una redacción diferente. Gracias por el aviso.
Sherlock9

Respuestas:

7

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

Pruébalo en línea!

Cómo funciona

Esto utiliza la serie de forma cerrada para P ( r , n ), y su derivada con respecto a r , reorganizada para la estabilidad numérica y la vectorización, para hacer una búsqueda del método de Newton para r tal que P ( r , n ) = 0.5, redondeando r a un número entero antes de cada paso, hasta que el paso se mueva r en menos de 3/4. Con una buena suposición inicial, esto generalmente toma solo una o dos iteraciones.

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 ix i + 1 = (-1) i + 1 ⋅ binom ( n , i + 1) ⋅ (( n - i - 1) / n ) r log (( n - i - 1) / n)
z ix i + 1 = d / d r P ( r , n )

Anders Kaseorg
fuente
1
Excelente trabajo en toda la respuesta! Primero, debería haberme dado cuenta de que 0.366512era logalgo de hace años. Lo usaré -log(log(2)en mi próxima iteración. En segundo lugar, la idea de utilizar el método de Newton también es muy inteligente y me alegra ver que esto funciona tan bien. Tercero, casi definitivamente voy a robar exp(log(binom(n, i+1)) + r * log((n-i-1)/n)): ¡P Kudos con una gran respuesta! : D
Sherlock9
1
¡He agregado el tiempo oficial! Buena respuesta por cierto :)
James
2
Estoy realmente confundido. Cambié la numpyimportación from numpy import *y, por alguna razón, el tiempo se redujo básicamente a 0 ... ¿ Probar en línea ?
notjagan
@notjagan caché hit tal vez?
NoOneIsHere
1
Me gustaría disculparme por varias cosas: 1) mi plagio de su respuesta cuando intenté encontrar mejoras; 2) no reconocerlo adecuadamente y solo tratar de arreglar mi respuesta; 3) que esta disculpa ha tardado tanto. Estaba tan mortificado que al principio simplemente abandoné este desafío. En un pequeño intento de reparación, supongo que es justo decirte que mi mejora principal en esta respuesta fue cambiar del método de Newton a incrementar r, ya que tu aproximación inicial ya es bastante buena. Espero verte en PPCG una vez más, y lo siento por todo.
Sherlock9