read(`code_finsubgames.txt`):

Help_aperiodic:=proc()
    print(` SearchForAperiodic(S, depth, stopfound:=false, prd:=false) `):
    print(` Wordize(period) , GeneralWordize(str, blocks) , Seq1(n) `):
    print(` Seq2(n) , Sn(n) , SnIndex(n, i) , XSn(n, valsonly:=false) `):
    print(` SGvalsOld(n, offset:=0) , SGvals(n, offset:=0) `):
    print(` SGvalsDecOld(n, offset:=0) , SGvalsSplit(n, offset:=0) `):
    print(` SGvalsSplitNested(n, offset:=0) , SGvalsSplitCS(n, offset:=0) `):
    print(` G(n, i, trace:=false) , GDiff(n, i, c:=0) `):
    print(` LCommonSuffix(S1, S2) , LCommonSuffixFormula(i) `):
    print(` SearchForExpansions(n, cap) , BasicGraph(n:=1) `):
    print(` BasicMexGraph(n:=1) , MexFixGraph(G, Gmex) , Mex1Seq(n) `):
    print(` SeqStart(n, i, len) , SeqComp(seq1, seq2) , AdjSeqComp(n, i) `):
    print(` Sn2(n) , Akn(k, n) , Fibk(k, n) , Fibk2(k, n) , Beatty(alpha, n) `):
    print(` SortDec(S) , SeqLE(f, n) , DigitExpansion(S, n) `):
    print(` kZeckendorf(k, n) , kZeckendorf2(k, n) , Zeckendorf(n) `):
    print(` Aknp1DE(k, n) , Seq1p1DE(n) , Wn(k, n) , ValSeq(S, val) `):
    print(` WValSeq(k, val, cap) , WordComp(seeds, rule) `):
    print(` AnalyzeWordComp(WC, n, suppress:=false) , TMWC(w, i) `):
    print(` Traceables(n) , TraceablesList(n) , TraceablesEL(n) `):
    print(` TraceablesELList(n) , TraceablesUL(n) `):
    print(` RandTraceable(n, k:=2, p:=0.5) `):
    print(` AnalyzeTraceable(trseq, n, suppress:=false) `):
    print(` BuildPWord(P, k:=2) , AnalyzePWord(P, suppress:=false) `):
    print(` BuildSet(n, k:=2, p:=0.5) , SGvalsGeneral(S, offset:=0) `):
    print(` BuildTraceable(n, k:=2, p:=0.5) `):
    print(` SGvalsGeneralizeBT(n, k:=2, p:=0.5) , Greedy2(n, k:=2) `):
    print(` Changek(w, k, S:=false) , Traceable(trseq, n:=false) `):
    print(` Promotion(sq, j:=1) , PromotionFn(n, sq, j:=1) , Increasek(w) `):
    print(` SGvalsGeneralizeBT2(n, k:=2, p:=0.5) , RandSeq(n, p:=0.3) `):
    print(` Promotest(n, p:=0.3) , TraceableReps(trseq, n:=false) `):
    print(` SubComp(L, i, suppress:=false) , SetFrom(L) `):
    print(` SearchForPotentialSubCompFail(n, k:=2, test:=0, suppress:=true) `):
    print(` SearchForOne(S) , PromoSubtract(L) , SkipSeq(L, skip:=2) `):
    print(` RepSkip(L, cap, skip:=2) , CompRep(Ls, cap) `):
    print(` ExtendRepNimSeqTest(L, cap, n:=3) `):
end:

#Try to find an aperiodic game with bounded Nim sequence
SearchForAperiodic:=proc(S, depth, stopfound:=false, prd:=false)
local S2, perd, period, prefix, i, cap, found, chkamt, j,
        chkseq, SGsub, tp, nval, ok:
    SGsub:=proc(per, n, st, top, extra) local s, st2:
        st2:={}:
        for s in st do
            if n - s > top then
                st2:=st2 union {extra[n - s - top]}:
            elif n - s >= 0 then
                st2:=st2 union {per[((n-s) mod nops(per)) + 1]}:
            fi:
        od:
        return mex(st2):
    end:
    if depth = 0 then
        printf("Check out %a further\n", S):
        return true:
    fi:
    if prd = false then
        perd:=GetPeriodAndPrefix(S)[1]:
    else
        perd:=prd:
    fi:
    #cap:=max(perd):
    cap:=2:
    S2:=S:
    found:=false:
    chkamt:=5:
    for i from 1 + max(S) to 7 * max(S) do:
        chkseq:=[]:
        tp:=i - 1:
        S2:=S union {i}:
        ok:=true:
        for j from i to i + chkamt - 1 do
            nval:=SGsub(perd, i, S2, tp, chkseq):
            if nval > cap then
                ok:=false:
                break:
            fi:
            chkseq:=[op(chkseq), nval]:
        od:
        if ok then
            #TODO add some quick check
            period, prefix:=GetPeriodAndPrefix(S2):
            if nops(prefix) = 0 and nops(period) > nops(perd)#2 * nops(perd)
                    and max(period) = cap then
                found:=SearchForAperiodic(S2 union {i},
                    depth - 1, stopfound, period) or found:
            fi:
        fi:
        if stopfound and found then
            break:
        fi:
        S2:=S minus {i}:
    od:
    if not found then
        printf("Failed to extend %a\n", S):
    fi:
    return found:
end:

#Convert 0,1 to A and 0,1,2 to B
Wordize:=proc(period) local i, ret:
    ret:="":
    i:=2:
    while i <= nops(period) do
        if period[i] = 0 then
            if period[i-1] = 1 then
                ret:=cat(ret, "A"):
            elif period[i-1] = 2 then
                ret:=cat(ret, "B"):
            else
                ret:=cat(ret, "X"):
            fi:
            i:=i + 1:
        elif period[i] = 2 and period[i-1] = 2 then
            ret:=cat(ret, "X"):
        fi:
        i:=i + 1:
    od:
    if period[-1] = 1 then
        ret:=cat(ret, "A"):
    elif period[-1] = 2 then
        ret:=cat(ret, "B"):
    else
        ret:=cat(ret, "X"):
    end:
    return ret:
end:

GeneralWordize:=proc(str, blocks)
local str2, i, j, current, newcurrent, oks, block, theblock, ret:
    ret:="":
    current:="":
    oks:=blocks:
    str2:=cat(str, "-"):
    for i from 1 to length(str2) do
        newcurrent:=cat(current, str2[i]):
        theblock:=false:
        for block in oks do
            if block[1] = current then
                theblock:=block:
            fi:
            if length(block[1]) < length(newcurrent) or
                    block[1][1..length(newcurrent)] <>
                    newcurrent then
                oks:=oks minus {block}:
            fi:
        od:
        if nops(oks) = 0 then
            if theblock <> false then
                ret:=cat(ret, theblock[2]):
            else
                ret:=cat(ret, "X"):
            fi:
            oks:=blocks:
            current:=str2[i]:
            for block in oks do
                if block[1][1] <> current then
                    oks:=oks minus {block}:
                fi:
            od:
        else
            current:=newcurrent:
        fi:
    od:
    return ret:
end:

#1, 4, 12, 28, 73, 163, 343, 867, 1915, 4011, 8203, 16588?
#2, 4, 4, 17, 17, 17, 181, 181, 181, 181, 181?
#4: AB
#12: AB AB B
#28: ABABB ABABB B
#73: ABABBABABBB ABABBABABBB ABABBB
#163: ABABBABABBBABABBABABBBABABBB ABABBABABBBABABBABABBBABABBB ABABBB
#343: ABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBB ABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBB ABABBB
#867: ABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBBABABBB ABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBBABABBB ABABBABABBBABABBABABBBABABBBABABBABABBBABABBABABBBABABBBABABBBABABBB

#ONE MORE WORDIZE STEP: A<-ABABB, B<-B
#12: A
#28: A A B
#73: AAB AAB AB
#163: AABAABAB AABAABAB AB
#343: AABAABABAABAABABAB AABAABABAABAABABAB AB
#867: AABAABABAABAABABABAABAABABAABAABABABAB AABAABABAABAABABABAABAABABAABAABABABAB AABAABABAABAABABABAB

#ONE MORE WORDIZE STEP: A<-AAB, B<-AB
#28: A
#73: A A B
#163: AAB AAB B
#343: AABAABB AABAABB B
#867: AABAABBAABAABBB AABAABBAABAABBB AABAABBB
#1915: AABAABBAABAABBBAABAABBAABAABBBAABAABBB AABAABBAABAABBBAABAABBAABAABBBAABAABBB AABAABBB

#ONE MORE WORDIZE STEP: A<-AAB, B<-B
#73: A
#163: A A B
#343: AAB AAB B
#867: AABAABB AABAABB AABB
#1915: AABAABBAABAABBAABB AABAABBAABAABBAABB AABB
#4011: AABAABBAABAABBAABBAABAABBAABAABBAABBAABB AABAABBAABAABBAABBAABAABBAABAABBAABBAABB AABB

#ONE MORE WORDIZE STEP: A<-AAB, B<-B
#163: A
#343: A A B
#867: AAB AAB AB
#1915: AABAABAB AABAABAB AB
#4011: AABAABABAABAABABAB AABAABABAABAABABAB AB
#8203: AABAABABAABAABABABAABAABABAABAABABABAB AABAABABAABAABABABAABAABABAABAABABAB ABAB

#Proposed sequence that has aperiodic bounded Nim sequence
#Seq1(i)=4Seq1(i-1)-2Seq1(i-2)+1
Seq1:=proc(n) local i:
option remember:
    #if n=-1 then
    #    return 1:
    #elif n=0 then
    #    return 4:
    #elif n=1 then
    #    return 12:
    #else
    #    return 3*Seq1(n-1) + n + 3 + add(Seq1(i), i=1..n-2):
    #fi:
    if n=-1 then
        return 1:
    elif n=0 then
        return 4:
    elif n=1 then
        return 12:
    elif n=2 then
        return 41:
    else
        return 4*Seq1(n-1) - 2*Seq1(n-2) + 1:
    fi:
end:

#Length of wordize of periods of Seq1
#Equals A161941
#Seq2(i)=4Seq2(i-1)-2Seq2(i-2)
Seq2:=proc(n) local i:
option remember:
    if n=1 then
        return 5:
    else
        return 3*Seq2(n-1) + 1 + add(Seq2(i), i=1..n-2):
    fi:
end:

#Nim sequence period n in our proposed sequence
#Returned as a ternary word
Sn:=proc(n) local i, ret:
option remember:
    if n=-1 then
        return "01":
    elif n=0 then
        return "01012":
    elif n=1 then
        return "0101201012012":
    else
        ret:=cat(Sn(n-1), Sn(n-1), Sn(n-1)):
        for i from n-2 by -1 to 1 do
            ret:=cat(ret, Sn(i)):
        od:
        ret:=cat(ret, "012"):
        return ret:
    fi:
end:

#ith index (starting from 1) in Sn
SnIndex:=proc(n, i) local pos, val, j, sn, k:
    val:=Seq1(n) + 1:
    j:=MapleMod(i, val):
    if n <= 1 then
        sn:=Sn(n):
        return parse(sn[j]):
    elif j = val - 2 then
        return 0:
    elif j = val - 1 then
        return 1:
    elif j = val then
        return 2:
    fi:
    val:=Seq1(n-1) + 1:
    pos:=val * 3:
    if j <= pos then
        return SnIndex(n-1, MapleMod(j, val)):
    fi:
    for k from n-2 by -1 to 1 do
        val:=pos:
        pos:=pos + (Seq1(k) + 1):
        if j <= pos then
            return SnIndex(k, j - val):
        fi:
    od:
    #Should never get here
    return FAIL:
end:

#Sn, but every index where we need to check for the mex is an X
XSn:=proc(n, valsonly:=false) local sn, i, val:
    sn:=Sn(n):
    for i from -1 to n do
        val:=Seq1(i)+1:
        sn:=cat(sn[1..-val-1], "X", sn[-val+1..-1]):
    od:
    return sn:
end:

#Vals you need to take the mex of to find SG(Seq1(n) + offset)
#in Sn
#OLD VERSION, DOES NOT WORK, NATHAN IS IDIOT
SGvalsOld:=proc(n, offset:=0) local i, val, vals:
    vals:=[]:
    for i from -1 to n do
        val:=Seq1(i):
        vals:=[op(vals), SnIndex(n, val + offset)]:
    od:
    return vals:
end:

#Vals you need to take the mex of to find SG(Seq1(n) + offset)
#in Sn
SGvals:=proc(n, offset:=0) local i, val, vals, seq1n:
    vals:=[]:
    seq1n:=Seq1(n):
    for i from -1 to n do
        val:=Seq1(i):
        vals:=[SnIndex(n, seq1n - val + offset + 1), op(vals)]:
    od:
    return vals:
end:

#Vals you need to take the mex of to find SG(Seq1(n) + offset)
#in Sn
#Decorated with DOESwhich type of pointer they are
#OLD THING, MAY NOT WORK, NATHAN IS IDIOT
SGvalsDecOld:=proc(n, offset:=0) local vals, i, val, ind:
    vals:=SGvals(n, offset):
    for i from 1 to nops(vals) do
        ind:=nops(vals)-i+1:
        val:=convert(vals[ind], string):
        if i=1 then
            vals[ind]:=cat(val, "I"):
        elif i=2 then
            vals[ind]:=cat(val, "F"):
        elif i=nops(vals) then
            vals[ind]:=cat(val, "Z"):
        elif offset >= Seq1(i-2) then
            vals[ind]:=cat(val, "E"):
        else
            vals[ind]:=cat(val, "S"):
        fi:
    od:
    return vals:
end:

#Vals you need to take the mex of to find SG(Seq1(n) + offset)
#in Sn
#Split by which pointers are inductive
SGvalsSplit:=proc(n, offset:=0) local i, val, vals, seq1n, bar:
    vals:=[]:
    seq1n:=Seq1(n):
    bar:=false:
    for i from -1 to n do
        val:=Seq1(i):
        if offset - 1 < val and not bar then
            vals:=["|", op(vals)]:
            bar:=true:
        fi:
        vals:=[SnIndex(n, seq1n - val + offset + 1), op(vals)]:
    od:
    return vals:
end:

#Vals you need to take the mex of to find SG(Seq1(n) + offset)
#in Sn
#Split by which pointers are inductive
#Split those recursively
SGvalsSplitNested:=proc(n, offset:=0)
local i, val, vals, seq1n, bar:
    vals:=[]:
    seq1n:=Seq1(n):
    bar:=false:
    for i from n by -1 to -1 do
        val:=Seq1(i):
        if offset > val then
            bar:=true:
            break:
        fi:
        vals:=[op(vals), SnIndex(n, seq1n - val + offset + 1)]:
    od:
    if bar then
        #TODO fix this
        vals:=[op(vals), "|",
            op(SGvalsSplitNested(i, offset - val))]:
    fi:
    return vals:
end:

#Vals you need to take the mex of to find SG(Seq1(n) + offset)
#in Sn
#Split by which pointers are in the common suffix
SGvalsSplitCS:=proc(n, offset:=0) local i, val, vals, seq1n, bar:
    vals:=[]:
    seq1n:=Seq1(n):
    bar:=false:
    for i from -1 to n do
        val:=Seq1(i):
        if offset - 1 < val - LCommonSuffixFormula(n)
                and not bar then
            vals:=["|", op(vals)]:
            bar:=true:
        fi:
        vals:=[SnIndex(n, seq1n - val + offset + 1), op(vals)]:
    od:
    return vals:
end:

#G function (indices in periods)
G:=proc(n, i, trace:=false)
local imod, an, anm1, sval, oldsval, k:
    if trace then
        printf("G(%d, %d)\n", n, i):
    fi:
    if n = -1 then
        return i mod 2:
    elif n = 0 then
        return parse("01012"[(i mod 5) + 1]):
    elif n = 1 then
        return parse("0101201012012"[(i mod 13) + 1]):
    else
        an:=Seq1(n):
        imod:=i mod (an + 1):
        if i >= an - 2 then
            return 2 + i - an:
        fi:
        anm1:=Seq1(n-1):
        sval:=3 * (anm1 + 1):
        if imod < sval then
            return G(n-1, imod mod (anm1 + 1), trace):
        fi:
        for k from 2 to n-1 do
            oldsval:=sval:
            sval:=sval + Seq1(n-k) + 1:
            if imod < sval then
                return G(n-k, imod - oldsval, trace):
            fi:
        od:
    fi:
    #SHOULD NEVER GET HERE!!!
    return FAIL:
end:

#Differences involving G functions
#G(n, Seq1(n-1) - Seq1(n-i) + c)
GDiff:=proc(n, i, c:=0)
    return G(n, Seq1(n-1) - Seq1(n-i) + c, true):
end:

#Length of common suffix
LCommonSuffix:=proc(S1, S2) local i,j,l:
    i:=length(S1):
    j:=length(S2):
    l:=0:
    while i>0 and j>0 and S1[i]=S2[j] do
        l:=l+1:
        i:=i-1:
        j:=j-1:
    od:
    return l:
end:

#Length of common suffix formula
#LCommonSuffix(Sn(i), Sn(i-1))
LCommonSuffixFormula:=proc(i) local j:
    if i=1 then
        return 3:
    elif i=2 then
        return 6:
    else
        #return -2 + Seq1(i) - 3*Seq1(i-1) +
        #    LCommonSuffixFormula(i-1):
        #return 1 + i + LCommonSuffixFormula(i-1) +
        #    add(Seq1(j), j=1..i-2):
        return 4*LCommonSuffixFormula(i-1) -
            2*LCommonSuffixFormula(i-2) + 4:
    fi:
end:

#Search for expansions of Seq1
SearchForExpansions:=proc(n, cap) local S, i, S2, E, perd, ppf:
    S:={seq(Seq1(i), i=-1..n)}:
    perd:=GetPeriodAndPrefix(S)[1]:
    E:={}:
    for i from 2 to cap do
        if not(i in S) then
            S2:=S union {i}:
            #TODO add some quick check
            ppf:=GetPeriodAndPrefix(S2):
            if ppf[1] = perd and nops(ppf[2]) = 0 then
                E:=E union {i}:
            fi:
        fi:
    od:
    return E:
end:

#Graph sending 0 to 1, 1 to 0 or 2, and 2 to 1
BasicGraph:=proc(n:=1) local i:
    if n=1 then
        return Digraph(["[0]","[1]","[2]"],
            {["[0]","[1]"],["[1]","[0]"],["[1]","[2]"],
            ["[2]","[0]"]}):
    else
        return TensorProductD(seq(BasicGraph(1),i=1..n)):
    fi:
end:

#Graph sending 0 to 1, 1 to 0 or 2, 2 to 1 or 3,
#and 3 to 0 or 2
BasicMexGraph:=proc(n:=1) local i:
    if n=1 then
        return Digraph(["[0]","[1]","[2]","[3]"],
            {["[0]","[1]"],["[1]","[0]"],["[1]","[2]"],
            ["[2]","[0]"],["[2]","[3]"],["[3]","[0]"],
            ["[3]","[2]"]}):
    else
        return TensorProductD(seq(BasicGraph(1),i=1..n)):
    fi:
end:

#Remove edges based on mexes, including possibly adding
#additional vertices for repeated mex values in fixing graph
MexFixGraph:=proc(G, Gmex) local V, E, T, j, vtx, e, vtx2,
        vertexmexchecks, vertexfix, edgechecks:
    vertexmexchecks:=proc(vtx) local S, vp, i:
        vp:=vtx[1]:
        S:={}:
        for i from 1 to nops(vp) do
            S:=S union {vp[i][1]}:
        od:
        return evalb(mex(S) = vtx[2][1]):
    end:
    vertexfix:=proc(vtx):
        if type(vtx[1][1], integer) then
            return [[op(vtx[..-2])], vtx[-1]]:
        else
            return vtx:
        fi:
    end:
    edgechecks:=proc(edge)
        return evalb(vertexmexchecks(vertexfix(Parse(edge[1])))
            and vertexmexchecks(vertexfix(Parse(edge[2])))):
    end:
    V:=[]:
    T:=cartprod([[seq(Parse(Vertices(G)[j]),
        j=1..nops(Vertices(G)))], [seq(Parse(Vertices(Gmex)[j]),
        j=1..nops(Vertices(Gmex)))]]):
    while not T[finished] do
        vtx:=T[nextvalue]():
        vtx2:=vertexfix(vtx):
        if vertexmexchecks(vtx2) then
            V:=[op(V), VertexEncode(vtx)]:
        fi:
    od:
    E:={}:
    T:=cartprod([Edges(G), Edges(Gmex)]):
    while not T[finished] do
        e:=T[nextvalue]():
        e:=EdgeMerge(e):
        if edgechecks(e) then
            E:=E union {e}:
        fi:
    od:
    #print(V, E):
    return Graph(V, E):
end:

#Encoded sequence of mex 1 guys
Mex1Seq:=proc(n) local sq, s, ret, i, val:
    sq:=[seq(SGvals(n, i), i=1..Seq1(n))]:
    ret:=[]:
    for s in sq do
        if mex(s) = 1 then
            if 2 in s then
                val:=0:
                for i from 2 to nops(s) - 1 do
                    val:=val * 2:
                    val:=val + s[i] / 2:
                od:
                ret:=[op(ret), val+1]:
            else
                ret:=[op(ret), 0]:
            fi:
        fi:
    od:
    return ret:
end:

#Sequence starting from position i, length len
SeqStart:=proc(n, i, len) local s, sn:
    sn:=Sn(n):
    s:=sn[((i mod Seq1(n)) + 1)..]:
    while length(s) < len do
        s:=cat(s, sn):
    od:
    return s[1..len]:
end:

#Compare sequences to see how different they are
#Retuns a list of block lengths same, different, same, etc
SeqComp:=proc(seq1, seq2) local ret, i, ct, same:
    ct:=0:
    ret:=[]:
    same:=true:
    for i from 1 to min(length(seq1), length(seq2)) do
        if (same and seq1[i] <> seq2[i]) or
                (not same and seq1[i] = seq2[i]) then
            same:=(not same):
            ret:=[op(ret), ct]:
            ct:=0:
        fi:
        ct:=ct + 1:
    od:
    return [op(ret), ct]:
end:

#Seqcomp on SeqStart(n, Seq1(n) - Seq1(i), Seq1(n) - Seq1(i) + 2)
#and same with n+1 and i+1
AdjSeqComp:=proc(n, i)
    return SeqComp(SeqStart(n, Seq1(n) - Seq1(i),
        Seq1(n) - Seq1(i) + 2), SeqStart(n+1,
        Seq1(n+1) - Seq1(i+1), Seq1(n+1) - Seq1(i+1) + 2)):
end:
#For difference 1, a_n=4a_(n-1)+2a_(n-2)+5
#In all cases, first Seq1(i) terms match, then a mismatch

#SeqStart(n+1, Seq1(n+1) - Seq1(n), Seq1(n))
Sn2:=proc(n)
    return SeqStart(n+1, Seq1(n+1) - Seq1(n), Seq1(n)):
end:

####STUFF BELOW HERE AFTER PROVING SPECIAL TERNARY CASE####

#Legal move sequence for (k+1)-ary Nim sequence
Akn:=proc(k, n)
option remember:
    if n < k then
        return n:
    fi:
    return 3*Akn(k, n-1) - Akn(k, n-2) + 1:
end:

#Fibonacci-like sequence starting with k, 1
Fibk:=proc(k, n)
option remember:
    if n = 1 then
        return k:
    elif n = 2 then
        return 1:
    else
        return Fibk(k, n-1) + Fibk(k, n-2):
    fi:
end:

#Fibonacci-like sequence starting with 1,2,3,...,k
Fibk2:=proc(k, n)
option remember:
    if n <= k then
        return n:
    else
        return Fibk2(k, n-1) + Fibk2(k, n-2):
    fi:
end:

#Beatty sequence
Beatty:=proc(alpha, n)
    return floor(n * alpha):
end:

#Sort decreasing
SortDec:=proc(S)
    return sort(S, `>`):
end:

#Values in a sequence less than or equal to n
#Returned in decreasing order
#f is a function of n>=1 that is strictly increasing
SeqLE:=proc(f, n) local S, cap, i, high, low:
    if n=0 then
        return []:
    fi:
    cap:=1:
    low:=0:
    high:=infinity:
    while true do
        if f(cap + 1) > n and f(cap) <= n then
            return SortDec([seq(f(i), i=1..cap)]):
        elif f(cap + 1) > n then
            high:=cap:
            cap:=floor((low + high) / 2):
        elif high = infinity then
            low:=cap:
            cap:=cap * 2:
        else
            low:=cap:
            cap:=ceil((low + high) / 2):
        fi:
    od:
end:

#Digital expansion
#S list in decreasing order
DigitExpansion:=proc(S, n) local ret, i, val:
    ret:=[0]:
    if n = 0 then
        return ret:
    end:
    val:=n:
    i:=1:
    while i <= nops(S) do
        if val < S[i] then
            if i < nops(S) then
                ret:=[op(ret), 0]:
            fi:
            i:=i + 1:
        else
            ret[i]:=ret[i] + 1:
            val:=val - S[i]:
        fi:
    od:
    return ret:
end:

#Zeckendorf-like expansion for Fibk(k, n) instead of Fibk(1, n)
kZeckendorf:=proc(k, n) local f:
    f:=proc(n)
        return Fibk(k, n+1):
    end:
    return DigitExpansion(SeqLE(f, n), n):
end:

#Zeckendorf-like expansion for Fibk2(k, n) instead of Fibk2(3, n)
#Note: k here equivalent to k-2 in kZeckendorf
kZeckendorf2:=proc(k, n) local f:
    f:=proc(n)
        return Fibk2(k, n):
    end:
    return DigitExpansion(SeqLE(f, n), n):
end:

#Zeckendorf expansion
Zeckendorf:=proc(n):
    return kZeckendorf2(3, n):
end:

#Akn+1 expansion
Aknp1DE:=proc(k, n) local f:
    f:=proc(n)
        if n = 1 then
            return 1:
        else
            return Akn(k, n-1) + 1:
        fi:
    end:
    return DigitExpansion(SeqLE(f, n), n):
end:

#Seq1 expansion
Seq1p1DE:=proc(n) local f:
    f:=proc(n)
        if n = 1 then
            return 1:
        else
            return Seq1(n-3) + 1:
        fi:
    end:
    return DigitExpansion(SeqLE(f, n), n):
end:

#w_n, where k is the max value
Wn:=proc(k, n) local i:
option remember:
    if n = 1 then
        return [seq(i, i=0..(k-1))]:
    else
        return [op(Wn(k, n-1)),
            seq(op(Wn(k, n-i)), i=1..(n-1)), k]:
    fi:
end:

#Sequence of positions with a given value
#Index from 0
ValSeq:=proc(S, val) local i, ret:
    ret:=[]:
    for i from 1 to nops(S) do
        if S[i] = val then
            ret:=[op(ret), i-1]:
        fi:
    od:
    return ret:
end:

#ValSeq on Wn, with appropriate n for cap
#May include values greater than cap, but not too much greater
WValSeq:=proc(k, val, cap) local w, i:
    w:=[]:
    i:=1:
    while nops(w) < cap do
        w:=Wn(k, i):
        i:=i + 1:
    od:
    return ValSeq(w, val):
end:

#Testing arbitrary word composition rules
#positive numbers in list are that wi
#negative numbers are offsets from current
#zeroes mean put mex of w1 in that position
#ranges mean do all of the range, in order
#names mean apply that function to next argument
#(all pieces if range)
#function can take 3 args (must take at least 1):
#1. string to operate on
#2. wi from which it is being called
#3. wi upon which it is being called
#strings mean special things:
#"shield": negative numbers that would be 1 or less are ignored
#(only works on numbers, not ranges)
#seeds is a list of w1 through wk (k=nops(seeds))
#rule is applied beyond k
#returns a procedure to make the words
#If pass 0 to return, gives description
WordComp:=proc(seeds, rule) local do_work, mexguy, fixthing:
    
    mexguy:=convert(mex({op(map(parse,
            convert(seeds[1], list)))}), string):
    
    fixthing:=proc(i, thing):
        if thing < 0 then
            return i + thing:
        else
            return thing:
        fi:
    end:

    do_work:=proc(i) local w, thing, rg, j, f, ft, shield:
    option remember:
        if i = 0 then
            return seeds, rule:
        elif i <= nops(seeds) then
            return seeds[i]:
        fi:
        w:="":
        f:=ident:
        shield:=false:
        for thing in rule do
            if type(thing, name) then
                f:=thing:
            elif type(thing, string) then
                if thing = "shield" then
                    shield:=true:
                fi:
            elif type(thing, range) then
                rg:=op(thing):
                rg:=[fixthing(i, rg[1]), fixthing(i, rg[2])]:
                for j from rg[1] by sign(rg[2]-rg[1]) to rg[2] do
                    w:=cat(w, f(do_work(j), i, j)):
                od:
                f:=ident:
            elif thing = 0 then
                w:=cat(w, f(mexguy)):
                f:=ident:
            else
                ft:=fixthing(i, thing):
                if not(shield) or thing > 0 or ft > 1 then
                    w:=cat(w, f(do_work(ft), i, ft)):
                fi:
                f:=ident:
            fi:
        od:
        return w:
    end:
    
    return do_work:
end:

#Analyze word composition rule
#WC = the rule
#Go through iteration n
AnalyzeWordComp:=proc(WC, n, suppress:=false) local i, L, w,
        zp, tp, p7, zpf, npf, npfs, tpf, msf, p7f,
        val, dexp, li, mx, L2:
    L:=[1]:
    for i from 1 to n-1 do
        L:=[length(WC(i)), op(L)]:
    od:
    w:=WC(n):
    L2:=[length(w), op(L)]:
    zp:=true:
    tp:=true:
    p7:=true:
    npfs:=[]:
    msf:=-1:
    li:=-1:
    mx:=parse(w[-1]):
    for i from 1 to length(w) do
        val:=parse(w[i]):
        if val > msf then
            msf:=val:
            tp:=true:
            npfs:=[op(npfs), -1]:
        fi:
        if i > L2[li] then
            li:=li-1:
        fi:
        dexp:=DigitExpansion(L, i-1):
        if zp and val = 0 and dexp[-1] <> 0 then
            zp:=false:
            zpf:=i-1:
        fi:
        if tp and val = msf and dexp[-1] <> 0 then
            tp:=false:
            tpf:=i-1:
        fi:
        if dexp[-1] = 0 and npfs[val+1] = -1 then
            npfs[val+1] := i-1:
        fi:
        if p7 then
            #if li = -1 then
            #    print(i-1, w[i], L2[li]):
            #else
            #    print(i-1, w[i], L2[li], L[li+1]):
            #fi:
            if i <= mx then
                if w[i] <> convert(i-1, string) then
                    p7:=false:
                    p7f:=i-1:
                fi:
            elif li < -2 and i = L2[li] then
                if w[i] <> convert(mx, string) then
                    p7:=false:
                    p7f:=i-1:
                fi:
            elif w[i - L2[li+1]] <> w[i] then
                p7:=false:
                p7f:=i-1:
            fi:
        fi:
        if not suppress then
            printf("Value: %d\tIndex: %d\tRep: %a\n", val,
                i-1, dexp):
        fi:
    od:
    if zp then
        printf("Zero property holds\n"):
    else
        printf("Zero property first fails in position %d\n", zpf):
    fi:
    if {op(2..-2, npfs)} = {-1} then
        printf("One property holds\n"):
    else
        npf:=infinity:
        for i from 2 to nops(npfs) - 1 do
            if npfs[i] > -1 and npfs[i] < npf then
                npf:=npfs[i]:
            fi:
        od:
        printf("One property first fails in position %d\n", npf):
    fi:
    if tp then
        printf("Two property holds\n"):
    else
        printf("Two property first fails in position %d\n", tpf):
    fi:
    if p7 then
        printf("Proposition 7 holds\n"):
    else
        printf("Proposition 7 first fails in position %d\n", p7f):
    fi:
end:

#Thue-Morse word composition thingy
TMWC:=proc(w, i) local tm:
    tm:=ThueMorse(i, "2", "3"):
    if tm[i-2] = "2" then
        return cat(w, w):
    else
        return cat(w, w, w):
    fi:
end:

#All (labeled) traceable ternary words of length n
Traceables:=proc(n) local L, w, ret, c:
option remember:
    if n = 1 then
        return {["0", 2]}:
    elif n = 2 then
        return {["01", 2]}:
    else
        L:=Traceables(n-1):
        ret:={}:
        for w in L do
            if w[1][-1] = "0" then
                ret:=ret union {[cat(w[1], "1"), w[2]]}:
            elif w[1][-1] = "2" then
                ret:=ret union {[cat(w[1], "0"), w[2]]}:
            else
                c:=w[1][-w[2]]:
                ret:=ret union {[cat(w[1], c), w[2]]}:
                ret:=ret union {[cat(w[1], "2"), n]}:
            fi:
        od:
        return ret:
    fi:
end:

#All (labeled) traceable ternary words of length n
#list form
TraceablesList:=proc(n) local L, w, ret, c:
option remember:
    if n = 1 then
        return [["0", 2]]:
    elif n = 2 then
        return [["01", 2]]:
    else
        L:=TraceablesList(n-1):
        ret:=[]:
        for w in L do
            if w[1][-1] = "0" then
                ret:=[op(ret), [cat(w[1], "1"), w[2]]]:
            elif w[1][-1] = "2" then
                ret:=[op(ret), [cat(w[1], "0"), w[2]]]:
            else
                c:=w[1][-w[2]]:
                ret:=[op(ret), [cat(w[1], c), w[2]]]:
                ret:=[op(ret), [cat(w[1], "2"), n]]:
            fi:
        od:
        return ret:
    fi:
end:

#All (essentially labeled) traceable ternary words of length n
#A traceable is essentially labeled if it is not the same word
#with a label that is a multiple of another
TraceablesEL:=proc(n) local L, w, ret, c, i, rett, s, oks, ok, j:
option remember:
    if n = 1 then
        return {["0", 2]}:
    elif n = 2 then
        return {["01", 2]}:
    else
        L:=TraceablesEL(n-1):
        ret:=[]:
        for w in L do
            if w[1][-1] = "0" then
                ret:=[op(ret), [cat(w[1], "1"), w[2]]]:
            elif w[1][-1] = "2" then
                ret:=[op(ret), [cat(w[1], "0"), w[2]]]:
            else
                c:=w[1][-w[2]]:
                ret:=[op(ret), [cat(w[1], c), w[2]]]:
                ret:=[op(ret), [cat(w[1], "2"), n]]:
            fi:
        od:
        rett:={}:
        s:="":
        oks:=[]:
        for i from 1 to nops(ret) do
            if ret[i][1] <> s then
                s:=ret[i][1]:
                oks:=[]:
            fi:
            ok:=true:
            for j in oks do
                if (ret[i][2] mod j) = 0 then
                    ok:=false:
                    break:
                fi:
            od:
            if ok then
                oks:=[op(oks), ret[i][2]]:
                rett:=rett union {ret[i]}:
            fi:
        od:
        return rett:
    fi:
end:

#All (essentially labeled) traceable ternary words of length n
#A traceable is essentially labeled if it is not the same word
#with a label that is a multiple of another
TraceablesELList:=proc(n) local L, w, ret, c, i, rett, s, oks, ok, j:
option remember:
    if n = 1 then
        return [["0", 2]]:
    elif n = 2 then
        return [["01", 2]]:
    else
        L:=TraceablesELList(n-1):
        ret:=[]:
        for w in L do
            if w[1][-1] = "0" then
                ret:=[op(ret), [cat(w[1], "1"), w[2]]]:
            elif w[1][-1] = "2" then
                ret:=[op(ret), [cat(w[1], "0"), w[2]]]:
            else
                c:=w[1][-w[2]]:
                ret:=[op(ret), [cat(w[1], c), w[2]]]:
                ret:=[op(ret), [cat(w[1], "2"), n]]:
            fi:
        od:
        rett:=[]:
        s:="":
        oks:=[]:
        for i from 1 to nops(ret) do
            if ret[i][1] <> s then
                s:=ret[i][1]:
                oks:=[]:
            fi:
            ok:=true:
            for j in oks do
                if (ret[i][2] mod j) = 0 then
                    ok:=false:
                    break:
                fi:
            od:
            if ok then
                oks:=[op(oks), ret[i][2]]:
                rett:=[op(rett), ret[i]]:
            fi:
        od:
        return rett:
    fi:
end:

#All (unlabeled) traceable ternary words of length n
TraceablesUL:=proc(n) local L, ret, w:
    L:=TraceablesEL(n):
    ret:={}:
    for w in L do
        ret:=ret union {w[1]}:
    od:
    return ret:
end:

#Random traceable strongly Fergusonian (not uniform)
#Never return the periodic 012..k guy
#p is the probability of accepting a number that works
RandTraceable:=proc(n, k:=2, p:=0.5) local r, w, c, i:
    if n <= k then
        w:="":
        for i from 0 to n-1 do
            w:=cat(w, convert(i, string)):
        od:
        return [w, [k, 1]]:
    elif n <= 2*k then
        w:="":
        for i from 0 to k-1 do
            w:=cat(w, convert(i, string)):
        od:
        for i from 0 to n-1-k do
            w:=cat(w, convert(i, string)):
        od:
        return [w, [k, 1]]:
    else
        w:=RandTraceable(n-1, k):
        if parse(w[1][-1]) < k-1 then
            return [cat(w[1],
                convert(parse(w[1][-1])+1, string)), w[2]]:
        elif w[1][-1] = convert(k, string) then
            return [cat(w[1], "0"), w[2]]:
        else
            c:=w[1][-w[2][1]]:
            r:=rand():
            if r/1e12 > p then
                return [cat(w[1], c), w[2]]:
            else
                return [cat(w[1], convert(k, string)),
                    [n, op(w[2])]]:
            fi:
        fi:
    fi:
end:

#Analyze traceable with given tracing sequence
#Go through length n+1
AnalyzeTraceable:=proc(trseq, n, suppress:=false) local i, L, w,
        zp, tp, p7, zpf, npf, npfs, tpf, msf, p7f,
        val, dexp, li, mx, L2, j:
    w:="01":
    j:=1:
    for i from 3 to n+1 do
        if j < nops(trseq) and i-1 = trseq[j+1]-1 then
            w:=cat(w, "2"):
            j:=j+1:
        else
            w:=cat(w, w[i-trseq[j]]):
        fi:
    od:
    L:=[1, op(trseq)]:
    L:=ListTools[Reverse](L):
    L2:=[length(w)+1, op(L)]:
    zp:=true:
    tp:=true:
    p7:=true:
    npfs:=[]:
    msf:=-1:
    li:=-1:
    mx:=2:
    for i from 1 to length(w) do
        val:=parse(w[i]):
        if val > msf then
            msf:=val:
            tp:=true:
            npfs:=[op(npfs), -1]:
        fi:
        if i > L2[li] then
            li:=li-1:
        fi:
        dexp:=DigitExpansion(L, i-1):
        if zp and val = 0 and dexp[-1] <> 0 then
            zp:=false:
            zpf:=i-1:
        fi:
        if tp and val = msf and dexp[-1] <> 0 then
            tp:=false:
            tpf:=i-1:
        fi:
        if dexp[-1] = 0 and npfs[val+1] = -1 then
            npfs[val+1] := i-1:
        fi:
        if p7 then
            #if li = -1 then
            #    print(i-1, w[i], L2[li]):
            #else
            #    print(i-1, w[i], L2[li], L[li+1]):
            #fi:
            if i <= mx then
                if w[i] <> convert(i-1, string) then
                    p7:=false:
                    p7f:=i-1:
                fi:
            elif li < -2 and i = L2[li] then
                if w[i] <> convert(mx, string) then
                    p7:=false:
                    p7f:=i-1:
                fi:
            elif w[i - L2[li+1]] <> w[i] then
                p7:=false:
                p7f:=i-1:
            fi:
        fi:
        if not suppress then
            printf("Value: %d\tIndex: %d\tRep: %a\n", val,
                i-1, dexp):
        fi:
    od:
    if zp then
        printf("Zero property holds\n"):
    else
        printf("Zero property first fails in position %d\n", zpf):
    fi:
    if {op(2..-2, npfs)} = {-1} then
        printf("One property holds\n"):
    else
        npf:=infinity:
        for i from 2 to nops(npfs) - 1 do
            if npfs[i] > -1 and npfs[i] < npf then
                npf:=npfs[i]:
            fi:
        od:
        printf("One property first fails in position %d\n", npf):
    fi:
    if tp then
        printf("Two property holds\n"):
    else
        printf("Two property first fails in position %d\n", tpf):
    fi:
    if p7 then
        printf("Proposition 7 holds\n"):
    else
        printf("Proposition 7 first fails in position %d\n", p7f):
    fi:
end:

#Build word with p values
#P is an array where index [i][j] is p_{i-1, i-j+1}
BuildPWord:=proc(P, k:=2) local W, i, j, l:
    W:=[""]:
    for i from 0 to k-1 do
        W[1]:=cat(W[1], convert(i, string)):
    od:
    for i from 1 to nops(P) do
        W:=[op(W), ""]:
        for j from 1 to nops(P[i]) do
            for l from 1 to P[i][j] do
                W[-1]:=cat(W[-1], W[i-j+1]):
            od:
        od:
        W[-1]:=cat(W[-1], convert(k, string)):
    od:
    return [seq(length(W[i]), i=1..nops(W))], W:
end:

#BuildPWord and then check if traceable
#Only does k=2
#Returns true if traceable, false if not
AnalyzePWord:=proc(P, suppress:=false) local sq, i, j, w:
    sq, w:=BuildPWord(P):
    w:=w[-1]:
    printf("Sequence: %a\n", sq):
    if not suppress then
        printf("Word: %s\n", w):
    fi:
    #AnalyzeTraceable(sq, sq[-1]*2, suppress):
    j:=1:
    for i from 3 to sq[-1] do
        if j < nops(sq) and i = sq[j+1] then
            j:=j+1:
            if w[i] <> "2" then
                printf("Traceability first fails (case 2) in position %d\n", i-1):
                return false:
            fi:
        elif w[i] <> w[i - sq[j]] then
            printf("Traceability first fails (case 3) in position %d\n", i-1):
            return false:
        fi:
    od:
    printf("Traceable\n"):
    return true:
end:

#Testing conjecture that (N+S)^N=empty => (N+V)^N=empty
#Build a set S of size n with (N+S)^N=empty
#S will contain 1 through k-1, but not k
#p is the probability of accepting a number into S that works
BuildSet:=proc(n, k:=2, p:=0.5)
local S, A, V, i, j, s, ok, broken:
    S:={seq(i, i=1..k-1)}:
    A:=[k, 1]:
    i:=k+1:
    while nops(S) < n do
        if DigitExpansion(A, i)[-1] <> 0 then
            i:=i+1:
            next:
        fi:
        S:=S union {i}:
        A:=[i+1, op(A)]:
        ok:=true:
        for j from i to 2*i+1 do
            if DigitExpansion(A, j)[-1] = 0 and
                    DigitExpansion(A, j+1)[-1] = 1 then
                for s in S do
                    if DigitExpansion(A, j+s)[-1] = 0 and
                            DigitExpansion(A, j+s+1)[-1] = 1 then
                        ok:=false:
                        break:
                    fi:
                od:
                if not ok then
                    break:
                fi:
            fi:
        od:
        if rand()/1e12 > p then
            ok:=false:
        fi:
        if not ok then
            S:=S minus {i}:
            A:=[op(2..,A)]:
        fi:
        i:=i+1:
    od:
    broken:=false:
    V:={1}:
    for j from 3 to A[1] do
        if DigitExpansion(A, j)[-1] = 0 then
            if DigitExpansion(A, j+1)[-1] = 0 then
                V:=V union {j}:
            elif broken = false then
                for s in V do
                    if DigitExpansion(A, j+s)[-1] = 0 and
                            DigitExpansion(A, j+s+1)[-1] = 1 then
                        broken:=[j, s]:
                        break:
                    fi:
                od:
            fi:
        fi:
    od:
    return S, A, V, broken:
end:
#Conjecture is false!

#Vals you need to take the mex of to find SG(S[n] + offset)
SGvalsGeneral:=proc(S, offset:=0) local i, vals, seqn, N:
    vals:=[]:
    seqn:=GetPeriodAndPrefix(S)[1]:
    N:=S[-1]+1:
    for i from 1 to nops(S) do
        vals:=[seqn[MapleMod(N + offset - S[i], N)], op(vals)]:
    od:
    return vals:
end:

#Build a set S of size n
#S will contain 1 through k-1, but not k
#A will be tracing sequence of Nim sequence of S
#p is the probability of accepting a number into S that works
BuildTraceable:=proc(n, k:=2, p:=0.5)
local r, w, c, i, per, pref:
    if n <= k then
        w:="":
        for i from 0 to n-1 do
            w:=cat(w, convert(i, string)):
        od:
        return [w, [k, 1], {seq(k-i, i=1..k-1)}]:
    elif n <= 2*k then
        w:="":
        for i from 0 to k-1 do
            w:=cat(w, convert(i, string)):
        od:
        for i from 0 to n-1-k do
            w:=cat(w, convert(i, string)):
        od:
        return [w, [k, 1], {seq(k-i, i=1..k-1)}]:
    else
        w:=BuildTraceable(n-1, k, p):
        if parse(w[1][-1]) < k-1 then
            return [cat(w[1],
                convert(parse(w[1][-1])+1, string)), w[2], w[3]]:
        elif w[1][-1] = convert(k, string) then
            return [cat(w[1], "0"), w[2], w[3]]:
        else
            c:=w[1][-w[2][1]]:
            r:=rand():
            if r/1e12 > p then
                return [cat(w[1], c), w[2], w[3]]:
            else
                per, pref:=GetPeriodAndPrefix({n-1} union w[3]):
                #print(n, w[2], w[3], per, pref):
                if nops(per) = n and max(per) <= k 
                        and nops(pref) = 0 then
                    return [cat(w[1], convert(k, string)),
                        [n, op(w[2])], {n-1} union w[3]]:
                else
                    return [cat(w[1], c), w[2], w[3]]:
                fi:
            fi:
        fi:
    fi:
end:

#Invoke SGvalsGeneral on BuildTraceable
#check if everybody is well-behaved
SGvalsGeneralizeBT:=proc(n, k:=2, p:=0.5)
local w, i, j, v, L, L2, LS, ok:
    w:=BuildTraceable(n, k, p):
    L:=[seq(SGvalsGeneral(w[-1], i), i=0..n)]:
    ok:=true:
    for i from 1 to nops(L) do
        L2:=L[i][1..-(k-1)]:
        LS:={seq(j, j in L2)}:
        v:=mex(L[i]):
        if v in {0, k} and not(LS subset
                {(v+1) mod (k+1), (v+2) mod (k+1)}) then
            ok:=false:
            printf("First failure at position %d, guys: %a\n",
                i-1, L[i]):
            break:
        fi:
    od:
    return w, L, ok:
end:
#Conjecture is false: {1,2,9,13,17}

#Greedy Algorithm redux
#Maybe faster?
Greedy2:=proc(n, k:=2)
local w, A, S, c, i, per, pref:
    w:="":
    for i from 0 to k-1 do
        w:=cat(w, convert(i, string)):
    od:
    A:=[k,1]:
    S:={seq(k-i, i=1..k-1)}:
    if n < k then
        return w[1..n], A, S:
    fi:
    for i from k+1 to n+1 do
        if parse(w[-1]) < k-1 then
            w:=cat(w, convert(parse(w[-1])+1, string))
        elif w[-1] = convert(k, string) then
            w:=cat(w, "0"):
        else
            c:=w[-A[1]]:
            if i <= A[1]*2 then
                w:=cat(w, c):
            else
                per, pref:=GetPeriodAndPrefix({i-1} union S):
                #print(n, w[2], w[3], per, pref):
                if nops(per) = i and max(per) <= k 
                        and nops(pref) = 0 then
                    w:=cat(w, convert(k, string)):
                    A:=[i, op(A)]:
                    S:={i-1} union S:
                else
                    w:=cat(w, c):
                fi:
            fi:
        fi:
    od:
    return w, A, S:
end:

#Change word or list w from its k value to value k
#Modify S accordingly, if S is passed
Changek:=proc(w, k, S:=false)
local lst, ww, i, j, retw, retS, L, oldk, li, rw:
    if type(w, string) then
        ww:=[]:
        for i from 1 to length(w) do
            ww:=[op(ww), parse(w[i])]:
        od:
    else
        ww:=w:
    fi:
    oldk:=max(ww):
    retw:=[]:
    if S <> false then
        L:=[]:
        for i in S do
            if i > oldk then
                L:=[op(L), i]:
            fi:
        od:
        L:=sort(L):
        L:=[op(L), infinity]:
        li:=1:
        retS:={seq(i, i=1..k-1)}:
    fi:
    for i from 1 to nops(ww) do
        if oldk < k then
            if ww[i] < oldk - 1 then
                retw:=[op(retw), ww[i]]:
                if S <> false and i = L[li]+1 then
                    retS:=retS union {nops(retw) - 1}:
                    li:=li + 1:
                fi:
            elif ww[i] = oldk - 1 then
                retw:=[op(retw), seq(j, j=ww[i]..k-1)]:
                if S <> false and i = L[li]+1 then
                    retS:=retS union {seq(j,
                        j=nops(retw)-(k-ww[i])..nops(retw)-1)}:
                    li:=li + 1:
                fi:
            else
                retw:=[op(retw), k]:
                if S <> false and i = L[li]+1 then
                    retS:=retS union {nops(retw) - 1}:
                    li:=li + 1:
                fi:
            fi:
        else
            if ww[i] < k then
                retw:=[op(retw), ww[i]]:
                if S <> false and i = L[li]+1 then
                    retS:=retS union {nops(retw) - 1}:
                    li:=li + 1:
                fi:
            elif ww[i] = oldk then
                retw:=[op(retw), k]:
                if S <> false and i = L[li]+1 then
                    retS:=retS union {nops(retw) - 1}:
                    li:=li + 1:
                fi:
            elif S <> false and i = L[li]+1 then
                li:=li + 1:
            fi:
        fi:
    od:
    if type(w, string) then
        rw:="":
        for i from 1 to nops(retw) do
            rw:=cat(rw, convert(retw[i], string)):
        od:
    else
        rw:=retw:
    fi:
    if S = false then
        return rw:
    else
        return rw, retS:
    fi:
end:

#Make a traceable word with given tracing sequence
#Length n (or length based on tracing sequence if omitted)
Traceable:=proc(trseq, n:=false) local len, w, i, k, ti:
    if n = false then
        len:=max(trseq):
    else
        len:=n:
    fi:
    k:=trseq[2]:
    ti:=2:
    w:=[]:
    for i from 1 to len do
        if i <= k then
            #w:=cat(w, convert(i-1, string)):
            w:=[op(w), i-1]:
        elif ti < nops(trseq) and i = trseq[ti+1] then
            #w:=cat(w, convert(k, string)):
            w:=[op(w), k]:
            ti:=ti + 1:
        else
            #w:=cat(w, w[i - trseq[ti]]):
            w:=[op(w), w[i - trseq[ti]]]:
        fi:
    od:
    return w:
end:

#j-Promotion of sequence
Promotion:=proc(sq, j:=1) local rsq, ret, i, l, drep:
option remember:
    ret:=[]:
    rsq:=ListTools[Reverse](sq, opt=descending):
    for i from 1 to nops(sq) do
        if i <= j then
            ret:=[op(ret), sq[i]]:
        else
            drep:=DigitExpansion(rsq, sq[i] - 1):
            drep:=ListTools[Reverse](drep):
            #print(i, sq[i]-1, drep, ret):
            if drep[j] > 0 then
                drep[j] := drep[j] + 1:
            fi:
            ret:=[op(ret), 1+add(drep[l]*ret[l], l=1..i-1)]:
        fi:
    od:
    return ret:
    #return ListTools[Reverse](ret):
end:

#j-Promotion function of value wrt sequence
PromotionFn:=proc(n, sq, j:=1) local promo, rsq, drep:
    promo:=Promotion(sq):
    rsq:=ListTools[Reverse](sq, opt=descending):
    drep:=DigitExpansion(rsq, n):
    drep:=ListTools[Reverse](drep):
    if drep[j] > 0 then
        drep[j] := drep[j] + 1:
    fi:
    return add(drep[l]*promo[l], l=1..nops(drep)):
end:

#Change k, using new morphism format
#Always increase k by 1
Increasek:=proc(w) local ww, rw, k, i, retw:
    if type(w, string) then
        ww:=[]:
        for i from 1 to length(w) do
            ww:=[op(ww), parse(w[i])]:
        od:
    else
        ww:=w:
    fi:
    k:=max(ww):
    retw:=[]:
    for i from 1 to nops(ww) do
        if ww[i] = 0 then
            retw:=[op(retw), 0]:
        fi:
        retw:=[op(retw), ww[i]+1]:
    od:
    if type(w, string) then
        rw:="":
        for i from 1 to nops(retw) do
            rw:=cat(rw, convert(retw[i], string)):
        od:
    else
        rw:=retw:
    fi:
    return rw:
end:

#Invoke SGvalsGeneral on BuildTraceable
#check if everybody is well-behaved
#This is trying to prove Theorem 7 (or at least not disprove)
SGvalsGeneralizeBT2:=proc(n, k:=2, p:=0.5)
local w, v, vw, m, i2, ok, i, wL, mS, s:
    w:=BuildTraceable(n, k, p):
    wL:=ListTools[Reverse](w[2]):
    v:=Promotion(wL):
    #vS:={1, seq(PromotionFn(i, wL), i in w[-1])}:
    vw:=GetPeriodAndPrefix(vS)[1]:
    #L:=[seq(SGvalsGeneral(w[-1], i), i=0..n)]:
    #L2:=[seq(SGvalsGeneral(vS, i), i=0..length(Increasek(w[1])))]:
    #print(wL, v):
    ok:=true:
    for i from 0 to length(w[1]) do
        mS:={}:
        for s in w[-1] do
            if s <= i then
                mS:=mS union {s}:
            fi:
        od:
        m:=mex(mS):
        i2:=Promotion(i, wL):
        for s in w[-1] do
            if s <= i then
                if parse(w[1][i+1-s]) <> vw[i2+1-PromotionFn(s, wL)] then
                    ok:=false:
                    printf("First failure at position %d/%d\n",
                        i-1, i2-1):
                    break:
                fi:
            fi:
        od:
    od:
    return wL, v, ok:
end:

#Random sequence of length n, starting with 1
RandSeq:=proc(n, p:=0.3) local sq, i:
    sq:=[1]:
    i:=2:
    while true do
        if rand()/1e12 < p then
            sq:=[op(sq), i]:
            if nops(sq) = n then
                return sq:
            fi:
        fi:
        i:=i+1:
    od:
end:

#Test if Increasek(Traceable(sq)) = Traceable(Promotion(sq))
#Length n, prob p (params for RandSeq)
Promotest:=proc(n, p:=0.3) local sq:
    sq:=RandSeq(n, p):
    return sq, Promotion(sq), evalb(Increasek(Traceable(sq)) =
        Traceable(Promotion(sq))):
end:

#Representations, values, indices in traceable
TraceableReps:=proc(trseq, n:=false) local len, w, i, rsq:
    if n <> false then
        len:=n:
    else
        len:=trseq[-1]-1:
    fi:
    w:=Traceable(trseq, n):
    rsq:=ListTools[Reverse](trseq):
    for i from 1 to nops(w) do
        printf("Val: %d\tInd: %d\tRep: %a\n", w[i], i-1,
            DigitExpansion(rsq, i-1)):
    od:
end:

#Comparing w[i-s] with v[f(i)-f(s)]
SubComp:=proc(L, i, suppress:=false) local v, w, s, ip, sp, F:
    w:=Traceable(L):
    v:=Increasek(w):
    ip:=PromotionFn(i,L):
    F:={}:
    for s from 1 to i do
        sp:=PromotionFn(s, L):
        if not(suppress) then
            printf("w[%d-%d]=%d,\tv[f(%d)-f(%d)]=v[%d-%d]=%d\n",
                i, s, w[i+1-s], i, s, ip, sp, v[ip+1-sp]):
        fi:
        if w[i+1-s] <> v[ip+1-sp] then
            F:=F union {s}:
            if not(suppress) and w[i+1-s] <> 0 then
                printf("*******\n"):
            fi:
        fi:
    od:
    return F:
end:

#Get S from L
SetFrom:=proc(L) local i:
    return {seq(i,i=1..L[2]-1),seq(L[i]-1, i=3..nops(L))}:
end:

#Look for places where w[i]=k, w[i-s]<>v[f(i)-f(s)],
#And we can add s to our subtraction set without breaking it
SearchForPotentialSubCompFail:=proc(n, k:=2, test:=0, suppress:=true)
local w, L, S, i, F, ok, tS:
    w:=BuildTraceable(n, k):
    L:=ListTools[Reverse](w[2]):
    S:=w[3]:
    w:=Traceable(L):
    F:={}:
    tS:={}:
    if type(test, integer) then
        tS:={test}:
    else
        tS:=test:
    fi:
    for i from 1 to nops(w) do
        if w[i] in tS then
            F:=F union SubComp(L, i-1, suppress):
            if not(suppress) then
                printf("\n"):
            fi:
        fi:
    od:
    ok:=true:
    for i in F do
        if GetPeriodAndPrefix(S union {i}) = w then
            ok:=false:
            printf("Check out %d further\n", i):
        fi:
    od:
    return L, F, ok:
end:

#Look for s with w[s]=1 to add to S
SearchForOne:=proc(S) local st, w, i:
    if type(S, list) then
        if nops(S) = 1 then
            st:=BuildTraceable(S[1])[-1]:
        else
            st:=BuildTraceable(S[1], S[2])[-1]:
        fi:
        print(st):
    else
        st:=S:
    fi:
    w:=GetPeriodAndPrefix(st)[1]:
    for i from 3 to nops(w) do
        if w[i] = 1 and
                GetPeriodAndPrefix(st union {i-1})[1] = w then
            return i-1:
        fi:
    od:
    return FAIL:
end:

#Check promotion and subtraction
#Do this for all zends
PromoSubtract:=proc(L)
local Lp, rL, rLp, i, ip, j, jp, repi, repj, repij,
        repip, repjp, repijp, S:
    Lp:=Promotion(L):
    rL:=ListTools[Reverse](L):
    rLp:=ListTools[Reverse](Lp):
    S:={seq(i, i=1..L[-1])}:
    for i from 1 to L[-1] do
        repi:=DigitExpansion(rL, i):
        if repi[-1] = 0 then
            ip:=PromotionFn(i, L):
            repip:=DigitExpansion(rLp, ip):
            for j from 1 to i do
                repj:=DigitExpansion(rL, j):
                repij:=DigitExpansion(rL, i-j):
                jp:=PromotionFn(j, L):
                repjp:=DigitExpansion(rLp, jp):
                repijp:=DigitExpansion(rLp, ip-jp):
                printf("i=%d=%a,\tj=%d=%a,\t", i, repi, j, repj):
                printf("i-j=%d=%a\n", i-j, repij):
                printf("ip=%d=%a,\tjp=%d=%a,\t",
                    ip, repip, jp, repjp):
                printf("ip-jp=%d=%a\n\n", ip-jp, repijp):
                if repij <> repijp then
                    S:=S minus {j}:
                fi:
            od:
        fi:
    od:
    return S:
end:

#Sequence obtained from skipping
SkipSeq:=proc(L, skip:=2) local i:
    return ListTools[Reverse]([L[-1],
            seq(L[-skip*i], i=1..floor(nops(L)/skip))]):
end:

#Compare representations when skipping elements of L
RepSkip:=proc(L, cap, skip:=2) local skips, i, s, LS:
    if type(skip, integer) then
        skips:=[skip]:
    else
        skips:=skip:
    fi:
    LS:=table():
    for s in skips do
        LS[s]:=SkipSeq(L, s):
        printf("Skip %d: %a\n", s, LS[s]):
    od:
    printf("\n"):
    for i from 0 to cap do
        printf("%d:\tSkip 1: %a", i, DigitExpansion(L, i)):
        for s in skips do
            printf("\tSkip %d: %a", s, DigitExpansion(LS[s], i)):
        od:
        printf("\n"):
    od:
end:

#Compare representations
CompRep:=proc(Ls, cap) local i, j, L:
    for i from 0 to cap do
        printf("%d:\t", i):
        for j from 1 to nops(Ls) do
            printf("\tLs[%d]: %a", j, DigitExpansion(Ls[j], i)):
        od:
        printf("\n"):
    od:
end:
#Non-volatile zend in skip vs end in even # zeroes no skip
#No relation visible

#Given a representation word with Nim sequence property
#can it always be extended by a(n)=3a(n-1)-a(n-2)
#L=rep. sequence (must grow at least double each term)
#cap=number of sums to test
ExtendRepNimSeqTest:=proc(L, cap, n:=3) local i, j, L2, S:
    L2:=L:
    for i from 1 to cap do
        L2:=[n*L2[1]-(n-2)*L2[2], op(L2)]:
        S:={seq(L2[j]-1, j=1..nops(L2)-1)} union
            {seq(j, j=1..L2[-2]-1)}:
        print(S):
        if nops(GetPeriodAndPrefix(S)[2]) <> 0 then
            return L2, S, false:
        fi:
    od:
    return true:
end: