Tutkitaan :n pituisia jonoja symboleista ja . 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 lasketaan vaihteluväli erotuksena , missä on kävelyn korkeus askeleen jälkeen.

Kysymys: kuinka monta kävelyistä on vaihteluväliltään korkeintaan , missä 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 kaava

brw(2, n) =

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

Koska emme kuitenkaan keksi tästä yleistä kaikille toimivaa kaavaa, ryhdytään kehittämään tehokkaampaa algoritmia. Kun jonoa luetaan alusta lähtien ja tarkastellaan korkeutta, jolla ollaan askeleen jälkeen, niin mieleen tulee tietenkin Markovin ketju suunnattu verkko, jossa tilana solmuina on . 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 . 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 (), 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)

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 -putkeen, johon päädyttiin. Huom. ”putki” on itse asiassa verkkona (suuntaamaton) polku.

Keskitytään aluksi yhteen -putkeen ja sen kävelyihin eli yllä mainittuun ”seilaamiseen”. Merkitään :lla ensimmäisestä päästä () lähtevien kävelyjen määrää pituisessa putkessa, jotka päätyvät putken solmuun . 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 -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ö

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 (). Koitetaan saada molempien maailmojen paremmat puolet ja käytetään keskikokoisille :n ja :n arvoille rekursiivista metodia ja pienille :n arvoille matriisimenetelmää. Arvolle 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 . 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 (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 kaava .
  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 . Onko näillä joku yhteys (onhan verkon osat polkuja)?
• Keksi kaava yleiselle .
• Tutki vaihteluvälin jakaumaa :n pituisille -kävelyille. Kuinka sen keskiarvo käyttäytyy?

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 osajonojen lukumäärä ? Oletetaan että :ssä esiintyvien merkkien lukumäärä on jokin ennalta kiinnitetty pienehkö luku. Halutaan löytää lineaariaikainen algoritmi.

Ratkaisu

Olkoon sana, jonka osajonot halutaan laskea. Olkoon :n pituus ja siinä esiintyvät merkit . Kun samaistetaan jokainen :n osajono siihen indeksien osajonoon, josta löytyy kaikkein aikaisimmin, vastaa jokainen sanan osajono solmusta lähtevää kävelyä suunnatulla graafilla , jonka solmut ovat ja solmusta lähtee kaaret missä on indeksi, jossa merkki esiintyy :ssä seuraavan kerran lukien indeksistä lähtien (ja jos ei enää esiinny, niin kaarta ei ole) ja tämä vastaavuus on bijektiivinen.

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

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

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 pituinen ja , joten summan suppenemisen eikä :n kääntyvyyden kanssa tule ongelmia.

Merkitään . Mehän tarvitsemme vain sen ensimmäisen rivin summan, joten koko matriisia ei tarvitse muodostaa, eikä summaa laskea, vaan riittää ratkaista

,

missä 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ä on alakolmiomatriisi, jossa on ykkösiä diagonaalilla ja jokaisessa sarakkeessa korkeintaan :ssä paikassa (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)