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ä.

Kaksi karsinaa

Kaikkihan tietävät, että ympyrä on paras muoto kun halutaan aidata tietyn suuruinen ala käyttäen mahdollisimman vähän aitaa (ks. isoperimetrinen epäyhtälö). Mutta entäpä, jos onkin tehtävä kaksi karsinaa, jotka voivat käyttää yhteistä aitaa? Entä yleisemmin n karsinaa?

Oletetaan että yhden karsinan ala on A = 1 ja yritetään löytää mahdollisimman pieni piirinen aitaus.

Ensimmäinen vaihtoehto, joka tulee mieleen on kaiketi ”tavallinen” suorakaiteen mallinen, jossa on yksi väliaita:

Optimoidaan piiri ja saadaan x = \frac{\sqrt{3}}{2} = 0.866 ja p=2\sqrt{3} + \frac{6}{\sqrt 3} = 6.9282.

Mutta onko tämä paras mahdollinen? Koska ympyrä on paras ratkaisu yhdelle karsinalle, niin kokeillaanpa jakaa ympyrä kahtia:

Saadaan r = \sqrt{\frac{2}{\pi}} ja p = 2\pi r + 2r = 2(\pi+1) \sqrt{\frac{2}{\pi}}  = 6.609. Siis parempi kuin suorakaide-ratkaisu.

Mutta entäpä jos yhdistettäisiin hieman kumpaakin: laitetaan keskelle suorakaide ja päihin puoliympyrät (ja tietenkin keskelle jakoaita) ja optimoidaan tämä muoto:

Optimiksi saadaan r = \sqrt{\frac{2}{2+\pi}}=0.6237 ja p=2\sqrt{4+2\pi}=6.413. Löydettiin siis halkaistua ympyrääkin parempi ratkaisu.

Toinen tapa yrittää käyttää ympyrää olisi ollut tehdä kaksi täysin erillistä ympyräkarsinaa. Tällöin yhdellä on pienempi säde (r=\sqrt{\frac{1}{\pi}}) mutta piiri on p = 2\times 2\pi r = 4\sqrt \pi = 7.089 > 7 ja p=7 saadaan jo siitä (optimoimattomasta) ratkaisusta, että laitetaan kaksi neliötä vierekkäin.

Olemme tähän asti ajatelleet, että aitauksen kuuluisi olla konveksi (tietenkin kummankin karsinan pitää olla konveksi, mutta että koko rakennelman?!?!?). Eihän sen täydykään: otetaan vain kaksi sopivan kokoista ympyrää ja läimäytetään ne yhteen. Ei keskikohdasta, vaan optimoidaan tämä kohta.

Tämä onkin niinsanottu tuplakupla ratkaisu ja myös todistetusti (tai itseasiassa 2d-versio on tämä) ongelman optimi ratkaisu, joten optimi on p = 6.3591.

Huomioita

  • Jos karsinat eivät olisikaan samankokoiset, niin tuplakupla ratkaisussa jakoaita ei olekaan jana, vaan ympyrän kaari (ks. tuplakuplan määritelmä), joten karsinatkaan eivät välttämättä ole konvekseja. Jakokaari määräytyy siitä, että kaarien täytyy kohdata toisensa 120 asteen kulmassa eli niiden tangentit jakavat täysikulman kolmeen yhtäsuureen osaan.
  • Tarkka ratkaisu, joka tästä 120 asteen vaatimuksesta nähdään on \alpha = \frac{\pi}{3}, r = \frac{1}{\sqrt{\frac{2\pi}{3} + \frac{\sqrt 3}{4}}} ja p=\frac{\frac{8\pi}{3} +\sqrt 3}{ \sqrt { \frac{2\pi}{3} + \frac{\sqrt 3}{4} }}

Suorakaiteen neliöiminen

Suorakaide halutaan täyttää päälekkäinmenemättömillä neliöillä. Oletetaan, että suorakaiteen sivujen pituudet ja myös neliöiden sivut ovat positiivisia kokonaislukuja. (Jos ne ovat rationaalilukuja, niin voimme skaalata kaikkien nimittäjien pienimmällä yhteisellä jaettavalla, jolloin päästään kokonaislukutilanteeseen.)

Tämä ”neliöinti” on tietenkin helppoa, kun lähdetään vain erottamaan (esim. mahdollisimman suuria) neliöitä (lopulta voidaan aina jakaa 1×1-neliöihin). Neliöinnissä saattaa tapahtua niin, että se koostuu kahdesta osasta eli suorakaide voidaan jakaa joko pysty- tai vaakaviivalla kahteen suorakaiteeseen siten, että ei kuljeta minkään neliön läpi. Toisin sanottuna neliöiden sivut ovat jossain kohtaa kaikki ”linjassa” koko suorakaiteen matkalta. Kutsutaan tällaisia neliöintejä jaettavissa oleviksi. Jos neliöinti ei ole jaettavissa oleva, on se jaoton. Jokainen neliöinti voidaan muodostaa jaottomista neliöinneistä niitä yhdistelemällä.

Jaettavissa oleva neliöinti (vasemmalla) voidaan jakaa keltaisen viivan kohdalta. Jaottomassa neliöinnistä (oikealla) tällaista jakokohtaa ei löydy.

Toisaalta erityisen kiinnostavia olisi myös sellaiset neliöinnit, joissa jokainen käytetty neliö on erikokoinen. Näitä kutsutaan täydellisiksi. Myös täydellinen neliöinti voi olla jaettavissa oleva, mutta tällöin sen osissa käytettyjen neliöiden täytyy tietenkin olla myös erikokoisia, joten mitä tahansa kahta täydellistä neliöintiä ei voi yhdistää täydelliseksi neliöinniksi.

Entistä rajaavampi oletus neliöinnille on ns. ”yksinkertaisuus”: neliöinti ei pidä sisällään mitään aidosti pienempää neliöintiä eli toisin sanottuna neliöt eivät muodosta sisälle pienempää vähintään kaksi neliötä sisältävää suorakaidetta. Jos neliöinti ei ole yksinkertainen, sen sanotaan olevan yhdistetty. Kuvan 1 molemmat esimerkit ovat yhdistettyjä (oikean puoleisesta löytyy monta sisä-suorakaidetta vaikkei koko suorakaiteen halkaisevaa viivaa löydykään). Sen sijaan kuvan 2 (ks. alla) neliöinti on yksinkertainen.

Yhteys tasoverkkoihin

Olkoon S neliöinti. Merkitään siinä olevien vaakasuorien janojen joukkoa V. Verkon kaaret vastaavat neliöinnin neliöitä (neliö yhdistää janat, joille sen vastakkaiset sivut kuuluvat). Suorakaiteen ylä- ja ala-sivuja vastaavia solmuja kutsutaan verkon pohjois- ja etelä-navoiksi. Esimerkiksi:

Neliön väri vastaa kaaren väriä.

Vastaavasti tämä konstruktio voitaisiin tehdä käyttämällä neliöinnin pystysuoria janoja (tällöin vasen ja oikea sivu ovat napoja). Tämä vastaa duaaliverkkoa. Tason duaaliverkko muodostetaan siten, että solmuiksi otetaan edellisen tahkot (sivujen rajoittamat tason alueet). Nyt tehdään kuitenkin vielä yksi lisäys: edellisen navat yhdistetään kaarella, joka jakaa ulkopuolen kahteen alueeseen ja nämä ovat duaaliverkon navat.

Huomataan, että neliöiden sivujen suuruudet liittyvät verkon virtauksiin. Asetetaan verkon sivulle sitä vastaavan neliön suuruinen virtaus. Virtauksen suunta on ilmeinen, kun katsotaan kuinka sivu suhtautuu suorakaiteeseen (jossa virtaus on alaspäin). Huomataan, että tällöin nämä virtaukset totetuttavat Kirchhoffin lait: jokaiseen solmuun saapuva virtaus on yhtä suuri kuin siitä lähtevä virtaus. Perustelu: solmu vastaa neliöinnin vaakaviivaa. Viivan yläpuolella on saapuvat neliöt (eli kaaret) ja alapuolella lähtevät. Kummatkin ovat yhteensä viivan pituus. Tämä toteutuu myös pohjois- ja etelä-navoilla, kun lisätään verkkoon vielä kaari etelä-navalta pohjois-navalle. Myös jokaisen tahkon ympäri kulkevan virtauksen summa on nolla. Tämä onkin edellistä vastaava vaatimus duaaliverkolle.

Nyt voidaankin miettiä toimisiko temppu myös toiseen suuntaan. Tehdään tasoverkolle seuraavaa. Valitaan kaksi solmua navoiksi, laitetaan pohjois-navalle saapuvaksi virta ja ratkaistaan virtaukset, siten että Kirchhoffin säännöt toteutuvat. Etelä-navalta lähtee saman suuruinen virta kuin mikä saapui pohjois-navalle ja sen arvolla voidaan ratkaisu skaalata kokonaislukuvektoriksi. Tästä ratkaisusta voidaan koota suorakaiteen neliöinti.

Katsotaan pientä esimerkkiä:

Tasoverkko, sen Kirchhoffin säännöt koodaava matriisi K, sen ydinvektori (eli yhtälön Kx = 0 ratkaisu, huom. voidaan skaalata kokonaislukuvektoriksi, sillä ratkaisu tuottaa aina rationaalilukuja) ja siitä koottu neliöinti. Huom. napoja ja ulkotahkoja ei ole otettu mukaan.

Tässä esimerkissä huomataan eräs erikoisuus: neliöintiin tuli vierekkäin kaksi saman kokoista (4-sivuista) neliötä. Ne tulivat kaarista 0-1 ja 1-2. Mikä näissä kaarissa on vikana? Jos muodostamme duaaliverkon, huomamme että se onkin multiverkko! Vihreän kolmion 012 ja violetin nelikulmion 0123 välillä on kaksi eri sivua. Tämä verkko ei ole 3-kaari-yhtenäinen, vaan solmu 1 voidaan irroittaa verkosta leikkaamalla kaksi kaarta. Kolme-kaari-yhtenäisyys onkin hyvä vaatimus verkolle, jotta saadaan ”hyviä” neliöintejä.

Jos palataan alussa mainittuun neliöinnin jaottomuuden käsitteeseen, niin nyt huomataan, että se vastaa verkon (tai sen duaalin) 1-solmu-epä-yhtenäisyyttä: jos koko suorakaiteen läpi menevää viivaa vastaava solmu poistetaan, tulee verkosta epäyhtenäinen.

Kysymyksiä

  • Jos tehtävää yleistetään kolmiulotteiseksi eli täytetään laatikko [0, a]\times[0, b]\times[0,c] kuutioilla, kuinka muuttuu vastaavasti muodostettu verkko, kun vaihdetaan akselia, jonka suhteen se muodostetaan? Verkon solmut ovat tasojen osia, joille kuutoiden tahkot asettuvat.

Linkkejä

Binääriblokit

Väritetään n:n peräkkäisen ruudun rivistä osa mustiksi ja lasketaan millaiset putket väritettyjä tuli. Esim:

n=10, värityksestä muodostuu putket [2, 1, 3]

Yhteensä erilaisia värityksiä on 2^n, sillä jokainen ruutu voidaan joko värittää tai jättää tyhjäksi (vrt. binääriluvut < 2^n), mutta kysymykset koskevat värityksestä muodostuvia putkia. Sanotaan putkien pituuksista koottua vektoria putkityypiksi. Yo. esimerkissä putkityyppi on [2, 1, 3].

  • a) Kuinka monta erilaista putkityyppiä voidaan n:n ruudun rivistä muodostaa?
  • b) Olkoon rivin pituus n. Kuinka moni eri väritys tuottaa putkityypin [k_1, k_2, \dots k_m]?

Mahdolliset väritykset ja putkityypit, n = 1, 2, 3, 4

Ratkaisu

a)

Putkityyppi [k_1, k_2, \dots, k_m] on mahdollinen muodostaa, mikäli se mahtuu n-riviin. Jokaisen putken välissä on oltava vähintään 1 tyhjä ruutu ja välejä on m-1, joten se mahtuu tasan silloin kun

\begin{aligned} m - 1 + \sum_{j=1}^m k_j \leq n \end{aligned}

Määritellään lukujonot a_n ja b_n seuraavasti

\begin{aligned} a_n = \#\{(k_1, k_2, \dots, k_m)   \hspace{.3cm}  |  \hspace{.3cm}   m - 1 + \sum_{j=1}^m k_j \leq n \}   \\ b_n =   \#\{(k_1, k_2, \dots, k_m)  \hspace{.3cm} |  \hspace{.3cm}   m - 1 + \sum_{j=1}^m k_j = n \}  \end{aligned}

eli b_n laskee sellaiset putkityypit, joita ”ei enää voi pidentää” ja a_n on kysytty jono.

Huomataan, että b_n itse asiassa laskee n-pituisia binäärijonoja, joissa ei ole yli yhden pituisia nollaputkia ja jotka alkavat ja päättyvät ykköseen. Näiden lukumäärä voidaan laskea seuraavan graafin n-1 -kävelyjen solmusta A solmuun A määränä.

Solmuja voidaan ajatella tiloina seuraavasti: A=”viimeksi on tullut ykkönen”, B=”viimeksi on tullut nolla”, C=”on tullut kahden nollan putki” . Tila C on hylkäävä, siitä ei pääse enää mihinkään. (Voi myös ajatella, että C:stä on siirtymä vain itseensä; todennäköisyydellä 1, jos graafi mielletään Markovin ketjuksi.) Aloitetaan tilasta A, sillä jonon täytyy alkaa ykkösellä, joten voidaan olettaa, että yksi ykkönen on jo tullut. Kävelyn täytyy myös päättyä tilaan A, sillä viimeinen numero halutaan olevan ykkönen.

Se, että näin todella saadaan haluttu lukumäärä, nähdään esimerkiksi ajattelemalla graafia Markovin ketjuna, jossa tilojen A ja B siirtymät tapahtuvat molemmat todennäköisyydellä \frac{1}{2} ja tilasta C on itsesiirtymä todennäköisyydellä 1. Näin saadaan todennäköisyys, että binäärijono on kysytynlainen (ks. tilojen selitykset kuvatekstissä). Tämä muutetaan lukumääräksi kertomalla kaikkien n-1-pituisten binäärijonojen lukumäärällä eli 2^{n-1}:llä.

Lasketaan siis vierekkäisyysmatriisin n-1:s potenssi diagonalisoimalla

\begin{aligned}     \left[ \begin{array}{ccc} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 0 & 0 \end{array} \right] ^{n-1}  =      \left[ \begin{array}{ccc} -1 & \frac{1-\sqrt 5}{2} &  \frac{1+\sqrt 5}{2}  \\ * & * & * \\ * & * & * \end{array} \right]   \left[ \begin{array}{ccc} 0 &  &    \\  &  \left( \frac{1-\sqrt 5}{2} \right)^{n-1}  &  \\  &  &  \left( \frac{1+\sqrt 5}{2} \right)^{n-1}  \end{array} \right]    \left[ \begin{array}{ccc} 0 & * & * \\ -\frac{1}{\sqrt 5} & * & * \\   \frac{1}{\sqrt 5} & * & * \end{array} \right]           \end{aligned}

Tämän alkio paikassa (0, 0) antaa luvun b_n. Saadaan

\begin{aligned}  b_n = \frac{1}{\sqrt 5} \left (  \frac{1+\sqrt 5}{2} \right)^n -     \frac{1}{\sqrt 5} \left (  \frac{1-\sqrt 5}{2} \right)^n  \end{aligned}

joten b_n on kuuluisa Fibonaccin lukujono f_n =f_{n-1} + f_{n-2}, f_1=1, f_2 = 1.

Ratkaistaan sitten a_n. Koska n+1-pituiseen riviin mahtuvalle putkityypille (k_1, k_2, \dots, k_m) pätee joko m - 1 + \sum_{j=1}^m k_j \leq n tai m - 1 + \sum_{j=1}^m k_j  = n+1 , niin nähdään rekursio

\begin{aligned}  a_{n+1} = a_n + b_{n+1}  \end{aligned}.

Osoitetaan nyt induktiolla, että a_n = f_{n+2}. Perustapaukset a_0 = 1 ja a_1 = 2 nähdään helposti luetteloimalla. Induktioaskel saadaan suoraan edellä olleesta rekursiosta, sillä, jos oletetaan a_n = f_{n+2}, niin

\begin{aligned}  a_{n+1} = a_n + b_{n+1}  = f_{n+2} + f_{n+1} = f_{n+3}. \end{aligned}

b)

Olkoon putkityyppi [k_1, k_2, \dots k_m]. Merkitään väritettyjen ruutujen määrää S = \sum_{j=1}^m k_j. Kysymys voidaan muotoilla toisin: Asetetaan n-S tyhjää ruutua jonoon ja kysytään kuinka monella tavalla näiden väliköistä (ja päädyistä) voidaan valita m kappaletta, joihin putkityypin putket (järjestyksessä) asetellaan. Välikköjä ja päätyjä on n-S-1 + 2 = n-S+1 kappaletta, joten vastaus on

\begin{aligned}   {n-S+1} \choose {m} \end{aligned}.

Huomioita

  • Jonon b_n voi todeta Fibonacci jonoksi myös huomaamalla rekursioyhtälön b_n = b_{n-1} + b_{n-2} (joko suoraan graafista tai matriisien kertolaskusta).