[To jest pytanie partnera, aby dokładnie obliczyć prawdopodobieństwo ]

To zadanie dotyczy pisania kodu w celu dokładnego i szybkiego obliczenia prawdopodobieństwa . Wynik powinien być precyzyjnym prawdopodobieństwem zapisanym jako ułamek w najbardziej zredukowanej formie. Oznacza to, że nigdy nie powinien generować, 4/8ale raczej 1/2.

Dla pewnej dodatniej liczby całkowitej n, rozważ jednolicie losowy ciąg 1s i -1s długości ni nazwij ją A. Teraz połączmy się do Ajej pierwszej wartości. To znaczy, A[1] = A[n+1]jeśli indeksowanie od 1. Ama teraz długość n+1. Teraz rozważ także drugi losowy ciąg długości, nktórego pierwsze nwartości to -1, 0 lub 1 z prawdopodobieństwem 1 / 4,1 / 2, 1/4 każdy i nazwij go B.

Rozważmy teraz wewnętrzną produkt A[1,...,n]oraz Bi wewnętrzną produkt A[2,...,n+1]i B.

Rozważmy na przykład n=3. Możliwe wartości Ai Bmogą być A = [-1,1,1,-1]i B=[0,1,-1]. W tym przypadku dwoma produktami wewnętrznymi są 0i 2.

Twój kod musi przedstawiać prawdopodobieństwo, że oba produkty wewnętrzne są równe zero.

Kopiując tabelę wyprodukowaną przez Martina Büttnera otrzymaliśmy następujące przykładowe wyniki.

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

Języki i biblioteki

Możesz korzystać z dowolnego darmowego języka i bibliotek, które ci się podobają. Muszę być w stanie uruchomić Twój kod, więc proszę podać pełne wyjaśnienie, jak uruchomić / skompilować kod w systemie Linux, jeśli to w ogóle możliwe.

Zadanie

Twój kod musi zaczynać się od n=1i dawać poprawne dane wyjściowe dla każdego rosnącego nw osobnej linii. Powinno się zatrzymać po 10 sekundach.

Wynik

Wynik jest po prostu najwyższy nosiągnięty, zanim kod przestanie działać po 10 sekundach od uruchomienia na moim komputerze. W przypadku remisu zwycięzca najszybciej zdobędzie najwyższy wynik.

Tabela wpisów

• n = 64w Pythonie . Wersja 1 autorstwa Mitcha Schwartza
• n = 106w Pythonie . Wersja z 11 czerwca 2015 r. Autor: Mitch Schwartz
• n = 151w C ++ . Odpowiedź Port of Mitch Schwartz autorstwa kirbyfan64sos
• n = 165w Pythonie . Wersja 11 czerwca 2015 r. W wersji „przycinanie” autorstwa Mitcha Schwartza z N_MAX = 165.
• n = 945w Pythonie przez Min_25 przy użyciu dokładnej formuły. Niesamowity!
• n = 1228w Pythonie autorstwa Mitcha Schwartza przy użyciu innej dokładnej formuły (na podstawie poprzedniej odpowiedzi Min_25).
• n = 2761w Pythonie autorstwa Mitcha Schwartza przy użyciu szybszej implementacji tej samej dokładnej formuły.
• n = 3250w Pythonie przy użyciu Pypy autorstwa Mitcha Schwartza przy użyciu tej samej implementacji. Ten wynik musi pypy MitchSchwartz-faster.py |tailunikać przewijania konsoli nad głową.

Pyton

Formuła zamknięta p(n)to

Wykładnicza funkcja generująca p(n)jest

gdzie I_0(x)jest zmodyfikowana funkcja Bessela pierwszego rodzaju.

Edytuj w dniu 2015-06-11:
- zaktualizowałem kod Python.

Edytuj w dniu 13.06.2015 r .:
- dodano dowód powyższej formuły.
- naprawiono time_limit.
- dodano kod PARI / GP.

Pyton

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!;

Dowód:
ten problem jest podobny do problemu dwuwymiarowego (ograniczonego) losowego chodzenia.

Jeśli A[i] = A[i+1]można przejść od (x, y)się (x+1, y+1)[1], sposobem (x, y)[2] lub sposoby (x-1, y-1)[1] sposobem.

Jeśli A[i] != A[i+1]można przejść od (x, y)się (x-1, y+1)[1], sposobem (x, y)[2] lub sposoby (x+1, y-1)[1] sposobem.

Pozwól a(n, m) = [x^m]((x+1)^n + (x-1)^n), b(n) = [x^n](1+x)^{2n}i c(n)jest wiele sposobów, aby przejść od (0, 0)do (0, 0)z netapów.

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

Ponieważ p(n) = c(n) / 8^nmożemy uzyskać powyższą formułę zamkniętą.

Pyton

Uwaga: gratulacje dla Min_25 za znalezienie rozwiązania w formie zamkniętej!

Dzięki za interesujący problem! Można to rozwiązać za pomocą DP, chociaż nie mam obecnie dużej motywacji do optymalizacji prędkości, aby uzyskać wyższy wynik. Przydałoby się golfa.

Kod osiągnął N=39w ciągu 10 sekund na starym laptopie z Pythonem 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

Dla krotki (a,b,s,t): ajest pierwszym elementem A, bjest ostatnim elementem B, sjest wewnętrznym produktem A[:-1]i B, i tjest wewnętrznym produktem A[1:-1]i B[:-1]przy użyciu notacji plastra Python. Mój kod nie przechowuje tablic Aani Bnigdzie, więc używam tych liter, aby odwoływać się odpowiednio do kolejnych elementów, które zostaną dodane do Ai B. Ten wybór nazw zmiennych sprawia, że ​​wyjaśnienie jest trochę niewygodne, ale pozwala na ładny wygląd A*b+a*Bsamego kodu. Zauważ, że dodawany element Ajest przedostatni, ponieważ ostatni element jest zawsze taki sam jak pierwszy. Wykorzystałem sztuczkę Martina Büttnera polegającą na włączeniu się 0dwukrotnieBkandydaci w celu uzyskania właściwego rozkładu prawdopodobieństwa. Słownik X(który został nazwany Yprzez N+1) śledzi liczby wszystkich możliwych tablic według wartości krotki. Zmienne ni doznaczają licznik i mianownik, dlatego przemianowałem nazwę nzadania na as N.

Kluczową częścią logiki jest to, że można zaktualizować z Ndo N+1korzystania tylko wartości w krotce. Dwa wewnętrzne produkty określone w pytaniu są podane przez s+A*Bi t+A*b+a*B. Jest to jasne, jeśli przyjrzysz się nieco definicjom; pamiętać, że [A,a]i [b,B]są dwa ostatnie elementy tablic Ai Bodpowiednio.

Zauważ, że si tsą małe i ograniczone zgodnie z N, a dla szybkiej implementacji w szybkim języku możemy uniknąć słowników na korzyść tablic.

Możliwe jest wykorzystanie symetrii, biorąc pod uwagę wartości różniące się tylko znakiem; Nie przyjrzałem się temu.

Uwaga 1 : Rozmiar słownika rośnie kwadratowo N, gdzie rozmiar oznacza liczbę par klucz-wartość.

Uwaga 2 : Jeśli ustawimy górną granicę N, możemy przycinać krotki dla których N_MAX - N <= |s|i podobnie dla t. Można tego dokonać, określając stan pochłaniania lub pośrednio za pomocą zmiennej, która przechowuje liczbę przyciętych stanów (które trzeba będzie pomnożyć przez 8 przy każdej iteracji).

Aktualizacja : ta wersja jest szybsza:

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

Wdrożone optymalizacje:

• włóż wszystko do main()- dostęp do zmiennych lokalnych jest szybszy niż globalny
• N=1jawnie obsługuj, aby uniknąć sprawdzania (1,-1) if N else [a](co wymusza spójność pierwszego elementu w krotce podczas dodawania elementów do Arozpoczynania od pustej listy)
• rozwinąć wewnętrzne pętle, co również eliminuje mnożenie
• podwoić liczbę cdodawania 0do Bzamiast wykonywania tych operacji dwukrotnie
• mianownik jest zawsze, 8^Nwięc nie musimy go śledzić
• teraz symetrię biorąc pod uwagę: możemy rozwiązać ten pierwszy element A, jak 1i podzielić przez mianownik 2, ponieważ prawidłowych połączeń (A,B)z A[1]=1i te A[1]=-1mogą być wprowadzane do korespondencji jeden-do-jednego przez zaprzeczenie A. Podobnie możemy naprawić pierwszy element Bjako nieujemny.
• teraz z przycinaniem. Będziesz musiał się bawić, N_MAXaby zobaczyć, jaki wynik może uzyskać na twoim komputerze. Może być przepisany, aby znaleźć odpowiedni N_MAXautomatycznie przez wyszukiwanie binarne, ale wydaje się niepotrzebny? Uwaga: nie musimy sprawdzać stanu przycinania przed dotarciem dookoła N_MAX / 2, więc możemy uzyskać niewielkie przyspieszenie poprzez iterację w dwóch fazach, ale zdecydowałem się tego nie robić dla uproszczenia i czystości kodu.

Pyton

Korzystając z pomysłu losowego marszu Min_25, udało mi się dojść do innej formuły:

Oto implementacja Pythona oparta na 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()

Wyjaśnienie / dowód:

Najpierw rozwiązujemy powiązany problem liczenia, na który zezwalamy A[n+1] = -A[1]; to znaczy dodatkowy element połączony Amoże być 1lub -1niezależnie od pierwszego elementu. Nie musimy więc śledzić, ile razy A[i] = A[i+1]ma miejsce. Mamy następujący losowy spacer:

Z (x,y)możemy przejść do (x+1,y+1)[1 sposób], (x+1,y-1)[1 sposób], (x-1,y+1)[1 sposób], (x-1,y-1)[1 sposób], (x,y)[4 sposoby]

