Steenbrink/resolucion.ipynb

In [1]:

# g = Graph(17)
# a = graph_editor(g)
# a.show()

In [2]:

# A =g.adjacency_matrix()
# A[0, 0] = -13
# for i in (1, 3, 5, 6, 7):
#     A[i, i] = -2
#     A[i + 8, i + 8] = -2
# for i in (2,):
#     A[i, i] = -4
#     A[i + 8, i + 8] = -4
# for i in (4, 8):
#     A[i, i] = -1
#     A[i + 8, i + 8] = -1
# # A = A.change_ring(QQ)

In [3]:

A = Matrix(QQ, 17, [-13, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -4, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -4, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1])
Vf = vector(4 * [0] + [1] + 7 * [0] + [1] + 4 * [0])
Vg = vector(8 * [0] + [1] + 7 * [0] + [1])
Vd = vector([13] + 7 * [0] + [1] + 7 * [0] + [1])
R.<x, y, z> = QQ[]
S.<s> = QQ[]
T.<t> = QQ[]

In [4]:

def trinomiales(A, Vf, Vg, Vd, total=True):
    n = A.nrows()
    A1 = -A^-1
    Nf, Ng, Nd = [A1 * N0 for N0 in (Vf, Vg, Vd)]
    diag = vector(A.diagonal())
    unos = vector(n * [1])
    if total:
        VV = Vf + Vg
    else:
        VV = Vf.parent()(0)
    euler = (-A + 2 * identity_matrix(n)) * unos + diag - VV - Vd
    nus = A1 * (diag + 2 * unos) + unos
    expos = {}
    for j in range(n):
        a, b = (Nf[j], Ng[j])
        eq = 3 * Nf[j] == 2 * Ng[j]
        nnd = min(3 * Nf[j], 2 * Ng[j]) == Nd[j]
        if nnd:
            expos[j] = ((a, b), (a, b), nus[j], eq, nnd)
        else:
            expos[j] = ((a, b), (Nf[j], Nd[j] - Ng[j]), nus[j], eq, nnd)
    cadenas = ('d',)
    if total:
        expos['f'] = ((1, 0), (0, 0), 1, False, True)
        expos['g'] = ((0, 1), (0, 0), 1, False, True)
        cadenas = ('f', 'g') + cadenas
    expos['d'] = ((0, 0), (0, 1), 1, True, False)
    pols = {}
    for j in range(n):
        if expos[j][-1]:
            a, b = expos[j][0]
            pols[j] = (z^3 + x^a * z + x^b, (expos[j][2], 1, 1), euler[j])
        else:
            a, b = expos[j][1]
            pols[j] = (z^3 + x^a * z^2 + x^b, (expos[j][2], 1, 1), euler[j])
    interdiv = [(i, j) for i, j in Combinations(n, 2) if A[i, j] != 0]
    interstrict = [(i, j) for i in range(n) for j in cadenas if eval('V' + j)[i] != 0]
    for i,j in interdiv + interstrict:
        if j in ZZ:
            ep = 1
        else:
            ep = eval('V' + j)[i]
        eqi, nndi = expos[i][-2:]
        eqj, nndj = expos[j][-2:]
        nd = eqi and eqj and not (nndi and nndj)
        if not nd:
            ai, bi = expos[i][0]
            aj, bj = expos[j][0]
            pols[i, j] = (z^3 + x^ai * y^aj * z + x^bi * y^bj,
                          (expos[i][2], expos[j][2], 1), ep)
        else:
            ai, bi = expos[i][1]
            aj, bj = expos[j][1]
            pols[i, j] = (z^3 + x^ai * y^aj * z^2 + x^bi * y^bj, 
                          (expos[i][2], expos[j][2], 1), ep)
    return n, Nf, Ng, Nd, euler, nus, expos, pols

def estrato(pol, pesos, euler):
    Z1 = ZetaFunctions(pol)
    zetatop = euler * Z1.topological_zeta(local=True, weights=pesos)
    zetamon = Z1.monodromy_zeta()^euler
    return zetatop, zetamon

def zetas(n, pols):
    zetatop=0
    zetamon=1
    for k in pols:
        zt, zm = estrato(*pols[k])
        zetatop += zt
        zetamon *= zm
    return zetatop, zetamon / (1 - t^3)^n

In [5]:

n, Nf, Ng, Nd, euler, nus, expos, pols = trinomiales(A, Vf, Vg, Vd)
zt0, zm0 = zetas(n, pols)

In [6]:

zt0.factor()

Out[6]:

(-7/39) * (s + 19/39)^-1 * (s + 1/2)^-1 * (s + 1)^-1 * (s^2 - 1207/546*s - 19/14)

In [7]:

for p, j in zm0.factor():
    print(p.is_cyclotomic(certificate=True), j)

2 4
1 5
6 4
3 16
26 4
13 4
78 4
39 15

In [8]:

zt0.partial_fraction_decomposition()

Out[8]:

(0, [76/65/(s + 19/39), -2/39/(s + 1/2), -253/195/(s + 1)])

In [9]:

n1, Nf1, Ng1, Nd1, euler1, nus1, expos1, pols1 = trinomiales(A, Vf, Vg, Vd, total=False)
zt1, zm1 = zetas(n1, pols1)

In [10]:

zt0 - zt1, zm0 / zm1

Out[10]:

(0, 1)

In [11]:

B = A.delete_columns([3, 4, 11, 12]).delete_rows([3, 4, 11, 12])
B[2, 2] = -2
B[8, 8] = -2
Wf = vector(Vf.list()[:2] + [2] + Vf.list()[5:10] + [2] + Vf.list()[13:])
Wg = vector(Vg.list()[:3] + Vg.list()[5:11] + Vg.list()[13:])
Wd = vector(Vd.list()[:3] + Vd.list()[5:11] + Vd.list()[13:])

In [12]:

n2, Nf2, Ng2, Nd2, euler2, nus2, expos2, pols2 = trinomiales(B, Wf, Wg, Wd, total=False)
zt2, zm2 = zetas(n2, pols2)

In [13]:

zt0 - zt2, zm1 / zm2

Out[13]:

(0, 1)

In [ ]: