Kruunaneliöt

Pöydällä on n\times n-neliöön aseteltuna n^2 kolikkoa. Jokainen kolikko on joko kruuna (merk. 1) tai klaava (merk. 0). Kuviosta löytyvä alineliö on kruunaneliö, jos sen kaikki kolikot ovat kruunoja.

Esimerkki, jossa n=5 ja suurin löytyvä kruunaneliön sivu on 2.

Kuviosta etsitään suurin mahdollinen kruunaneliö. Kysymys kuuluu: kuinka monta kolikoiden konfiguraatiota on, joissa suurin löytyvä kruunaneliö on sivultaan k? Merkitään tätä lukumäärää h(n, k).

Ratkaisu

Tutkitaan kahta laskutapaa: Inkluusio-ekskluusio ja Markovin ketju. Näistä ensimmäinen toimii kun k on suuri ja jälkimmäinen, kun k on pieni (2 tai 3). Kumpikaan ei toki kovin suurille n toimi. Suurin tapaus, jonka saamme näillä keinoin kokonaan ratkaistua on n=8.

Inkluusio-ekskluusio

Lasketaan n\times n-kolikkoneliöiden lukumäärä, joista löytyy (ainakin) k\times k-kruunaneliö. Merkitään tätä lukumäärää f(n, k). (Kysytty h(n, k) saadaan sitten h(n, k) = f(n, k)-f(n, k+1).) Tätä varten tehdään määritelmä (jossa indeksointi aloitetaan nollasta)

A_{ij} = niiden neliöiden joukko, joissa on k-sivuinen kruunaneliö lähtien indeksistä (i, j), i, j  = 0, 1, \dots , n-k

Esimerkki joukosta A_{ij}. Tässä n=6 ja k=3.

Nythän f(n, k) = | \cup_{i,j} A_{ij}| ja inkluusio-ekskluusio kaavan mukaan saadaan

f(n, k) = | \bigcup_{i,j} A_{ij}| = \sum_{\emptyset \neq J \subset \{A_{ij}\} } (-1)^{|J|+1} |\bigcap_{A_{ij} \in J} A_{ij}|

Leikkauksen \bigcap_{A_{ij} \in J} A_{ij} koko saadaan laskettua selvittämällä vastaavan kolikkojen yhdisteen koko U_J = |\bigcup_{A_{ij} \in J} S_{ij}|, missä S_{ij} on kyseisen alineliön (sen, jonka vaaditaan olevan kruunaneliö) kolikot. Kaikkien tämän yhdisteen kolikkojen on siis oltava kruunoja, kun taas loput n^2 - U_J saavat olla mitä vaan, joten

|\bigcap_{A_{ij} \in J} A_{ij}| = 2^{n^2-U_J}

Joukkojen A_{ij} lukumäärän (n-k+1)^2 kasvaessa läpikäytävien joukkojen J lukumäärä 2^{(n-k+1)^2} - 1 kasvaa hyvin nopeasti. Eri kokoisia yhdisteitä U_J on kuitenkin maksimissaan n^2 kappaletta. Jos saman yhdistekoon tuottavien J-joukkojen lukumäärän laskulle olisi tehokas keino, paranisi tämä algoritmi huomattavasti. Kun yhdisteessä on r joukkoa, niin suurille r yhdisteen koko näyttäisi olevan suurinpiirtein normaalijakautunut (paitsi loppupäässä, jossa koko on rajoitettu n^2:lla; kaikkien (r = (n-k+1)^2) yhdisteen kokohan on tietysti n^2). Tätä tutkintalinjaa ei kuitenkaan jatkettu pitemmälle.

Markovin ketju

Ajatellaan että kolikkoneliö muodostetaan rivi kerrallaan. Pidetään jokaiselle sarakkeelle kirjaa, kuinka pitkä putki kruunoja siinä on tullut. Riittää rajoittaa putken pituudeksi k-1 eli sen yli ei mennä vaikka tulisikin kruuna. Tiloina on siis kaikki n-vektorit (r_j), missä r_j \in \{ 0,1,\dots, k-1 \} (k^n kappaletta) ja siirtymät näiden välillä tehdään mahdollisien rivien (2^n kappaletta) avulla. Lisäksi on yksi absorboiva tila (kutsutaan lopputilaksi), johon siirrytään, jos k-sivuinen kruunaneliö on muodustunut. Kruunaneliön muodustuminen pystytään päättelemään sarakkeiden kruunaputkista: kun uusi rivi lisätään, kasvatetaan putkia ja nyt annetaan putken kasvaa aina k:hon asti. Kruunaneliö syntyy, jos tästä putkivektorista löytyy k-putki, jossa jokainen arvo on k. Jos kruunaneliö syntyi, siirrytään lopputilaan. Jos taas ei, niin rajoitetaan putket taas korkeintaan k-1:een ja siirrytään syntyneeseen tilaan.

Tämähän ei oikeastaan ole Markovin ketju, sillä siirtymämatriisiin ei laiteta todennäköisyyksiä, vaan tämän siirtymän tuottavan rivien lukumäärä. Tämän takia esim. lopputilan itsesiirtymäänkin laitetaan arvo 2^n. Lähtötila on nollavektori ja koko kolikkoneliö saadaan lisäämällä n riviä, joten lukumäärä f(n, k) saadaan siirtymämatriisin n:nnen potenssin paikasta, joka vastaa tiloja (alkutila, lopputila).

Otetaan esimerkiksi n=3, k=2. Tällöin tilat ovat (merkitään lopputilaa (-1, -1, -1))

[(0, 0, 0), (0, 0, 1), (0, 1, 0), (0, 1, 1), (1, 0, 0), (1, 0, 1), (1, 1, 0), (1, 1, 1), (-1, -1, -1)]

Siirtymämatriisi ja sen kolmas potenssi ovat

”Siirtymä”-matriisi, joka laskee 3×3 kolikkoneliöiden, joista löytyy 2×2 kruunaneliö, lukumäärän

Koodit

Seuraavassa yllä kuvaillut algoritmit toteutettuna SageMath-koodina. Funktio f(n, k) valitsee sopivamman algoritmin riippuen parametreistä n ja k.

import itertools

def fMarkov(n, k):
    
    def hasRun(s):
        run = 0
        for x in s:
            if x==k:
                run += 1
                if run==k: return True
            else:
                run = 0
        return False
    
    rows = list(itertools.product(*[[0,1]]*n))
    states = list(itertools.product(*[range(k) for _ in range(n)]))
    startState = (0,)*n
    endState = (-1,)*n
    states.append(endState)
    stateToMatInd = {s: i for i,s in enumerate(states)}
    startI = stateToMatInd[startState]
    endI = stateToMatInd[endState]
    mat = matrix(ZZ, len(states), sparse=False)
    for i1, s1 in enumerate(states[:-1]):
        for r in rows:
            s2 = tuple( min(k, s1[j]+1) if r[j]==1 else 0  for j in range(n))
            if hasRun(s2):
                mat.add_to_entry(i1, endI, 1)
            else:
                s2 = tuple(min(k-1, x) for x in s2)
                mat.add_to_entry(i1, stateToMatInd[s2], 1)
    mat.add_to_entry(endI, endI, len(rows))
    return (mat**n)[startI][endI]

def fIE(n, k):
    squares = [(i,j) for i in range(n-k+1) for j in range(n-k+1)]
    squareSets = {(i,j): set((i+i1, j+j2) for i1 in range(k)
                             for j2 in range(k)) for i,j in squares}
    sizeForUSize = {u: 2**(n**2-u) for u in range(n**2+1) } 
    tot = 0
    for r in range(1, len(squares)+1):
        for J in itertools.combinations(squares, r):
            U = set().union(*[squareSets[t] for t in J])
            tot += (-1)**(r+1) * sizeForUSize[len(U)]
    return tot

def f(n, k):
    if k>n: return 0
    if k==n: return 1
    if k==0: return 2**(n**2)
    if k==1: return 2**(n**2)-1
    if k==2 and n<14 : return fMarkov(n, k)
    if k==3 and n<9: return fMarkov(n, k)
    return fIE(n, k)
    

import time
startTime = time.time()
n = 6
fVals = [f(n, k) for k in range(n+1)]
#print (fVals)
print ([fVals[k]-(fVals[k+1] if k<n else 0) for k in range(n+1)])
print ("Time: ", time.time()-startTime)

Tulokset

nk=0k=1k=2k=3k=4k=5k=6k=7k=8
31416941
4142175229134461
51157012721736490348633719181
6121418970800455750229271716842688863238579341
7110722441743615943411992984822521481265058911124196411584144634177322541
811968910887183791358154770750168649491709944381160397146396311672705366016837855028243223653179531300461
Lukumäärät h(n, k)

