[Ceci est une question partenaire pour calculer une probabilité exactement ]

Cette tâche consiste à écrire du code pour calculer une probabilité exactement et rapidement . La sortie doit être une probabilité précise écrite sous forme de fraction dans sa forme la plus réduite. C'est-à-dire qu'il ne devrait jamais sortir 4/8mais plutôt 1/2.

Pour un entier positif n, considérez une chaîne uniformément aléatoire de 1 et -1 de longueur net appelez-la A. Maintenant, concaténez à Asa première valeur. C'est A[1] = A[n+1]si l'indexation à partir de 1. a Amaintenant une longueur n+1. Considérons maintenant également une deuxième chaîne aléatoire de longueur ndont les premières nvaleurs sont -1, 0 ou 1 avec probabilité 1 / 4,1 / 2, 1/4 chacune et appelons-la B.

Considérons maintenant le produit intérieur de A[1,...,n]et Bet le produit intérieur de A[2,...,n+1]et B.

Par exemple, réfléchissez n=3. Les valeurs possibles pour Aet Bpourraient être A = [-1,1,1,-1]et B=[0,1,-1]. Dans ce cas, les deux produits intérieurs sont 0et 2.

Votre code doit générer la probabilité que les deux produits internes soient nuls.

En copiant le tableau produit par Martin Büttner, nous avons les exemples de résultats suivants.

n   P(n)
1   1/2
2   3/8
3   7/32
4   89/512
5   269/2048
6   903/8192
7   3035/32768
8   169801/2097152

Langues et bibliothèques

Vous pouvez utiliser n'importe quelle langue et bibliothèque librement disponibles. Je dois être en mesure d'exécuter votre code, veuillez donc inclure une explication complète sur la façon d'exécuter / compiler votre code sous Linux si possible.

La tâche

Votre code doit commencer par n=1et donner la sortie correcte pour chaque n croissant sur une ligne distincte. Il devrait s'arrêter après 10 secondes.

Le score

Le score est simplement le plus élevé natteint avant que votre code ne s'arrête après 10 secondes lorsqu'il est exécuté sur mon ordinateur. S'il y a égalité, le vainqueur sera celui qui obtiendra le plus rapidement le score le plus élevé.

Tableau des entrées

• n = 64en Python . Version 1 par Mitch Schwartz
• n = 106en Python . Version 11 juin 2015 par Mitch Schwartz
• n = 151en C ++ . Port de la réponse de Mitch Schwartz par kirbyfan64sos
• n = 165en Python . Version 11 juin 2015 la version "élagage" de Mitch Schwartz avec N_MAX = 165.
• n = 945en Python par Min_25 en utilisant une formule exacte. Incroyable!
• n = 1228en Python par Mitch Schwartz en utilisant une autre formule exacte (basée sur la réponse précédente de Min_25).
• n = 2761en Python par Mitch Schwartz en utilisant une implémentation plus rapide de la même formule exacte.
• n = 3250en Python en utilisant Pypy par Mitch Schwartz en utilisant la même implémentation. Ce score doit pypy MitchSchwartz-faster.py |tailéviter les frais généraux de défilement de la console.

Python

Une formule sous forme fermée de p(n)est

Une fonction de génération exponentielle de p(n)est

où I_0(x)est la fonction de Bessel modifiée du premier type.

Edit on 2015-06-11:
- mise à jour du code Python.

Modifier le 13/06/2015:
- ajout d'une preuve de la formule ci-dessus.
- fixe le time_limit.
- ajout d'un code PARI / GP.

Python

