tsscpp_interpolation.sage

import os


############################################################
####        Reading in K Matrices
############################################################


def retrieve_mat(name, n):
    if name == 'M' or name == 'M_inv' or name == 'K':
        mat_file = open('tsscpp_data/tsscpp_{0}/{0}_{1}'.format(name, n), 'r')
        mat_data = sage_eval(mat_file.readline())
        mat = Matrix(QQ, mat_data)
        mat_file.close()
        return mat
    else:
        print("Error: first argument must be 'M', 'M_inv', or 'K'")
        return

def make_K_mats(N):
    K_mats = []
    for n in [2 .. N]:
        K_mats.append(retrieve_mat('K', n))
    return K_mats


############################################################
####        Interpolating the Kernel K_n
############################################################


def make_denom(n, x, y):
    x_max = ceil(x/2)
    y_max = ceil(y/2)
    return prod(4*n^2 - (2*m - 1)^2 for m in [1 .. x_max]) * prod(4*n^2 - (2*m - 1)^2 for m in [1 .. y_max])

# Creates the denominator of the conjectured formula for K_n, as a symbolic expression type

def make_Q_interpolated(x, y, alpha, beta, K_mats):
    i = 2*(x-1) + alpha
    j = 2*(y-1) + beta

    # K_n(x,y)_(alpha, beta) corresponds to entry (i,j) of K_n

    m = max(i,j)
    n_min = ceil(1 + m/4)
    N = len(K_mats) + 1

    # n_min is the lowest index n for which K_n contains entry (i,j)

    entries = []
    for n in [n_min .. N]:
        entries.append(K_mats[n-2][i-1, j-1])

    # creates a sequence of entries K_n(x,y)_(alpha, beta) in n

    interp_pts = [(n, 2^(2*(x+y)) * entries[n-n_min] * make_denom(n, x, y)) for n in [n_min .. N]]
    R = QQ['x']
    Q_expr = R.lagrange_polynomial(interp_pts)
    return Q_expr

# Creates the polynomial Q_(alpha, beta, x, y) of the conjectured form of K_n, as a polynomial type in x

def make_rat_fn(x, y, alpha, beta, K_mats):
    i = 2*(x-1) + alpha
    j = 2*(y-1) + beta

    var('n')
    Q = make_Q_interpolated(x, y, alpha, beta, K_mats)
    denom(n) = make_denom(n, x, y)
    rat_fn(n) = 2^(-2*(x + y)) * Q(n) / denom(n)

    return rat_fn.simplify_rational()

def make_K_interpolated(K_mats):
    N = len(K_mats) + 1
    x_max = floor((2*N - 9) / 5)

    # x_max is the greatest index x for which K_n(x,x) can be interpolated, using the data from K_mats

    K_interpolated = Matrix(SR, 2*x_max)
    for x in [1 .. x_max]:
        for y in [x .. x_max]:
            if x == y:
                alpha = 1
                beta = 2
                i = 2*(x-1) + alpha
                j = 2*(y-1) + beta
                K_interpolated[i-1, j-1] = make_rat_fn(x, y, alpha, beta, K_mats)
            else:
                for alpha in [1,2]:
                    for beta in [1,2]:
                        i = 2*(x-1) + alpha
                        j = 2*(y-1) + beta
                        K_interpolated[i-1, j-1] = make_rat_fn(x, y, alpha, beta, K_mats)
        if x % 5 == 0:
            print('x={0} completed'.format(x))

    return K_interpolated - K_interpolated.transpose()


############################################################
####        Limits of Rational Functions
############################################################


def leading_degree(poly):
    num_terms = len(poly.coefficients())
    return(poly.coefficients()[num_terms-1][1])

def leading_coeff(poly):
    num_terms = len(poly.coefficients())
    return(poly.coefficients()[num_terms-1][0])

def rat_limit(rat):
    num = rat.numerator()
    den = rat.denominator()
    if leading_degree(num) < leading_degree(den):
        return 0
    elif leading_degree(num) == leading_degree(den):
        return leading_coeff(num) / leading_coeff(den)
    else:
        pm = sgn(leading_coeff(num) / leading_coeff(den))
        return(pm * oo)


############################################################
####        FIXME
############################################################


def write_M_infinity(N, path):
    """
    write_M_infinity(N, path)

    Takes a positive integer N (large in practice) and creates the matrix
    M_N, which is written to a file

    Arguments
    N: size of the matrix M_N to be calculated. Note that if N is odd,
       then N is rounded to the next even number so that the indexing
       of M_N starts at delta = 0
    path: directory path in which to write
    """

    N_next_even = 2 * ceil(N/2)
    M_infinity = make_Mn(N_next_even)

    file_name = 'M_infinity_trunc_{0}'.format(N_next_even)
    text_file = open(os.path.join(path, file_name), 'w')
    text_file.write(str(mat_to_list(M_infinity)))
    text_file.close()

    print('M_{0} successfully written'.format(N_next_even))
    return()

def read_M_infinity(N, path):
    """
    read_M_infinity(N, path)

    Reads in the matrix M_N (created by read_M_infinity(N, path))
    as a list of lists

    Arguments
    N: size of the matrix M_N to read in. Specifically, this function
       reads the file created by make_M_infinity(N, path)
    path: directory path containing the file to read
    """

    N_next_even = 2 * ceil(N/2)
    file_name = 'M_infinity_trunc_{0}'.format(N_next_even)
    text_file = open(os.path.join(path, file_name), 'r')
    return(sage_eval(text_file.readline()))

def make_M_infinity_submat(n, L):
    """
    make_M_infinity_submat(n, L)

    Creates the matrix M_n from a list L of the entries of
    M_infinity

    Arguments:
    n: order of the matrix M_n
    L: list of lists containing the data of M_infinity (part of
       this data, more precisely)
    """
    delta = n % 2
    return(Matrix(ZZ, L)[delta:n, delta:n])