read(`code_finsubgames.txt`):

Help_wythoffsums:=proc()
    print(` ExtForm(S) , DirSum(Glist) `):
    print(` WythSum(Glist, LegalMoves, LegalPileCounts) `):
    print(` LegalMovesESG(G, counts, restrict:=false) `):
    print(` SGESG(G, counts) , SGSeqDS2(S, cap) `):
    print(` SGSeqWS2Legal(S, cap) , SGSeqWS2All(S, cap) `):
    print(` SG2Seq(G, cap) , SGFSGgenseq(S,cap) `):
    print(` SGFSGgenseqext(S,cap,L) , SGSubarray(L, vals, x) `):
    print(` SGTableToStr(tbl, full:=false) `):
    print(` IsPossibleEventualPeriodAssumptionsGen(per,startp,S,piles,x) `):
    print(` GetPeriodAndPrefixGen(S, piles, x) `):
    print(` SGPP(S, counts, index, x) `):
    print(` GetFormattedDataGen(S,maxp) `):
    print(` GenerateSelfWythPeriodData(cap,maxp,out,append:=true) `):
end:

#Extended subtraction game format: [Slist, metadata]
#Slist is a list of subtraction sets, one for each pile
#metadata encodes the Wythoff couplings
#metadata is a list of lists, each of which has the format
#[piles, moves, counts]
#piles is a set of sets denoting the sets of piles under
#consideration, moves denotes
#the legal Wythoff-type moves for those piles, and counts
#denotes the legal pile counts for those moves
#moves can instead by one of the following:
#"LEGALII", "LEGALI", "LEGALIIM", "LEGALIM", "LEGALUM", "LEGALU",
#"LEGALUUM", "LEGALUU", "ALL"
#"LEGALII" will allow a move iff it is legal in ALL piles
#"LEGALIIM" will allow a move iff it is legal in ALL piles and, if
#a compound Wythoff-type move is being played, it is also legal
#in all sub-moves (no longer a true "Wythoff sum")
#"LEGALI" will allow a move iff it is legal in ALL piles it is
#applied to
#"LEGALIM" will allow a move iff it is legal in ALL piles it is
#made in (no longer a true "Wythoff sum")
#"LEGALUM" will allow a move iff it is legal in SOME pile it is
#made in (no longer a true "Wythoff sum")
#"LEGALU" will allow a move iff it is legal in SOME pile it is
#applied to
#"LEGALUUM" will allow a move iff it is legal in SOME pile or, if
#a compound Wythoff-type move is being played, it is legal
#in some sub-moves (no longer a true "Wythoff sum")
#"LEGALUU" will allow a move iff it is legal in SOME pile
#"ALL" will allow any move
#counts can instead be one of the following:
#"ALL"
#"ALL" will allow any number of piles

#Convert a subtraction set into extended subtraction game format
#An extended subtraction game is an element of the closure of the
#set of subtraction games under direct and Wythoff sums
#S is a subtraction set
ExtForm:=proc(S)
    return [[S],[]]:
end:

#Direct sum of subtraction games
#Glist is a list of extended subtraction games
DirSum:=proc(Glist) local Slist, i, j, k, l, entry, piles,
        pls, offset, metadata, tmdata:
    Slist:=[]:
    metadata:=[]:
    for i from 1 to nops(Glist) do
        tmdata:=Glist[i][2]:
        offset:=nops(Slist):
        Slist:=[op(Slist), op(Glist[i][1])]:
        for j from 1 to nops(tmdata) do
            entry:=tmdata[j]:
            piles:={}:
            for pls in entry[1] do
                piles:=piles union {{seq(l + offset, l in pls)}}:
            od:
            metadata:=[op(metadata), [piles, entry[2], entry[3]]]:
        od:
    od:
    return [Slist, metadata]:
end:

#Wythoff sum of subtraction games
#Glist is a list of extended subtraction games
#LegalMoves is a set of positive integers denoting the numbers
#of counters that can be removed in a Wythoff-type move
#LegalPileCounts is a set of integers greater than 1 denoting
#the number of piles that can be encompassed by a Wythoff-type move
WythSum:=proc(Glist, LegalMoves, LegalPileCounts) local Slist, i,
        j, k, l, entry, piles, pils, pls, offset, metadata,
        tmdata:
    Slist:=[]:
    metadata:=[]:
    piles:={}:
    for i from 1 to nops(Glist) do
        tmdata:=Glist[i][2]:
        offset:=nops(Slist):
        piles:=piles union {{seq(j, j = offset + 1..offset +
            nops(Glist[i][1]))}}:
        Slist:=[op(Slist), op(Glist[i][1])]:
        for j from 1 to nops(tmdata) do
            entry:=tmdata[j]:
            pils:={}:
            for pls in entry[1] do
                pils:=pils union {{seq(l + offset, l in pls)}}:
            od:
            metadata:=[op(metadata), [pils, entry[2], entry[3]]]:
        od:
    od:
    metadata:=[op(metadata), [piles, LegalMoves,
        LegalPileCounts]]:
    return [Slist, metadata]:
end:

#Legal moves for an extended subtraction game
#counts is the number of counters in each pile
#counts can be an integer, in which case the return values are the
#legal subtractions, assuming all counts are at most counts
#counts can also be false, in which case the return value is the
#set of sets of coupled piles
#restrict can be a set of piles, in which case
#legal moves are only computed for moving in those piles
LegalMovesESG:=proc(G, counts, restrict:=false)
local ret, i, j, k, s, legalii, legaluu, piles, moves, cts, st,
        pls, all, restmoves, ok, ct, comb, noor, mv, pl, ctr,
        ctrval, l, m, movesr, msupp:
option remember:
    ret:={}:
    legalii:=false:
    legaluu:=false:
    all:=false:
    movesr:=moves:
    #Regular moves
    for i from 1 to nops(G[1]) do
        if restrict <> false and not(i in restrict) then
            next:
        fi:
        for s in G[1][i] do
            if counts = false then
                ret:=ret union {{i}}:
            elif type(counts, integer) and counts >= s then
                ret:=ret union {[seq(0, j=1..i-1), s,
                    seq(0, j=i+1..nops(G[1]))]}:
            elif counts[i] >= s then
                ret:=ret union {[op(1..i-1, counts),
                    counts[i] - s, op(i+1..nops(counts),
                    counts)]}:
            fi:
        od:
    od:
    #Wythoff-type moves
    for i from 1 to nops(G[2]) do
        piles:=G[2][i][1]:
        if restrict <> false then
            ok:=true:
            for pls in piles do
                for j in pls do
                    if not(j in restrict) then
                        ok:=false:
                        break:
                    fi:
                od:
                if not ok then
                    break:
                fi:
            od:
            if not ok then
                next:
            fi:
        fi:
        moves:=G[2][i][2]:
        cts:=G[2][i][3]:
        if moves = "LEGALII" or moves = "LEGALIIM" then
            if legalii = false then
                legalii:=G[1][1]:
                for j from 2 to nops(G[1]) do
                    legalii:=legalii intersect G[1][j]:
                od:
            fi:
            moves:=legalii:
        elif moves = "LEGALUU" or moves = "LEGALUUM" then
            if legaluu = false then
                legaluu:=G[1][1]:
                for j from 2 to nops(G[1]) do
                    legaluu:=legaluu union G[1][j]:
                od:
            fi:
            moves:=legaluu:
        elif moves = "LEGALI" or moves = "LEGALIM" then
            st:=false:
            for pls in piles do
                for j in pls do
                    if st = false then
                        st:=G[1][j]:
                    else
                        st:=st intersect G[1][j]:
                    fi:
                od:
            od:
            moves:=st:
        elif moves = "LEGALU" or moves = "LEGALUM" then
            st:={}:
            for pls in piles do
                for j in pls do
                    st:=st union G[1][j]:
                od:
            od:
            moves:=st:
        elif moves = "ALL" then
            if all = false then
                if type(counts, list) then
                    all:={seq(j, j=1..max(counts))}:
                else
                    all:={seq(j, j=1..counts)}:
                fi:
            fi:
            moves:=all:
        fi:
        restmoves:=[]:
        for pls in piles do
            restmoves:=[op(restmoves), LegalMovesESG(G, false,
                pls)]:
        od:
        if cts="ALL" then
            cts:={seq(j, j=2..nops(G[1]))}:
        fi:
        for ct in cts do
            noor:= nops(restmoves):
            if ct > noor then
                break:
            fi:
            comb:=firstcomb(noor, ct):
            while comb <> FAIL do
                pl:=[]:
                for j in comb do
                    pl:=[op(pl), restmoves[j]]:
                od:
                ctr:=InitCompCounter([seq(nops(pl[j]),
                    j=1..nops(pl))]):
                while ctr <> false do
                    ctrval:=ValCompCounter(ctr):
                    if counts = false then
                        mv:={}:
                        for l from 1 to nops(ctrval) do
                            for m in pl[l][ctrval[l]] do
                                mv:=mv union {m}:
                            od:
                        od:
                        ret:=ret union {mv}:
                    else
                        msupp:={}:
                        if movesr = "LEGALUUM" or
                                movesr = "LEGALUM" then
                            for l from 1 to nops(ctrval) do
                                for m in pl[l][ctrval[l]] do
                                    msupp:=msupp union G[1][m]:
                                od:
                            od:
                        fi:
                        for k in moves union msupp do
                            if type(counts, integer) and
                                    counts < k then
                                next:
                            elif type(counts, integer) then
                                mv:=[seq(0, l=1..nops(G[1]))]:
                            else
                                mv:=counts:
                            fi:
                            ok:=true:
                            for l from 1 to nops(ctrval) do
                                for m in pl[l][ctrval[l]] do
                                    if mv[m] < k and type(counts,
                                            list) then
                                        ok:=false:
                                        break:
                                    elif (movesr = "LEGALIIM" or
                                            movesr = "LEGALIM")
                                            and not(k in G[1][m])
                                            then
                                        ok:=false:
                                        break:
                                    elif type(counts, list)
                                            then
                                        mv[m]:=mv[m] - k:
                                    else
                                        mv[m]:=k:
                                    fi:
                                od:
                                if not ok then
                                    break:
                                fi:
                            od:
                            if ok then
                                ret:=ret union {mv}:
                            fi:
                        od:
                    fi:
                    ctr:=IncCompCounter(ctr):
                od:
                comb:=nextcomb(comb, noor);
            od:
        od:
    od:
    return ret:
end:

SGESG:=proc(G, counts) local i:
option remember:
    return mex({seq(SGESG(G, i), i in LegalMovesESG(G, counts))}):
end:

SGSeqDS2:=proc(S, cap) local i, j, G, SE:
    SE:=ExtForm(S):
    G:=DirSum([SE, SE]):
    return [seq([seq(SGESG(G, [i, j]), j=0..i)], i=0..cap)]:
end:

SGSeqWS2Legal:=proc(S, cap) local i, j, G, SE:
    SE:=ExtForm(S):
    G:=WythSum([SE, SE], "LEGALI", "ALL"):
    return [seq([seq(SGESG(G, [i, j]), j=0..i)], i=0..cap)]:
end:

SGSeqWS2All:=proc(S, cap) local i, j, G, SE:
    SE:=ExtForm(S):
    G:=WythSum([SE, SE], "ALL", "ALL"):
    return [seq([seq(SGESG(G, [i, j]), j=0..i)], i=0..cap)]:
end:

SG2Seq:=proc(G, cap) local i, j:
    return [seq([seq(SGESG(G, [i, j]), j=0..cap)], i=0..cap)]:
end:

#Sprague-Grundy sequence for a given generalized subtraction set S
#(cap is a list with the number of piles of entries,
#each entry is the cap for that dimension of the sequence)
#The return value is a multidimensional list encoding the
#Sprague-Grundy values
SGFSGgenseq:=proc(S,cap)
    return SGFSGgenseqext(S,cap,[]):
end:

#Sprague-Grundy sequence for a given generalized subtraction set S
#(cap is a list with the number of piles of entries,
#each entry is the cap for that dimension of the sequence)
#The return value is a multidimensional list encoding the
#Sprague-Grundy values
#The argument L is an upper left corner of the sequence of values
#you don't want to recompute
SGFSGgenseqext:=proc(S,cap,L) local ret,vals,ctr,s,vctr,rng,idex,i:
    rng:=(n) -> 0..n:
    ret:=Array(op(map(rng, cap)), L, fill=-1):
    ctr:=InitCompCounter(cap, 0):
    while ctr <> false do
        vctr:=ValCompCounter(ctr):
        if ret[op(vctr)] = -1 then
            vals:={}:
            for s in S do
                idex:=vctr - s:
                if {seq(evalb(idex[i] >= 0),i=1..nops(idex))} =
                        {true} then
                    vals:=vals union {ret[op(idex)]}:
                fi:
            od:
            ret[op(vctr)]:=mex(vals):
        fi:
        ctr:=IncCompCounter(ctr):
    od:
    return ret:
end:

#Take a sub-array where some piles are fixed
#vals is the list of pile counts,
#but x is used for a variable pile
SGSubarray:=proc(L, vals, x) local dims, i:
    dims:=[ArrayDims(L)]:
    for i from 1 to nops(dims) do
        if vals[i] <> x then
            dims[i]:=vals[i]:
        fi:
    od:
    return L[op(dims)]:
end:

SGTableToStr:=proc(tbl, full:=false) local f, tbl2, i, j:
    tbl2:=convert(tbl, list, nested=true):
    if full then
        for i from 1 to nops(tbl) do
            for j from nops(tbl[i])+1 to nops(tbl) do
                tbl2[i]:=[op(tbl2[i]), tbl2[j][i]]:
            od:
        od:
    fi:
    f:=(lst) -> Join(map(convert, lst, string), " "):
    return Join(map(f, tbl2), "\n"):
end:

#Is is possible that per is an eventual period of the subtraction
#game with subtraction set S?
#Assumptions: period already known to occur as an eventual period
#at least twice for some value
#If this is not true, this may not work!
#piles are the values of the piles.  There should be one field
#that is x (the variable field)
IsPossibleEventualPeriodAssumptionsGen:=proc(per,startp,S,piles,x)
local chk,i,vals,seed,s,index,LL,pstart,plcm,
        period,prefix,newpiles,unx,addme:
    if S={} then
        return evalb({seq(per[i],i=1..nops(per))}={0}):
    fi:
    
    for i from 1 to nops(piles) do
        if piles[i] = x then
            index:=i:
            break:
        fi:
    od:
    
    unx:=proc(arg)
        if arg = x then
            return 0:
        else
            return arg:
        fi:
    end:
    
    seed:=piles:
    seed[index]:=0:
    pstart:=0:
    plcm:=1:
    for s in S do
        seed[index]:=max(seed[index], s[index]):
        
        newpiles:=seed - s:
        newpiles[index]:=x:
        if newpiles <> piles and {seq(evalb(
                unx(newpiles[i]) >= 0),
                i=1..nops(newpiles))} = {true} then
            #Fetch other periods/prefixes
            period,prefix:=GetPeriodAndPrefixGen(S, newpiles, x):
            #Get the period lengths/where they start
            pstart:=max(pstart, nops(prefix)):
            plcm:=ilcm(plcm, nops(period)):
        fi:
    od:
    
    if seed[index]=0 then
        return evalb({seq(per[i],i=1..nops(per))}={0}):
    elif nops(chk)>=seed[index]+1 then
        return true:
    fi:
    
    chk:=per:
    while nops(chk)<=2*seed[index]*plcm+pstart do
        chk:=[op(chk),op(per)]:
    od:
    
    addme:=(startp - 1) mod nops(per):
    if addme > 0 then
        chk:=[op(per[-addme..-1]),op(chk)]:
    fi:
    
    #print("chk:", chk):
    #print("stt:", startp):
    #print("piles:", piles):
    #print("per:", per):
    #print("start:", seed[index]*plcm+pstart+1):
    
    for i from seed[index]*plcm+pstart+1 to
            2*seed[index]*plcm+pstart do
        vals:={}:
        #print("i now equals:", i):
        for s in S do
            newpiles:=seed - s:
            newpiles[index]:=x:
            #print(s, "s"):
            #print(newpiles):
            if newpiles = piles then
                #print(i - s[index], "inn"):
                vals:=vals union {chk[i-s[index]]}:
                #print(vals):
            else
                newpiles[index]:=i-s[index]:
                #print(newpiles, "in"):
                if {seq(evalb(newpiles[i] >=0),
                        i=1..nops(newpiles))} = {true} then
                    vals:=vals union {SGPP(S, newpiles,
                        index, x)}:
                fi:
                #print(vals):
            fi:
        od:
        #print(vals):
        if mex(vals)<>chk[i] then
            return false:
        fi:
    od:
    return true:
end:

#Given a subtraction set S of a generalized subtraction game,
#along with a list of pile counts with one unknown, x
#determine its (normalized) period and its prefix
GetPeriodAndPrefixGen:=proc(S, piles, x)
local L,cap,periods,per,k,
        capconst,perconst,#Always good to avoid magic numbers!
        seed,s,i,index,LL,startp,perdata:
option remember:#So we can call this multiple times without
                #needing to recompute
    
    capconst:=2:
    perconst:=2:
    
    if S={} then
        return [0],[]:
    fi:
    
    for i from 1 to nops(piles) do
        if piles[i] = x then
            index:=i:
            break:
        fi:
    od:
    
    seed:=piles:
    seed[index]:=0:
    for s in S do
        seed[index]:=max(seed[index], s[index]):
    od:
    
    if seed[index] = 0 then
        return [0],[]:
    fi:
    
    cap:=capconst*seed:
    L:=SGFSGgenseq(S,cap):
    while true do
        LL:=convert(SGSubarray(L, piles, x), list):
        periods:=GetEventualPeriods(LL,perconst,true):
        for perdata in periods do
            per,startp:=op(perdata):
            if IsPossibleEventualPeriodAssumptionsGen(per,startp,
                    S,piles,x) then
                k:=max(2,ceil(2*seed[index]/nops(per))):
                if HasEventualPeriodT(LL,nops(per),k) then
                    return NormalizeEventualPeriod(LL,per):
                fi:
            fi:
        od:
        
        cap:=2*cap:
        L:=SGFSGgenseqext(S,cap,L):
    od:
end:

#Use period and prefix to get SG value
SGPP:=proc(S, counts, index, x) local period, prefix, piles, n:
    piles:=counts:
    n:=piles[index]:
    piles[index]:=x:
    period,prefix:=GetPeriodAndPrefixGen(S, piles, x):
    if n < nops(prefix) then
        return prefix[n+1]:
    else
        return period[((n - nops(prefix)) mod nops(period)) + 1]:
    fi:
end:

#Given a subtraction set S, return its double self Wythoff
#period data in the form:
#S [numbers in set in increasing order]
#F# [prefix for # in other pile]
#P# [period for # in other pile]
#G# [prefix length for # in other pile]
#L# [period length for # in other pile]
#E [minimal contraction 1]
#E [minimal contraction 2]
#etc.
#(here, brackets do not denote lists, but rather data
#field descriptions)
GetFormattedDataGen:=proc(S,maxp) local E,get,s,i,x:
    get:=[seq([GetPeriodAndPrefixGen(
        {seq(op([[s,0],[0,s],[s,s]]), s in S)},
        [i,x], x)], i=0..maxp)]:
    E:=GetMinimalContractions(S):
    return EncapsulateDataGen(S,get,E):
end:

#Generate period data on all self double Wythoff sums of
#subtraction sets whose elements are <=cap
#Output the data to file out
#If append is true and file out exists already,
#start where you left off
#maxp=max other pile count to consider
GenerateSelfWythPeriodData:=proc(cap,maxp,out,append:=true)
local ctr,i,fd,S,
        appval,line,spline,spl,chrd:
    appval:=false:
    if out<>terminal then
        if append and Exists(out) then
            fd:=fopen(out,READ):
            try:
                Position(fd,infinity):
                while true do
                    Position(fd,Position(fd)-2):
                    chrd:=ReadCharacter(fd):
                    if chrd="S" then
                        Position(fd,Position(fd)+1):
                        line:=ReadLine(fd):
                        spline:=Split(line):
                        appval:=parse(sprintf("%a",spline)):
                        appval:=[seq(parse(appval[i]),
                            i=1..nops(appval))]:
                        break:
                    fi:
                od:
            finally:
                fclose(fd):
            end try:
            fd:=fopen(out,APPEND):
        else
            fd:=fopen(out,WRITE):
        fi:
    else
        fd:=out:
    fi:
    
    try:
        ctr:=InitCounter(cap):
        if appval<>false then
            for i in appval do
                ctr[i]:=1:
            od:
            
            i:=1:
            ctr:=IncCounter(ctr):
        fi:
        while ctr<>false do
            S:=CounterToSet(ctr):
            
            fprintf(fd,"%s\n\n", GetFormattedDataGen(S,maxp)):
            
            ctr:=IncCounter(ctr):
        od:
    finally:
        if out<>terminal then
            fclose(fd):
        fi:
    end try:
end: