# torhom.mpl # # Copyright (C) 2000-2016 Matthias Franz # # Torhom is a Maple package that provides functions to study # the Borel-Moore homology of real and complex toric varieties. # See http://www-fourier.ujf-grenoble.fr/~franz/maple/torhom.html # for more information. # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation; either version 2 of the License, or # (at your option) any later version. # macros for functions from the convex package macro(basiclists = CONVEX['_idmat']); macro(matmatmul = CONVEX['_matmatmul']); macro(det = CONVEX['_det']); macro(hermite3 = CONVEX['_hermite3']); macro(transpose = CONVEX['_transpose']); macro(Dotprod = `convex/dotprod`); macro(MODZ_0 = MODZ['_zero']); macro(ormap = `std/exists`); macro(parsemod = CONVEX['_parsemod']); macro(ERRORnargs = ERROR("too many arguments")); torhom := module() # description "functions to compute the Borel-Moore homology of toric varieties"; option package, # copyright = "2000-2016 Matthias Franz", load = init; export E2hom, homreal, printarray, total; local prmul, B, sortind, extmul, make_d1, is_hom_manifold, prmul_real, make_d_real, Z2map, init; init := proc() printf("Torhom version 1.3.1, Copyright (C) 2000-2016 Matthias Franz\n\ This package is distributed under the GNU General Public License\n"); end; # # COMPLEX # prmul := proc(A, B) local C, c, i, j; C := table(); for i to nops(A) do C[i] := table(sparse); for j to nops(B) do c := Dotprod(A[i], B[j]); if c <> 0 then C[i][j] := c fi od od; eval(C) end; # # exterior multiplication # B := proc(p, q) option remember; # making bases of exterior powers if q = 0 then [[]] elif q = p+1 then [] else [op(B(p-1, q)), op(map((x, y) -> [op(x), y], B(p-1, q-1), p))] fi end; # multiplication of exterior powers sortind := proc(ind2, ind1) option remember; local i1; for i1 to nops(ind1) while ind1[i1] < ind2 do od; if i1 <= nops(ind1) and ind1[i1] = ind2 then RETURN(FAIL, 0) else RETURN([op(1..i1-1, ind1), ind2, op(i1..-1, ind1)], -(-1)^i1) fi; end; extmul := proc(P1, P2) local P, ind1, ind2, ind, s; P := table(sparse); for ind1 in [indices(P1)] do for ind2 in [indices(P2)] do ind, s := sortind(op(ind1), ind2); if s <> 0 then #lprint(ind1, ind2, ind, s); P[op(ind)] := P[op(ind)]+s*P1[op(ind1)]*P2[op(ind2)] fi od od; eval(P) end; # # composing differential # make_d1 := proc(pr_pq, Ci, Cj, Bi, Bj, o) local i, j, Ci_inv, Cj_inv, Bi_inv, Bj_inv, A, d1, c, i2, j2, i_off, j_off; Ci_inv := table([seq(Ci[i] = i, i = 1..nops(Ci))]); Cj_inv := table([seq(Cj[i] = i, i = 1..nops(Cj))]); Bi_inv := table([seq(Bi[i] = i, i = 1..nops(Bi))]); Bj_inv := table([seq(Bj[i] = i, i = 1..nops(Bj))]); d1 := array(1..nops(Ci)*nops(Bi), 1..nops(Cj)*nops(Bj), sparse); if Bi = [[]] then for i in Ci do for j in Cj do c := o[i, j]; if c <> 0 then d1[Ci_inv[i], Cj_inv[j]] := c fi od od else for i in Ci do for j in Cj do c := o[i, j]; if c = 0 then next fi; i_off := (Ci_inv[i]-1)*nops(Bi); j_off := (Cj_inv[j]-1)*nops(Bj); A := eval(pr_pq[i, j]); for i2 in [indices(A)] do for j2 in [indices(A[op(i2)])] do d1[i_off+Bi_inv[i2], j_off+Bj_inv[j2]] := c*A[op(i2)][op(j2)] od od od od fi; eval(d1) end; is_hom_manifold := proc(F, n, indC, o) # tries to guess whether the support of F is a homology manifold # "true" means "yes", "false" means "maybe not" local MC, i; MC := `if`(type(F, CONE), [F], convex['maximal'](F, 'list')); i := nops([indices(o[n])]); if convex['codim'](MC[-1]) <> 0 then false elif nops(MC) = 1 then userinfo(4, torhom, '`fan is affine`'); true elif (nops(MC) = 2 and nops(indC[n-1]) = i-1) then userinfo(4, torhom, '`fan consists of 2 full-dimensional cones with common facet`'); true elif 2*nops(indC[n-1]) = i then userinfo(4, torhom, '`fan is complete`'); true else false fi end; # # E2hom # E2hom := proc(F::{FAN0, CONE0}, F0) local na, COEFF, n, o, pr, bases, indC, reg_dim, support_h_mf, extmatmul, p, q, ind, inds, E2, A, stop_rs, orient; na, COEFF := parsemod(args); n := convex['ambientdim'](F); if na = 1 then stop_rs := {} elif na > 2 then ERRORnargs elif type(F0, FAN0) and convex['ambientdim'](F0) = n then stop_rs := convert(map(x -> convert(convex['rays'](x), set), convex['maximal'](F0, 'list')), set) else ERROR("second argument must be a subfan of the first") fi; # userinfo(2, torhom, `if`(COEFF['_char'] = 2, '`computing lattice bases and projections ...`', '`computing lattice bases, orientations, and projections ...`')); # indC := Array(0..n, [{}$n+1]); o := Array(0..n, [seq(table(sparse), p = 0..n)]); pr := Array(0..n, 0..n); bases := proc(d::nonnegint, rs::set) option remember; local A, U, V, i; if d = nops(rs[1]) # full-dimensional cone? # it is not necessary to check whether it is regular then RETURN([], [], basiclists(d)) fi; A := hermite3(transpose(convert(rs, list)), false, U, V); if d <= reg_dim then if d = nops(rs) then for i to d do if abs(A[i, i]) <> 1 then reg_dim := d-1; break fi od else reg_dim := d-1 fi fi; [V[d+1..-1], U[d+1..-1], V[1..d]] end; bases(0, {}) := [basiclists(n)$2, []]; reg_dim := max(n-1, 0); if not stop_rs subset convex['traverse2'](F, proc(f, rs) local d; if ormap(proc(x) rs subset x end, stop_rs) then false else d := CFACE['dim'](f); indC[d] := indC[d] union {rs}; true fi end, `if`(COEFF['_char'] = 2, proc(f2, f1, rs2, rs1) # f2 is a facet of f1 local d; d := CFACE['dim'](f1); o[d][rs1, rs2] := 1; pr[d-1, 1][rs1, rs2] := prmul(bases(d, rs1)[2], bases(d-1, rs2)[1]) end, proc(f2, f1, rs2, rs1) # f2 is a facet of f1 local d, B1, B2; d := CFACE['dim'](f1); B1 := bases(d, rs1); B2 := bases(d-1, rs2); o[d][rs1, rs2] := sign(det(matmatmul(B1[3], [op(B2[3]), (rs1 minus rs2)[1]]))); pr[d-1, 1][rs1, rs2] := prmul(B1[2], B2[1]) end)) #if not stop_rs subset (result of traverse2) then ERROR("second argument must be a subfan of the first") fi; bases := NULL; # to free memory # preparing shortcuts for computation (following Brion) support_h_mf := stop_rs = {} and is_hom_manifold(F, n, indC, o); if support_h_mf then userinfo(2, torhom, `support is full-dim homology manifold with boundary`); if reg_dim = n-1 and convex['isregular'](F) then reg_dim := n fi; userinfo(3, torhom, `fan is regular in dimension`, reg_dim) else reg_dim := 0 fi; if reg_dim <> n then # userinfo(2, torhom, '`computing exterior powers of projections ...`'); # # multiplication of maps between exterior powers extmatmul := proc(P1, P2, n, k) local P, ind; for ind in B(n, k) do P[op(ind)] := extmul(P1[ind[1]], P2[op(2..-1, ind)]) od; eval(P) end; # computing powers recursively for p from 0 to n-1 do if pr[p, 1] = 0 then next fi; # already defined in transversal for ind in [indices(pr[p, 1])] do inds := op(ind); pr[p, 0][inds] := 1; for q from 2 to n-p-1 do pr[p, q][inds] := extmatmul(pr[p, 1][inds], pr[p, q-1][inds], n-p-1, q); od od od; fi; # userinfo(2, torhom, '`composing differentials and computing homology ...`'); # # computing ranks and torsions E2 := Array(0..n, 0..n, `if`(COEFF = 'MODZ', () -> MODZ_0, () -> 0)); A := Array(0..n); for q from 0 to n do userinfo(4, torhom, `q =`, q); for p from q to n do A[p] := nops(indC[n-p])*nops(B(p, q)) od; # Brion's reduction for p from n to q+n-reg_dim+1 by -1 do A[p-1] := A[p-1]-A[p]; A[p] := 0; od; COEFF['_homology'](q, n, A, proc(p) convert( make_d1(eval(pr[n-p, q]), indC[n-p+1], indC[n-p], B(p-1, q), B(p, q), o[n-p+1]), listlist) end); for p from q to n do E2[p, q] := A[p] od od; E2 end; # # REAL # prmul_real := proc(A, B) # returns a listlist whose rows are the images of the basis vectors # of the quotient by the small cone in the quotient by the big cone # # WARNING: different from complex case! local v, w; [seq([seq(Dotprod(v, w), v = A)], w = B)] end; Z2map := proc(A::listlist) local v, l; l := [[0$nops(A[1])]]; for v in A do l := [op(l), op(map(`+`, l, v))] od; map(modp, l, 2) end; # # composing differential # make_d_real := proc(pr, Ci, Cj, di, dj, o) local ci, cj, d, m, n, i, j, oij, l, ki, kj, v, k; ci := 2^di; cj := 2^dj; d := array(1..nops(Ci)*ci, 1..nops(Cj)*cj, sparse); m := 0; n := 0; for i in Ci do for j in Cj do oij := o[i, j]; if oij = 0 then n := n+cj; next fi; l := Z2map(pr[i, j]); for kj from 1 to cj do v := l[kj]; ki := `+`(seq(v[k]*2^k, k = 1..di)); d[m+ki/2+1, n+kj] := oij; od; n := n+cj; od; m := m+ci; n := 0; od; eval(d) end; # # homreal # homreal := proc(F::{FAN0, CONE0}, F0) local na, COEFF, n, o, pr, bases, indC, # indC[d] is the set of raysets of cones of dim d p, A, stop_rs, orient; na, COEFF := parsemod(args); n := convex['ambientdim'](F); if na = 1 then stop_rs := {} elif na > 2 then ERRORnargs elif type(F0, FAN0) and convex['ambientdim'](F0) = n then stop_rs := convert(map(x -> convert(convex['rays'](x), set), convex['maximal'](F0, 'list')), set) else ERROR("second argument must be a subfan of the first") fi; # userinfo(2, torhom, `if`(COEFF['_char'] = 2, '`computing lattice bases and projections ...`', '`computing lattice bases, orientations, and projections ...`')); # indC := Array(0..n, [{}$n+1]); o := Array(0..n, [seq(table(sparse), p = 0..n)]); pr := Array(0..n, 0..n); bases := proc(d::nonnegint, rs::set) option remember; local U, V; if d = nops(rs[1]) # full-dimensional cone? then RETURN([], [], basiclists(d)) fi; hermite3(transpose(convert(rs, list)), false, U, V); [V[d+1..-1], U[d+1..-1], V[1..d]] end; bases(0, {}) := [basiclists(n)$2, []]; if not stop_rs subset convex['traverse2'](F, proc(f, rs) local d; if ormap(proc(x) rs subset x end, stop_rs) then false else d := CFACE['dim'](f); indC[d] := indC[d] union {rs}; true fi end, `if`(COEFF['_char'] = 2, proc(f2, f1, rs2, rs1) local d; d := CFACE['dim'](f1); o[d][rs1, rs2] := 1; pr[d-1, 1][rs1, rs2] := prmul_real(bases(d, rs1)[2], bases(d-1, rs2)[1]) end, proc(f2, f1, rs2, rs1) local d, B1, B2; d := CFACE['dim'](f1); B1 := bases(d, rs1); B2 := bases(d-1, rs2); o[d][rs1, rs2] := sign(det(matmatmul(B1[3], [op(B2[3]), (rs1 minus rs2)[1]]))); pr[d-1, 1][rs1, rs2] := prmul_real(B1[2], B2[1]) end)) #if not stop_rs subset (result of traverse2) then ERROR("second argument must be a subfan of the first") fi; bases := NULL; # to free memory # userinfo(2, torhom, '`composing differentials and computing homology ...`'); # # computing ranks and torsions A := Array(0..n, [seq(nops(indC[n-p])*2^p, p = 0..n)]); COEFF['_homology'](0, n, A, proc(p) local d; d := make_d_real(eval(pr[n-p, 1]), indC[n-p+1], indC[n-p], p-1, p, o[n-p+1]); # bugfix for rational case: CONVEX[_gausselim] doesn't accept array if COEFF = 'MODQ' then convert(d, listlist) elif COEFF <> 'MODZ' # bugfix for mod p case: LinearAlgebra:-Modular:-Mod doesn't # accept array together with integer[] then convert(d, Matrix) else d fi end); A end; # # MISCELLANEOUS # printarray := proc(E0::{array, Array}) local E, B, i, j, n; if nargs >1 then ERRORnargs fi; E := convert(E0, array); B := [op(2, eval(E))]; if nops(B) = 1 then array([seq(E[i], i = op(B))]) elif nops(B) = 2 then n := op([2, 2], B)+op([2, 1], B); array([seq([seq(E[j, n-i], j = B[1])], i = B[2])]) else ERROR("illegal array") fi # n := op([2, 1, 2], eval(E)); # array(1..n+1, 1..n+1, [seq([seq(E[j, n-i], j = 0..n)], i = 0..n)]) end; total := proc(E::{Array(nonnegint), Array(MODZ)}) local s, n, k, j; if nargs >1 then ERRORnargs fi; s := `if`(type(E, 'Array'(MODZ)), MODZ['directsum'], `+`); n := op([2, 1, 2], E); Array(0..2*n, [seq(s(seq(E[j, k-j], j = max(0, k-n)..min(k, n))), k = 0..2*n)]) end; #bettisum2 := proc(A::Array(MODZ)) #local b, i; # if nargs >1 then ERRORnargs fi; # # b := 0; # for i from 0 to op([2, 2], A) do # b := b+rank(A[i])+2*nops(select(x -> evalb(x mod 2 = 0), torsion(A[i]))) # od; # b #end; #$include "cellhom.mt" end module; # # saving the module to a repository # $define torhomrep "torhom.lib" march('create', torhomrep, 16); savelibname := currentdir(); savelib('torhom'); map(march, ['gc', 'reindex', 'pack'], torhomrep); march('setattribute', torhomrep, mode="READONLY");