Arvo f(n, 2) (eli nxn-matriisien lukumäärä, joista löytyy (ainakin) 2×2 kruunaneliö) on OEIS-jono A255938. Funktio fMarkov() kykenee laskemaan vielä termin

f(13, 2) = 747083368584493829509302940089141575261663191900672

Tehtäviä

  1. Tutkitaan nyt 3 \times n kolikkosuorakaidetta. Olkoon lisäksi jokainen kolikko painotettu, niin että kruunan todennäköisyys on p.
    a) Olkoon p = \frac{1}{2}. Mikä on ensimmäinen n, jolle todennäköisyys, että suorakaiteesta löytyy 3 \times 3 kruunaneliö, on yli 0.5?
    b) Mille arvolle p todennäköisyys 3 \times 3 kruunaneliön löytymiselle 3 \times 5 suorakaiteesta on 0.25?
    c) Oletetaan, että saat aluksi valita suorakaiteen pituuden n. Mutta tämän myötä (kaikkien kolikoiden kruuna-tn) p muuttuu kaavan p = \max(0, \frac{2}{3} - \frac{n}{1000}) mukaan. Miten n kannattaa valita, jotta 3 \times 3 kruunaneliön löytymisen todennäköisyys maksimoituu?
  2. Inkluusio-ekskluusio koodissa yhdiste lasketaan jokaiselle J-joukolle uudestaan. Tätä voisi parannella: generoidaan J-joukot itse ja pidetään aina siihen mennessä syntynyt yhdiste muistissa. Ks. myös OEIS-linkin kaavaa (eikö se tosin ole hyvin epätehokas, koska siinä on vielä tulo J-joukon potenssijoukon yli??)

Suorakaide ja ympyrä kolmiossa

Olkoon annettu kolmio, jonka kantakulmat ovat terävät. Kolmion sisälle on asetettava akselien suuntainen suorakaide, jonka yläpuolelle (ja myös kolmion sisälle) ympyrä. Tehtävänä on maksimoida suorakaiteen ja ympyrän alojen summa.

Asetetaan kolmio koordinaatistoon seuraavalla tavalla ja tehdään kuvassa olevat merkinnät.

Merkitään lisäksi kolmion kannan OB pituutta b:llä ja sivujen pituuksia OA = s_1, AB = s_2 sekä kolmion piiriä p = b+s_1+s_2. Merkitään vielä kärjen A x-koordinaattia a:lla (sen y-koordinaattihan on korkeus h).

Ratkaistaan suorakaiteen x-rajat x_1 ja x_2. Yhdenmuotoisista kolmioista saadaan

\frac{x_1}{y} = \frac{a}{h} ja \frac{b-x_2}{y} = \frac{b-a}{h}

Näin ollen suorakulmion vaakasivu on x_2-x_1 = b(1-\frac{y}{h}).

Ratkaistaan sitten ympyrän säde. Maksimitapauksessahan ympyrä on alkuperäisen kolmion sivujen ja suorakaiteen yläsivun rajoittaman kolmion K_y, sisään piirretty ympyrä, joten käytetään sille löytyvää kaavaa

r^2 = \frac{(s-a_1)(s-b_1)(s-c_1)}{s},

missä s on K_y:n piirin puolikas ja a_1, b_1, c_1 sen sivut. Yksi sivu on suorakaiteen vaakasivu ja toiset saadaan alkuperäisen kolmion sivuista kertomalla suhdeluvulla \frac{h-y}{h} (käytetään taas yhdenmuotoisia kolmioita). Jokaisessa termissä on siis kertoimena tämä suhdeluku. Nimittäjään jää sen neliö. Siis

r^2 = \frac{(p/2-s_1)(p/2-s_2)(p/2-b)(1-y/h)^2}{p/2}

= \frac{p/2(p/2-s_1)(p/2-s_2)(p/2-b)1/h^2(h-y)^2}{(p/2)^2}

Toisaalta Heronin kaavan mukaan kolmion ala on \sqrt{\frac{p}{2} (\frac{p}{2}-s_1)(\frac{p}{2}-s_2)(\frac{p}{2}-b)}, joten saamme

r^2 =  \frac{b^2}{p^2} (h-y)^2

Näin saamme suorakaiteen ja ympyrän yhteiselle alalle funktion

f(y) = b(1-\frac{y}{h})y + \pi \frac{b^2}{p^2} (h-y)^2

Toisen asteen kerroin on \frac{\pi b^2}{p^2} - \frac{b}{h} <0, (ks. Tehtävä 1). Näin ollen paraabeli f on alaspäin aukeava. Derivaatan nollakohdasta löydämme maksimikohdan

y_{max} = \frac{h}{2} \left(  1- \frac{\pi bh}{p^2-\pi bh} \right).

Tehtäviä

  1. Osoita, että kolmiolle, jonka kanta on b, korkeus h ja piiri p pätee \pi \frac{b^2}{p^2} < \frac{b}{h} .
  2. Osoita että edellä vasen puoli voidaan vielä kertoa 2:lla ja epäyhtälö säilyy. Kenties vielä kolmella… Mikä on suurin luku, jolla epäyhtälö vielä on voimassa?

Rajoitettu vaihteluväliset kävelyt

Tutkitaan n:n pituisia jonoja symboleista -1 ja 1. Kun tällainen mielletään kävelyksi, määräytyy sille kunkin askeleen jälkeen korkeus, jolla kävely on (lähdetään nollasta). Kävelylle w lasketaan vaihteluväli h erotuksena h(w) = \max_j h_j(w) - \min_j h_j(w), missä h_j on kävelyn korkeus j askeleen jälkeen.

Kaksi esimerkkikävelyä ja niiden vaihteluvälin lasku. Tässä n=7.

Kysymys: kuinka monta kävelyistä on vaihteluväliltään korkeintaan r, missä r \geq 0 on annettu parametri.

Ratkaisu

Tehdään ensin kokeilua varten algoritmi, joka käy kaikki kävelyt (s.o. -1, 1 -jonot) läpi. Kielenä SageMath

def brw(r, n):
    count = 0
    for w in cartesian_product([(-1,1)]*n):
        ss = [sum(w[:j]) for j in range(n+1)]
        if max(ss)-min(ss) <= r: count += 1
    return count

for r in range(7): print ([brw(r, n) for n in range(10)])

Huomataan, että arvoille r = 2, r = 3, r = 4 ja r = 5 jono löytyy OEIS:sta. Sieltä löytyy generoivat funktiot ja tapaukselle r=2 kaava

brw(2, n) = 2^{\lfloor \frac{n+2}{2} \rfloor} + 2^{\lfloor \frac{n+1}{2} \rfloor} - 2

Myös tapauksille r=3,4 on OEIS:ssa lähes yhtä siistit kaavat (niissä esiintyy Fibonacci-luvut).

Koska emme kuitenkaan keksi tästä yleistä kaikille r toimivaa kaavaa, ryhdytään kehittämään tehokkaampaa algoritmia. Kun jonoa w luetaan alusta lähtien ja tarkastellaan korkeutta, jolla ollaan askeleen j jälkeen, niin mieleen tulee tietenkin Markovin ketju suunnattu verkko, jossa tilana solmuina on h_j. Tämä ei kuitenkaan itsessään riitä, vaan tilaan tarvitaan myös tieto siitä missä maksimi- ja minimikorkeuksissa on käyty siihen mennessä. Siirtymät tietenkin tehdään vain sallituissa rajoissa. Seuraava Sage-koodi tekee tämän ja laskee kävelyjen määrän siirtymämatriisin n. potenssista (summataan lähtötilaa vastaava rivi).

def brw(r, n):
    bds = [(i,j) for i in range(-r, 1) for j in range(r+1) if j-i<=r]
    states = [] #triples (i, j, k) where (i,j) are the min, max -heights so far and k is current height
    for (i,j) in bds:
        states += [(i,j,k) for k in range(i, j+1)]
    matN = len(states)
    startState = states.index((0,0,0))
    arr = [[0]*matN for _ in range(matN)]
    for ind1, (i,j,k) in enumerate(states):
        i2 = i-1 if k-1<i else i
        if j-i2<=r:
            s2 = (i2, j, k-1)
            arr[ind1][states.index(s2)] += 1
        j2 = j+1 if k+1>j else j
        if j2-i<=r:
            s2 = (i, j2, k+1)
            arr[ind1][states.index(s2)] += 1
    return sum((Matrix(arr)**n)[startState])

Tilojen määrä kuitenkin kasvaa nopeasti, kun r kasvaa (O(r^2)), joten tutkitaan verkon rakennetta ja huomataan siitä laskemista helpottavia rekursioita.

Piirretään verkko koodilla

def getBrwMat(r):
    bds = [(i,j) for i in range(-r, 1) for j in range(r+1) if j-i<=r]
    states = []
    for (i,j) in bds:
        states += [(i,j,k) for k in range(i, j+1)]
    matN = len(states)
    ret = [[0]*matN for _ in range(matN)]
    for ind1, (i,j,k) in enumerate(states):
        i2 = i-1 if k-1<i else i
        if j-i2<=r:
            s2 = (i2, j, k-1)
            ret[ind1][states.index(s2)] += 1
        j2 = j+1 if k+1>j else j
        if j2-i<=r:
            s2 = (i, j2, k+1)
            ret[ind1][states.index(s2)] += 1
    return Matrix(ret), states


myR = 3
a, states = getBrwMat(myR)
startState = states.index((0,0,0))
g = DiGraph(a)
vColors = {}
heights = {}
for r,color in enumerate(['blue', 'purple', 'pink', 'orange', 'red']):
    vsOfH = [j for j,s in enumerate(states) if s[1]-s[0]==r]
    vColors[color] = vsOfH
    heights[r] = vsOfH
g.show(figsize=8, vertex_size=70, vertex_labels=None,
       vertex_colors=vColors, iterations=700,
      layout='ranked', heights=heights)
Tutkittava verkko, kun r=3. Lähtötila on alhaalla.

Verkon rakenteellahan on selvä merkitys: ”putket” ovat sen hetkisiä (min, max) -rajoja. Putkessa voi ”seilata” edes takaisin ja päästä siirtyä aina seuraavaan putkeen eli siirtää (min, max) -rajoja yhdellä. Kunnes max-min = r, jolloin vaihteluväli ei saa enää kasvaa, ja jäädään seilaamaan siihen r+1-putkeen, johon päädyttiin. Huom. ”putki” on itse asiassa verkkona (suuntaamaton) polku.

”Putki” eli polkuverkko

Keskitytään aluksi yhteen m-putkeen ja sen kävelyihin eli yllä mainittuun ”seilaamiseen”. Merkitään p(m, j, n):lla ensimmäisestä päästä (j=0) lähtevien kävelyjen määrää m pituisessa putkessa, jotka päätyvät putken solmuun j. Huomaa, että tämä on tietysti symmetrinen, jos putki käännetään.

Tutkitaan sitten koko verkkoa ja tehdään hieman yleisempi määritelmä: brwRec(m, r, n) = kävelyiden määrä verkolla, jossa ensimmäinen putki on m-pituinen ja lähdetään tämän putken (jommasta kummasta) päästä. Huomaa, että haluttu lukumäärä on tällöin erikoistapaus brw(r, n) = brwRec(1, r, n). Nyt kävely voi seilata ensimmäisessä putkessa jonkun aikaa ja siirtyä sitten, päätyessään jompaan kumpaan päätyyn, seuraavalle tasolle, jolle jää loput askeleet (-1 sillä myös tasolle siirtymä kuluttaa yhden askeleen). Kävely voi myös olla siirtymättä seuraavalle tasolle ja päätyä lopuksi mihin tahansa putken somuun. Saadaan siis rekursioyhtälö

\text{brwRec}(r, n) \\= \sum_{k=0}^{n-1} (p(m, 0, k) + p(m, m-1, k)) \text{brwRec}(m+1, r-1, n-1-k) \\+ \sum_{j=0}^m p(m, j, n)

Koodataan nämä rekursiot ja käytetään memoisatiota:

pMemo = {}
def p(m, j, n):
    if n==0: return 1if j==0 else 0
    memoKey = (m,j,n)
    if memoKey in pMemo:
        return pMemo[memoKey]
    ret = 0
    if j>0: ret += p(m, j-1, n-1)
    if j<m-1: ret += p(m, j+1, n-1)
    pMemo[memoKey] = ret
    return ret

brwMemo = {}
def brwRec(m, r, n):
    if n==0: return 1
    memoKey = (m,r,n)
    if memoKey in brwMemo:
        return brwMemo[memoKey]
    ret = 0
    if r>0:
        for j in range(n):
            ret += (p(m, 0, j) + p(m, m-1, j)) * brwRec(m+1, r-1, n-1-j)
    ret += sum(p(m, j, n) for j in range(m)) #stay on the first path
    brwMemo[memoKey] = ret
    return ret

brw = lambda r,n: brwRec(1, r, n)

Tämä algoritmi toimii jo kohtuullisessa ajassa suurehkoillekin r:n arvoille, mutta tässä ongelmaksi tulee suuret n:n arvot. Matriisimenetelmässä taas suuret n:n arvot eivät tuottaneet ongelmaa, koska matriisin potenssiinkorotus on nopeaa (O(log n)). Koitetaan saada molempien maailmojen paremmat puolet ja käytetään keskikokoisille r:n ja n:n arvoille rekursiivista metodia ja pienille r:n arvoille matriisimenetelmää. Arvolle r=2 voidaan käyttää myös valmista kaavaa. Koska ratkaisujen lukuarvot kasvavat nopeasti hyvin suuriksi, annetaan testin vuoksi vastaukset jossain suuressa modulossa, esim 7347249713. Tehdään laskenta siis renkaassa \mathbb{Z} / 7347249713\mathbb{Z}. Tämä onnistuu Sage-koodilla

MOD = 7347249713
MY_RING = Integers(MOD)
ZERO = MY_RING(0)
ONE = MY_RING(1)
TWO = MY_RING(2)

pMemo = {}
def p(m, j, n):
    if n==0: return ONE if j==0 else ZERO
    memoKey = (m,j,n)
    if memoKey in pMemo:
        return pMemo[memoKey]
    ret = ZERO
    if j>0: ret += p(m, j-1, n-1)
    if j<m-1: ret += p(m, j+1, n-1)
    pMemo[memoKey] = ret
    return ret

brwMemo = {}
def brwRec(m, r, n):
    if n==0: return ONE
    memoKey = (m,r,n)
    if memoKey in brwMemo:
        return brwMemo[memoKey]
    ret = ZERO
    if r>0:
        for j in range(n):
            ret += (p(m, 0, j) + p(m, m-1, j)) * brwRec(m+1, r-1, n-1-j)
    ret += sum(p(m, j, n) for j in range(m)) #stay on the first path
    brwMemo[memoKey] = ret
    return ret

def brwMatrixMethod(r, n):
    bds = [(i,j) for i in range(-r, 1) for j in range(r+1) if j-i<=r]
    states = []
    for (i,j) in bds:
        states += [(i,j,k) for k in range(i, j+1)]
    matN = len(states)
    startState = states.index((0,0,0))
    arr = [[0]*matN for _ in range(matN)]
    for ind1, (i,j,k) in enumerate(states):
        i2 = i-1 if k-1<i else i
        if j-i2<=r:
            s2 = (i2, j, k-1)
            arr[ind1][states.index(s2)] += 1
        j2 = j+1 if k+1>j else j
        if j2-i<=r:
            s2 = (i, j2, k+1)
            arr[ind1][states.index(s2)] += 1
    return sum( (Matrix(MY_RING, arr)**n)[startState])


def brw(r, n):
    if r==2: return (TWO**((n+3)/2) if n%2 else 3*TWO**(n/2)) - 2
    #if r==3: TODO
    if r<5: return brwMatrixMethod(r, n)
    return brwRec(1, r, n)

testCases = [
    (1, 3),
    (2, 3),
    (2, 4),
    (2, 7),
    (3, 5),
    (3, 12),
    (4, 15),
    (5, 16),
    (2, 20),
    (2, 109),
    (3, 67),
    (3, 174),
    (11, 23),
    (12, 43),
    (15, 38),
    (21, 45),
    (3, 923),
    (2, 242398),
    (4, 923847),
    (2, 12345678910),
]

for (r, n) in testCases:
    print (r, n, brw(r,n))

Tämä tulostaa testitapauksiin vastaukset

"""
1 3 2
2 3 6
2 4 10
2 7 30
3 5 26
3 12 1028
4 15 12328
5 16 37144
2 20 3070
2 109 868658622
3 67 1447052816
3 174 1595158214
11 23 8099490
12 43 1135120293
15 38 2365381958
21 45 5576619597
3 923 1918313630
2 242398 3965505254
4 923847 2261069692
2 12345678910 4322229703
"""

Tehtäviä

  • Yllä, arvolle r=2 (kakkosen potensseja sisältävä) kaava toimi myös kun laskenta tehtiin modulossa. Pohdi kuinka se onnistuisi muille kaavoille (ts. jos potenssiin korotettava luku ei ole kokonaisluku).
  • Perustele tapauksen r=2 kaava 2^{\lfloor \frac{n+2}{2} \rfloor} + 2^{\lfloor \frac{n+1}{2} \rfloor} - 2.
    Vinkki: Käytä edellä ollutta rekursiota. Huomaa, että 3-putken kävelyt ovat muotoa askel keskelle ja takaisin jompaan kumpaan päätyyn jne. Käytä geometrisen summan kaavaa.
  • Laske verkon spektri. On tunnettua, että polun spektri on \{2\cos \frac{j\pi}{n+1}, j=1,\dots, n\}. Onko näillä joku yhteys (onhan verkon osat polkuja)?
  • Keksi kaava yleiselle r.
  • Tutki vaihteluvälin jakaumaa n:n pituisille {-1, +1}-kävelyille. Kuinka sen keskiarvo käyttäytyy?

Etäisyyksien tulo on 1

Olkoon pisteet A=(1,0), B=(\cos \alpha, \sin \alpha) ja C=(\cos \alpha, -\sin \alpha). Määritellään funktio f:\mathbb{R}^2 \to [0, \infty), f(Z) = tulo pisteen Z etäisyyksistä pisteisiin A, B ja C.

Funktion f tasa-arvo käyrä {f=1} sivuaa ympyrää pisteessä D.

Selvitä parametrin \alpha arvo, jolla yo. kuvan tilanne toteutuu eli f:n tasa-arvo käyrä \{f=1\} sivuaa yksikköympyrää oikeassa puolitasossa.

Ratkaisu

Kompleksitulkinnat

Tulkitaan taso kompleksitasoksi, jolloin A=1, B=e^{i\alpha} ja C=\bar B = e^{-i \alpha}. Tutkitaan funktiota g(t) = |f(e^{it})|^2-1, t\in (0, \alpha). (Huom: etäisyyden neliö on 1 joss etäisyys on 1.) Jotta sivuaminen tapahtuu pisteessä e^{it_0}, täytyy arvon t_0 \in (0, \alpha) olla funktion g lokaali maksimi ja g(t_0)=0 sekä g'(t_0)=0. Selvitetään ensin funktion g lauseke:

g(t) + 1\\= |e^{it}-1|^2 |(e^{it}-e^{i\alpha})(e^{it}-e^{-i\alpha})|^2 \\= (\sin^2 t + (1-\cos t)^2)|e^{2it}-2e^{it}(e^{i\alpha}+e^{-i\alpha})+e^{i\alpha}e^{-i\alpha}|^2 \\= 2(1-\cos t)|e^{2it}-2e^{it}\cos(\alpha)+1|^2  \\= 2(1-\cos t)|e^{it}(e^{it}-2\cos(\alpha)+e^{-it})|^2  \\= 2(1-\cos t)|2\Re(e^{it})-2\cos(\alpha)|^2  \\= 8(1-\cos t)(\cos t - \cos(\alpha))^2

Merkitään a = \cos \alpha. Nyt voidaan derivoida

g'(t) = \sin t (\cos t - a)(\cos t - a - 2(1-\cos t))

Koska t_0 \in (0, \alpha), on a = 3\cos t_0 - 2, joka sijoitettuna yhtälöön g(t_0) = 0 antaa

(1-\cos t_0)(2-2\cos t_0)^2 - \frac{1}{8} = 0

jonka ratkaisuna saadaan \cos t_0 = 1-2^{-5/3} ja näin ollen a = 1-\frac{3}{2  \sqrt[^3]{4} }.

Siis \cos \alpha = 1-\frac{3}{2 \sqrt[^3]{4} } \approx 0.055 ja \sin \alpha = \sqrt {\frac{3}{\sqrt[^3]{4}} - \frac{9}{8\sqrt[^3]{2}}} \approx 0.998.

Tehtäviä

  • Yleistys useammalle pisteelle. Kuinka parametri-pisteet täytyy valita? Tasavälein valinta ei toimi, kuten ao. kuvasta neljälle pisteelle nähdään.
  • undefined
  • Osoita, että valittiinpa kuinka monta pistettä tahansa (äärellisen monta kuitenkin) ja asetettiinpa ne kuinka tahansa, niin ympyrän kehältä löytyy aina piste, jossa f>1. (Määritelmässä otetaan tietenkin tulo kaikkien parametri-pisteiden yli.)
  • Kuinka suuri osa ympyrän kehästä voi korkeintaan olla joukon \{f \leq 1\} sisällä ja millä parametri-pisteiden (n kpl) asettelulla tämä saavutetaan? Onko kolmelle pisteelle alkuperäisessä tehtävässä ratkaistu tapaus optimaalinen?

Integraali 1/(1+x^a)

Tehtävä

Olkoon a>1. Laske integraali

I_a = \int_0^{\infty}  \frac{1}{1+x^a} dx

Ratkaisu

Rationaaliluku-tapaus

Lasketaan ensin I_{\frac{m}{n}} rationaaliluvulle \frac{m}{n}. Oletetaan, että \gcd (m, n) = 1 ja m>n\geq 1. Käytetään apuna kompleksista integrointia ja määritellään joukossa \{\Re z>0\} funktio

f(z) = \frac{1}{1+z^{\frac{m}{n}}}

Selvitetään f:n navat. Niitä ovat pisteet joille z^{\frac{m}{n}} = -1. On oltava |z|=1, joten z=e^{i\theta} jollekin reaaliselle \theta. Jotta z^{\frac{m}{n}} = -1, on oltava i\theta\frac{m}{n} = \pi i (2j+1) eli yhtäpitävästi \theta = \pi(2j+1)\frac{n}{m}. Näin ollen navat ovat z_j = e^{i\pi(2j+1)\frac{n}{m}}, j=0,1,\dots,m-1.

Positiivisella reaaliakselilla funktio f on juuri haluttu integrandi. Lisäksi huomataan, että puolisuoralla te^{2 \pi i \frac{n}{m}}, t \geq 0 pätee f(z) = \frac{1}{1+t^{\frac{m}{n}}}. Nämä huomiot johtavat siihen että integrointi poluksi kannattaa valita (ks. Huomioita-osion tarkennus):

Polku koostuuu kolmesta osasta: jana origosta R:ään, ympyrän kaari R:stä Re^{2 \pi i n/m}:ään, josta lopuksi jana takaisin origoon. Interaktiivinen versio: https://www.desmos.com/calculator/nwrl7z0311

Navoista yksi, z_0, on polun \gamma sisällä ja polku kiertää sen kerran vastapäivään. Näin ollen Residylauseen mukaan

\int_{\gamma} f(z) dz = 2\pi i Res(f, z_0)

Koska kaarella \gamma_2, jossa z=Re^{it}, pätee |f(z)| = \frac{1}{|1+(Re^{it})^{\frac{m}{n}}|} \leq \frac{1}{|R^{\frac{m}{n}}-1|}, ja kaaren pituus on korkeintaan 2\pi R, niin \int_{\gamma_2} f(z)dz \to 0, kun R \to \infty, sillä \frac{m}{n} > 1.

Kuten tuli jo mainittua, \int_{\gamma_1} f(z)dz = I_{\frac{m}{n}}.

Tehdään merkintöjä helpottava huomio: e^{2\pi i \frac{n}{m}} = z_0^2. Nyt, koska polun \gamma_3 vastapolulle \hat \gamma_3 [0, R]->\mathbb{C}, \hat \gamma_3(t) = tz_0^2 derivaatta on vakio z_0^2, saamme aiemman huomion, että f tuottaa halutun integrandin puolisuoralla tz_0^2, nojalla että \int_{\gamma_3} f(z) dz = -z_0^2 I_{\frac{m}{n}}.

Siis \lim_{R\to \infty} \int _{\gamma} f(z) dz = (1-z_0^2) I_{\frac{m}{n}}. Nyt tarvitsee enää laskea Res(f, z_0). Tätä varten kirjoitetaan f(z) muodossa

f(z) = \frac{1}{z-z_0} \frac{z-z_0}{1+z^{\frac{m}{n}}}.

Nyt l’Hopitalin säännön nojalla ((1+z^{\frac{m}{n}})' ei ole 0 z_0:ssa, joten voidaan solvetaa myös kompleksifunktioille)

\lim_{z\to z_0} \frac{z-z_0}{1+z^{\frac{m}{n}}} = \frac{1}{\frac{m}{n}z_0^{\frac{m}{n}-1}} = \frac{n}{m}e^{-i\pi\frac{n}{m}\frac{m-n}{n}} = \frac{n}{m}e^{i\pi\frac{n-m}{m}} = \frac{n}{m}e^{i\pi\frac{n}{m}}e^{-i \pi} = -\frac{n}{m}z_0.

Siispä Res(f, z_0) = -\frac{n}{m}z_0. Olemme täten ratkaisseet

I_{\frac{m}{n}} = -2\pi i \frac{n}{m}z_0 \frac{1}{1-z_0^2} = \frac{n \pi}{m} \frac{2i}{(z_0^2-1)\bar z_0} = \frac{n \pi}{m} \frac{1}{\Im z_0} = \frac{n \pi}{m \sin (\frac{n\pi }{m})}.

Yleinen tapaus

Tapauksen a=\frac{m}{n} ratkaisun muoto I_a = \frac{\pi}{a \sin (\frac{\pi}{a})} antaa vihjeen, että kaava pätisi myös kaikilla reaaliluvuilla a>1. Tämä pitääkin paikkansa. Tehdään tätä varten itse asiassa vielä yleisemmin kaikille kompleksiluvuille w=\sigma + it, kun \sigma > 1 (olkoon tämä joukko U), määritelmä

\phi(w) = \int_0^\infty \frac{1}{1+x^w} dx

Merkitään g(x, w) = \frac{1}{1+x^w}. Määritelmä on järkevä, sillä integraali suppenee: |g(x,w)| \leq (|x^w| - 1)^{-1} \leq 2|x^w|^{-1} = 2x^{-\sigma} , kun x on suuri. Lisäksi w\mapsto g(x,w)=\frac{1}{1+ \exp(\log(x)w)} on holomorfinen U:ssa kaikilla x, sillä \exp(\log(x)w) = x^\sigma \exp(it\log(x)) \neq -1, kun \sigma > 1.

Nyt kaikille kompakteille K \subset U löytyy \sigma_0 > 1, siten että \Re w \geq \sigma_0 kaikilla w\in K. Näin ollen \int_K \int_0^\infty |g(x,w)| dx dw < \infty ja Fubini-Tonellin lauseen mukaan integrointijärjestyksen voi vaihtaa. Sovelletaan tätä suljetulle suoristuvalle käyrälle \gamma: [0,1] \to U, jolloin

\int_\gamma \phi(w)dw = \int_\gamma \int _0^\infty g(x,w) dx dw = \int _0^\infty \int_\gamma g(x,w) dw dx = \int _0^\infty 0 dx = 0

Tässä \int_\gamma g(x,w) dw = 0 Cauchyn-Goursaut’n lauseen mukaan, sillä U on yhdesti yhtenäinen. Nyt Moreran lauseen mukaan \phi on holomorfinen U:ssa (\phi:n jatkuvuus on ilmeistä).

Siis \phi: U \to \mathbb{C} on holomorfinen funktio. Mutta osoitimme rationaaliluku-tapauksessa, että \phi(w) = \frac{\pi}{w\sin (\frac{\pi}{w})}, kun w \in \mathbb{Q}. Koska \frac{\pi}{w\sin (\frac{\pi}{w})} on myös U:ssa holomorfinen (ks. Huomiot) ja \mathbb{Q}:lla on U:ssa kasautumispiste, niin Identtisyyslauseen mukaan funktiot ovat yhtäsuuret kaikilla w \in U, erityisesti kaikilla reaaliluvuilla a>1.

Huomioita

  • Ei-kokonaisluku eksponentti on hieman ongelmallinen origossa, koska haaroittuminen. Tästä selvitään, kun rajoitetaan integraali ensin alhaalta ja otetaan sitten raja-arvo. Polussa tehdään origon lähellä pieni \epsilon-kaari.
  • Kaikki funktion \sin juuret ovat reaalisia: \sin(z) = \frac{\exp(iz)-\exp(-iz)}{2i}=0 \implies e^{iz} = e^{-iz} \implies iz = -iz+2k \pi i \implies \Re (iz)=0 \implies \Im z=0.
  • Edellisen perusteella ratkaisun kaavassa nimittäjässä oleva \sin(\frac{\pi}{w}) on holomorfinen U:ssa (kun w \in (1, \infty), niin \frac{\pi}{w} \in (0, \pi), jossa sinillä ei ole juurta).
  • Koska \mathbb{Q} \subset \mathbb{R} on tiheä, niin yleinen tapaus seuraa jo pelkästä \phi:n jatkuvuudesta. Mutta yo. tavalla tulos saatiin osoitettua kaikille kompleksiluvuille puolitasossa \Re w > 1.

Tehtäviä

  • Laske integraali \int_0^\infty \frac{x^b}{1+x^a}dx, kun 0 \leq b<a-1. Vinkki: tee sopiva muuttujanvaihdos, joka johtaa ratkaistuun tapaukseen.
  • Perustele \phi:n jatkuvuus. Vinkki: määrittele \phi_n(w) = \int_0^\infty \frac{dx}{1+x^{w+h_n}}, missä h_n \to 0 ja käytä Dominoidun konvergenssin lausetta.

Sanan osajonojen lukumäärä

Sanan osajono on sanasta järjestyksessä poimituista merkeistä muodostettu uusi sana. Esimerkiksi sanalla ’ABBA’ on osajonot (11 kpl) [”, ’A’, ’AA’, ’AB’, ’ABA’, ’ABB’, ’ABBA’, ’B’, ’BA’, ’BB’, ’BBA’]. Huomaa, että esim. ’AB’ saadaan kahdella eri tavalla, mutta lasketaan vain kerran.

Kuinka laskea sanan w osajonojen lukumäärä f(w)? Oletetaan että w:ssä esiintyvien merkkien lukumäärä on jokin ennalta kiinnitetty pienehkö luku. Halutaan löytää lineaariaikainen algoritmi.

Ratkaisu

Olkoon w sana, jonka osajonot halutaan laskea. Olkoon w:n pituus n ja siinä esiintyvät merkit \{c_1, c_2, \dots, c_m\}. Kun samaistetaan jokainen w:n osajono v siihen indeksien osajonoon, josta v löytyy kaikkein aikaisimmin, vastaa jokainen sanan osajono solmusta -1 lähtevää kävelyä suunnatulla graafilla G, jonka solmut ovat V(G) = \{-1, 0,1,\dots, n-1 \} ja solmusta j lähtee kaaret (j, k_c) missä k_c on indeksi, jossa merkki c esiintyy w:ssä seuraavan kerran lukien indeksistä j+1 lähtien (ja jos c ei enää esiinny, niin kaarta ei ole) ja tämä vastaavuus on bijektiivinen.

Esimerkiksi sanalle w = ”ABBA” graafi näyttää seuraavalta

Sanan ”ABBA” graafi. theFolls on dicti, jonka ao Python-koodi antaa, koodissa indeksit ovat siirretty yhdellä ja itse asiassa theFolls antaa seuraavan etsimisen alkaa itse siitä indeksistä missä ollaan.

Olkoon A edellä kuvatun suunnatun graafin vierekkäisyysmatriisi. Osajonojen lukumäärä on tällöin ensimmäisen rivin summa seuraavasta matriisista

\begin{aligned}  \sum_{j=0}^\infty A^j  = \left(I - A \right)^{-1}   \end{aligned}

Huomaa, että edellä summa on itse asiassa vain äärellinen, sillä kaaret vievät aina eteen päin, joten graafin kävely voi olla korkeintaan n:n pituinen ja A^{n+1}=0, joten summan suppenemisen eikä I-A:n kääntyvyyden kanssa tule ongelmia.

Merkitään B=I-A. Mehän tarvitsemme vain sen ensimmäisen rivin summan, joten koko matriisia A ei tarvitse muodostaa, eikä summaa \sum_{j=0}^{n}A^j laskea, vaan riittää ratkaista

B^{T}x = e_1,

missä e_1 on sarakevektori, jonka ensimmäinen komponentti on 1 ja loput nollia. Kun tämä kirjoitetaan Gauss-Jordanin menetelmän vaatimaan muotoon (jossa ei ratkaista koko käänteismatriisia, vain sen ensimmäinen rivi), niin huomataan, että ratkaiseminen on erittäin helppoa, sillä B on alakolmiomatriisi, jossa on ykkösiä diagonaalilla ja jokaisessa sarakkeessa korkeintaan m:ssä paikassa -1 (ja muuten nollia). Näin ollen saadaan seuraava lineaariaikainen algoritmi (Python):

def theFolls(w):
    """map c -> list a where a[j] is the next index of c in w
    when start searching from the index j"""
    alph = list(set(c for c in w))
    n = len(w)
    ret = {c: [n]*n for c in alph}
    prevs = {c: -1 for c in alph}
    for i, c in enumerate(w):
        for j in range(prevs[c]+1, i+1):
            ret[c][j] = i
        prevs[c] = i
    return ret

def countSubsequences(w):
    """How many words can be obtained from w as subsequences"""
    n = len(w)
    folls = theFolls(w)
    edges = [
        [folls[c][i]+1 for c in folls if folls[c][i]<n]
        for i in range(n)
    ]
    ret = 0
    a = [0 for _ in range(n+1)]
    a[0] = 1
    for i in range(n):
        for j in edges[i]:
            a[j] += a[i]
    return sum(a)

Junarata-tehtävä 2

Tiedosto fenwickTree.py

class FenwickTree:

    def __init__(self, n):
        self.n = n+1
        self.arr = [0]*self.n

    def incr(self, ind, val):
        """Increment the value at index ind by val"""
        if ind<0: raise Exception("Index out of bounds")
        curr = ind+1
        while curr<self.n:
            self.arr[curr] += val
            curr += curr&(-curr)

    def getSum(self, upto):
        """Get the prefix sum up to index upto"""
        s = 0
        ind = upto+1
        while ind>0:
            s += self.arr[ind]
            ind -= ind&(-ind)
        return s

Tiedosto junarata6.py

from fenwickTree import FenwickTree

class IntvalSet:
    """Set where the included natural numbers
        are given as intervals [a, b]"""
    def __init__(self):
        self.intvals = []

    def has(self, x):
        for (a,b) in self.intvals:
            if a<=x and x<=b: return True
        return False

    def add(self, a, b):
        self.intvals.append((a, b))

    def getAll(self, lowerBdd=0):
        """generates all numbers in the set.
            If lowerBdd is given,
            generates only numbers >= lowerBdd"""
        for (a,b) in self.intvals:
            for x in range(max(a, lowerBdd), b+1):
                yield x

    def __str__(self):
        return "[]" if len(self.intvals)==0 else " U ".join(["[%d, %d]" %x for x in self.intvals])


def getData(num):
    ret = []
    with open("junar"+str(num)+".in") as f:
        ret = f.read().split("\n")
    return [int(x) for x in ret[1:] if len(x)>0]

def getPoss(p):
    """Get list of positions of elements in the
        permutation p of [1..n]
        e.g. getPoss([2,5,3,1,4]) = [None, 3, 0, 2, 4, 1]
    """
    ret = [None]*(len(p)+1)
    for i,x in enumerate(p):
        ret[x] = i
    return ret

"""
Let's denote types of orderings of k-1, k and k+1
(or k and k+1; k-1 and k in cases k=1; k=n  like this:

123: 0
132: 1
213: 2
231: 3
312: 4
321: 5
_12: 6
_21: 7
12_: 8
21_: 9
"""
def getTyyppi(a, b, c):
    if a<b:
        if b<c: return 0
        if a<c: return 1
        return 3
    else: #b<a
        if b>c: return 5
        if a<c: return 2
        return 4


def getNumOfRounds(poss):
    """How many rounds does it take to gather the numbers.
        @param poss: index-positions of the permutation
    """
    curr = 1
    prevInd = poss[1]
    r = 1
    for curr in range(1, len(poss)):
        currInd = poss[curr]
        if currInd<prevInd: r += 1
        prevInd = currInd
    return r


def countScore2Pairs(arr):
    n = len(arr)
    poss = getPoss(arr)
    toConsider = []
    t = None
    tempSet = None
    for i, x in enumerate(arr):
        if x==1: t = 6 if poss[x+1] > i else 7
        elif x==n: t = 8 if poss[x-1] < i else 9
        else: t = getTyyppi(poss[x-1], i, poss[x+1])

        tempSet = IntvalSet()
        
        if t==2:
            tempSet.add(poss[x-1], poss[x+1]-1)
        elif t==5:
            tempSet.add(0, poss[x+1])
            tempSet.add(poss[x-1], n-1)
        elif t==1:
            tempSet.add(poss[x-1]+1, poss[x+1])
        elif t==7:
            tempSet.add(0, poss[x+1])
        elif t==9:
            tempSet.add(poss[x-1], n-1)
        
        toConsider.append(tempSet)


    """
    mat = [[1 if toConsider[i].has(j) else 0 for j in range(n)]
           for i in range(n)]

    print ("To Consider Matrix:")
    #print (mat)
    for row in mat: print row
    """

    tree = FenwickTree(n)

    events = [[] for _ in range(n+1)] #grouped by time
    #print ("Events len is "+str(len(events)))
    for i, c in enumerate(toConsider):
        for (j1, j2) in c.intvals:
            #print ("Pushing events for interval [%d, %d]" %(j1, j2))
            events[j1].append([j1, i, 1])
            events[j2+1].append([j2+1, i, -1])
    
    numPairs = 0
    for k, evs in enumerate(events[:-1]):
        for (j, i, d) in evs:
            tree.incr(i, d)
        for (j1, j2) in toConsider[k].intvals:
            numPairs += tree.getSum(j2) - tree.getSum(j1-1)
        #subtract the adjacent ones
        x = arr[k]
        if (x>1 and toConsider[k].has(poss[x-1])
            and toConsider[poss[x-1]].has(k) ): numPairs -= 1
        if (x<n and toConsider[k].has(poss[x+1])
            and toConsider[poss[x+1]].has(k) ): numPairs -= 1
    
    return numPairs/2


def junarata(arr):
    n = len(arr)
    poss = getPoss(arr)
    origRounds = getNumOfRounds(poss)
    numPairs = countScore2Pairs(arr)
    return "%d %d %d" %(n, origRounds-2, numPairs)



import random

#a = range(1, 10); random.shuffle(a); print(junarata(a));

#Solve putkaposti
for i in range(1, 10):
    print(junarata(getData(i)))



#Test execution time
"""
import time
startT = time.time()
print (junarata(getData(9)))
print ("Time: "+str(time.time()-startT))
"""

Junarata-tehtävä

class IntvalSet:
    """Set where the included natural numbers
        are given as intervals [a, b]"""
    def __init__(self):
        self.intvals = []

    def has(self, x):
        for (a,b) in self.intvals:
            if a<=x and x<=b: return True
        return False

    def add(self, a, b):
        self.intvals.append((a, b))

    def getAll(self, lowerBdd=0):
        """generates all numbers in the set.
            If lowerBdd is given,
            generates only numbers >= lowerBdd"""
        for (a,b) in self.intvals:
            for x in range(max(a, lowerBdd), b+1):
                yield x

    def __str__(self):
        return ", ".join([str(x) for x in self.getAll()])
    


def getData(num):
    ret = []
    with open("junar"+str(num)+".in") as f:
        ret = f.read().split("\n")
    return [int(x) for x in ret[1:] if len(x)>0]


def getPoss(p):
    """Get list of positions of elements in the
        permutation p of [1..n]
        e.g. getPoss([2,5,3,1,4]) = [None, 3, 0, 2, 4, 1]
    """
    ret = [None]*(len(p)+1)
    for i,x in enumerate(p):
        ret[x] = i
    return ret

"""
Let's denote types of orderings of k-1, k and k+1
(or k and k+1; k-1 and k in cases k=1; k=n  like this:

123: 0
132: 1
213: 2
231: 3
312: 4
321: 5
_12: 6
_21: 7
12_: 8
21_: 9
"""
def getTyyppi(a, b, c):
    if a<b:
        if b<c: return 0
        if a<c: return 1
        return 3
    else: #b<a
        if b>c: return 5
        if a<c: return 2
        return 4


def getNumOfRounds(poss):
    """How many rounds does it take to gather the numbers.
        @param poss: index-positions of the permutation
    """
    curr = 1
    prevInd = poss[1]
    r = 1
    for curr in range(1, len(poss)):
        currInd = poss[curr]
        if currInd<prevInd: r += 1
        prevInd = currInd
    return r

"""
As the above algorithm shows, the number of rounds for a permutation is
1 + #{k | poss[k+1]<poss[k]}
A swap can decrease the rounds by at most 2.
This can be seen from the following:
Let's consider for each k in [1..n] the positions that into which it is
moved, decrease the rounds. This can only ever be at most 1.
But a swap moves two, so if the other one also is a decreasing move,
we get 2, except in the case when we swap some k and k+1,
since that only decreases the rounds by 1.

(Look at the types, and consider for each how it decreases rounds.)

Let's hope the answer will always be 'score 2' (decreases rounds by 2) as
the existence of such swap is highly probably when we have a large random
permutation.
"""

def countScore2Pairs(arr):
    n = len(arr)
    poss = getPoss(arr)
    toConsider = []
    t = None
    tempSet = None
    for i, x in enumerate(arr):
        if x==1: t = 6 if poss[x+1] > i else 7
        elif x==n: t = 8 if poss[x-1] < i else 9
        else: t = getTyyppi(poss[x-1], i, poss[x+1])

        tempSet = IntvalSet()
        
        if t==2:
            tempSet.add(poss[x-1], poss[x+1]-1)
        elif t==5:
            tempSet.add(0, poss[x+1])
            tempSet.add(poss[x-1], n-1)
        elif t==1:
            tempSet.add(poss[x-1]+1, poss[x+1])
        elif t==7:
            tempSet.add(0, poss[x+1])
        elif t==9:
            tempSet.add(poss[x-1], n-1)
        
        toConsider.append(tempSet)
    
    numPairs = 0
    for i, x in enumerate(arr):
        #if i%1000==0: print ("i="+str(i))
        for j in toConsider[i].getAll(i+1):
            if toConsider[j].has(i) and abs(x-arr[j])>1:
                numPairs += 1
    return numPairs


def junarata(arr):
    n = len(arr)
    poss = getPoss(arr)
    origRounds = getNumOfRounds(poss)
    numPairs = countScore2Pairs(arr)
    return "%d %d %d" %(n, origRounds-2, numPairs)


##---------------------------------------------------
## Random tests

import random
def randomJuna(n):
    a = range(1, n+1)
    random.shuffle(a)
    return [int(x) for x in junarata(a).split()]

from collections import Counter
def genJunaData(n, simuN):
    d = [randomJuna(n) for _ in range(simuN)]
    #print ([x[2] for x in d])
    c1 = Counter([x[1] for x in d])
    c2 = Counter([x[2] for x in d])
    c = c2
    return [ [k, c[k]] for k in sorted(c.keys()) ]
    
##-------------------------------------------------------


#Solve putkaposti
for i in range(1, 10):
    print(junarata(getData(i)))


""" Test execution time
import time
startT = time.time()
print (junarata(getData(9)))
print ("Time: "+str(time.time()-startT))
"""

Kahden vääntyneen kolikon peli

Pelissä on kaksi kolikkoa, kutsutaan niitä ensimmäinen ja toinen kolikko. Todennäköisyydet että näillä tulee kruuna (H) ovat vastaavasti x ja y.

Peli etenee seuraavasti: Pelaaja heittää kolikkoja toistuvasti ja pistää merkille minkä tuloksen hän sai (HT, TH, HH tai TT). Jos hän jossain vaiheessa saa kaksi kertaa putkeen tuloksen HH tai kaksi kertaa putkeen tuloksen TT, niin peli päättyy. Jos taas hänen viime tuloksensa on HH ja hän saa heti tämän jälkeen tuloksen TT, niin hän pääsee bonuspeliin. Bonuspelissä hän heittää molempia kolikkoja (erikseen) kaksi kertaa. Jos ensimmäisellä kolikolla tulee ensimmäisellä heitolla H, pelaaja voittaa 5 pistettä, mutta jos tämän jälkeen toisella heitolla tulee H, pelaaja menettää 3 pistettä (eli voittaa vain 2 pistettä yhteensä). Mikäli ensimmäisellä heitolla tulee T, pelaaja ei voita eikä häviä mitään. Vastaavat luvut toiselle kolikolle ovat 3 ja -4 eli HT antaa 3 pistettä ja HH 3-4 = -1 pistettä. Bonuspelin jälkeen peli jatkuu normaalisti (bonuspelin heittoja ei rekisteröidä, eli pelissä viimeksi on tullut TT).

Miten vääntyneet kolikot peliin kannattaa ottaa (eli kuinka muuttujat x ja y väliltä [0, 1] valita), kun halutaan maksimoida pelaajan voittojen odotusarvo?

Ratkaisu

Lasketaan todennäköisyydet

  • a = \mathbb{P}(HT \wedge TH) = x+y-2xy
  • b = \mathbb{P}(HH) = xy
  • c = \mathbb{P}(TT) = 1-x-y+xy

Muodostetaan Markovin ketju tiloina ”HT tai TH”, ”HH”, ”TT” ja ”End”.

Pelin Markovin ketju, sen siirtymätodennäköisyysmatriisi ja fundamentaalimatriisi. Olemme kiinnostuneet fundamentaalimatriisin ensimmäisen rivin toisesta alkiosta, joka kertoo odotusarvon kuinka monta kertaa ketju käy tilassa HH ennen absorboitumista, kun se alkaa tilasta ”HT tai TH”. (Aloitetaan tästä tilasta, koska silloin ei ole väliä mitä on tullut viimeksi.)

Siirtymätodennäköisyydet ovat muuten tilasta riippumattomat, paitsi tiloista ”HH” ja ”TT” ei mennäkkään vastaavilla heitoilla niihin itseihinsä vaan tilaan ”End”. Tila ”End” on absorboiva.

Pelin voiton odotusarvo saadaan (iteroidun odotusarvon kaavan nojalla) kertomalla yhden bonuspelin voiton odotusarvo (5x -3x^2 + 3y - 4y^2) sillä kuinka monta kertaa siihen odotusarvoisesti päästään. Tämä odotusarvo taas saadaan kertomalla tilassa ”HH” vierailujen odotusarvo siirtymätodennäköisyydellä ”HH” -> ”TT”, joka on c. Näin ollen pelin voiton odotusarvo on

\begin{aligned} \frac{( 5x -3x^2 + 3y - 4y^2 )cb(c+1)}{1-a(b+1)(c+1)-bc} = \frac{ (5x -3x^2 + 3y - 4y^2)xy(1-x-y+xy)(2-x-y+xy) }{1-(x+y-2xy)(1+xy)(2-x-y+xy)-xy(1-x-y+xy)} \end{aligned}

Kuvasta nähdään, että tämän maksimikohta joukossa [0,1]^2 on jossain (0.7, 0.3):n kohdilla:

Odotusarvo riippuu muuttujista x ja y. Kuvaaja ja tasa-arvo käyrä-kuvaaja.

Tässä javascript-funktiot odotusarvon ja sen gradientin laskemiseksi:

/*
(x(5x - 3x^2 + 3y-4y^2)y(1 - x - y + xy)(2 - x - y + xy))
/
(1 - xy(1 - x - y + xy)
- (x + y - 2xy) (1 + xy) (2 - x - y + xy))
, x\in[0,1], y\in [0,1]
*/

const f = (x, y) => {
    let g = 5*x - 3*x**2 + 3*y-4*y**2;
    let num = g*x*y*(1 - x - y + x*y)*(2 - x - y + x*y);
    let denom = 1 - x*y*(1 - x - y + x*y)
        - (x + y - 2*x*y)*(1 + x*y)*(2 - x - y + x*y);
    
    return num/denom;
};

const gradF = (x, y) => {
    let p1 = 5*x - 3*x**2 + 3*y-4*y**2;
    let p2 = x*y*(1-x-y+x*y); //= xy - x^2y - xy^2 + x^2y^2
    let p3 = 2 - x - y + x*y;
    let num = p1*p2*p3;
    let denom = 1 - x*y*(1 - x - y + x*y)
            - (x + y - 2*x*y)*(1 + x*y)*(2 - x - y + x*y);
    
    let dxp1 = 5 - 6*x;
    let dxp2 = y - 2*x*y - y**2 + 2*x*y**2;
    let dxp3 = -1+y;
    let dxnum = dxp1*p2*p3 + p1*dxp2*p3 + p1*p2*dxp3;
    let dxdenom = (y-2)*(y-1)**2 + 3*x**2*y*(1-3*y+2*y**2)
                + x*(2-8*y+14*y**2-6*y**3);
    
    let dyp1 = 3-8*y;
    let dyp2 = x - x**2 - 2*x*y + 2*x**2*y;
    let dyp3 = -1 + x;
    let dynum = dyp1*p2*p3 + p1*dyp2*p3 + p1*p2*dyp3;
    let dydenom = 2*(y-1) + x**2*(-4+14*y-9*y**2)
                + x*(5-8*y+3*y**2) + x**3*(1-6*y+6*y**2)
    
    return [
        (dxnum*denom-num*dxdenom)/denom**2,
        (dynum*denom-num*dydenom)/denom**2
    ];
    
};

Jyrkimmän nousun menetelmällä (gradient descent) tai muuten löydetään maksimi

(x,y) =  (0.71046, 0.28314)

ja odotusarvo

1.2863.

Eräs noppapeli

Pelaaja saa aluksi yhden nopan. Joka kierros pelaaja heittää kaikki noppansa ja menettää ne joiden osoittamaa arvoa tuli useampi kuin yksi. Kierroksen päätteeksi pelaaja saa aina yhden uuden nopan. Tavoite on saada seitsemän noppaa (tai yleisemmin d+1 noppaa, kun nopassa on d tahkoa). Toisin sanoen peli päättyy, kun pelaajalla on kuusi noppaa ja niillä tulee jokaisella eri arvo. Mikä on yhden pelin kierrosten määrän (N) odotusarvo?

Ratkaisu

Peli on Markovin ketju, jossa tilat ovat noppien määrät 1, 2, \dots, d+1. Meidän täytyy selvittää siirtymätodennäköisyydet näiden välillä. Pienille arvoille, kuten d = 6, tämä voidaan tehdä tutkimalla kaikki mahdolliset heitot ja kirjaamalla kuinka moni niistä säilyttää kuinkakin monta noppaa. Seuraava Sage-koodi tekee tämän:

from itertools import product
from collections import Counter

def getTransitionsProbsBrute1(diceAmounts, dieFaces):
    n = dieFaces+1
    ret = [0]*n
    for dice in product(range(dieFaces), repeat=diceAmounts):
        keep = sum(x for x in Counter(dice).values() if x==1)
        #we're adding 1 die after round, but 1 is the
        #first state, so keep is the correct index
        ret[keep] += 1
    tot = dieFaces**diceAmounts
    return [x/tot for x in ret]

Laskentaa voidaan hieman nopeuttaa käyttämällä kombinaatioita (takaisin panolla):

from itertools import combinations_with_replacement as cmb_rep
from collections import Counter

def getType(dice):
    return list(Counter(dice).values())

def getTransitionsProbsBrute2(diceAmount, dieFaces):
    n = dieFaces+1
    ret = [0]*n
    for dice in cmb_rep(range(dieFaces), r=diceAmount):
        keep = sum(x for x in Counter(dice).values() if x==1)
        #multinomial of the type gives the amount
        #of orderings this corresponds to when account for all
        weight = multinomial(getType(dice))
        ret[keep] += weight
    tot = dieFaces**diceAmount
    return [x/tot for x in ret]

Tässä tarvittiin tietää kombinaatiosta sen tyyppi eli siinä esiintyvien lukujen esiintymien lukumäärät (koottuna vektoriksi).

Nyt voidaan muodostaa siirtymämatriisi ja käyttää tavanomaisia keinoja \mathbb{E}[N]:n selvittämiseen:

def getTransitionMatrixBrute(dieFaces):
    arr = [getTransitionsProbsBrute2(i, dieFaces)
           for i in range(1, dieFaces+2)]
    return Matrix(arr)

def getFundamentalMatrix(pMat):
    n = pMat.dimensions()[0]-1
    arr = [[(1 if i==j else 0)-pMat[i][j]
            for j in range(n)]
           for i in range(n)]
    return Matrix(arr)**-1

def getExpectedGameLen(dieFaces):
    P = getTransitionMatrixBrute(dieFaces)
    fundMat = getFundamentalMatrix(P)
    return sum(fundMat[0])


print (getExpectedGameLen(6))

Ja näin saadaan tulos

\begin{aligned} \mathbb{E}[N] = \frac{133191861}{12500} = 10655.34888 \end{aligned}

Kombinaatioiden (eritoteen takaisinpanon kanssa) määrä kuitenkin kasvaa melko nopeasti, kun d kasvaa. Edellä jo mainittua tyyppi-ideaa voidaan jatkojalostaa: se kuinka monta noppaa heittotuloksesta menetetään riippuu vain tyypistä (kaikki, joiden arvoa yli yksi kappale menetetään). Tyypit ovat noppien lukumäärän ositukset (paitsi viimeisen tilan ositus [1,1,\dots,1] ei kelpaa, sillä se tarkoittaa d+1 eri arvoa, mutta nopassa on vain d eri arvoa. (Tästä vielä huomio-osioissa.) Nyt täytyy enää laskea kuinka moni noppien heitto vastaa kutakin ositusta. Tätä on selvitetty koodin kommenteissa:

def getFactorOfType(t, dieFaces):
    #first choose a value for each number in the type
    ret = factorial(dieFaces)/factorial(dieFaces-len(t))
    #divide away the orderings in the type
    #(we need the type's type)
    ret /= factorialProd(getType(t))
    #now we can multiply by the total number of orderings
    ret *= multinomial(t)
    return ret

def getTransitionMatrix(dieFaces):
    n = dieFaces+1
    arr = [[0]*n for i in range(n)]
    for i in range(1, n+1):
        for t in Partitions(i):
            if (len(t)>dieFaces): continue
            toState = i-sum(x for x in t if x>1)
            weight = getFactorOfType(t, dieFaces)
            arr[i-1][toState] += weight
    for i in range(n):
        tot = dieFaces**(i+1)
        for j in range(n):
            arr[i][j] /= tot
    return Matrix(arr)

Tässä vielä lopullinen tyypeillä laskeva koodi kokonaisuudessaan ja taulukko tuloksesta muutamilla d:n arvoilla:

from collections import Counter

def getType(dice):
    return list(Counter(dice).values())

def factorialProd(arr):
    return reduce(lambda acc,x: acc*factorial(x), arr, 1)

def getFactorOfType(t, dieFaces):
    ret = factorial(dieFaces)/factorial(dieFaces-len(t))
    ret /= factorialProd(getType(t))
    ret *= multinomial(t)
    return ret

def getTransitionMatrix(dieFaces):
    n = dieFaces+1
    arr = [[0]*n for i in range(n)]
    for i in range(1, n+1):
        for t in Partitions(i):
            if (len(t)>dieFaces): continue
            toState = i-sum(x for x in t if x>1)
            weight = getFactorOfType(t, dieFaces)
            arr[i-1][toState] += weight
    for i in range(n):
        tot = dieFaces**(i+1)
        for j in range(n):
            arr[i][j] /= tot
    return Matrix(arr)

def getFundamentalMatrix(pMat):
    n = pMat.dimensions()[0]-1
    arr = [[(1 if i==j else 0)-pMat[i][j]
            for j in range(n)]
           for i in range(n)]
    return Matrix(arr)**-1

def getExpectedGameLen(dieFaces):
    P = getTransitionMatrix(dieFaces)
    fundMat = getFundamentalMatrix(P)
    return sum(fundMat[0])


res = [[d, float(getExpectedGameLen(Integer(d)))]
      for d in range(1, 21)]
print (res)
dE[N]
11
24
315
480.96296
5718.61364
610655.34888
7261903.04770
810636087.48690
9713183947.85983
107.89467e+10
111.44259e+13
124.35102e+15
132.16594e+18
141.77942e+21
152.41246e+24
165.39710e+27
171.99232e+31
181.21348e+35
191.21944e+39
202.02174e+43

Huomioita

  • Viimeisen tilan siirtymiä ei tarvitsisi laskea (siinä mielessä, että peli jatkuisi) vaan se asetetaan absorboivaksi eikä sitä tarvita odotusarvon laskussa ollenkaan. Itse asiassa sen laskeminen on juuri kaikista matriisin riveistä vaativin, sillä tarvittavia objekteja (permutaatioita takaisinpanolla, kombinaatioita takaisin panolla tai osituksia) on eniten.
  • Myös ensimmäinen tila on käytännössä turha: yhtä noppaa on turha edes heittää, koska sitä ei voi menettää; ottaa vaan suoraan toisen nopan. Mutta pidetään se mukana yksinkertaisuuden vuoksi ja lasketaan pelin pituuteen mukaan myös nämä deterministiset askeleet.

Tehtäviä

  • Muokkaa koodia, niin että turhia viimeisen tilan siirtymiä ei lasketa. Aseta viimeiseen tilaan itsesiirtymä todennäköisyydellä 1 tai unohda koko viimeinen tila, sillä sitä ei laskuissa tarvita.
  • Laske N:n varianssi. Vinkki: katso Wikipedian artikkeli.
  • Entä jos sovitaan, että pelaaja häviää pelin, mikäli hän jollain kierroksella menettää kaikki noppansa. Mikä on häviämisen todennäköisyys?
  • Voidaan sopia mitä erilaisimpia sääntöjä siitä miten pelaaja saa lisänopan /-noppia tai menettää niitä riippuen hänen heitostaan. Tutki jotain tällaista muokattua peliä.