gdzie xi ystanąć na dwóch iloczyn skalarny i liczymy ilość sposobów, aby przejść od (0,0)do (0,0)w nkrokach. Liczba ta zostanie następnie pomnożona przez, 2aby uwzględnić fakt, że Amożna zacząć od 1lub -1.

Odnosimy się do pobytu w (x,y)postaci zerowym ruchu .

Powtarzamy liczbę niezerowych ruchów i, które muszą być równe, aby wrócić do (0,0). Ruchy poziome i pionowe tworzą dwa niezależne jednowymiarowe losowe spacery, które można policzyć C(i,i/2)^2, gdzie C(n,k)jest współczynnikiem dwumianowym. (W przypadku marszu z kkrokami w lewo i kkrokami w prawo istnieją C(2k,k)sposoby wyboru kolejności kroków.) Ponadto istnieją C(n,i)sposoby umieszczania ruchów i 4^(n-i)sposoby wybierania ruchów zerowych. Otrzymujemy więc:

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

Teraz musimy wrócić do pierwotnego problemu. Zdefiniuj dopuszczalną parę, (A,B)która ma być zamieniana, jeśli Bzawiera zero. Zdefiniuj parę, (A,B)która będzie prawie dopuszczalna, jeśli A[n+1] = -A[1]i dwa produkty kropkowe są równe zero.

Lemat: W danym nprzypadku prawie dopuszczalne pary są w korespondencji jeden-do-jednego z parami wymiennymi.

Możemy (odwracalnie) przekonwertować parę konwertowalną (A,B)na prawie dopuszczalną parę (A',B'), negując A[m+1:]i B[m+1:], gdzie mjest indeks ostatniego zera w B. Sprawdzenie tego jest proste: jeśli ostatnim elementem Bjest zero, nie musimy nic robić. W przeciwnym razie, gdy negujemy ostatni element A, możemy zanegować ostatni element B, aby zachować ostatni składnik iloczynu z przesuniętą kropką. Ale to neguje ostatnią wartość nie przesuniętego produktu kropkowego, więc naprawiamy to, negując element przedostatni A. Ale wtedy to odrzuca ostatnią wartość przesuniętego produktu, więc negujemy element przedostatni B. I tak dalej, aż do osiągnięcia zerowego elementu w B.

Teraz musimy tylko pokazać, że nie ma prawie dopuszczalnych par, dla których Bnie ma zera. Aby iloczyn punktowy był równy zeru, musimy mieć taką samą liczbę 1i -1warunki do anulowania. Każdy -1termin składa się z (1,-1)lub (-1,1). Zatem parzystość liczby takich -1zdarzeń jest ustalona zgodnie z n. Jeśli pierwszy i ostatni element Amają różne znaki, zmieniamy parzystość, więc jest to niemożliwe.

Więc rozumiemy

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

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

co daje powyższą formułę (ponowne indeksowanie za pomocą i' = i/2).

Aktualizacja: Oto szybsza wersja korzystająca z tej samej formuły:

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

Wdrożone optymalizacje:

• funkcja inline p(n)
• użyj rekurencji dla współczynników dwumianowych C(n,k)zk <= n/2
• użyj rekurencji dla centralnych współczynników dwumianowych

Wyjaśnienie wzoru Min_25

Min_25 opublikował świetny dowód, ale jego wykonanie zajęło trochę czasu. To jest trochę wyjaśnienia do wypełnienia między wierszami.

a (n, m) reprezentuje liczbę sposobów wyboru A, tak że A [i] = A [i + 1] m razy. Wzór na a (n, m) jest równoważny a (n, m) = {2 * (n wybierz m) dla nm nawet; 0 dla nm nieparzyste.} Dopuszczalna jest tylko jedna parzystość, ponieważ A [i]! = A [i + 1] musi wystąpić parzystą liczbę razy, aby A [0] = A [n]. Współczynnik 2 wynika z początkowego wyboru A [0] = 1 lub A [0] = -1.

Po ustaleniu liczby (A [i]! = A [i + 1]) na q (nazwanej i we wzorze c (n)), dzieli się ona na dwa losowe spacery 1D o długości q i nq. b (m) to liczba sposobów na wykonanie jednowymiarowego losowego spaceru m kroków, który kończy się w tym samym miejscu, w którym zaczął, i ma 25% szansy na ruch w lewo, 50% szansy na pozostanie w bezruchu i 25% szansy na poruszać się w prawo. Bardziej oczywistym sposobem określenia funkcji generowania jest [x ^ m] (1 + 2x + x ^ 2) ^ n, gdzie 1, 2x i x ^ 2 oznaczają odpowiednio lewą, brak ruchu i prawą. Ale potem 1 + 2x + x ^ 2 = (x + 1) ^ 2.

C ++

Tylko część (doskonałej) odpowiedzi Pythona autorstwa Mitcha Schwartza. Podstawową różnicą jest to, że użyłem 2do reprezentowania -1dla azmiennej i zrobił coś podobnego do b, który pozwolił mi używać tablicy. -O3Używam Intel C ++ z , mam N=141! Moja pierwsza wersja ma N=140.

To używa wzmocnienia. Próbowałem wersji równoległej, ale miałem problemy.

#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;
}