def solve():
  # straightforward implementation

  from time import time
  from itertools import count

  def binom(n, r):
    return facts[n] // (facts[r] * facts[n - r])

  def p(N):
    ans = 0
    for i in range(1 + N // 2):
      t = binom(2 * (N - 2 * i), N - 2 * i)
      t *= binom(N, 2 * i)
      t *= binom(4 * i, 2 * i)
      ans += t
    e = (ans & -ans).bit_length() - 1
    numer = ans >> e
    denom = 1 << (3 * N - 1 - e)
    return numer, denom

  facts = [1]
  time_limit = 10.0 + time()

  for i in count(1):
    facts.append(facts[-1] * (2 * i - 1))
    facts.append(facts[-1] * (2 * i))

    n, d = p(i)

    if time() > time_limit:
      break

    print("%d %d/%d" % (i, n, d))

solve()

PARI / GP

p(n) = polcoeff( (exp(x/2) + 1) * besseli(0, x/4) ^ 2, n) * n!;

Preuve:
ce problème est similaire à un problème de marche aléatoire en 2 dimensions (restreint).

Si A[i] = A[i+1], nous pouvons passer de (x, y)la (x+1, y+1)[1 voie], (x, y)[2] ou moyens (x-1, y-1)[1 voie].

Si A[i] != A[i+1], nous pouvons passer de (x, y)la (x-1, y+1)[1 voie], (x, y)[2] ou moyens (x+1, y-1)[1 voie].

Laissez a(n, m) = [x^m]((x+1)^n + (x-1)^n), b(n) = [x^n](1+x)^{2n}et c(n)le nombre de façons de passer de (0, 0)la (0, 0)avec nétapes.

Alors, c(n) = \sum_{i=0}^n a(n, i) * b(i) * b(n-i).

Depuis p(n) = c(n) / 8^n, nous pouvons obtenir la formule sous forme fermée ci-dessus.

Python

Remarque: Félicitations à Min_25 pour avoir trouvé une solution sous forme fermée!

Merci pour le problème intéressant! Cela peut être résolu en utilisant DP, bien que je ne me sente pas très motivé actuellement pour optimiser la vitesse pour obtenir un score plus élevé. Ça pourrait être bien pour le golf.

Le code a atteint N=39en 10 secondes sur cet ancien ordinateur portable exécutant Python 2.7.5.

from time import*
from fractions import*
from collections import*

X={(1,0,0,0):1,(-1,0,0,0):1}

T=time()
N=0

while 1:
    Y=defaultdict(lambda:0)
    n=d=0
    for a,b,s,t in X:
        c=X[(a,b,s,t)]
        for A in ( (1,-1) if N else [a] ):
            for B in 1,0,0,-1:
                n+=c*(s+A*B==0==t+A*b+a*B)
                d+=c
                Y[(a,B,s+A*B,t+A*b)]+=c
    if time()>T+10: break
    N+=1
    print N,Fraction(n,d)
    X=Y

Pour tuple (a,b,s,t): aest le premier élément de A, best le dernier élément de B, sest le produit interne de A[:-1]et B, et test le produit interne de A[1:-1]et B[:-1], en utilisant la notation de tranche Python. Mon code ne stocke pas les tableaux Aou Bn'importe où, donc j'utilise ces lettres pour faire référence aux éléments suivants à ajouter à Aet B, respectivement. Ce choix de nom de variable rend l'explication un peu gênante, mais permet une belle apparence A*b+a*Bdans le code lui-même. Notez que l'élément ajouté Aest l'avant-dernier, car le dernier élément est toujours le même que le premier. J'ai utilisé l'astuce de Martin Büttner consistant à inclure 0deux foisBcandidats afin d'obtenir la distribution de probabilité appropriée. Le dictionnaire X(qui est nommé Ypour N+1) garde la trace du nombre de tous les tableaux possibles en fonction de la valeur du tuple. Les variables net dreprésentent le numérateur et le dénominateur, c'est pourquoi j'ai renommé l' nénoncé du problème comme N.

L'élément clé de la logique est que vous pouvez mettre à jour à partir Nde l' N+1utilisation des seules valeurs du tuple. Les deux produits intérieurs spécifiés dans la question sont donnés par s+A*Bet t+A*b+a*B. C'est clair si vous examinez un peu les définitions; notez que [A,a]et [b,B]sont les deux derniers éléments des tableaux Aet Brespectivement.

Notez que set tsont petits et bornés selon N, et pour une implémentation rapide dans un langage rapide, nous pourrions éviter les dictionnaires en faveur des tableaux.

Il peut être possible de profiter de la symétrie en considérant des valeurs ne différant que par le signe; Je n'ai pas examiné cela.

Remarque 1 : La taille du dictionnaire augmente de façon quadratique N, où taille signifie le nombre de paires clé-valeur.

Remarque 2 : Si nous définissons une limite supérieure N, nous pouvons élaguer les tuples pour lesquels N_MAX - N <= |s|et de même pour t. Cela pourrait être fait en spécifiant un état absorbant, ou implicitement avec une variable pour contenir le nombre d'états élagués (qui devraient être multipliés par 8 à chaque itération).

Mise à jour : Cette version est plus rapide:

from time import*
from fractions import*
from collections import*

N_MAX=115

def main():
    T=time()

    N=1
    Y={(1,0,0,0):1,(1,1,1,0):1}
    n=1
    thresh=N_MAX

    while time() <= T+10:
        print('%d %s'%(N,Fraction(n,8**N/4)))

        N+=1
        X=Y
        Y=defaultdict(lambda:0)
        n=0

        if thresh<2:
            print('reached MAX_N with %.2f seconds remaining'%(T+10-time()))
            return

        for a,b,s,t in X:
            if not abs(s)<thresh>=abs(t):
                continue

            c=X[(a,b,s,t)]

            # 1,1

            if not s+1 and not t+b+a: n+=c
            Y[(a,1,s+1,t+b)]+=c

            # -1,1

            if not s-1 and not t-b+a: n+=c
            Y[(a,1,s-1,t-b)]+=c

            # 1,-1

            if not s-1 and not t+b-a: n+=c
            Y[(a,-1,s-1,t+b)]+=c

            # -1,-1

            if not s+1 and not t-b-a: n+=c
            Y[(a,-1,s+1,t-b)]+=c

            # 1,0

            c+=c

            if not s and not t+b: n+=c
            Y[(a,0,s,t+b)]+=c

            # -1,0

            if not s and not t-b: n+=c
            Y[(a,0,s,t-b)]+=c

        thresh-=1

main()

Optimisations mises en œuvre:

• tout mettre en place main()- l'accès aux variables locales est plus rapide que global
• gérer N=1explicitement pour éviter la vérification (1,-1) if N else [a](qui impose que le premier élément du tuple soit cohérent, lors de l'ajout d'éléments à Apartir de la liste vide)
• dérouler les boucles intérieures, ce qui élimine également la multiplication
• doubler le nombre cpour ajouter un 0à Bau lieu de faire ces opérations deux fois
• le dénominateur est toujours 8^Nainsi nous n'avons pas besoin de le suivre
• maintenant en tenant compte de la symétrie: nous pouvons fixer le premier élément de Aas 1et diviser le dénominateur par 2, puisque les paires valides (A,B)avec A[1]=1et celles avec A[1]=-1peuvent être mises en correspondance biunivoque par négation A. De même, nous pouvons fixer le premier élément de Bcomme non négatif.
• maintenant avec l'élagage. Vous devrez jouer avec N_MAXpour voir quel score il peut obtenir sur votre machine. Il pourrait être réécrit pour trouver N_MAXautomatiquement un approprié par recherche binaire, mais semble inutile? Remarque: Nous n'avons pas besoin de vérifier l'état de l'élagage avant d'atteindre N_MAX / 2, donc nous pouvons obtenir un peu d'accélération en itérant en deux phases, mais j'ai décidé de ne pas le faire pour la simplicité et la propreté du code.

Python

En utilisant l'idée de marche aléatoire de Min_25, j'ai pu arriver à une formule différente:

Voici une implémentation Python basée sur Min_25:

from time import*
from itertools import*

def main():
    def binom(n, k):
        return facts[n]/(facts[k]*facts[n-k])

    def p(n):
        numer=0
        for i in range(n/2+1):
            t=binom(2*i,i)
            t*=t
            t*=binom(n,2*i)
            t<<=2*(n-2*i)
            numer+=t
        if not n&1:
            numer+=t
        e=(numer&-numer).bit_length()-1
        numer>>=e
        denom=1<<(3*n-e)
        return numer, denom

    facts=[1]
    time_limit=time()+10

    for i in count(1):
        facts.append(facts[-1]*i)

        n,d=p(i)

        if time()>time_limit:
            break

        print("%d %d/%d"%(i,n,d))

main()

Explication / preuve:

D'abord, nous résolvons un problème de comptage connexe, où nous le permettons A[n+1] = -A[1]; c'est-à-dire que l'élément supplémentaire concaténé en Apeut être 1ou -1indépendamment du premier élément. Nous n'avons donc pas besoin de garder une trace du nombre de fois A[i] = A[i+1]. Nous avons la marche aléatoire suivante:

De (x,y)nous pouvons passer à (x+1,y+1)[1 voie], (x+1,y-1)[1 voie], (x-1,y+1)[1 voie], (x-1,y-1)[1 voie], (x,y)[4 voies]

où xet ypour les deux produits scalaires, et nous comptons le nombre de façons de passer de (0,0)à (0,0)par nétapes. Ce nombre serait ensuite multiplié par 2pour tenir compte du fait que cela Apeut commencer par 1ou -1.

Nous nous référons à rester à (x,y)un mouvement zéro .

Nous itérons sur le nombre de mouvements non nuls i, qui doit être égal pour revenir (0,0). Les mouvements horizontaux et verticaux constituent deux marches aléatoires unidimensionnelles indépendantes, qui peuvent être comptées par C(i,i/2)^2, où C(n,k)est le coefficient binomial. (Pour une marche avec des kmarches à gauche et des kmarches à droite, il existe des C(2k,k)façons de choisir l'ordre des marches.) De plus, il existe des C(n,i)façons de placer les mouvements et des 4^(n-i)façons de choisir les mouvements zéro. On obtient donc:

a(n) = 2 * sum_{i in (0,2,4,...,n)} C(i/2,i)^2 * C(n,i) * 4^(n-i)

Maintenant, nous devons revenir au problème d'origine. Définissez une paire autorisée (A,B)à convertir si elle Bcontient un zéro. Définissez une paire (A,B)comme étant presque admissible si A[n+1] = -A[1]et les deux produits scalaires sont tous deux nuls.

Lemme: Pour une donnée n, les paires presque autorisées sont en correspondance biunivoque avec les paires convertibles.

Nous pouvons (de façon réversible) convertir une paire convertible (A,B)en une paire presque admissible (A',B')en niant A[m+1:]et B[m+1:], où mest l'indice du dernier zéro B. La vérification est simple: si le dernier élément de Best nul, nous n'avons rien à faire. Sinon, lorsque nous nions le dernier élément de A, nous pouvons nier le dernier élément de Bafin de conserver le dernier terme du produit scalaire décalé. Mais cela annule la dernière valeur du produit scalaire non décalé, nous corrigeons donc cela en annulant l'avant-dernier élément de A. Mais cela annule l'avant-dernière valeur du produit déplacé, nous annulons donc l'avant-dernier élément de B. Et ainsi de suite, jusqu'à atteindre un élément zéro B.

Maintenant, nous devons simplement montrer qu'il n'y a pas de paires presque autorisées pour lesquelles Bne contient pas de zéro. Pour qu'un produit scalaire soit égal à zéro, nous devons avoir un nombre 1et des -1termes égaux à annuler. Chaque -1terme est composé de (1,-1)ou (-1,1). Ainsi, la parité du nombre de ceux -1qui se produisent est fixée en fonction de n. Si les premier et dernier éléments de Aont des signes différents, nous changeons la parité, c'est donc impossible.

Nous obtenons donc

c(n) = a(n)/2 if n is odd, else a(n)/2 + C(n,n/2)^2

p(n) = c(n) / 8^n

qui donne la formule ci-dessus (réindexation avec i' = i/2).

Mise à jour: voici une version plus rapide utilisant la même formule:

from time import*
from itertools import*

def main():
    time_limit=time()+10

    binoms=[1]
    cb2s=[1]
    cb=1

    for n in count(1):
        if n&1:
            binoms=[a+b for a,b in zip([0]+binoms,binoms)]
        else:
            binoms=[a+b for a,b in zip([0]+binoms,binoms+[binoms[-1]])]
            cb=(cb<<2)-(cb+cb)/(n/2)
            cb2s.append(cb*cb)

        numer=0
        for i in xrange(n/2+1):
            t=cb2s[i]*binoms[min(2*i,n-2*i)]
            t<<=2*(n-2*i)
            numer+=t
        if not n&1:
            numer+=t
        e=(numer&-numer).bit_length()-1
        numer>>=e
        denom=1<<(3*n-e)

        if time()>time_limit:
            break

        print("%d %d/%d"%(n,numer,denom))

main()

Optimisations mises en œuvre:

• fonction en ligne p(n)
• utiliser la récurrence pour les coefficients binomiaux C(n,k)aveck <= n/2
• utiliser la récurrence pour les coefficients binomiaux centraux

Explication de la formule de Min_25

Min_25 a publié une excellente preuve, mais il a fallu un certain temps pour la suivre. C'est un peu d'explication à remplir entre les lignes.

a (n, m) représente le nombre de façons de choisir A telles que A [i] = A [i + 1] m fois. La formule pour a (n, m) est équivalente à a (n, m) = {2 * (n choisissez m) pour nm pair; 0 pour nm impair.} Une seule parité est autorisée car A [i]! = A [i + 1] doit se produire un nombre pair de fois pour que A [0] = A [n]. Le facteur 2 est dû au choix initial A [0] = 1 ou A [0] = -1.

Une fois que le nombre de (A [i]! = A [i + 1]) est fixé à q (nommé i dans la formule c (n)), il se sépare en deux marches aléatoires 1D de longueur q et nq. b (m) est le nombre de façons de faire une marche aléatoire unidimensionnelle de m pas qui se termine au même endroit qu'il a commencé, et a 25% de chances de se déplacer vers la gauche, 50% de chances de rester immobile et 25% de chances de se déplaçant à droite. Une façon plus évidente de déclarer la fonction de génération est [x ^ m] (1 + 2x + x ^ 2) ^ n, où 1, 2x et x ^ 2 représentent respectivement gauche, pas de mouvement et droite. Mais alors 1 + 2x + x ^ 2 = (x + 1) ^ 2.

C ++

Juste un portage de la (excellente) réponse Python de Mitch Schwartz. La principale différence est que je 2représentais -1la avariable et que je faisais quelque chose de similaire pour b, ce qui me permettait d'utiliser un tableau. En utilisant Intel C ++ avec -O3, j'ai compris N=141! Ma première version a eu N=140.

Cela utilise Boost. J'ai essayé une version parallèle mais j'ai rencontré des problèmes.

#include <boost/multiprecision/gmp.hpp>
#include <boost/typeof/typeof.hpp>
#include <boost/rational.hpp>
#include <boost/chrono.hpp>
#include <boost/array.hpp>
#include <iostream>
#include <utility>
#include <map>

typedef boost::multiprecision::mpz_int integer;
typedef boost::array<boost::array<std::map<int, std::map<int, integer> >, 3>, 2> array;
typedef boost::rational<integer> rational;

int main() {
    BOOST_AUTO(T, boost::chrono::high_resolution_clock::now());

    int N = 1;
    integer n = 1;
    array* Y = new array, *X = NULL;
    (*Y)[1][0][0][0] = 1;
    (*Y)[1][1][1][0] = 1;

    while (boost::chrono::high_resolution_clock::now() < T+boost::chrono::seconds(10)) {
        std::cout << N << " " << rational(n, boost::multiprecision::pow(integer(8), N)/4) << std::endl;
        ++N;
        delete X;
        X = Y;
        Y = new array;
        n = 0;

        for (int a=0; a<2; ++a)
            for (int b=0; b<3; ++b)
                for (BOOST_AUTO(s, (*X)[a][b].begin()); s != (*X)[a][b].end(); ++s)
                    for (BOOST_AUTO(t, s->second.begin()); t != s->second.end(); ++t) {
                        integer c = t->second;
                        int d = b&2 ? -1 : b, e = a == 0 ? -1 : a;

                        if (s->first == -1 && t->first+d+e == 0) n += c;
                        (*Y)[a][1][s->first+1][t->first+d] += c;

                        if (s->first == 1 && t->first-d+e == 0) n += c;
                        (*Y)[a][1][s->first-1][t->first-d] += c;

                        if (s->first == 1 && t->first+d-e == 0) n += c;
                        (*Y)[a][2][s->first-1][t->first+d] += c;

                        if (s->first == -1 && t->first-d-e == 0) n += c;
                        (*Y)[a][2][s->first+1][t->first-d] += c;

                        c *= 2;

                        if (s->first == 0 && t->first+d == 0) n += c;
                        (*Y)[a][0][s->first][t->first+d] += c;

                        if (s->first == 0 && t->first-d == 0) n += c;
                        (*Y)[a][0][s->first][t->first-d] += c;
                    }
    }

    delete X;
    delete Y;
}