Code

#!/usr/bin/python
#-*- coding: utf-8 -*-
import numpy as np
import operator
import matplotlib.pyplot as plt
import pandas as panda
import time
import math

'''=================  RULES  =================== '''

def rules(a,b):                                                                  
    i=1
    tmp1 = [[a-(1+i),b] for i in range(a)]
    tmp2 = [[a,b-(1+i)] for i in range(b)]
    
    if a<b:
        tmp3 = [[a-(1+i),b-(1+i)] for i in range(a)]
    else:
        tmp3 = [[a-(1+i),b-(1+i)] for i in range(b)]
        
    return tmp1 + tmp2 + tmp3


'''================  CALCULS  ================== '''
def mex(arr):
    arr = np.unique(arr)
    for i in range(arr.size):
        if arr[i]!=i:
            return i
    return arr.size
           
def losing(k, symmetrical=False):#return the losing positions
    x=time.time()
    if k <= 0:
        return None
    sortie = [[0,0]]
    a=0
    b=1
    while a < k:
        while b < k:
            losing = True
            for option in rules(a,b):
                if option in sortie:
                    losing = False
                    break
            if losing:
                sortie.append([a,b])
                if symmetrical:
                    sortie.append([b,a])
            b=b+1
        a=a+1
        if symmetrical:
            b=a
        else:
            b=0
    print 'time : '+str(time.time()-x)
    return sortie

def Grundy(k, symmetrical=False):#returning nimbers for positions from 0 to k
    x=time.time()
    if k <= 0:
        return None
    a=0
    dico = [[0 for i in range(k)] for j in range(k)]
    b=1
    while a < k:
        while b < k:
            tmp = [dico[option[0]][option[1]] for option in rules(a,b)]
            dico[a][b] = mex(tmp)
            if symmetrical:
                dico[b][a]=dico[a][b]
            b=b+1
        a=a+1
        if symmetrical:
            b=a
        else:
            b=0
    print 'elapsed time : '+str(time.time() - x)
    return dico

def GrundyValeur(x,y, symmetrical=False):#Return the nimber of a given position
    elapsedTime = time.time()
    if x < 0 or y < 0:
        return None
    if symmetrical:
        if x>y:
            dico = [[-1 for a in range(x+1)] for b in range(x+1)]
        else:
            dico = [[-1 for a in range(y+1)] for b in range(y+1)]
    else:
        dico = [[-1 for a in range(x+1)] for b in range(y+1)]
    dico[0][0] = 0
    if dico[y][x] != -1:
        return dico[y][x]
    stack = [[y,x]]
    stack.extend(rules(y,x))
#Game initialisation
    while dico[y][x]==-1:
        IN = True
        if dico[stack[-1][0]][stack[-1][1]] != -1:
            stack.pop()
        else:
            tmp = rules(stack[-1][0],stack[-1][1])
            for option in tmp:
                if dico[option[0]][option[1]] == -1:
                    IN = False
                    break
            if IN: # if all the options are in dico
                Nmin = mex([dico[option[0]][option[1]] for option in tmp])
                dico[stack[-1][0]][stack[-1][1]] = Nmin
                if symmetrical:
                    dico[stack[-1][1]][stack[-1][0]] = Nmin
                stack.pop()
            else:
                stack.extend(tmp)
    print 'elapsed time : '+str(time.time() - elapsedTime)
    return dico[y][x]

def GrundyGameSum(k):#Find nimbers if the rule is compatible with game sum
    t = time.time()
    tableau = [0 for i in range(k)]
    for i in range(k):
        tableau[i] = mex([tableau[option[1]] for option in rules(0,i)])
    sortie = [[tableau[i] ^ tableau[j] for j in range(k)] for i in range(k)]
    print str(time.time() - t)
    return sortie

def ultraFastWythoffLosing(k):#find easily all the losing positions for a wythoff game
    t = time.time()
    gold = (1+math.sqrt(5))/2
    tab = [[0,0]]
    n=1
    while tab[-1][0] < k:
        tab.extend([[int(n*gold), int(n*gold*gold)],[int(n*gold*gold), int(n*gold)]])
        n+=1
    del tab[-1]
    del tab[-1]
    print str(time.time() - t)
    return tab


'''================  DISPLAY  ================== '''

def displayGrundy(tab):#display all the grundy values of an array
    fig, ax = plt.subplots()
    ax.imshow(tab, cmap=plt.cm.magma, interpolation='nearest')
    ax.set_title('tableau des valeurs de Grundy')
    ax.invert_yaxis()
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    ax.set_xticks(np.arange(0, len(tab), 1))
    ax.set_yticks(np.arange(0, len(tab), 1))
    #plt.grid()
    plt.show()  

def displayLosing(arr, taille, Grille = False):#display all the losing values of an array
    tableau = np.zeros((taille,taille), dtype = np.int)
    for losing in arr:
        tableau[losing[0]][losing[1]] = 1

    fig, ax = plt.subplots()
    ax.imshow(tableau, cmap=plt.cm.Greys, interpolation='nearest')
    ax.set_title('losing positions')
    ax.invert_yaxis()
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)
    ax.yaxis.set_ticks_position('left')
    ax.xaxis.set_ticks_position('bottom')
    ax.set_xticks(np.arange(0, taille, 1))
    ax.set_yticks(np.arange(0, taille, 1))
    if Grille:
        plt.grid()
    plt.show()  

def sortGrundy(tab):#sort position by nimbers
    maximum = max(max(tab))
    sortie = [['nimber : ' +str(i)] for i in range(maximum+1)]
    for i in range(len(tab)):
        for j in range(len(tab)):
            sortie[tab[i][j]].append([i,j])
    return sortie
    
def saveCSV(tab):#save an array in a CSV file
    df = panda.DataFrame(tab)
    df.to_csv("save.csv")

def cleanDisplay(tab):#display an array of nimbers legibly
    for i in range(len(tab[0])):
        for j in range(len(tab[0])):
            print "(",i,",",j,") =",tab[i][j]

    
'''================  MAIN  ================== '''
    
if __name__=='__main__':


    print " Menu :"
    print " 1 - Search for losing positions"
    print " 2 - Search for nimbers"
    print " 3 - Search a nimber for a given position"
    choice1 = input("Your choice : ")
    if choice1 != 3:
        print " Menu 2 :"
        print " 1 - Display results (list)"
        print " 2 - Save results"
        print " 3 - Display results (table)"
        choice2 = input("Your 2nd choice : ")
    else:
        choice2 = 1

    boolean = False
    retour = None
    if choice1 == 1:
        print " Is your game symmetrical"
        print " 1 - Yes"
        print " 2 - No"
        symetrie = input("Your choice : ")
        if symetrie == 1 :
            boolean = True
        choice3 = input("Upper bound : ")
        retour = losing(choice3, boolean)
    elif choice1 == 2:
        print "Menu 3 :"
        print " Is your game symmetrical"
        print " 1 - Yes"
        print " 2 - No"
        symetrie = input("Your choice : ")
        if symetrie == 1:
            boolean = True
        choice3 = input("Upper bound : ")
        retour = Grundy(choice3,boolean)
    elif choice1 == 3:
        print "Is your game symmetrical?"
        print " 0 - no"
        print " 1 - yes"
        symetrie = input("Your choice : ")
        if symetrie == 1 :
            boolean = True
        choice3x = input("Position to compute (first heap) : ")
        choice3y = input("Position to compute (second heap) : ")
        retour = GrundyValeur(choice3x, choice3y, boolean)
        
    if choice2 == 1:
        if choice1==2:
            cleanDisplay(retour)
        else:
            print retour
    elif choice2 == 2:
        if choice1==2:
            saveCSV(sortGrundy(retour))
        else:
            saveCSV(retour)
    elif choice2 == 3:
        if choice1 == 1:
            displayLosing(retour, choice3)
        if choice1 == 2:
            displayGrundy(retour)

Example of rules

Rules allows us to define what you can play for a given position in a game. It takes two integers as an input. This is the size of each heap. And it return a python list. In this list you will find all the allowed move in the current context. We recommend to use list comprehension, because it's faster.
For example, in the code above, the rules are those of Wythoff. It allows to pick how many token you desire in one heap, or the same number of token in both heap.

#!/usr/bin/python

import numpy as np
import operator
import matplotlib.pyplot as plt 


def rules(a,b):                                                                  
    i=1
    tmp1 = [[a-(1+i),b] for i in range(a)]
    tmp2 = [[a,b-(1+i)] for i in range(b)]
    
    if a<b:
        tmp3 = [[a-(1+i),b-(1+i)] for i in range(a)]
    else:
        tmp3 = [[a-(1+i),b-(1+i)] for i in range(b)]
        
    return tmp1 + tmp2 + tmp3