From b4efbf5b9e0b6a2499870b04bd542bda70418384 Mon Sep 17 00:00:00 2001 From: Vincent Delecroix Date: Fri, 4 Dec 2020 09:35:45 +0100 Subject: [PATCH 1/5] Reorganize and expand generalized multiple zeta values --- .../__init__.py | 4 + .../generalized_multiple_zeta_values.py | 1322 +++++++++++++++++ .../gmzv_generic_reduction.py | 812 ++++++++++ .../gmzv_three_variables.py | 217 +++ .../gmzv_two_variables.py | 59 + .../options.py | 4 + ...additive_multivariate_generating_series.py | 2 +- .../misc/divergent_multiple_zeta_values.py | 128 ++ surface_dynamics/misc/factored_denominator.py | 17 + .../misc/factored_denominator.py.orig | 1121 ++++++++++++++ .../misc/generalized_multiple_zeta_values.py | 985 ------------ surface_dynamics/misc/groups.py | 26 + surface_dynamics/misc/linalg.py | 144 +- surface_dynamics/misc/linalg.py.orig | 458 ++++++ .../test_generalized_multiple_zeta_values.py | 3 + 15 files changed, 4315 insertions(+), 987 deletions(-) create mode 100644 surface_dynamics/generalized_multiple_zeta_values/__init__.py create mode 100644 surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py create mode 100644 surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py create mode 100644 surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py create mode 100644 surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py create mode 100644 surface_dynamics/generalized_multiple_zeta_values/options.py create mode 100644 surface_dynamics/misc/divergent_multiple_zeta_values.py create mode 100644 surface_dynamics/misc/factored_denominator.py.orig delete mode 100644 surface_dynamics/misc/generalized_multiple_zeta_values.py create mode 100644 surface_dynamics/misc/groups.py create mode 100644 surface_dynamics/misc/linalg.py.orig diff --git a/surface_dynamics/generalized_multiple_zeta_values/__init__.py b/surface_dynamics/generalized_multiple_zeta_values/__init__.py new file mode 100644 index 00000000..f111f2b0 --- /dev/null +++ b/surface_dynamics/generalized_multiple_zeta_values/__init__.py @@ -0,0 +1,4 @@ +from .generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, convergent_multizeta, generalized_multiple_zeta_functions +from .gmzv_two_variables import is_Z2_convergent, Z2 +from .gmzv_three_variables import is_Z3_convergent, Z3 +from .gmzv_generic_reduction import gmzv_solvers diff --git a/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py b/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py new file mode 100644 index 00000000..d6fe26f4 --- /dev/null +++ b/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py @@ -0,0 +1,1322 @@ +r""" +Generalized multiple zeta values + +Most functions in this module have a pair ``(n, den_tuple)`` as +arguments. ``n`` stands for the dimension (number of variables) +and ``den_tuple`` is a list of pairs ``(v, p)`` where ``v`` +is a vector of length ``n`` with entries `\{0, 1\}` corresponding +to a linear form and ``p`` is a positive integer corresponding +to a power. For example, the standard mzv + +.. MATH:: + + \sum_{x,y,z \ge 1} \frac{1}{x (x + y) (x + y + z)^2} + +would be encoded as ``[(1, (1, 0, 0)), (1, (1, 1, 0)), (2, (1, 1, 1))]``. +The ordering of the pairs ``(v, p)`` is irrelevant. +""" +#***************************************************************************** +# Copyright (C) 2019-2023 Vincent Delecroix <20100.delecroix@gmail.com> +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +#***************************************************************************** + +import itertools +from collections import defaultdict +import cypari2 + +from sage.misc.cachefunc import cached_method, cached_function +from sage.all import ZZ, QQ, matrix, bernoulli_polynomial, prod, FreeModule +from sage.arith.all import binomial, factorial +from sage.sets.disjoint_set import DisjointSet +from sage.combinat.permutation import Permutations + + +from surface_dynamics.misc.linalg import linearly_independent_vectors +from .options import VERBOSE, DIVERGENT_MZV + +try: + from sage.modular.multiple_zeta import Multizetas +except (ImportError, cypari2.handle_error.PariError): + def Multizetas(*args, **kwds): + raise ValueError('your sage version does not support multiple zeta values') + +from surface_dynamics.misc.linalg import disjoint_vectors +from .gmzv_two_variables import Z2 +from .gmzv_three_variables import Z3, Z3_sort_abc + +# NOTE: ordering of variables (= inverse base 2) +# 1 2 3 4 5 6 7 ... +# x y x+y z x+z y+z x+y+z t x+t (y+t) (x+y+t) (z+t) ... +# (1000) (0100) (1100) (0010) (1010) (0110) (1110) (0001) (1001) +def linear_forms(d): + r""" + Return the `2^d - 1` linear forms with coefficients `\{0, 1\}` in ``d`` variables. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import linear_forms + + sage: list(linear_forms(1)) + [(1)] + + sage: list(linear_forms(2)) + [(1, 0), (0, 1), (1, 1)] + + sage: L1 = list(linear_forms(3)) + sage: L2 = L1[:] + sage: L2.sort(key = lambda x: x[::-1]) + sage: assert L1 == L2 + """ + F = FreeModule(ZZ,d) + for n in range(1, 2**d): + v = F(ZZ(n).digits(2, padto=d)) + v.set_immutable() + yield v + + +class DivergentZetaError(Exception): + pass + + +def handle_term(n, den_tuple): + r""" + Return a linear combination of standard mzv equivalent to the gmzv ``(n, den_tuple)``. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import handle_term, is_convergent + + sage: M = Multizetas(QQ) # optional: mzv + + sage: V1 = FreeModule(ZZ, 1) + sage: v = V1((1,)); v.set_immutable() + sage: dt = ((v,3),) + sage: assert is_convergent(1, dt) and handle_term(1, dt) == M((3,)) # optional: mzv + + + sage: V2 = FreeModule(ZZ, 2) + sage: va = V2((1,0)); va.set_immutable() + sage: vb = V2((0,1)); vb.set_immutable() + sage: vc = V2((1,1)); vc.set_immutable() + sage: dt = ((va,2), (vc,3)) + sage: assert is_convergent(2, dt) and handle_term(2, dt) == M((2,3)) # optional: mzv + sage: dt1 = ((va,2),(vb,3)) + sage: dt2 = ((va,3),(vb,2)) + sage: assert is_convergent(2,dt1) and is_convergent(2,dt2) # optional: mzv + sage: assert handle_term(2, ((va,2), (vb,3))) == handle_term(2, ((va,3), (vb,2))) == M((2,)) * M((3,)) # optional: mzv + + sage: V3 = FreeModule(ZZ, 3) + sage: va = V3((1,0,0)); va.set_immutable() + sage: vb = V3((0,1,0)); vb.set_immutable() + sage: vc = V3((0,0,1)); vc.set_immutable() + sage: vd = V3((1,1,0)); vd.set_immutable() + sage: ve = V3((1,0,1)); ve.set_immutable() + sage: vf = V3((0,1,1)); vf.set_immutable() + sage: vg = V3((1,1,1)); vg.set_immutable() + sage: assert handle_term(3, ((va,2), (vd,3), (vg,4))) == M((2,3,4)) # optional: mzv + sage: assert handle_term(3, ((va,2), (vb,3), (vc,4))) == handle_term(3, ((va,3), (vb,2), (vc,4))) # optional mzv + sage: assert handle_term(3, ((va,2), (vb,3), (vc,4))) == M((2,)) * M((3,)) * M((4,)) # optional: mzv + sage: assert handle_term(3, ((va,1), (vc,2), (vd,3))) == handle_term(3, ((va,1), (vb,2), (ve,3))) # optional: mzv + sage: assert handle_term(3, ((va,1), (vc,2), (vd,3))) == handle_term(3, ((va,2), (vb,1), (vf,3))) # optional: mzv + sage: assert handle_term(3, ((va,1), (vc,2), (vd,3))) == M((2,)) * M((1,3)) # optional: mzv + """ + if n == 1: + M = Multizetas(QQ) + ans = M.zero() + for v, p in den_tuple: + # sum (v[0]x)^p + ans += QQ.one() / v[0] **p * M((p,)) + return ans + + elif n == 2: + dat = to_Z2(den_tuple) + if dat is None: + raise NotImplementedError("generalized mzv {}".format(den_tuple)) + return Z2(*dat) + + elif n == 3: + dat = to_Z3(den_tuple) + if dat is None: + raise NotImplementedError("generalized mzv {}".format(den_tuple)) + return Z3(*dat) + + return _handle_term(n, den_tuple) + + +@cached_function +def _handle_term(n, den_tuple): + if VERBOSE: + print("handle term({}, {})".format(n, den_tuple)) + assert all(len(v) == n for v,p in den_tuple), (n, den_tuple) + + if any(x < 0 or x > 1 for v,p in den_tuple for x in v): + raise ValueError("unhandled zeta values {}".format(den_tuple)) + + # 0. check convergence + if not DIVERGENT_MZV and not is_convergent(n, den_tuple): + raise DivergentZetaError("{} Z({})".format(n, den_tuple)) + + # 1. multizeta + Z = is_multizeta(n, den_tuple) + if Z is not None: + if VERBOSE: + print("standard multizeta") + return Z + + # 2. apply stuffle + P = is_stufflisable(n, den_tuple) + if P is not None: + if VERBOSE: + print("stuffle") + return sum(handle_term(nn, dd) for nn, dd in P) + + # 3. equal rows + # HOW? + + # 4. diminish the number of linear forms without creating convergence problem + data = has_term_sum_of_smaller_terms(n, den_tuple) + if data is not None: + if VERBOSE: + print("relation between linear forms: L_i = sum L_j with L_i > L_j") + return sum(coeff * handle_term(n, new_den_tuple) for coeff, new_den_tuple in kill_relation(n, den_tuple, data[0], data[1])) + + # 5. make "big terms" + data = is_reducible(n, den_tuple) + if data is not None: + if VERBOSE: + print("reduction") + return sum(coeff * handle_term(n, new_den_tuple) for coeff, new_den_tuple in data) + + # 5. 3-variables + if n == 3: + dat = to_Z3(den_tuple) + if dat is None: + raise NotImplementedError("generalized mzv {}".format(den_tuple)) + return Z3(*dat) + + raise NotImplementedError("unhnandled generalized multiple zeta value {}".format(den_tuple)) + + +def clean_term(n, den_tuple): + D = {} + for den, p in den_tuple: + if den in D: + D[den] += p + if not D[den]: + del D[den] + else: + D[den] = p + return tuple(sorted(D.items())) + + +def to_Z2(den_tuple): + r""" + Converts ``den_tuple`` to arguments ``(a, b, c)`` for the function `Z2`. + """ + a = b = c = 0 + for v,p in den_tuple: + v = tuple(v) + if len(v) != 2: + raise ValueError + if v == (1,0): + a += p + elif v == (0,1): + b += p + elif v == (1,1): + c += p + else: + return + return a,b,c + + +def to_Z3(den_tuple, sort=True): + r""" + Converts ``den_tuple`` to arguments ``(a, b, c, d, e, f, g)`` for the function `Z3`. + """ + a = b = c = d = e = f = g = 0 + for v,p in den_tuple: + v = tuple(v) + if len(v) != 3: + raise ValueError + if v == (1,0,0): + a += p + elif v == (0,1,0): + b += p + elif v == (0,0,1): + c += p + elif v == (1,1,0): + d += p + elif v == (1,0,1): + e += p + elif v == (0,1,1): + f += p + elif v == (1,1,1): + g += p + else: + return + # now normalize + if sort: + a,b,c,d,e,f,g = Z3_sort_abc(a,b,c,d,e,f,g) + if b == c: + if a == b: + d,e,f = sorted([d,e,f], reverse=True) + else: + d,e = sorted([d,e], reverse=True) + return a,b,c,d,e,f,g + + +def linear_form(R, v): + return sum(R.gen(i) for i,j in enumerate(v) if j) + + +def negative_rays(n): + l = [] + v = [0]*n + for i in range(n): + v[i] = -1 + l.append(v[:]) + v[i] = 0 + return l + + +def is_convergent(n, den_tuple): + r""" + Test whether the GMZV is convergent. + + TESTS:: + + sage: import itertools + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import is_convergent + sage: V = FreeModule(ZZ, 3) + sage: va = V((1,0,0)); va.set_immutable() + sage: vb = V((0,1,0)); vb.set_immutable() + sage: vc = V((0,0,1)); vc.set_immutable() + sage: vd = V((1,1,0)); vd.set_immutable() + sage: ve = V((1,0,1)); ve.set_immutable() + sage: vf = V((0,1,1)); vf.set_immutable() + sage: vg = V((1,1,1)); vg.set_immutable() + sage: gens = [va,vb,vc,vd,ve,vf,vg] + sage: N = 0 + sage: for p in itertools.product([0,1,2], repeat=7): # optional: mzv + ....: if sum(map(bool,p)) == 3 and is_convergent(3, list(zip(gens,p))): + ....: print(p) + ....: N += 1 + (0, 0, 0, 0, 1, 1, 2) + (0, 0, 0, 0, 1, 2, 1) + (0, 0, 0, 0, 1, 2, 2) + (0, 0, 0, 0, 2, 1, 1) + (0, 0, 0, 0, 2, 1, 2) + (0, 0, 0, 0, 2, 2, 1) + (0, 0, 0, 0, 2, 2, 2) + (0, 0, 0, 1, 0, 1, 2) + ... + (2, 0, 2, 0, 0, 2, 0) + (2, 0, 2, 2, 0, 0, 0) + (2, 1, 0, 0, 0, 0, 2) + (2, 1, 0, 0, 0, 2, 0) + (2, 2, 0, 0, 0, 0, 2) + (2, 2, 0, 0, 0, 2, 0) + (2, 2, 0, 0, 2, 0, 0) + (2, 2, 2, 0, 0, 0, 0) + sage: print(N) # optional: mzv + 125 + """ + from sage.geometry.polyhedron.constructor import Polyhedron + from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing + + assert all(len(v) == n for v, p in den_tuple), (n, den_tuple) + + # TODO: fast code path + + R = PolynomialRing(QQ, 'x', n) + x = R.gens() + den_poly = prod(linear_form(R, v)**p for v, p in den_tuple) + newton_polytope = Polyhedron(vertices=den_poly.exponents(), rays=negative_rays(n)) + V = newton_polytope.intersection(Polyhedron(rays=[[1]*n])).vertices() + r = max(max(v.vector()) for v in V) + return r > 1 + + +def convergent_multizeta(t): + r""" + Multizeta value at a convergent index ``t``. + + TESTS:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import convergent_multizeta + sage: assert all(convergent_multizeta(t) == Multizeta(*t) for t in [(2,),(3,),(1,2),(3,2),(1,1,2)]) # optional: mzv + + sage: convergent_multizeta((0,3)) # optional: mzv + ζ(2) - ζ(3) + sage: convergent_multizeta((0,2,2)) # optional: mzv + ζ(1,2) - ζ(2,2) + sage: convergent_multizeta((1,0,3)) # optional: mzv + ζ(1,2) - ζ(1,3) - ζ(2) + ζ(3) + sage: convergent_multizeta((0,1,3)) # optional: mzv + -ζ(1,3) + ζ(2) - ζ(3) + + sage: convergent_multizeta((0, 4)) # optional: mzv + ζ(3) - ζ(4) + sage: convergent_multizeta((-1, 5)) # optional: mzv + 1/2*ζ(3) - 1/2*ζ(4) + sage: convergent_multizeta((-2, 5)) # optional: mzv + 1/3*ζ(2) - 1/2*ζ(3) + 1/6*ζ(4) + + sage: convergent_multizeta((-1, 3, 4)) # optional: mzv + 1/2*ζ(1,4) - 1/2*ζ(2,4) + sage: convergent_multizeta((-1, -1, 8)) # optional: mzv + 1/8*ζ(4) - 5/12*ζ(5) + 3/8*ζ(6) - 1/12*ζ(7) + sage: convergent_multizeta((-2,-2,10)) # optional: mzv + 1/18*ζ(4) - 4/15*ζ(5) + 31/72*ζ(6) - 1/4*ζ(7) + 1/72*ζ(8) + 1/60*ζ(9) + sage: convergent_multizeta((-1, -2, 10)) # optional: mzv + 1/10*ζ(5) - 3/8*ζ(6) + 5/12*ζ(7) - 1/8*ζ(8) - 1/60*ζ(9) + + sage: convergent_multizeta((4,-4,10)) # optional: mzv + -1/3*ζ(1,10) + 1/30*ζ(3,10) + 1/5*ζ(4,5) - 1/2*ζ(4,6) + 1/3*ζ(4,7) - 1/30*ζ(4,9) - 1/10*ζ(8) - 2/5*ζ(9) + 1/2*ζ(10) + sage: convergent_multizeta((2,-1,8)) # optional: mzv + -1/2*ζ(1,8) + 1/2*ζ(2,6) - 1/2*ζ(2,7) - 1/2*ζ(7) + 1/2*ζ(8) + + sage: convergent_multizeta((0,3,2)) # optional: mzv + ζ(2,2) - ζ(3,2) + + Divergent cases:: + + sage: convergent_multizeta((0,2)) # optional: mzv + Traceback (most recent call last): + ... + DivergentZetaError: divergent multizeta value (0, 2) + sage: convergent_multizeta((0,0,3)) # optional: mzv + Traceback (most recent call last): + ... + DivergentZetaError: divergent multizeta value (0, 0, 3) + sage: convergent_multizeta((1,0,2)) # optional: mzv + Traceback (most recent call last): + ... + DivergentZetaError: divergent multizeta value (1, 0, 2) + sage: convergent_multizeta((0,1,2)) # optional: mzv + Traceback (most recent call last): + ... + DivergentZetaError: divergent multizeta value (0, 1, 2) + """ + if VERBOSE: + print("convergent_multizeta({})".format(t)) + n = len(t) + if not DIVERGENT_MZV: + for i in range(1,n+1): + if sum(t[-i:]) <= i: + raise DivergentZetaError("divergent multizeta value {}".format(t)) + if all(x > 0 for x in t): + M = Multizetas(QQ) + W = M.basis().keys() + return M.term(W(t)) + else: + # find the first non-positive index and drop the summation on the + # corresponding variable using Faulhaber's formula + i = 0 + while t[i] > 0: + i += 1 + t = list(t) + tt = t[:i] + t[i+1:] + x = QQ['x'].gen() + if t[i] == 0: + # bug: bernoulli_polynomial(x, 0) is an integer + P = x - 1 + else: + P = bernoulli_polynomial(QQ['x'].gen(), -t[i]).integral() + s = Multizetas(QQ).zero() + for c,e in zip(P.coefficients(), P.exponents()): + tt[i] -= e + s += c * convergent_multizeta(tt[:]) + tt[i] += e + if i > 0: + Q = P(x+1) + for c,e in zip(Q.coefficients(), Q.exponents()): + tt[i-1] -= e + s -= c * convergent_multizeta(tt[:]) + tt[i-1] += e + return s + + +def is_multizeta(n, den_tuple): + r""" + Return the corresponding zeta values if it is one. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import linear_forms, is_multizeta + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + + sage: is_multizeta(3, ((vb, 2), (vf, 5), (vg, 2))) # optional: mzv + ζ(2,5,2) + + sage: is_multizeta(3, [(vg, 5)]) # optional: mzv + 1/2*ζ(3) - 3/2*ζ(4) + ζ(5) + + sage: assert is_multizeta(3, ((va,3), (vb,3), (vc,3))) is None # optional: mzv + sage: assert is_multizeta(3, ((vb,2), (ve,5), (vg,1))) is None # optional: mzv + """ + assert all(len(v) == n for v,p in den_tuple), (n, den_tuple) + + by_support = [None for _ in range(n+1)] + for vp in den_tuple: + v,p = vp + if any(x > 1 for x in v): + return + m = sum(vp[0]) + if by_support[m] is None: + by_support[m] = vp + else: + return + + k = 0 + linear_order = [None] * n + for step, vp in enumerate(by_support): + if vp is None: + continue + v,p = vp + # discover new variables + for i in range(n): + if v[i] and linear_order[i] is None: + linear_order[i] = k + k += 1 + if k > step: + return + + if all(v is not None for v in linear_order[1:]): + return convergent_multizeta([0 if vp is None else vp[1] for vp in by_support[1:]]) + + +def is_stufflisable(n, den_tuple): + m = len(den_tuple) + possibilities = [] + for i in range(n): + for j in range(i): + if all(not (v[i] and v[j]) for v,p in den_tuple): + # see if we can perform a stuffle that will lower rank in next step + for l in range(n): + if all(v[i] + v[j] == v[l] for v,p in den_tuple): + return stuffle(n, den_tuple, i, j) + possibilities.append((i,j)) + + # pick the first candidate ? + # multiple columns shuffling ? + if possibilities: + i, j = possibilities[0] + return stuffle(n, den_tuple, i, j) + + +def stuffle(n, den_tuple, i, j): + r""" + Apply a stuffle relation for the variables i and j. + """ + m = matrix([v for v,p in den_tuple]).transpose() + + # big subsimplex 1 + mm = m.__copy__() + mm[i] += mm[j] + D1 = clean_term(n, tuple((c,p) for c,(v,p) in zip(mm.columns(), den_tuple))) + + # big subsimplex 2 + mm = m.__copy__() + mm[j] += mm[i] + D2 = clean_term(n, tuple((c,p) for c,(v,p) in zip(mm.columns(), den_tuple))) + + # small subsimplex + mm = m.__copy__() + mm[i] += mm[j] + mm = mm.delete_rows([j]) + D3 = clean_term(n-1, tuple((c,p) for c,(v,p) in zip(mm.columns(), den_tuple))) + + return [(n,D1),(n,D2),(n-1,D3)] + + +def is_product(n, den_tuple): + r""" + Test whether the GMZV is a product of GMZV of smaller dimensions. If so, return + a pair of linear combinations whose product is equivalent. + + INPUT: + + - ``n`` - number of variables + + - ``den_tuple`` - tuple of pairs ``(vector, power)`` + + TESTS:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import is_product + sage: term1, term2 = is_product(3, [((1,0,0),2), ((0,1,0),5), ((1,1,0),1), ((0,0,1),3)]) + sage: term1 + (2, (((0, 1), 5), ((1, 0), 2), ((1, 1), 1))) + sage: term2 + (1, (((1), 3),)) + + sage: is_product(3, [((1,0,0),2), ((0,1,0),3), ((1,0,1),1), ((0,0,1),5)]) + [(2, (((0, 1), 5), ((1, 0), 2), ((1, 1), 1))), (1, (((1), 3),))] + + sage: is_product(3, [((1,1,1),3)]) is None + True + """ + D = DisjointSet(n) + assert all(len(v) == n for v,p in den_tuple), (n, den_tuple) + + # 1. product structure + for v,_ in den_tuple: + i0 = 0 + while not v[i0]: + i0 += 1 + i = i0 + 1 + while i < n: + if v[i]: + D.union(i0, i) + i += 1 + if D.number_of_subsets() == 1: + # no way to split variables + return + + # split variables + Rdict = D.root_to_elements_dict() + keys = sorted(Rdict.keys()) + key_indices = {k: i for i,k in enumerate(keys)} + values = [Rdict[k] for k in keys] + values_indices = [{v:i for i,v in enumerate(v)} for v in values] + n_list = [len(J) for J in values] + F = [FreeModule(ZZ, nn) for nn in n_list] + new_terms = [[] for _ in range(len(Rdict))] + for v,p in den_tuple: + i0 = 0 + while not v[i0]: + i0 += 1 + i0 = D.find(i0) + assert all(D.find(i) == i0 for i in range(n) if v[i]), (i0, [D.find(i) for i in range(n) if v[i]]) + k = key_indices[i0] + vv = F[k]() + for i in range(n): + if v[i]: + vv[values_indices[k][i]] = v[i] + vv.set_immutable() + new_terms[k].append((vv,p)) + + return list(zip(n_list, [tuple(sorted(terms)) for terms in new_terms])) + +# TODO: apply a power of a relation at once with multinomial coefficients +def apply_relation(n, den_tuple, i, relation): + # iterator through the pairs (coeff, den_tuples) obtained by applying relation on the i-th term + ans = [] + for j in range(n): + if relation[j] and j != i: + term = list(den_tuple) + v, p = term[i] + term[i] = (v, p+1) + v, p = term[j] + term[j] = (v, p-1) + yield (relation[j], tuple(term)) + + +def kill_relation(n, den_tuple, i, relation): + r""" + Make calls to :func:`apply_relation` until we get rid of term. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import kill_relation + + sage: V = FreeModule(ZZ, 2) + sage: va = V((1,0)); va.set_immutable() + sage: vb = V((0,1)); vb.set_immutable() + sage: vc = V((1,1)); vc.set_immutable() + + sage: den_tuple = ((va,2),(vb,3),(vc,4)) + sage: kill_relation(2, den_tuple, 2, [1,1,0]) + [(1, (((0, 1), 3), ((1, 1), 6))), + (2, (((0, 1), 2), ((1, 1), 7))), + (1, (((1, 0), 2), ((1, 1), 7))), + (3, (((0, 1), 1), ((1, 1), 8))), + (3, (((1, 0), 1), ((1, 1), 8)))] + """ + assert len(relation) == len(den_tuple) + assert 0 <= i < len(den_tuple) + assert sum(relation[j] * den_tuple[j][0] for j in range(len(den_tuple))) == den_tuple[i][0] + D = {den_tuple: QQ.one()} + todo = [den_tuple] + while todo: + den_tuple = todo.pop(0) + coeff1 = D.pop(den_tuple) + for coeff, new_den_tuple in apply_relation(n, den_tuple, i, relation): + coeff *= coeff1 + if new_den_tuple in D: + D[new_den_tuple] += coeff + elif any(not new_den_tuple[j][1] for j in range(len(den_tuple))): + new_den_tuple = tuple(x for x in new_den_tuple if x[1]) + if new_den_tuple not in D: + D[new_den_tuple] = coeff + else: + D[new_den_tuple] += coeff + else: + todo.append(new_den_tuple) + D[new_den_tuple] = coeff + + return [(coeff, new_den_tuple) for new_den_tuple, coeff in D.items()] + + +def has_term_sum_of_smaller_terms(n, den_tuple): + r""" + Look for a vector v_i and {v_j} with sum(v_j) = v_i + + This is useful only useful when the linear forms have relations between them. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import linear_forms, has_term_sum_of_smaller_terms + + sage: va,vb,vc = linear_forms(2) + sage: has_term_sum_of_smaller_terms(2, ((va,1),(vb,1),(vc,1))) + [2, [1, 1, 0]] + sage: assert has_term_sum_of_smaller_terms(2, ((va,1),(vc,1))) is None + sage: assert has_term_sum_of_smaller_terms(2, ((vb,1),(vc,1))) is None + sage: assert has_term_sum_of_smaller_terms(2, ((va,1),(vb,1))) is None + + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + sage: has_term_sum_of_smaller_terms(3, ((va,1),(vb,1),(vc,1),(vg,1))) + [3, [1, 1, 1, 0]] + sage: has_term_sum_of_smaller_terms(3, ((va,1),(vb,1),(ve,1),(vf,1),(vg,1))) + [4, [0, 1, 1, 0, 0]] + sage: has_term_sum_of_smaller_terms(3, ((vd,1),(ve,1),(vf,1),(vg,1))) + [3, [1/2, 1/2, 1/2, 0]] + sage: has_term_sum_of_smaller_terms(3, ((va,1),(vd,1),(ve,1),(vg,1))) + [3, [-1, 1, 1, 0]] + """ + l = len(den_tuple) + for iv in range(l): + v = den_tuple[iv][0] + if sum(v) == 1: + continue + + candidates = [k for k,(u,_) in enumerate(den_tuple) if sum(u) < sum(v) and all(u[i] <= v[i] for i in range(n))] + M = matrix(QQ, [den_tuple[k][0] for k in candidates]) + try: + coeffs = M.solve_left(v) + except ValueError: + continue + full = [0] * l + for k,c in zip(candidates, coeffs): + full[k] = c + return [iv, full] + + +def write_array(level, explanations): + for i,key in enumerate(sorted(level, key=lambda x: 1000 if level[x] is None else level[x])): + print("%s & %s \\\\" % (lin_prod(key),'\\infty' if level[key] is None else level[key])) + + +def is_reducible(n, den_tuple): + r""" + Test whether the given GMZV is reducible and if so return an equivalent linear expression of simpler GMZV. + + We call a GMZV is reducible If (x1+x2+...+xd) is not present, use a linear relation to create it. Then + try to kill using other forms. Then try to write as P(x1,x2,...,x_{d-1}) (x1+x2+...+xd) and recurse. + + For example, this function expresses all Tornheim sum as linear combinations of MZV. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import linear_forms, is_reducible + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + sage: is_reducible(3, ((va, 3), (vb, 3), (vc, 3))) + [(1, (((0, 0, 1), 3), ((0, 1, 1), 3), ((1, 1, 1), 3))), + (1, (((0, 1, 0), 3), ((0, 1, 1), 3), ((1, 1, 1), 3))), + (3, (((0, 0, 1), 2), ((0, 1, 1), 4), ((1, 1, 1), 3))), + (3, (((0, 1, 0), 2), ((0, 1, 1), 4), ((1, 1, 1), 3))), + ... + (90, (((1, 0, 0), 1), ((1, 0, 1), 1), ((1, 1, 1), 7))), + (90, (((0, 1, 0), 1), ((1, 1, 0), 1), ((1, 1, 1), 7))), + (90, (((1, 0, 0), 1), ((1, 1, 0), 1), ((1, 1, 1), 7)))] + """ + if len(den_tuple) == 1: + return + + # force the max vector (1,1,...,1) to appear possibly with an exponent 0 + F = FreeModule(ZZ, n) + vmax = F([1] * n) + vmax.set_immutable() + imax = None + for i, (v, p) in enumerate(den_tuple): + if v == vmax: + imax = i + break + if imax is None: + imax = len(den_tuple) + den_tuple = den_tuple + ((vmax, 0),) + if imax != len(den_tuple) - 1: + den_tuple = list(den_tuple) + den_tuple.append(den_tuple.pop(imax)) + den_tuple = tuple(den_tuple) + imax = len(den_tuple) - 1 + + assert den_tuple[imax][0] == vmax + + M = matrix(QQ, [v for v,p in den_tuple if v != vmax]) + try: + relation = M.solve_left(vmax) + except ValueError: + if den_tuple[-1][1] == 0: + return + den_tuple2 = den_tuple[:-1] + variables = set().union(*[[i for i in range(n) if v[i]] for v,p in den_tuple2]) + if len(variables) == n: + return + killed = [(1, den_tuple)] + else: + killed = kill_relation(n, den_tuple, imax, list(relation) + [0]) + + ans = defaultdict(QQ) + for coeff, den_tuple2 in killed: + assert den_tuple2[-1][0] == vmax + pmax = den_tuple2[-1][1] + assert pmax + den_tuple2 = den_tuple2[:-1] + + variables = set().union(*[[i for i in range(n) if v[i]] for v,p in den_tuple2]) + if len(variables) == n: + # removing the maximal vector (1,1,...,1) is not enough to make a variable disappear + data = is_reducible(n, den_tuple2) + if data is None: + data = [(1, den_tuple2)] + else: + # less variables! + nn = len(variables) + variables = sorted(variables) + new_indices = {j:i for i,j in enumerate(variables)} + G = FreeModule(ZZ, nn) + new_den_tuple2 = [] + for v,p in den_tuple2: + vv = G([v[i] for i in variables]) + vv.set_immutable() + new_den_tuple2.append((vv,p)) + new_den_tuple2 = tuple(new_den_tuple2) + data = is_reducible(nn, new_den_tuple2) + if data is None: + data = [(1, new_den_tuple2)] + + # lift to the n variables version + new_data = [] + for coeff3, den_tuple3 in data: + den_tuple3 = [(F([v[new_indices[j]] if j in new_indices else 0 for j in range(n)]), p) for v,p in den_tuple3] + for v,p in den_tuple3: + v.set_immutable() + new_data.append((coeff3, tuple(sorted(den_tuple3)))) + data = new_data + + # update the answer + for coeff3, den_tuple3 in data: + imax = None + for i,(v,p) in enumerate(den_tuple3): + if v == vmax: + imax = i + break + if imax is None: + den_tuple3 = den_tuple3 + ((vmax,pmax),) + else: + den_tuple3 = den_tuple3[:imax] + ((vmax, den_tuple3[imax][1] + pmax),) + den_tuple3[imax+1:] + ans[den_tuple3] += coeff * coeff3 + + if len(ans) > 1: + return [(coeff, den_tuple) for den_tuple, coeff in ans.items()] + + +class GeneralizedMultipleZetaFunction: + r""" + Generalized multiple zeta function. + + A generalized multiple zeta function is a series of the form + + .. MATH:: + + Z(s_1, \ldots, s_r) = \sum_{x_1, x_2, \ldots, x_d \geq 1} L_1(x)^{-s_1} L_2(x)^{-s_2} \cdots L_r(x)^{-s_r} + + where the sum is over integers $x_i$ and the $L_i$ are linear forms in `x = + (x_1, x_2, \ldots, x_d)` with `\{0,1\}` coefficients. It could be thought of + as as generalized polylogarithm in the variable `s = (s_1, \ldots, s_d)` + evaluated at `z = (1, 1, \ldots, 1)`. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import linear_forms, GeneralizedMultipleZetaFunction + + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + sage: f = GeneralizedMultipleZetaFunction([va, vb, vc, vg]) + sage: f + GeneralizedMultipleZetaFunction('(0)(1)(2)(0+1+2)') + sage: f((1, 1, 1, 1)) + 6*ζ(1,1,2) + """ + def __init__(self, *args, as_rows=False, reduced=None): + if len(args) == 1 and isinstance(args[0], str): + s = args[0].replace(' ', '').replace(')*(', ')(') + if not s.startswith('(') or not s.endswith(')'): + raise ValueError + s = [lin.split('+') for lin in s[1:-1].split(')(')] + letters = set().union(*s) + indices = {letter:i for i,letter in enumerate(sorted(letters))} + F = FreeModule(ZZ, len(indices)) + args = ([sum(F.gen(indices[letter]) for letter in lin) for lin in s],) + + self._entries = matrix(ZZ, *args) + if not as_rows: + self._entries = self._entries.transpose() + self._reduced = reduced + + def children(self): + r""" + Compatibility with function with recursive calls. + """ + return [] + + @classmethod + def from_matrix(cls, m, canonicalize, sort_rows, set_immutable, reduced): + L = GeneralizedMultipleZetaFunction.__new__(GeneralizedMultipleZetaFunction) + L._reduced = reduced + if canonicalize: + L._entries = m.permutation_normal_form() + nr = m.nrows() + nc = m.ncols() + elif sort_rows: + L._entries = matrix(ZZ, sorted(m.rows(), reverse=True)) + else: + L._entries = matrix(ZZ, m) + + if set_immutable: + L._entries.set_immutable() + return L + + def __getitem__(self, key): + return self._entries[key] + + def __eq__(self, other): + if type(self) is not type(other): + raise TypeError + return self._entries == other._entries + + def __ne__(self, other): + if type(self) is not type(other): + raise TypeError + return self._entries != other._entries + + def __call__(self, e): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, convergent_multizeta + sage: GeneralizedMultipleZetaFunction([[1,0,0],[1,1,1],[1,0,1]])((2,4,3)) + ζ(2,3,4) + sage: GeneralizedMultipleZetaFunction([[1,0,0],[1,1,1]])((2,4)) == convergent_multizeta((2,0,4)) + True + + sage: G = GeneralizedMultipleZetaFunction([[1,1,1,1],[1,1,1,0],[1,1,0,0]]) + sage: G((2,2,2)) + ζ(1,2,2) - ζ(2,2,2) + """ + M = Multizetas(QQ) + d = self.nrows() + n = self.ncols() + if len(e) != n: + raise ValueError('invalid argument: generalized mzv with {} columns but got e={}'.format(n, e)) + if d == 1: + # pure zeta + return M((sum(e),)) + elif d == 2: + return Z2(*to_Z2(list(zip(self.columns(), e)))) + elif d == 3: + return Z3(*to_Z3(list(zip(self.columns(), e)))) + + dat = self.is_multizeta() + if dat: + s = [0] * self.nrows() + for x, c in zip(e, self.columns()): + s[sum(c)-1] += x + return convergent_multizeta(s) + + raise NotImplementedError('generalized mzv {}'.format(self.lin_prod_string())) + + def copy(self, immutable=False): + if immutable and self._entries.is_immutable(): + return self + L = GeneralizedMultipleZetaFunction.__new__(GeneralizedMultipleZetaFunction) + L._reduced = self._reduced + L._entries = self._entries.__copy__() + if immutable: + L._entries.set_immutable() + return L + + def __repr__(self): + return "GeneralizedMultipleZetaFunction({!r})".format(self.lin_prod_string()) + + def lin_prod_string(self): + nr = self.nrows() + return ''.join('(' + '+'.join(str(i) for i in range(nr) if c[i]) + ')' for c in self._entries.columns()) + + def subs(self, substitution_dict, canonicalize=False, set_immutable=False, reduced=False): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction + sage: L = GeneralizedMultipleZetaFunction([[1,0,0],[0,1,0],[0,0,1]], as_rows=True) + sage: L.subs({0: [(1,0),(1,1)], 1: [(1,0),(1,1),(1,2)]}) + GeneralizedMultipleZetaFunction('(0+1)(0+1)(1+2)') + """ + M = self._entries.__copy__() + for i, lin_comb in substitution_dict.items(): + M[i] = sum(c * self._entries[j] for c,j in lin_comb) + return GeneralizedMultipleZetaFunction.from_matrix(M, canonicalize, False, set_immutable, reduced) + + def set_immutable(self): + self._entries.set_immutable() + + @cached_method + def is_reduced(self): + if self._reduced is None: + self._reduced = self._entries.column_module().rank() == self._entries.ncols() + return self._reduced + + def __hash__(self): + return hash(self._entries) + + def __eq__(self, other): + if type(self) is not type(other): + raise TypeError + return self._entries == other._entries + def __ne__(self, other): + if type(self) is not type(other): + raise TypeError + return self._entries == other._entries + def nrows(self): + return self._entries.nrows() + def ncols(self): + return self._entries.ncols() + def rows(self): + return self._entries.rows() + def columns(self): + return self._entries.columns() + def delete_row(self, i): + self._entries = self._entries.delete_rows([i]) + def delete_column(self, i): + self._entries = self._entries.delete_columns([i]) + def symmetric(self): + r""" + Return all possible symmetric version obtained by permuting variables + """ + d = self.ncols() + for p in Permutations(self.columns()): + yield GeneralizedMultipleZetaFunction.from_matrix(matrix(ZZ, p).transpose(), False, True, True, self._reduced) + + def has_no_zero_row_or_column(self): + M = self._entries + return all(any(self._entries[i][j] for j in range(M.ncols())) for i in range(M.nrows())) and all(any(self._entries[i][j] for i in range(M.nrows())) for j in range(M.ncols())) + + def normalized_columns(self): + f = ZZ(self.nrows()).factorial() + nc = [f//sum(v) * v for v in self.columns()] + for v in nc: v.set_immutable() + return nc + + def is_product(self): + r""" + Return ``None`` or a list of matrices. + """ + d = self.nrows() + n = self.ncols() + M = self._entries + + # make group of rows + D = DisjointSet(d) + for c in self.columns(): + i0 = 0 + while i0 < d and not c[i0]: + i0 += 1 + if i0 == d: + raise ValueError('zero column') + i = i0 + 1 + while i < d: + if c[i]: + D.union(i0, i) + i += 1 + + if D.number_of_subsets() == 1: + # no way to split variables + return + + ans = [] + for rows in D: + cols = set().union(*[M.nonzero_positions_in_row(i) for i in rows]) + ans.append(GeneralizedMultipleZetaFunction(M.matrix_from_rows_and_columns(rows, sorted(cols)), as_rows=True)) + return ans + + def new_columns(self, forward_only=False): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, linear_forms + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + sage: L = GeneralizedMultipleZetaFunction([va, vd, vg]) + sage: [x[0] for x in L.new_columns()] == [vb, vc, ve, vf] + True + + sage: L = GeneralizedMultipleZetaFunction([va, vg]) + sage: [x[0] for x in L.new_columns()] == [vf] + True + + sage: L = GeneralizedMultipleZetaFunction([vd,ve,vf]) + sage: for c in L.new_columns(False): print(c) + ((1, 0, 0), (1/2, 1/2, -1/2)) + ((0, 1, 0), (1/2, -1/2, 1/2)) + ((0, 0, 1), (-1/2, 1/2, 1/2)) + ((1, 1, 1), (1/2, 1/2, 1/2)) + sage: for c in L.new_columns(True): print(c) + ((1, 1, 1), (1/2, 1/2, 1/2)) + """ + d = self._entries.nrows() + cols = set(self._entries.columns()) + V = FreeModule(ZZ, d) + for n in range(1, 2**d): + c = V(ZZ(n).digits(2,padto=d)) + c.set_immutable() + if c not in cols: + try: + s = self._entries.solve_right(c) + except ValueError: + continue + if forward_only and any(not c[i] and x[i] for j,x in enumerate(self.columns()) for i in range(d) if s[j]): + continue + yield c, s + + def is_multizeta(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, linear_forms + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + sage: GeneralizedMultipleZetaFunction([va,vd,vg]).is_multizeta() + [0, 1, 2] + sage: GeneralizedMultipleZetaFunction([vd,vb,vg]).is_multizeta() + [1, 0, 2] + + sage: GeneralizedMultipleZetaFunction([vb, vg]).is_multizeta() + [1, 0, 2] + sage: GeneralizedMultipleZetaFunction([vf, vg, vc]).is_multizeta() + [2, 1, 0] + sage: GeneralizedMultipleZetaFunction([va, vf, vg]).is_multizeta() + False + sage: GeneralizedMultipleZetaFunction([vf, vg]).is_multizeta() + [1, 2, 0] + """ + M = self._entries + d = M.nrows() + columns = M.columns() + columns.sort(key = lambda x: sum(x)) + variables = [False] * M.nrows() + order = [] + for c in columns: + for j in range(d): + if variables[j]: + if not c[j]: + return False + elif c[j]: + variables[j] = True + order.append(j) + return order + + def canonicalize(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction + sage: for cols in [[[1,0,1,0],[0,0,1,0],[1,0,0,0],[0,1,0,1]], [[0,1,0],[1,0,1],[1,1,1]], [[1,0,1,0],[0,0,1,0],[0,1,0,1]], [[1,0,1],[1,1,1],[1,0,0]], [[1,0,1],[0,0,1],[0,1,1]], [[0,1],[1,0],[0,1],[1,1],[0,1]]]: + ....: L = GeneralizedMultipleZetaFunction(cols) + ....: LL = L.copy() + ....: rp, cp = LL.canonicalize() + ....: assert all(L[rp(i+1)-1,cp(j+1)-1] == LL[i,j] for i in range(L.nrows()) for j in range(L.ncols())) + """ + assert not self._entries.is_immutable() + self._entries, (pr, pc) = self._entries.permutation_normal_form(True) + nr = self._entries.nrows() + nc = self._entries.ncols() +# self._entries = self._entries.parent()([[self._entries[nr-i-1, nc-j-1] for j in range(nc)] for i in range(nr)]) +# Sr = pr.parent(); dr = Sr.degree() +# Sc = pc.parent(); dc = Sc.degree() +# pr *= Sr([dr-i for i in range(dr)]) +# pc *= Sc([dc-i for i in range(dc)]) + return (pr, pc) + + def is_convergent(self, e): + r""" + Return whether the given exponent is convergent + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import linear_forms, GeneralizedMultipleZetaFunction + sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) + sage: L = GeneralizedMultipleZetaFunction([va,vb,vc,vd,ve,vf,vg]) + sage: L + GeneralizedMultipleZetaFunction('(0)(1)(2)(0+1)(0+2)(1+2)(0+1+2)') + sage: L.is_convergent((0,0,0,0,2,1,1)) + True + """ + if self.ncols() != len(e): + raise ValueError("generalized mzv with {} variables but got e={}".format(self.ncols(), e)) + + from sage.geometry.polyhedron.constructor import Polyhedron + + m = self.ncols() # number of linear forms + d = self.nrows() # number of variables + ZZd = FreeModule(ZZ, d) + zero = ZZd.zero() + + U = set([zero]) + V = set() + W = set() + for i,c in enumerate(self.columns()): + V.clear() + V.add(zero) + for j,x in enumerate(c): + if x: + v = e[i] * ZZd.gen(j) + v.set_immutable() + V.add(v) + W.clear() + for u in U: + for v in V: + w = u + v + w.set_immutable() + W.add(w) + U, W = W, U + P = Polyhedron(vertices=U) + vertices = P.intersection(Polyhedron(rays=[[1]*d])).vertices() + r = max(max(v.vector()) for v in vertices) + return r > 1 + + def stuffle(self, subset, only_highest_weight=False, canonicalize=False, sort_rows=False): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction + + sage: L = GeneralizedMultipleZetaFunction([[1,0,0],[0,1,0],[0,0,1]]) + sage: L((2,3,4)) == sum(LL((2,3,4)) for LL in L.stuffle([0,1])) + True + sage: L((2,3,4)) == sum(LL((2,3,4)) for LL in L.stuffle([1,2])) + True + sage: L((2,3,4)) == sum(LL((2,3,4)) for LL in L.stuffle([0,1,2])) + True + """ + d = self.ncols() + n = self.nrows() + rows = self.rows() + ans = [] + s = sum(rows[i] for i in subset) + if any(x > 1 for x in s): + raise ValueError('invalid subset for stuffle') + for i in subset: + new_rows = rows[:] + new_rows[i] = s + ans.append(GeneralizedMultipleZetaFunction.from_matrix(matrix(ZZ, new_rows), canonicalize, sort_rows, True, self._reduced)) + if only_highest_weight: + return ans + for k in range(2, len(subset) + 1): + for ss in itertools.combinations(subset, k): + new_rows = [rows[i] for i in range(n) if i not in ss] + new_rows.append(s) + ans.append(GeneralizedMultipleZetaFunction.from_matrix(matrix(ZZ, new_rows), canonicalize, sort_rows, True, self._reduced)) + return ans + + def stuffles(self, only_highest_weight=False, canonicalize=False, sort_rows=False): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, generalized_multiple_zeta_functions + sage: L = GeneralizedMultipleZetaFunction([[1,0,0],[0,1,0],[0,0,1]]) + sage: for x in L.stuffles(): + ....: print(x[0], [vv.lin_prod_string() for vv in x[1]]) + [0, 1] ['(0)(0+1)(2)', '(0+1)(1)(2)', '(1)(1)(0)'] + [0, 1, 2] ['(0)(0+1)(0+2)', '(0+1)(1)(1+2)', '(0+2)(1+2)(2)', '(1)(1)(0+1)', '(1)(0+1)(1)', '(0+1)(1)(1)', '(0)(0)(0)'] + [0, 2] ['(0)(1)(0+2)', '(0+2)(1)(2)', '(1)(0)(1)'] + [1, 2] ['(0)(1)(1+2)', '(0)(1+2)(2)', '(0)(1)(1)'] + + sage: for L in generalized_multiple_zeta_functions(3, 3): + ....: for subset, stuffle in L.stuffles(): + ....: assert L((2,3,4)) == sum(LL((2,3,4)) for LL in stuffle), (L, stuffle) + ....: assert L((2,4,3)) == sum(LL((2,4,3)) for LL in stuffle), (L, stuffle) + ....: assert L((3,2,4)) == sum(LL((3,2,4)) for LL in stuffle), (L, stuffle) + """ + d = self.ncols() + rows = self.rows() + ans = [] + for subset, s in disjoint_vectors(rows, min_size=2, max_size=None): + yield subset, self.stuffle(subset, only_highest_weight, canonicalize, sort_rows) + + def bad_minimals(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction + sage: L = GeneralizedMultipleZetaFunction(matrix([[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1],[1,1,1]]).transpose()) + """ + from sage.geometry.polyhedron.constructor import Polyhedron + m = self.ncols() # number of linear forms + d = self.nrows() # number of variables + ZZd = FreeModule(ZZ, d) + ZZm = FreeModule(ZZ, m) + ZZmd = FreeModule(ZZ, d + m) + zero = ZZmd.zero() + + rays = [] + for i,c in enumerate(self.columns()): + for j,x in enumerate(c): + if x: + rays.append(ZZmd(list(ZZd.gen(j)) + list(ZZm.gen(i)))) + C = Polyhedron(rays=rays) + Q = Polyhedron(vertices=[[1]*d + [0]*m], rays=[[0]*d + list(b) for b in ZZm.basis()]) + bad_minimals = [v[d:] for v in C.intersection(Q).vertices_list()] + return bad_minimals + + +def generalized_multiple_zeta_functions(d, r): + r""" + Return the set of generalized multiple zeta functions up to the permutation + action on variables and linear forms. + + INPUT: + + - ``d`` -- dimension + + - ``r`` -- rank + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import generalized_multiple_zeta_functions + sage: for f in generalized_multiple_zeta_functions(3, 3): + ....: print(f) + GeneralizedMultipleZetaFunction('(0)(1)(2)') + GeneralizedMultipleZetaFunction('(0+1+2)(0)(1)') + GeneralizedMultipleZetaFunction('(0+1)(0+2)(1)') + GeneralizedMultipleZetaFunction('(0+1)(0+2)(0)') + GeneralizedMultipleZetaFunction('(0+1+2)(0+1)(0)') + GeneralizedMultipleZetaFunction('(0+1)(0)(2)') + GeneralizedMultipleZetaFunction('(0+1)(0+2)(1+2)') + GeneralizedMultipleZetaFunction('(0+1+2)(0+1)(0+2)') + """ + assert 1 <= r <= d + V = set() + vectors = list(linear_forms(d)) + + if d == r: + for _, lins in linearly_independent_vectors(vectors, min_size=d, max_size=d): + V.add(GeneralizedMultipleZetaFunction.from_matrix(lins, True, False, True, True)) + else: + for _, lins in linearly_independent_vectors(vectors, min_size=r, max_size=r): + v = GeneralizedMultipleZetaFunction.from_matrix(lins.transpose(), True, False, True, True) + if v.has_no_zero_row_or_column(): + V.add(v) + return V diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py new file mode 100644 index 00000000..72e28582 --- /dev/null +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py @@ -0,0 +1,812 @@ +r""" +Automatic reduction of generalized multiple zeta values + + +EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.gmzv_generic_reduction import GMZVSolver + sage: S = GMZVSolver(3, 3) + sage: S + Solve strategy d=3 r=3 (0 known and 8 unknown) + sage: S.update_multizetas() + sage: S + Solve strategy d=3 r=3 (1 known and 7 unknown) +""" +#***************************************************************************** +# Copyright (C) 2019-2023 Vincent Delecroix <20100.delecroix@gmail.com> +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +#***************************************************************************** + +import itertools +import collections + +from sage.misc.misc_c import prod +from sage.misc.prandom import random, choice +from sage.misc.cachefunc import cached_method + +from sage.rings.integer_ring import ZZ +from sage.rings.rational_field import QQ +from sage.groups.perm_gps.permgroup_named import SymmetricGroup +from sage.modules.free_module_element import vector +from sage.modules.free_module import VectorSpace, FreeModule +from sage.matrix.constructor import matrix + +from surface_dynamics.misc.linalg import linearly_independent_vectors + +from .gmzv_two_variables import Z2 +from .gmzv_three_variables import Z3 +from .generalized_multiple_zeta_values import linear_forms, convergent_multizeta, DivergentZetaError, GeneralizedMultipleZetaFunction, to_Z2, to_Z3, generalized_multiple_zeta_functions +from .options import DIVERGENT_MZV + + +class EvalAlgebraic: + r""" + """ + def __init__(self, S, v, lin_comb): + assert v.ncols() == len(lin_comb) + self.v = v + self.lin_comb = lin_comb + self.support = [i for i in range(v.ncols()) if lin_comb[i]] + columns = v.columns() + c = sum(x * col for x, col in zip(lin_comb, columns)) + assert all(x == 0 or x == 1 for x in c), c + self.c = c = vector(ZZ, c) + self.S = S + self.new_forms = {} + for i in self.support: + new_columns = columns[:] + new_columns.pop(i) + new_columns.append(c) + vv = GeneralizedMultipleZetaFunction(new_columns) + _, cp = vv.canonicalize() + vv.set_immutable() + assert S.strategy[vv] is not None + self.new_forms[i] = (vv, cp) + + def children(self): + r""" + Return the list of zeta used + """ + return [nf[0] for nf in self.new_forms.values()] + + def __repr__(self): + return 'algebraic {} -> {}'.format(self.lin_comb, self.c) + + def __call__(self, exponents): + # TODO: implement multiplicative version + e = tuple(exponents) + (0,) + D = {e: ZZ.one()} + ans = 0 + todo = [e] + while todo: + e = todo.pop(0) + if all(e[j] for j in self.support): + coeff = D.pop(e) + if not coeff: + continue + for j in self.support: + ee = list(e) + ee[j] -= 1 + ee[-1] += 1 + ee = tuple(ee) + if ee not in D: + D[ee] = self.lin_comb[j] * coeff + todo.append(ee) + else: + assert ee in todo + D[ee] += self.lin_comb[j] * coeff + + for e, coeff in D.items(): + for j in self.support: + if e[j] == 0: + break + assert e[j] == 0 + e = e[:j] + e[j+1:] + vv, cp = self.new_forms[j] + e = [e[cp(i+1)-1] for i in range(self.v.ncols())] + ans += coeff * self.S.eval(vv, e) + + return ans + + +class FunctionWrap: + def __init__(self, f, groups, exp_position=0, extra_args=None): + self.f = f + self.groups = tuple(map(tuple, groups)) + if extra_args is None: + extra_args=() + self.extra_args = extra_args + self.pos = exp_position + def __hash__(self): + hash((self.f, self.groups, self.extra_args)) + def __repr__(self): + gstring = '(' + ','.join('+'.join(map(str, g)) for g in self.groups) + ')' + frepr = None + if isinstance(self.f, GeneralizedMultipleZetaFunction): + dat = self.f.is_multizeta() + frepr = self.f.lin_prod_string() + else: + try: + is_eval = self.f.__func__ == GMZVSolver.eval + except AttributeError: + is_eval = False + else: + if is_eval: + frepr = self.extra_args[0].lin_prod_string() + if frepr is None: + frepr = repr(self.f) + + return frepr + gstring + + def __call__(self, e): + ee = [sum(e[i] for i in g) for g in self.groups] + args = self.extra_args[:self.pos] + (ee,) + self.extra_args[self.pos:] + return self.f(*args) + + +def stuffle_to_lin_comb(v, rp, cp, stuffle): + r""" + Return a linear combination of functions taking list of functions + + INPUT: + + - ``v`` + + - ``rp`` -- row permutation + + - ``cp`` -- column permutation + + - ``stuffle`` -- + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, linear_forms + sage: from surface_dynamics.generalized_multiple_zeta_values.gmzv_generic_reduction import stuffle_to_lin_comb + sage: va, vb, vc, vd, ve, vf, vg = linear_forms(3) + sage: GeneralizedMultipleZetaFunction([vd, va, vb]) + sage: S = SymmetricGroup(3) + sage: p = S.one() + sage: stuffle_to_lin_comb + """ + vtmp = v.copy() + _, _ = vtmp.canonicalize() + assert all(v._entries[rp(i+1)-1, cp(j+1)-1] == vtmp._entries[i][j] for i in range(v.nrows()) for j in range(v.ncols())) + + if isinstance(stuffle, (tuple, list)): + # the subset is given with respect to the canonical form + subset, subtasks = stuffle + subset = [rp(i + 1) - 1 for i in subset] + stuffle = v.stuffle(subset, only_highest_weight=False) + + main = stuffle[:len(subset)] + assert len(main) == len(subtasks) + for vv, (rrp, ccp, task) in zip(main, subtasks): + assert vv.nrows() == v.nrows() + assert vv.ncols() == v.ncols() + for (vvv, rrrp, cccp) in stuffle_to_lin_comb(vv, rrp*rp, ccp*cp, task): + yield (vvv, rrrp , cccp) + remainder = stuffle[len(subset):] + for vv in remainder: + assert vv.nrows() < v.nrows() + assert vv.ncols() == v.ncols() + yield (vv, None, rp.parent().one()) + else: + # check that we got the correct element + assert stuffle == vtmp + yield (stuffle, rp, cp) + + +class LinearCombination: + def __init__(self, comb): + self.comb = [(coeff, f) for coeff, f in comb if coeff] + + def children(self): + ans = [] + for coeff, f in self.comb: + if isinstance(f, FunctionWrap): + try: + is_eval = isinstance(f.f.__self__, GMZVSolver) + except AttributeError: + is_eval = False + if is_eval: + ans.append(f.extra_args[0]) + else: + ans.append(f.f) + else: + ans.append(f) + return ans + + def __repr__(self): + return ' + '.join(str(f) if coeff == 1 else '%s %s' % (coeff, f) for coeff, f in self.comb) + + def __call__(self, e): + return sum(coeff * f(e) for coeff, f in self.comb) + + +class GMZVSolver: + r""" + A solver for finding relations between generalized multiple zeta functions and values. + + At any stage, the solver has some "known" gmzv that can be reduced to linear combinations + of standard mzv and some "unknown" forms. One can then try different options to look for + new relations and getting more "known" gmzv. + + Attributes: + + - ``d`` -- dimension + + - ``r`` -- rank + + - ``V`` -- list of multiple zeta functions to be handled + + - ``strategy`` -- dictionary + + - ``what`` -- dictionary + + - ``level`` -- dictionary + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.gmzv_generic_reduction import GMZVSolver + sage: S = GMZVSolver(3, 3) + sage: S + Solve strategy d=3 r=3 (0 known and 8 unknown) + """ + def __init__(self, d, r): + self.d = d + self.r = r + self.V = list(generalized_multiple_zeta_functions(d, r)) + + # a list v -> function to evaluate at v + self.strategy = {v: None for v in self.V} + self.what = {v: None for v in self.V} + self.level = {v: None for v in self.V} + + def graphviz_string(self, filename=None): + if filename is None: + import sys + f = sys.stdout + else: + f = open(filename, 'w') + + to_index = {v:i for i,v in enumerate(self.V)} + f.write("digraph G {\n") + f.write(" rankdir=TB\n") + for u, s in self.strategy.items(): + i = to_index[u] + if isinstance(s[0], EvalAlgebraic): + f.write(" %d [style=filled fillcolor=blue]\n" % i) + elif isinstance(s[0], LinearCombination): + f.write(" %d [style=filled fillcolor=red]\n" % i) + elif isinstance(s[0], GeneralizedMultipleZetaFunction): + f.write(" %d [style=filled fillcolor=yellow]\n" % i) + for v in s[0].children(): + if v not in to_index: + continue + f.write("%d -> %d\n" % (i, to_index[v])) + + f.write("}\n") + + if filename is not None: + f.close() + + def try_small_exponents(self, e_bound=3, s_bound=7, vlist=None): + if vlist is None: + vlist = sorted(self.V, key=lambda v: self.level[v]) + for v in vlist: + for i in IntegerListsLex(length=self.r, min_part=1, max_part=e_bound, max_sum=s_bound): + if v.is_convergent(i): + try: + self.eval(v, tuple(i)) + except DivergentZetaError: + print("failure at level {} with\n{}\nand exponent {}".format(self.level[v], v, i)) + except ValueError: + print("called for divergent zeta at level {} with\n{}\nand exponent {}".format(self.level[v], v, i)) + + def eval_abstract_sum(self, s): + return sum(QQ(coeff) * self.eval_den_tuple(den_tuple) for den_tuple, coeff in s._data.items()) + + def eval_den_tuple(self, den_tuple): + columns = [] + exponents = [] + for col, exp in den_tuple._dict.items(): + columns.append(col) + exponents.append(exp) + v = GeneralizedMultipleZetaFunction(columns) + rp, cp = v.canonicalize() + v.set_immutable() + exponents = [exponents[cp(i+1)-1] for i in range(len(columns))] + return self.eval(v, exponents) + + def eval(self, v, exponents): + if len(exponents) != v.ncols(): + raise ValueError("generalized multiple zeta values with {} linear forms but got e={}".format(v.ncols(), exponents)) + if v not in self.strategy: + raise NotImplementedError('v={} not in strategy'.format(v.lin_prod_string())) + if self.strategy[v] is None: + raise ValueError('no available strategy') + f, extra_args = self.strategy[v] + return f(exponents, *extra_args) + + def __repr__(self): + known = sum(x is not None for x in self.strategy.values()) + unknown = sum(x is None for x in self.strategy.values()) + return 'Solve strategy d={} r={} ({} known and {} unknown)'.format(self.d, self.r, known, unknown) + + def known(self): + return [v for v,status in self.strategy.items() if status is not None] + + def unknown(self): + return [v for v,status in self.strategy.items() if status is None] + + def update_multizetas(self): + r""" + Set trivial evaluation for multizetas and set their levels to 0 + """ + full_mzv = [[1] * k + [0] * (self.d - k) for k in range(self.d-1, 0, -1)] + for cols in itertools.combinations(full_mzv, self.r-1): + v = GeneralizedMultipleZetaFunction(([1]*self.d,) + cols) + v.set_immutable() + assert v in self.strategy, v + self.strategy[v] = (v, ()) + self.what[v] = ('multizeta', None) + self.level[v] = 0 + + def update_algebraic(self, v, lin_comb): + # take lin_comb of the current columns to generate a new one + f = EvalAlgebraic(self, v, lin_comb) + self.strategy[v] = (f, ()) + self.what[v] = ('algebraic', lin_comb) + self.level[v] = max(self.level[w] for w in f.children()) + 1 + + def algebraic(self, forward_only=False): + r""" + Iterate through available algebraic relations. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.gmzv_generic_reduction import GMZVSolver + sage: S = GMZVSolver(3, 3) + sage: S.update_multizetas() + sage: S + Solve strategy d=3 r=3 (1 known and 7 unknown) + sage: algebraic = list(S.algebraic()) + sage: print([ans[0].lin_prod_string() for ans in algebraic]) + ['(0+1+2)(0)(1)'] + + One can make the solver aware of this relation via:: + + sage: S.update_algebraic(*algebraic[0]) + sage: S + Solve strategy d=3 r=3 (2 known and 6 unknown) + """ + for v in self.strategy: + assert v.ncols() == self.r + if self.strategy[v] is not None: + continue + for c, lin_comb in v.new_columns(forward_only): + unknown = False + for i in range(self.r): + if not lin_comb[i]: + continue + m = v._entries.__copy__() + m[:,i] = c + vv = GeneralizedMultipleZetaFunction.from_matrix(m, True, False, True, True) + if self.strategy[vv] is None: + unknown = True + break + if unknown: + continue + yield v, lin_comb + + def update_forward_stuffle(self, v, stuffle): + d = v.nrows() # number of variables + n = v.ncols() + + # stuffle is a tree: each element is either a LFF or a list + # subset to stuffle, list of triples (rperm, cperm, new node) + comb = [] + Sn = SymmetricGroup(n) + level = 0 + for vv, rp, cp in stuffle_to_lin_comb(v, Sn.one(), Sn.one(), stuffle): + groups = [(cp(i+1) - 1,) for i in range(n)] + if vv.nrows() < d: + # assume .eval is working for less rows + comb.append((1, FunctionWrap(vv, groups))) + else: + level = max(level, self.level[vv]) + comb.append((1, FunctionWrap(self.eval, groups, exp_position=1, extra_args=(vv,)))) + + self.strategy[v] = (LinearCombination(comb), ()) + self.what[v] = ('forward stuffle', stuffle) + self.level[v] = level + 1 + + def _find_iterated_stuffle(self, v): + if self.strategy[v] is not None: + return v + for subset, stuffle in v.stuffles(only_highest_weight=True): + has = True + for i, x in enumerate(stuffle): + y = x.copy() + rperm, cperm = y.canonicalize() + y.set_immutable() + ans = self._find_iterated_stuffle(y) + if not ans: + has = False + break + stuffle[i] = (rperm, cperm, ans) + if has: + return [subset, stuffle] + return False + + def iterated_forward_stuffles(self): + r""" + Return the unknown gmzv that could be reduced via iteration of forward + stuffles to known gmzv. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.gmzv_generic_reduction import GMZVSolver + sage: S = GMZVSolver(3, 3) + sage: S.update_multizetas() + sage: stuffles = list(S.iterated_forward_stuffles()) + sage: list(ans[0].lin_prod_string() for ans in stuffles) + ['(0)(1)(2)', '(0+1)(0+2)(0)', '(0+1)(0)(2)'] + sage: S.update_forward_stuffle(*stuffles[1]) + sage: S + Solve strategy d=3 r=3 (2 known and 6 unknown) + """ + for v in self.strategy: + if self.strategy[v] is not None: + continue + stuffle = self._find_iterated_stuffle(v) + if stuffle: + yield v, stuffle + + def update_multistuffle(self, v, vnsym, lin_comb): + r""" + INPUT: + + - ``v`` + + - ``vnsym`` + + - ``lin_comb`` - list of ``(typ, coeff, vector)`` + """ + D = {} # our linear combination of full data + rem = {} # remainders + level = 0 + for t, coeff, data in lin_comb: + if t == 0: + # by design, the highest length values compensate and we ignore them + vv, subset = data + for vvv in vv.stuffle(subset, only_highest_weight=False, sort_rows=True): + if vvv.nrows() == v.nrows(): + continue + if vvv in D: + D[vvv] += coeff + else: + D[vvv] = coeff + elif t == 1: + # known value + vv = data + if vv in D: + D[vv] += coeff + else: + D[vv] = coeff + + for vv in list(D): + if D[vv] == 0: + del D[vv] + + comb = [] + for vv, coeff in D.items(): + vv = vv.copy() + rp, cp = vv.canonicalize() + vv.set_immutable() + groups = [(cp(i+1)-1,) for i in range(self.r)] + if vv.nrows() < self.d: + # assume .eval is working for less rows + comb.append((coeff, FunctionWrap(vv, groups))) + else: + comb.append((coeff, FunctionWrap(self.eval, groups, exp_position=1, extra_args=(vv,)))) + level = max(level, self.level[vv]) + + self.strategy[v] = (LinearCombination(comb), ()) + self.what[v] = ('multistuffle', lin_comb, D) + self.level[v] = level + 1 + + def multistuffles(self, verbose=True): + # TODO: we should make all matrices sparse + known, M_known = self.known_matrix() + relations, M_rel = self.stuffle_relation_matrix() + M = matrix(ZZ, M_rel + M_known) + M = M.sparse_matrix() + M.set_immutable() + Vnsym, nsym_to_can = self._all_forms() + indicesnsym = {v:i for i,v in enumerate(Vnsym)} + Fnsym = FreeModule(ZZ, len(indicesnsym)) + nrel = len(M_rel) + for i, v in enumerate(self.V): + if self.strategy[v] is not None: + continue + vnsym = next(v.symmetric()) # ? + try: + sol = M.solve_left(Fnsym.gen(indicesnsym[vnsym])) + except ValueError: + sol = None + if sol is not None: + yield v, vnsym, [(0, x, relations[i]) if i < nrel else (1, x, known[i-nrel]) for i,x in enumerate(sol) if x] + + @cached_method + def _all_forms(self): + d = self.d + r = self.r + assert 1 <= r <= d + vectors = list(linear_forms(d)) + + if d == r: + Vnsym = [] + for _, lins in linearly_independent_vectors(vectors, min_size=d, max_size=d): + Vnsym.append(GeneralizedMultipleZetaFunction.from_matrix(lins, False, True, True, True)) + else: + # columns are independent + # rows are sorted + Vnsym = set() + for _, lins in linearly_independent_vectors(vectors, min_size=r, max_size=r): + for p in itertools.permutations(range(r)): + v = GeneralizedMultipleZetaFunction.from_matrix(matrix(ZZ, [lins[p[i]] for i in range(r)]).transpose(), False, True, True, True) + if v.has_no_zero_row_or_column(): + Vnsym.add(v) + Vnsym = list(Vnsym) + + nsym_to_can = {} + for v in Vnsym: + vv = v.copy() + rp, cp = vv.canonicalize() + vv.set_immutable() + nsym_to_can[v] = (vv, rp, cp) + + return Vnsym, nsym_to_can + + @cached_method + def stuffle_relation_matrix(self): + r""" + return relations, mat + """ + Vnsym, nsym_to_can = self._all_forms() + indices = {v: i for i,v in enumerate(Vnsym)} + F = VectorSpace(QQ, len(Vnsym)) + assert len(Vnsym) == len(indices) + + relations = [] + mat = [] + for v in Vnsym: + for subset, stuffle in v.stuffles(only_highest_weight=True, sort_rows=True): + relations.append((v, subset[:])) + mat.append(F.gen(indices[v]) - sum(F.gen(indices[vv]) for vv in stuffle)) + + return relations, mat + + def known_matrix(self): + Vnsym, _ = self._all_forms() + indices = {v: i for i,v in enumerate(Vnsym)} + F = VectorSpace(QQ, len(Vnsym)) + elts = [] + mat = [] + for v, status in self.strategy.items(): + if status is not None: + for vv in v.symmetric(): + mat.append(F.gen(indices[vv])) + elts.append(vv) + return elts, mat + +def print_array(level, explanations): + width = max(len(L.lin_prod_string()) for L in level) + + line = '{L:<{width}} {level:<2} {handling:<10} {neighbors}' + for L in sorted(level, key=lambda x: 100000 if level[x] is None else level[x]): + Ls = '{L:<{width}}'.format(L=L.lin_prod_string(), width=width) + if level[L] is None: + print(line.format(L=Ls, width=width, level='oo', handling='x', neighbors='')) + continue + + handling = explanations[L][0] + N = explanations[L][2] + Ns = ', '.join(vv.lin_prod_string()+':'+str(level[vv]) for vv in N) + print(line.format(L=Ls, width=width, level=level[L], handling=handling, neighbors=Ns)) + + +def print_stuffle3(): + gens = V3gens() + simplices = [c for c in itertools.combinations(gens,3) if matrix(c).rank() == 3] + todo = [] + for s1,s2,s3 in simplices: + if s1[0] + s1[1] <= 1 and s2[0] + s2[1] <= 1 and s3[0] + s3[1] <= 1: + todo.append((s1,s2,s3,0,1,2)) + if s1[0] + s1[2] <= 1 and s2[0] + s2[2] <= 1 and s3[0] + s3[2] <= 1: + todo.append((s1,s2,s3,0,2,1)) + if s1[1] + s1[2] <= 1 and s2[1] + s2[2] <= 1 and s3[1] + s3[2] <= 1: + todo.append((s1,s2,s3,1,2,0)) + + a,b,c = ZZ['a,b,c'].gens() + for s1,s2,s3,i,j,k in todo: + D = ((s1,a),(s2,b),(s3,c)) + + # big subsimplex 1 + m = matrix([s1,s2,s3]).transpose() + m[i] += m[j] + t1,t2,t3 = m.columns() + D1 = clean_term(3, ((t1,a),(t2,b),(t3,c))) + + # big subsimplex 2 + m = matrix([s1,s2,s3]).transpose() + m[j] += m[i] + t1,t2,t3 = m.columns() + D2 = clean_term(3, ((t1,a),(t2,b),(t3,c))) + + # join + m = matrix([s1,s2,s3]).transpose() + m[i] += m[j] + m = m.delete_rows([j]) + t1,t2,t3 = m.columns() + D3 = clean_term(2, ((t1,a),(t2,b),(t3,c))) + + dat = to_Z3(D) + if dat[0] == 0: + continue + mzv = Z3_to_mzv(*dat) + if mzv is not None: + LHS = "ζ{}".fomat(mzv) + else: + LHS = "Z3{}".format(dat) + + dat1 = to_Z3(D1) + mzv1 = Z3_to_mzv(*dat1) + if mzv1 is not None: + RHS1 = "ζ{}".format(mzv1) + else: + RHS1 = "Z3{}".format(dat1) + + dat2 = to_Z3(D2) + mzv2 = Z3_to_mzv(*dat2) + if mzv2 is not None: + RHS2 = "ζ{}".format(mzv2) + else: + RHS2 = "Z3{}".format(dat2) + + dat3 = to_Z2(D3) + mzv3 = Z2_to_mzv(*dat3) + if mzv3 is not None: + RHS3 = "ζ{}".format(mzv3) + else: + RHS3 = "Z2{}".format(dat3) + + print("{} = {} + {} + {}".format(LHS, RHS1, RHS2, RHS3)) + + + +def try_relation(n, den_tuple): + r""" + Assuming that ``den_tuple`` has full rank, make it so that it has only ``n`` columns. + + EXAMPLES:: + + sage: from surface_dynamics.misc.generalized_multiple_zeta_values import try_relation + sage: V = FreeModule(ZZ, 2) + sage: va = V((1,0)); va.set_immutable() + sage: vb = V((0,1)); vb.set_immutable() + sage: vc = V((1,1)); vc.set_immutable() + + sage: den_tuple = ((va,2),(vb,2),(vc,2)) + sage: for r in try_relation(2, den_tuple): + ....: print(r) + """ + # assume it is full rank + if len(den_tuple) <= n: + return + M = matrix([v for v,p in den_tuple]) + for relation in sorted(M.left_kernel().basis(), key=lambda x: sum(bool(cc) for cc in x)): + if sum(x < 0 for x in relation) > sum(x > 0 for x in relation): + relation = -relation + for i in range(len(den_tuple)): + if relation[i] < 0: + yield kill_relation(n, den_tuple, i, relation) + + +def try_algebraic_forward(S): + L = list(S.algebraic(True)) + if not L: + return 0 + S.update_algebraic(*choice(L)) + return 1 + + +def try_algebraic(S): + L = list(S.algebraic(False)) + if not L: + return 0 + S.update_algebraic(*choice(L)) + return 1 + + +def try_forward_stuffle(S): + L = list(S.iterated_forward_stuffles()) + if not L: + return 0 + S.update_forward_stuffle(*choice(L)) + return 1 + +def try_multistuffle(S): + L = list(S.multistuffles()) + if not L: + return 0 + S.update_multistuffle(*choice(L)) + return 1 + + +class Solvers: + @staticmethod + def random(d, r): + import random + S = GMZVSolver(d, r) + S.update_multizetas() + update_functions = [try_algebraic, try_algebraic_forward, try_forward_stuffle, try_multistuffle] + while any(status is None for status in S.strategy.values()): + random.shuffle(update_functions) + for F in update_functions: + status = F(S) + if status: + break + if status == 0: + raise ValueError('not completable strategy') + return S + + @staticmethod + def multistuffle_first(d, r, verbose=False): + r""" + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import gmzv_solvers, GeneralizedMultipleZetaFunction + sage: S = gmzv_solvers.multistuffle_first(3, 3) + sage: l = GeneralizedMultipleZetaFunction([[1,0,0], [0,1,0], [0,0,1]]) + sage: l.set_immutable() + sage: assert S.eval(l, (2,2,2)) == Multizeta(2)**3 + """ + S = GMZVSolver(d, r) + S.update_multizetas() + update_functions = [try_multistuffle, try_algebraic_forward, try_algebraic] + while any(status is None for status in S.strategy.values()): + if verbose: + print(" {}/{}".format(len(S.known()), len(S.strategy))) + for F in update_functions: + res = F(S) + if res: + break + if not res: + print("NOT COMPLETABLE") + return S + return S + + @staticmethod + def reasonable(d, r, verbose=False): + S = GMZVSolver(d, r) + S.update_multizetas() + update_functions = [try_forward_stuffle, try_algebraic_forward, try_multistuffle, try_algebraic] + while any(status is None for status in S.strategy.values()): + if verbose: + print(" {}/{}".format(len(S.known()), len(S.strategy))) + for F in update_functions: + res = F(S) + if res: + break + if not res: + print("NOT COMPLETABLE") + return S + return S + +gmzv_solvers = Solvers() diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py new file mode 100644 index 00000000..29a03578 --- /dev/null +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py @@ -0,0 +1,217 @@ +r""" +Three variables generalized multiple zeta values + +This module is about the sums + +.. MATH:: + + Z3(a, b, c, d, e, f, g) + = + \sum_{x,y,z \geq 1} \frac{1}{x^a y^b z^c (x+y)^d (x+z)^e (y+z)^f (x+y+z)^g} +""" +#***************************************************************************** +# Copyright (C) 2019-2023 Vincent Delecroix <20100.delecroix@gmail.com> +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +#***************************************************************************** + +from sage.rings.rational_field import QQ +from sage.arith.misc import binomial +from sage.modular.multiple_zeta import Multizetas + +from .options import VERBOSE, DIVERGENT_MZV + +def Z3_sort_abc(a,b,c,d,e,f,g): + if b > a: + if c > b: + # c > b > a + a,b,c,d,e,f,g = c,b,a,f,e,d,g + elif c > a: + # b >= c > a + a,b,c,d,e,f,g = b,c,a,f,d,e,g + else: + # b > a >= c + a,b,c,d,e,f,g = b,a,c,d,f,e,g + elif c > a: + # c > a >= b + a,b,c,d,e,f,g = c,a,b,e,f,d,g + elif c > b: + # a >= c >= b + a,b,c,d,e,f,g = a,c,b,e,d,f,g + else: + pass + assert a >= b >= c + return a,b,c,d,e,f,g + + +def is_Z3_convergent(a, b, c, d, e, f, g): + r""" + Return whether `Z3(a,b,c,d,e,f,g)` is convergent. + + TESTS:: + + sage: from surface_dynamics.misc.generalized_multiple_zeta_values import is_Z3_convergent + + Convergent examples:: + + sage: assert is_Z3_convergent(2,0,2,2,0,0,0) + sage: assert is_Z3_convergent(0,0,0,0,0,0,4) + sage: assert is_Z3_convergent(1,0,0,1,0,0,2) + sage: assert is_Z3_convergent(0,1,0,1,0,0,2) + sage: assert is_Z3_convergent(0,1,0,0,0,1,2) + + Divergent examples:: + + sage: assert not is_Z3_convergent(0,0,0,1,1,1,0) + sage: assert not is_Z3_convergent(0,0,0,0,0,0,3) + + """ + from sage.geometry.polyhedron.constructor import Polyhedron + x, y, z = ZZ['x,y,z'].gens() + poly = x**a * y**b * z**c * (x+y)**d * (x+z)**e * (y+z)**f * (x+y+z)**g + newton_polytope = Polyhedron(vertices=poly.exponents(), rays=[(-1,0,0),(0,-1,0),(0,0,-1)]) + V = newton_polytope.intersection(Polyhedron(rays=[(1,1,1)])).vertices() + r = max(max(v.vector()) for v in V) + return r > 1 + + +def Z3(a, b, c, d, e, f, g, check_convergence=True): + r""" + The function ``Z3(a,b,c,d,e,f,g)``. + + The reduction algorithm was designed by Bill Allombert. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import Z3 + + sage: M = Multizetas(QQ) # optional: mzv + + sage: Z3(1,1,1,1,1,1,1) # optional: mzv + 21/2*ζ(1,1,5) + 9/2*ζ(1,2,4) - 3/2*ζ(1,3,3) - 3/2*ζ(1,4,2) + 9/2*ζ(1,6) + sage: Z3(3,0,0,0,0,3,0) # optional: mzv + 6*ζ(1,4) - 12*ζ(1,5) + 3*ζ(2,3) - 6*ζ(2,4) + ζ(3,2) - 2*ζ(3,3) + + sage: assert Z3(2,3,4,0,0,0,0) == M((2,)) * M((3,)) * M((4,)) # optional: mzv + sage: assert Z3(1,0,0,2,0,0,3) == M((1,2,3)) # optional: mzv + + sage: assert Z3(0,0,0,2,0,1,1) == M((4,)) / 2 # optional: mzv + sage: assert Z3(1,0,1,1,0,0,1) == 3 * M((1,1,2)) # optional: mzv + + sage: assert Z3(0,0,0,0,0,0,4) == 1/2 * M((2,)) - 3/2 * M((3,)) + M((4,)) # optional: mzv + + sage: assert Z3(0,0,0,2,0,1,1) == 2 * M((1,1,2)) - M((2,2)) - 3 *M((1,3)) # optional: mzv + """ + M = Multizetas(QQ) + CHECK_CONVERGENCE = False + + # x^a y^b z^c (x+y)^d (x+z)^e (y+z)^f (x+y+z)^g + if VERBOSE: + print("Z3({},{},{},{},{},{},{})".format(a,b,c,d,e,f,g)) + if a < 0 or b < 0 or c < 0 or d < 0 or e < 0 or f < 0 or g < 0: + raise ValueError("invalid exponents for Z3: a={} b={} c={} d={} e={} f={} g={}".format(a,b,c,d,e,f,g)) + + if not DIVERGENT_MZV and CHECK_CONVERGENCE and not is_Z3_convergent(a,b,c,d,e,f,g): + raise DivergentZetaError("divergent Z3({},{},{},{},{},{},{})".format(a,b,c,d,e,f,g)) + + # step 1: try to get rid of the terms (x+y), (x+z), (y+z) + if d and e and f: + if VERBOSE: + print("reduction (x+y+z) = ((x+y) + (x+z) + (y+z)) / 2") + return (Z3(a,b,c,d-1,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c,d,e-1,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c,d,e,f-1,g+1,CHECK_CONVERGENCE)) / 2 + if a and f: + # x^a (y+z)^f -> (x+y+z)^g + # additive version Z3(a-1,b,c,d,e,f,g+1) + Z3(a,b,c,d,e,f-1,g+1) + if VERBOSE: + print("reduction (x+y+z) = (x) + (y+z)") + return sum(binomial(a+k-1, k) * Z3(0,b,c,d,e,f-k,g+a+k,CHECK_CONVERGENCE) for k in range(f)) + \ + sum(binomial(f+k-1, k) * Z3(a-k,b,c,d,e,0,g+f+k,CHECK_CONVERGENCE) for k in range(a)) + if b and e: + # y^b (x+z)^e -> (x+y+z)^g + # additive version Z3(a,b-1,c,d,e,f,g+1) + Z3(a,b,c,d,e-1,f,g+1) + if VERBOSE: + print("reduction (x+y+z) = (y) + (x+z)") + return sum(binomial(b+k-1, k) * Z3(a,0,c,d,e-k,f,g+b+k,CHECK_CONVERGENCE) for k in range(e)) + \ + sum(binomial(e+k-1, k) * Z3(a,b-k,c,d,0,f,g+e+k,CHECK_CONVERGENCE) for k in range(b)) + if c and d: + # z^c (x+y)^d -> (x+y+z)^g + # additive version Z3(a,b,c-1,d,e,f,g+1) + Z3(a,b,c,d-1,e,f,g+1) + if VERBOSE: + print("reduction (x+y+z) = (z) + (x+y)") + return sum(binomial(c+k-1, k) * Z3(a,b,0,d-k,e,f,g+c+k,CHECK_CONVERGENCE) for k in range(d)) + \ + sum(binomial(d+k-1, k) * Z3(a,b,c-k,0,e,f,g+d+k,CHECK_CONVERGENCE) for k in range(c)) + + assert d*e*f == a*f == b*e == c*d == 0 + + a,b,c,d,e,f,g = Z3_sort_abc(a,b,c,d,e,f,g) + + # step 2: kill c + if c: + if VERBOSE: + print("reduction (x+y+z) = (x) + (y) + (z)") + return Z3(a-1,b,c,d,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b-1,c,d,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c-1,d,e,f,g+1,CHECK_CONVERGENCE) + + assert c == 0 + + if b: + # x^a y^b -> (x+y)^d + # additive version Z3(a-1,b,0,d+1,e,f,g) + Z3(a,b-1,0,d+1,e,f,g) + if VERBOSE: + print("reduction (x+y) = (x) + (y)") + return sum(binomial(a+k-1, k) * Z3(0,b-k,c,d+a+k,e,f,g,CHECK_CONVERGENCE) for k in range(b)) + \ + sum(binomial(b+k-1, k) * Z3(a-k,0,c,d+b+k,e,f,g,CHECK_CONVERGENCE) for k in range(a)) + + assert b == c == 0 + + assert b == c == d*e*f == 0 + + if a and f: + raise RuntimeError + + assert b == c == d*e*f == a*f == 0 + + if a == 0: + d,e,f = sorted([d,e,f], reverse=True) + + if a and d and e: + if VERBOSE: + print("reduction (x+y+z) = (x+y) + (x+z) - (x)") + return Z3(a,b,c,d-1,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c,d,e-1,f,g+1,CHECK_CONVERGENCE) - Z3(a-1,b,c,d,e,f,g+1,CHECK_CONVERGENCE) + + assert b == c == d*e*f == a*f == a*d*e == 0 + + # x^a (x+y)^d (x+y+z)^g + if e == f == 0: + from .generalized_multiple_zeta_values import convergent_multizeta + return convergent_multizeta((a,d,g)) + # x^a (x+z)^e (x+y+z)^g + if d == f == 0: + from .generalized_multiple_zeta_values import convergent_multizeta + return convergent_multizeta((a,e,g)) + if a == f == 0: + # (x+y)^d (x+z)^e (x+y+z)^g + if g == 0: + # (x+y)^d (x+z)^e + assert d and e + return M((d-1,e)) - M((d,e)) + M((e-1,d)) - M((e,d)) + M((d+e-1,)) - M((d+e,)) + elif d == e == 1: + if VERBOSE: + print("Bill formulas (x+y) (x+z) (x+y+z)^n") + return sum(M((1,g+1-i,i)) for i in range(2,g+1)) - 3*M((1,g+1)) + elif g == 1: + if VERBOSE: + print("[0,0,0,d,0,f,1] = [f,0,d-1,1,0,0,1] - [f,d,1] - zeta([d,1,f]) - zeta([1+d,f])") + from .gmzv_two_variables import Z2 + f,d = sorted([d,e]) + return Z3(f,0,d-1,1,0,0,1,CHECK_CONVERGENCE) - Z2(f,d,1,CHECK_CONVERGENCE) - M((f,1,d)) - M((f,d+1)) + else: + if VERBOSE: + print("[0,0,0,0,e,f,g] = [e,0,f,g,0,0,0] - [e,0,f,0,0,0,g] - Z(e,g,f) - Z(e,f+g) - Z(e+f,g)") + from .gmzv_two_variables import Z2 + f,d = sorted([d,e]) + return Z3(f,0,d,g,0,0,0,CHECK_CONVERGENCE) - Z3(f,d,0,0,0,0,g,CHECK_CONVERGENCE) - Z2(f,d,g,CHECK_CONVERGENCE) - M((f,g,d)) - M((f,g+d)) + + raise RuntimeError diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py new file mode 100644 index 00000000..f6c7b123 --- /dev/null +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py @@ -0,0 +1,59 @@ +r""" +Two variables generalized multiple zeta values + +This module is about the sums + +.. MATH:: + + Z2(a, b, c) = sum_{x, y} 1 / (x^a y^b (x + y)^c) +""" +#***************************************************************************** +# Copyright (C) 2019-2023 Vincent Delecroix <20100.delecroix@gmail.com> +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +#***************************************************************************** + + +from sage.rings.rational_field import QQ +from sage.arith.misc import binomial +from sage.modular.multiple_zeta import Multizetas + +from .options import VERBOSE, DIVERGENT_MZV + +def is_Z2_convergent(a, b, c): + r""" + Return whether Z2(a, b, c) is convergent. + + TESTS:: + + sage: from surface_dynamics.generalized_multiple_zeta_values import is_Z2_convergent + + Convergent examples:: + + sage: assert is_Z2_convergent(1,1,1) + sage: assert is_Z2_convergent(2,2,0) + sage: assert is_Z2_convergent(0,0,3) + + Divergent examples:: + + sage: assert not is_Z2_convergent(0,0,2) + sage: assert not is_Z2_convergent(1,2,0) + """ + from sage.geometry.polyhedron.constructor import Polyhedron + x, y = ZZ['x,y'].gens() + poly = x**a * y**b * (x+y)**c + newton_polytope = Polyhedron(vertices=poly.exponents(), rays=[(-1,0),(0,-1)]) + V = newton_polytope.intersection(Polyhedron(rays=[(1,1)])).vertices() + r = max(max(v.vector()) for v in V) + return r > 1 + + +def Z2(a, b, c, check_convergence=True): + M = Multizetas(QQ) + if a == 0 and b == 0: + return M((c-1,)) - M((c,)) + else: + return sum(binomial(a+i-1,i) * M((b-i,c+a+i)) for i in range(b)) + sum(binomial(b+i-1,i) * M((a-i,c+b+i)) for i in range(a)) diff --git a/surface_dynamics/generalized_multiple_zeta_values/options.py b/surface_dynamics/generalized_multiple_zeta_values/options.py new file mode 100644 index 00000000..633ecb79 --- /dev/null +++ b/surface_dynamics/generalized_multiple_zeta_values/options.py @@ -0,0 +1,4 @@ +VERBOSE = False +DIVERGENT_MZV = True + + diff --git a/surface_dynamics/misc/additive_multivariate_generating_series.py b/surface_dynamics/misc/additive_multivariate_generating_series.py index ab46cc51..2368d29d 100644 --- a/surface_dynamics/misc/additive_multivariate_generating_series.py +++ b/surface_dynamics/misc/additive_multivariate_generating_series.py @@ -58,7 +58,7 @@ def integral_sum_as_mzv(self): sage: f.integral_sum_as_mzv() # optional - mzv 36*ζ(1,1,4) + 24*ζ(1,2,3) + 12*ζ(1,3,2) + 12*ζ(2,1,3) + 6*ζ(2,2,2) """ - from .generalized_multiple_zeta_values import handle_term + from surface_dynamics.generalized_multiple_zeta_values.generaralized_multiple_zeta_values import handle_term n = self.parent().free_module().rank() return sum(QQ(num) * handle_term(n, den._tuple) for den, num in self._data.items()) diff --git a/surface_dynamics/misc/divergent_multiple_zeta_values.py b/surface_dynamics/misc/divergent_multiple_zeta_values.py new file mode 100644 index 00000000..96da9eb6 --- /dev/null +++ b/surface_dynamics/misc/divergent_multiple_zeta_values.py @@ -0,0 +1,128 @@ +from sage.modular.multiple_zeta import Multizetas + +M = Multizetas(QQ) + +# TODO: without this, it is not possible to create polynomials with Multizeta coefficients +M.is_prime_field = lambda: False + +MT = M['T'] # shuffle regularization +MX = M['X'] # stuffle regularization +T = MT.gen() +X = MX.gen() + +# concatenation product +def prefix(M, p, g): + p = M.basis().keys()(p) + ans = M.zero() + for u, coeff in g.monomial_coefficients().items(): + ans += coeff * M.term(p + u) + return ans + +def suffix(M, f, s): + s = M.basis().keys()(s) + ans = M.zero() + for u , coeff in f.monomial_coefficients().items(): + ans += coeff * M.term(u + s) + return ans + +def concatenate(M, left, right): + r""" + Concatenation product in the free algebra + + EXAMPLES:: + + sage: concatenate(M, M((2,)) + M((3,)), M((5,)) + M((7,))) + ζ(2,5) + ζ(2,7) + ζ(3,5) + ζ(3,7) + """ + ans = M.zero() + for u1, c1 in left.monomial_coefficients(): + for u2, c2 in right.monomial_coefficients(): + ans += c1 * c2 * M.term(u1 + u2) + return ans + +def stuffle_on_basis(M, w1, w2): + r""" + Return the stuffle product of two words + + EXAMPLES:: + + sage: stuffle_on_basis(M, Word([2]), Word([5])) + ζ(2,5) + ζ(5,2) + ζ(7) + sage: stuffle_on_basis(M, Word([2,2]), Word([3,4])) + ζ(2,2,3,4) + ζ(2,3,2,4) + ζ(2,3,4,2) + ζ(2,3,6) + ζ(2,5,4) + ζ(3,2,2,4) + ζ(3,2,4,2) + ζ(3,2,6) + ζ(3,4,2,2) + ζ(3,6,2) + ζ(5,2,4) + ζ(5,4,2) + ζ(5,6) + """ + W = M.basis().keys() + if not w1: + return M.term(W(w2)) + if not w2: + return M.term(W(w1)) + y1 = w1[0] + u1 = w1[1:] + y2 = w2[0] + u2 = w2[1:] + if y1 == y2: + return prefix(M, W([y1]), stuffle_on_basis(M, u1, w2) + stuffle_on_basis(M, w1, u2)) + \ + prefix(M, W([y1+y2]), stuffle_on_basis(M, u1, u2)) + else: + return prefix(M, W([y1]), stuffle_on_basis(M, u1, w2)) + \ + prefix(M, W([y2]), stuffle_on_basis(M, w1, u2)) + \ + prefix(M, W([y1+y2]), stuffle_on_basis(M, u1, u2)) + +def stuffle_product(M, left, right): + ans = M.zero() + for u1, c1 in left.monomial_coefficients().items(): + for u2, c2 in right.monomial_coefficients().items(): + ans += c1 * c2 * stuffle_on_basis(M, u1, u2) + return ans + +def stuffle_power(M, a, n): + if n < 0: + raise ValueError + if n == 0: + return M.one() + apow = a + while not (n & 1): + apow = stuffle_product(M, apow, apow) + n >>= 1 + + res = apow + n >>= 1 + while n: + apow = stuffle_product(M, apow, apow) + if n & 1: + res = stuffle_product(M, apow, res) + n >>= 1 + + return res + + +def stuffle_regularization(s): + r""" + Given a multizeta index s returns it as an expression in MX + + EXAMPLES:: + + sage: divergent_stuffle_regularization((1,)) + sage: divergent_stuffle_regularization((2,)) + sage: divergent_stuffle_regularization((2,1)) + sage: divergent_stuffle_regularization((1,1)) + sage: divergent_stuffle_regularization((2,1,1)) + sage: divergent_stuffle_regularization((1,1,1)) + """ + pass + +def shuffle_regularization(s): + r""" + Given a multizeta index s returns it as an expression in MT + """ + pass + +def shuffle_to_harmonic(n): + r""" + Image of X^n under the morphism rho: M[X] -> M[T] of Ihara-Kaneko-Zagier + """ + +def harmonic_to_shuffle(n): + r""" + Image of T^n under the morphism rho^{-1}: M[T] -> M[X] of Ihara-Kaneko-Zagier + """ diff --git a/surface_dynamics/misc/factored_denominator.py b/surface_dynamics/misc/factored_denominator.py index 6f5264ad..a1ba95fc 100644 --- a/surface_dynamics/misc/factored_denominator.py +++ b/surface_dynamics/misc/factored_denominator.py @@ -78,6 +78,23 @@ def vector_to_monomial_string(u, var_names): return '*'.join(s) if s else '1' +def vector_to_monomial_latex(u, var_names): + s = [] + for i,j in enumerate(u): + if j: + if isinstance(var_names, str): + var = '%s%d' % (var_names, i) + else: + var = var_names[i] + + if j == 1: + s.appen(var) + else: + s.append('%s^{%d}' % (var, j)) + + return ' '.join(s) if s else '1' + + def vector_to_linear_form_string(u, var_names): r""" EXAMPLES:: diff --git a/surface_dynamics/misc/factored_denominator.py.orig b/surface_dynamics/misc/factored_denominator.py.orig new file mode 100644 index 00000000..74d2abd9 --- /dev/null +++ b/surface_dynamics/misc/factored_denominator.py.orig @@ -0,0 +1,1121 @@ +r""" +Factored denominators and associated free modules +""" +# **************************************************************************** +# Copyright (C) 2019 Vincent Delecroix <20100.delecroix@gmail.com> +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + +from __future__ import absolute_import, print_function +from six.moves import range, zip + +from six import iteritems + +from sage.misc.cachefunc import cached_method +from sage.structure.element import Element +from sage.structure.parent import Parent +from sage.structure.richcmp import op_LT, op_EQ, op_NE, op_LE, op_GE, op_GT +from sage.structure.unique_representation import UniqueRepresentation +from sage.categories.all import Rings +from sage.rings.all import ZZ, QQ +from sage.modules.vector_integer_dense import Vector_integer_dense + + +def laurent_monomial(R, arg): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import laurent_monomial + sage: R = LaurentPolynomialRing(QQ, 'x', 2) + sage: laurent_monomial(R, (1,2)) + x0*x1^2 + sage: laurent_monomial(R, vector((-1r,3r))) + x0^-1*x1^3 + sage: laurent_monomial(R, (-1,3)) * laurent_monomial(R, (1,1)) == laurent_monomial(R, (0,4)) + True + """ + # only in beta 8.1 versions of Sage + # return R.monomial(*arg) + arg = tuple(arg) + if len(arg) != R.ngens(): + raise TypeError("tuple key must have same length as ngens") + + from sage.misc.misc_c import prod + return prod(x**int(i) for (x, i) in zip(R.gens(), arg)) + + +def vector_to_monomial_string(u, var_names): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import vector_to_monomial_string + sage: vector_to_monomial_string((0,4,3), 'x') + 'x1^4*x2^3' + sage: vector_to_monomial_string((0,4,3), 'hello') + 'hello1^4*hello2^3' + sage: vector_to_monomial_string((0,4,3), ('x','y','z')) + 'y^4*z^3' + sage: vector_to_monomial_string((1,0,-1), 'x') + 'x0*x2^-1' + """ + s = [] + for i, j in enumerate(u): + if j: + if isinstance(var_names, str): + var = '%s%d' % (var_names, i) + else: + var = var_names[i] + + if j == 1: + s.append(var) + else: + s.append('%s^%d' % (var, j)) + + return '*'.join(s) if s else '1' + +<<<<<<< HEAD +======= +def vector_to_monomial_latex(u, var_names): + s = [] + for i,j in enumerate(u): + if j: + if isinstance(var_names, str): + var = '%s%d' % (var_names, i) + else: + var = var_names[i] + + if j == 1: + s.appen(var) + else: + s.append('%s^{%d}' % (var, j)) + + return ' '.join(s) if s else '1' +>>>>>>> cd714fd (more latexization) + +def vector_to_linear_form_string(u, var_names): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import vector_to_linear_form_string + sage: vector_to_linear_form_string((0,4,3), 'x') + '4*x1 + 3*x2' + sage: vector_to_linear_form_string((0,4,3), 'hello') + '4*hello1 + 3*hello2' + sage: vector_to_linear_form_string((0,4,3), ('x','y','z')) + '4*y + 3*z' + sage: vector_to_linear_form_string((1,0,-1), 'x') + 'x0 - x2' + sage: vector_to_linear_form_string((2,0,-3), 'x') + '2*x0 - 3*x2' + sage: vector_to_linear_form_string((-1,1,0), 'x') + '-x0 + x1' + sage: vector_to_linear_form_string((-2,3,0), 'x') + '-2*x0 + 3*x1' + sage: vector_to_linear_form_string((0,-1,1,0), 'x') + '-x1 + x2' + """ + s = '' + first = True + for i, j in enumerate(u): + if j: + if isinstance(var_names, str): + var = '%s%d' % (var_names, i) + else: + var = var_names[i] + + if j == 1: + if first: + s += var + else: + s += ' + %s' % var + elif j == -1: + if first: + s += '-%s' % var + else: + s += ' - %s' % var + elif first: + s += '%d*%s' % (j, var) + elif j < 0: + s += ' - %d*%s' % (-j, var) + else: + s += ' + %d*%s' % (j, var) + first = False + + return '0' if not s else s + + +# NOTE: should this be an instance of Factorization? +class FactoredDenominator(object): + r""" + Factored denominator + + This class is a simple datastructure to handle factored denominator, that is + a list of pairs ``(v, m)`` where the ``v`` are integer vectors (of fixed + dimension) and ``m`` are multiplicities (i.e. positive integers). + + It is used for at least two purposes: + + - (Factored) product of polynomials of the form `(1 - m)^d` where `m` is a + monomial + + - generalized multiple zeta values, where denominators are products of + linear forms + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ**3 + sage: f1 = FactoredDenominator([((1,0,0), 2)], V) + sage: f2 = FactoredDenominator([((0,1,2), 3), ((1,1,1), 1)], V) + sage: f3 = FactoredDenominator([((0,-1,2), 1), ((1,0,0), 1), ((0,0,2), 1)], V) + sage: f1 + {(1, 0, 0): 2} + sage: f1 * f2 * f3 + {(1, 0, 0): 3, (0, 1, 2): 3, (1, 1, 1): 1, (0, -1, 2): 1, (0, 0, 2): 1} + sage: hash(f1) # random + 9164823176457928306 + + sage: FactoredDenominator(f1) + {(1, 0, 0): 2} + """ + __slots__ = ['_dict', '_tuple'] + + def __init__(self, data, V=None): + if isinstance(data, (tuple, list)): + if V is None: + raise ValueError("a ZZ-free module V must be provided") + d = V.dimension() + self._dict = {} + for mon, mult in data: + mon = V(mon) + if mon.is_zero(): + raise ValueError('zero in denominator') + for i in range(d): + if mon[i]: + break + mult = ZZ(mult) + mon.set_immutable() + self._dict[mon] = mult + + elif isinstance(data, dict): + if V is not None: + self._dict = {} + for k, v in data.items(): + k = V(k) + k.set_immutable() + self._dict[k] = v + else: + self._dict = data + + elif isinstance(data, FactoredDenominator): + self._dict = data._dict.copy() + self._tuple = data._tuple[:] + return + + elif isinstance(data, Vector_integer_dense): + data = V(data) + data.set_immutable() + self._dict = {data: ZZ.one()} + self._tuple = ((data, ZZ.one()),) + return + + else: + raise TypeError('invalid data of type {} to initialized a FactoredDenominator'.format(type(data))) + + self._tuple = tuple(sorted(self._dict.items())) + + def __len__(self): + r""" + Return the number of factors (without multiplicities). + """ + return len(self._tuple) + + def __iter__(self): + r""" + Iterates through the pairs ``(monomial, exponent)``. + """ + return iter(self._tuple) + + def subs(self, m): + r""" + Matrix substitution. + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + sage: V = ZZ**3 + sage: f = FactoredDenominator([((0,1,2), 3), ((1,1,1), 1)], V) + + sage: m = matrix(3, [1,1,0,0,1,1,1,0,1]) + sage: f.subs(m) == FactoredDenominator([((1,3,2), 3), ((2,2,2), 1)], V) + True + + sage: m = matrix(3, [1,1,1,1,1,1,1,1,1]) + sage: f.subs(m) + {(3, 3, 3): 4} + + sage: m = matrix(3, [1,0,0,1,0,0,1,0,0]) + sage: f.subs(m) + Traceback (most recent call last): + ... + ValueError: zero denominator + """ + if not self._tuple: + return self + new_dict = {} + for mon, mult in self._dict.items(): + new_mon = m * mon + new_mon.set_immutable() + if not new_mon: + raise ValueError('zero denominator') + if new_mon in new_dict: + new_dict[new_mon] += mult + else: + new_dict[new_mon] = mult + + return FactoredDenominator(new_dict) + + def to_multiplicative_polynomial(self, S, extra_var=False): + r""" + Return the product of the term in a given polynomial ring ``S``. + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ ** 3 + sage: f = FactoredDenominator([((1,0,0), 2)], V) + sage: g = FactoredDenominator([((1,0,0), 2), ((1,1,1),1)], V) + + sage: R1 = QQ['x,y,z'] + sage: f.to_multiplicative_polynomial(R1) + x^2 - 2*x + 1 + sage: g.to_multiplicative_polynomial(R1) + -x^3*y*z + 2*x^2*y*z - x*y*z + x^2 - 2*x + 1 + + sage: f.to_multiplicative_polynomial(R1['E'], extra_var=True) + x^2*E^2 - 2*x*E + 1 + sage: g.to_multiplicative_polynomial(R1['E'], extra_var=True) + -x^3*y*z*E^5 + 2*x^2*y*z*E^4 - x*y*z*E^3 + x^2*E^2 - 2*x*E + 1 + + sage: g.to_multiplicative_polynomial(QQ['t1,t2,t3']) + -t1^3*t2*t3 + 2*t1^2*t2*t3 - t1*t2*t3 + t1^2 - 2*t1 + 1 + """ + if extra_var: + R = S.base_ring() + T = S.gen() + else: + R = S + T = S.one() + + ans = S.one() + for a, i in self._tuple: + ans *= (S.one() - R.monomial(*a) * T ** sum(a)) ** i + return ans + + def to_additive_polynomial(self, S, extra_var=False): + r""" + Return the product of the term in a given polynomial ring ``S``. + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ ** 3 + sage: f = FactoredDenominator([((1,0,0), 2)], V) + sage: g = FactoredDenominator([((1,0,0), 2), ((1,1,1),1)], V) + + sage: R1 = QQ['x,y,z'] + sage: f.to_additive_polynomial(R1) + x^2 + sage: g.to_additive_polynomial(R1) + x^3 + x^2*y + x^2*z + + sage: f.to_additive_polynomial(R1['E'], extra_var=True) + x^2*E^2 + sage: g.to_additive_polynomial(R1['E'], extra_var=True) + (x^3 + x^2*y + x^2*z)*E^3 + + sage: g.to_additive_polynomial(QQ['t1,t2,t3']) + t1^3 + t1^2*t2 + t1^2*t3 + """ + if extra_var: + R = S.base_ring() + T = S.gen() + else: + R = S + T = S.one() + + ans = S.one() + for a, mult in self._tuple: + monomial = sum(coeff * R.gen(i) for i, coeff in enumerate(a)) + ans *= (monomial * T) ** mult + return ans + + def degree(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ**3 + sage: FactoredDenominator([], V).degree() + 0 + sage: FactoredDenominator([((1,0,0), 2)], V).degree() + 2 + sage: FactoredDenominator([((0,1,2), 3), ((1,1,1), 1)], V).degree() + 4 + sage: FactoredDenominator([((0,-1,2), 1), ((1,0,0), 1), ((0,0,2), 1)], V).degree() + 3 + + TESTS:: + + sage: parent(FactoredDenominator([], V).degree()) + Integer Ring + sage: parent(FactoredDenominator([((1,0,0), 2)], V).degree()) + Integer Ring + """ + if not self._tuple: + return ZZ.zero() + return sum(m for _, m in self._tuple) + + # TODO + # this method does not really make sense at this level of generality! + # this can be defined on denominators and extended by linearity + # we should just provide a method to do this in an abstract parent algebra + # + # _image_from_image_of_denominators(function, elt) + # return sum(num * function(den) for num,den in elt) + def logarithmic_minus_derivative(self, j): + r""" + Minus the logarithmic derivative -v' / v taken with respect to the ``j``-th variable + + Since this denominator has the form 1 / (1 - m1) (1 - m2) ... (1 - mk) the + logarithmic derivative is simply + + m1'/(1-m1) + m2'/(1-m2) + ... + mk'/(1-mk) + + INPUT: + + - ``j`` -- index of a variable + + OUTPUT: iterator of triples ``(mult, v, monomial)`` + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ**3 + + sage: f = FactoredDenominator([((1,0,0), 2)], V) + sage: f + {(1, 0, 0): 2} + sage: list(f.logarithmic_minus_derivative(0)) + [(2, (0, 0, 0), (1, 0, 0))] + sage: list(f.logarithmic_minus_derivative(1)) + [] + sage: list(f.logarithmic_minus_derivative(2)) + [] + + sage: f = FactoredDenominator([((1,1,1), 1)], V) + sage: list(f.logarithmic_minus_derivative(0)) + [(1, (0, 1, 1), (1, 1, 1))] + + sage: f = FactoredDenominator([((1,0,0), 1), ((0,1,0), 1), ((0,0,1), 1)], V) + sage: list(f.logarithmic_minus_derivative(0)) + [(1, (0, 0, 0), (1, 0, 0))] + + sage: f = FactoredDenominator([((1,0,0), 2), ((1,1,0), 3), ((1,1,1), 1)], V) + sage: f + {(1, 0, 0): 2, (1, 1, 0): 3, (1, 1, 1): 1} + sage: list(f.logarithmic_minus_derivative(0)) + [(2, (0, 0, 0), (1, 0, 0)), + (3, (0, 1, 0), (1, 1, 0)), + (1, (0, 1, 1), (1, 1, 1))] + """ + for mon, mult in self._tuple: + if mon[j]: + v = mon.__copy__() + mult *= v[j] + v[j] -= 1 + yield (mult, v, mon) + + def __eq__(self, other): + r""" + Equality test + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ**3 + + sage: f1 = FactoredDenominator([((1,0,0), 2)], V) + sage: f2 = FactoredDenominator([((1,0,0), 2)], V) + sage: f3 = FactoredDenominator([((1,0,0), 1)], V) + sage: f4 = FactoredDenominator([((1,0,1), 2)], V) + + sage: f1 == f1 and f1 == f2 + True + sage: f1 == f3 or f1 == f4 + False + """ + if type(self) is not type(other): + raise TypeError + return self._tuple == other._tuple + + def __ne__(self, other): + if type(self) is not type(other): + raise TypeError + return self._tuple != other._tuple + + def __lt__(self, other): + if type(self) is not type(other): + raise TypeError + return self._tuple < other._tuple + + def __le__(self, other): + if type(self) is not type(other): + raise TypeError + return self._tuple <= other._tuple + + def __gt__(self, other): + if type(self) is not type(other): + raise TypeError + return self._tuple > other._tuple + + def __ge__(self, other): + if type(self) is not type(other): + raise TypeError + return self._tuple >= other._tuple + + def copy(self): + res = FactoredDenominator.__new__(FactoredDenominator) + res._dict = self._dict.copy() + res._tuple = self._tuple + return res + + def __mul__(self, other): + r""" + Multiplication + + TESTS:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + sage: V = ZZ**3 + + sage: f = FactoredDenominator([((1,0,0), 2)], V) + """ + if type(self) is not type(other): + raise TypeError + + new_data = self._dict.copy() + for i, j in other._dict.items(): + if i in new_data: + new_data[i] += j + if new_data[i].is_zero(): + del new_data[i] + else: + new_data[i] = j + + return FactoredDenominator(new_data, None) + + def __truediv__(self, other): + r""" + (Partial) division. + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + sage: V = ZZ**3 + + sage: f1 = FactoredDenominator([((1,0,0), 3)], V) + sage: f2 = FactoredDenominator([((1,0,0), 1)], V) + sage: f1/f2 + {(1, 0, 0): 2} + """ + if type(self) != type(other): + raise TypeError + + sd = self._dict + od = other._dict + nd = sd.copy() + for i, j in od.items(): + jj = sd.get(i, -1) + if j > jj: + raise ArithmeticError + elif j == jj: + del nd[i] + else: + nd[i] -= j + + return FactoredDenominator(nd, None) + + __div__ = __truediv__ + + def __nonzero__(self): + return bool(self._dict) + + def is_one(self): + return not self._dict + + def str(self): + return str(self._dict) + + def __repr__(self): + return repr(self._dict) + + def __hash__(self): + return hash(self._tuple) + + def lcm_update(self, other): + if type(self) is not type(other): + raise TypeError + + sd = self._dict + od = other._dict + for i, j in od.items(): + if j > sd.get(i, -1): + sd[i] = j + + self._tuple = tuple(sorted(self._dict.items())) + + def lcm(self, other): + r""" + Return the lcm of two factored denominators. + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ**3 + + sage: f1 = FactoredDenominator([((1,0,0), 2),((0,1,0),1)], V) + sage: f2 = FactoredDenominator([((1,0,0), 1),((0,1,0),2)], V) + sage: f3 = FactoredDenominator([((1,0,0), 1),((0,1,0),1),((0,0,1),1)], V) + sage: f4 = FactoredDenominator([], V) + + sage: f1.lcm(f2) + {(1, 0, 0): 2, (0, 1, 0): 2} + sage: f1.lcm(f3) + {(1, 0, 0): 2, (0, 1, 0): 1, (0, 0, 1): 1} + sage: f2.lcm(f3) + {(1, 0, 0): 1, (0, 1, 0): 2, (0, 0, 1): 1} + + sage: f1.lcm(f4) == f1 and f4.lcm(f1) == f1 + True + """ + res = self.copy() + res.lcm_update(other) + return res + + def gcd_update(self, other): + if type(self) is not type(other): + raise TypeError + + sd = self._dict + od = other._dict + for i, j in od.items(): + if j < sd.get(i, j): + sd[i] = j + + self._tuple = tuple(sorted(self._dict.items())) + + def gcd(self, other): + r""" + Return the gcd of two factored denominator. + + EXAMPLES:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + + sage: V = ZZ**3 + + sage: f1 = FactoredDenominator([((1,0,0), 2),((0,1,0),1)], V) + sage: f2 = FactoredDenominator([((1,0,0), 1),((0,1,0),2)], V) + sage: f3 = FactoredDenominator([((1,0,0), 1),((0,1,0),1),((0,0,1),1)], V) + sage: f4 = FactoredDenominator([], V) + + sage: f1.gcd(f2) + {(1, 0, 0): 1, (0, 1, 0): 1} + sage: f1.gcd(f2) == f1.gcd(f3) == f2.gcd(f3) + True + + sage: f4.gcd(f1) == f4 and f1.gcd(f4) == f1 + True + """ + res = self.copy() + res.gcd_update(other) + return res + + # TODO + # this does not make much sense here + def str_monomials(self, var_names='x'): + terms = [] + for mon, mul in self._tuple: + term = '(1 - %s)' % vector_to_monomial_string(mon, var_names) + if not mul.is_one(): + term += '^%d' % mul + terms.append(term) + return '*'.join(terms) + + # TODO + # this does not make much sense here + def str_linear(self, var_names='h'): + terms = [] + for mon, mul in self._tuple: + term = '(' + vector_to_linear_form_string(mon, var_names) + ')' + if not mul.is_one(): + term += '^%d' % mul + terms.append(term) + + return '*'.join(terms) + + +# TODO: make it possible to use Laurent polynomials in denominator +# TODO: monomial substitution +# TODO: coefficient expansion Verdoolaege-Woods +# 1 - monomial versus linear form +class AbstractMSum(Element): + r""" + An abstract multiple sum + """ + def __init__(self, parent, data, allow_multiple=False, check=True): + r""" + _data is a dictionary: {denominator: polynomial} + """ + Element.__init__(self, parent) + V = self.parent().free_module() + R = self.parent().polynomial_ring() + self._data = {} + + if not data: + return + + if isinstance(data, (tuple, list)): + data = iter(data) + + elif isinstance(data, dict): + if not check: + self._data = data + return + data = list(data.items()) + + elif isinstance(data, AbstractMSum): + self._data = data.copy() + return + + for term in data: + if not isinstance(term, (tuple, list)) or len(term) != 2: + raise ValueError + den, num = term + + num = R(num) + if num.is_zero(): + continue + + if not isinstance(den, FactoredDenominator): + den = FactoredDenominator(den, V) + + if den in self._data: + if not allow_multiple: + raise ValueError('multiple times the same denominator in the input') + else: + self._data[den] += num + else: + self._data[den] = num + + def _den_str(self, den): + return str(den) + + def _repr_(self): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: R = M.polynomial_ring() + sage: x0, x1 = R.gens() + sage: M.term(-x0 + x1^3, [((1,1), 1), ((1,2),2)]) + (x1^3 - x0)/((1 - x0*x1)*(1 - x0*x1^2)^2) + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('t,u') + sage: t, u = M.polynomial_ring().gens() + sage: M.term(-t + u^3, [((1,1), 1), ((1,2),2)]) + (u^3 - t)/((1 - t*u)*(1 - t*u^2)^2) + """ + if not self._data: + return '0' + + terms = [] + for den in sorted(self._data): + num = self._data[den] + fraction = '(' + str(num) + ')' + if not den.is_one(): + fraction += '/' + fraction += '(' + self._den_str(den) + ')' + + terms.append(fraction) + + return ' + '.join(terms) + + def is_trivial_zero(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: M.term(0, []).is_trivial_zero() + True + sage: M.term(0, [((1,1),1)]).is_trivial_zero() + True + sage: M.term(1, []).is_trivial_zero() + False + """ + return not self._data + + def is_trivial_one(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: M.term(0, []).is_trivial_one() + False + sage: M.term(0, [((1,1),1)]).is_trivial_one() + False + sage: M.term(1, []).is_trivial_one() + True + sage: M.term(1, [((1,1),1)]).is_trivial_one() + False + + sage: f = M.term(1, [((1,0),1)]) + M.term(1, [((0,1),1)]) + sage: f.is_trivial_one() + False + """ + if len(self._data) != 1: + return False + (den, num) = next(iteritems(self._data)) + return num.is_one() and den.is_one() + + def _add_(self, other): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 3) + sage: m1 = M.term(1, [((1,1,0),1)]) + sage: m2 = M.term(-2, [((1,0,0),1),((0,1,0),2)]) + sage: m1 + m2 + (-2)/((1 - x1)^2*(1 - x0)) + (1)/((1 - x0*x1)) + """ + if self.is_trivial_zero(): + return other + if other.is_trivial_zero(): + return self + + ans_data = self._data.copy() + for den, num in other._data.items(): + if den in ans_data: + ans_data[den] += num + if ans_data[den].is_zero(): + del ans_data[den] + else: + ans_data[den] = num + + R = self.parent() + return R.element_class(R, ans_data) + + def _mul_(self, other): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 3) + sage: m1 = M.term(1, [((1,1,0),1)]) + sage: m2 = M.term(-2, [((1,0,0),1),((0,1,0),2)]) + sage: m1 * m2 + (-2)/((1 - x1)^2*(1 - x0)*(1 - x0*x1)) + + sage: R = M.polynomial_ring() + sage: M.one() * R.one() + (1) + + sage: R = M.polynomial_ring() + sage: x0, x1, x2 = R.gens() + sage: m1 = M.term(x0 - 2*x2 + 1, [((1,1,1),2)]) + sage: m2 = M.term(x0*x1 + 1, [((0,1,2),1)]) + sage: m3 = M.term(1 + x0 + x1, [((1,0,1),1)]) + sage: m4 = M.term(x0 - 1, [((1,0,0),2),((1,0,1),1)]) + sage: (m1 + m2) * (m3 + m4) - m1 * m3 - m2 * m3 - m1 * m4 - m2 * m4 + 0 + """ + ans_data = {} + for den1, num1 in self._data.items(): + for den2, num2 in other._data.items(): + den = den1 * den2 + num = num1 * num2 + if den in ans_data: + ans_data[den] += num + if not ans_data[den]: + del ans_data[den] + else: + ans_data[den] = num + + M = self.parent() + return M.element_class(M, ans_data) + + def __neg__(self): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 3) + sage: m1 = M.term(1, [((1,1,0),1)]) + sage: m1 + (1)/((1 - x0*x1)) + sage: -m1 + (-1)/((1 - x0*x1)) + sage: m1 - m1 + 0 + """ + ans_data = {} + for den, num in self._data.items(): + ans_data[den] = -num + M = self.parent() + return M.element_class(M, ans_data) + + def _richcmp_(self, other, op): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: x, y = M.polynomial_ring().gens() + sage: M.term(x*y+1, [((1,1),1)]) == M.term(x*y+1, [((1,1),1)]) + True + """ + if op != op_EQ and op != op_NE: + raise TypeError('no comparison available for multivariate series') + + if self._data == other._data: + return op == op_EQ + + difference = self - other + return bool(difference.factor()._data) == (op == op_NE) + + def degrees(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: from surface_dynamics.misc.additive_multivariate_generating_series import AdditiveMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: f = M.term(1, [((1,0), 1), ((1,1),2)]) + M.term(2, [((1,0), 2), ((1,2),3)]) + sage: f.degrees() + {3, 5} + + sage: A = AdditiveMultivariateGeneratingSeriesRing('x', 2) + sage: f = A.term(1, [((1,0), 1), ((1,1),2)]) + A.term(2, [((1,0), 2), ((1,2),3)]) + sage: f.degrees() + {3, 5} + """ + res = set() + for den, num in self._data.items(): + if num(*self.parent()._critical_point()).is_zero(): + raise ValueError('singular numerator') + res.add(den.degree()) + return res + + +# TODO: the UniqueRepresentation should normalize the Laurent polynomial +class AbstractMSumRing(UniqueRepresentation, Parent): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: MultiplicativeMultivariateGeneratingSeriesRing(5, 'x') is MultiplicativeMultivariateGeneratingSeriesRing(['x0','x1','x2','x3','x4']) + True + """ + def __init__(self, poly_ring): + from sage.modules.free_module import FreeModule + self._polynomial_ring = poly_ring + dim = ZZ(poly_ring.ngens()) + self._free_module = FreeModule(ZZ, dim) + + # univariate extension of the polynomial ring + # (needed in several algorithms) + self._polynomial_ring_extra_var = self._polynomial_ring['EXTRA_VAR'] + + Parent.__init__(self, category=Rings(), base=poly_ring.base_ring()) + + def _critical_point(self): + raise NotImplementedError + + def ngens(self): + r""" + Return the number of generators. + + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: from surface_dynamics.misc.additive_multivariate_generating_series import AdditiveMultivariateGeneratingSeriesRing + + sage: MultiplicativeMultivariateGeneratingSeriesRing('x', 3).ngens() + 3 + sage: AdditiveMultivariateGeneratingSeriesRing('x', 3).ngens() + 3 + """ + return self.polynomial_ring().ngens() + + def free_module(self): + return self._free_module + + def polynomial_ring(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: from surface_dynamics.misc.additive_multivariate_generating_series import AdditiveMultivariateGeneratingSeriesRing + + + sage: MultiplicativeMultivariateGeneratingSeriesRing('x', 3).polynomial_ring() + Multivariate Laurent Polynomial Ring in x0, x1, x2 over Rational Field + sage: AdditiveMultivariateGeneratingSeriesRing('x', 3).polynomial_ring() + Multivariate Polynomial Ring in x0, x1, x2 over Rational Field + """ + return self._polynomial_ring + + def polynomial_ring_extra_var(self): + return self._polynomial_ring_extra_var + + def with_extra_var(self): + return self['EXTRA_VAR'] + + @cached_method + def zero(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: M.zero() + 0 + sage: M.zero().parent() is M + True + sage: M.zero().is_zero() + True + """ + return self._element_constructor_(QQ.zero()) + + @cached_method + def one(self): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: M.zero() + 0 + sage: M.zero().parent() is M + True + sage: M.one().is_one() + True + """ + return self._element_constructor_(QQ.one()) + + def term(self, num, den): + r""" + Return the term ``num / den``. + + INPUT: + + - ``num`` - a Laurent polynomial + + - ``den`` - a list of pairs ``(vector, power)`` or a dictionary + whose keys are the vectors and the values the powers. The + vector ``v = (v_0, v_1, \ldots)`` with power ``n`` corresponds + to the factor `(1 - x_0^{v_0} x_1^{v_1} \ldots x_k^{v_k})^n`. + + EXAMPLES:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: from surface_dynamics.misc.additive_multivariate_generating_series import AdditiveMultivariateGeneratingSeriesRing + + + sage: A = AdditiveMultivariateGeneratingSeriesRing('x', 3) + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 3) + sage: M.term(1, [([1,1,0],1),([1,0,-1],2)]) + (1)/((1 - x0*x2^-1)^2*(1 - x0*x1)) + sage: M.term(1, {(1,1,0): 1, (1,0,2): 2}) + (1)/((1 - x0*x2^2)^2*(1 - x0*x1)) + + sage: A.term(1, [([1,1,0],1),([1,0,-1],2)]) + (1)/((x0 - x2)^2*(x0 + x1)) + sage: A.term(1, {(1,1,0): 1, (1,0,2): 2}) + (1)/((x0 + 2*x2)^2*(x0 + x1)) + + Also works if the denominator is already a factored denominator:: + + sage: from surface_dynamics.misc.factored_denominator import FactoredDenominator + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 3) + sage: den = FactoredDenominator({(1,0,0): 2, (1,1,1):1}) + sage: M.term(1, den) + (1)/((1 - x0)^2*(1 - x0*x1*x2)) + + sage: A.term(1, den) + (1)/((x0)^2*(x0 + x1 + x2)) + + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: R = M.polynomial_ring() + sage: x0, x1 = R.gens() + sage: d = {} + sage: v0 = vector(QQ, (1,0)); v0.set_immutable() + sage: v1 = vector(QQ, (1,1)); v1.set_immutable() + sage: d[v1] = 1 + sage: d[v0] = 1 + sage: M.term(1*x0^2*x1, d) + (x0^2*x1)/((1 - x0)*(1 - x0*x1)) + """ + return self.element_class(self, [(den, num)], self.free_module()) + + def _element_constructor_(self, arg): + r""" + TESTS:: + + sage: from surface_dynamics.misc.multiplicative_multivariate_generating_series import MultiplicativeMultivariateGeneratingSeriesRing + sage: from surface_dynamics.misc.additive_multivariate_generating_series import AdditiveMultivariateGeneratingSeriesRing + + sage: M = MultiplicativeMultivariateGeneratingSeriesRing('x', 2) + sage: A = AdditiveMultivariateGeneratingSeriesRing('x', 2) + sage: M(1) + (1) + sage: A(1) + (1) + + sage: R = M.polynomial_ring() + sage: M(R.0) + (x0) + sage: A(R.0) + (x0) + """ + num = self._polynomial_ring(arg) + return self.element_class(self, [([], num)], self.free_module()) diff --git a/surface_dynamics/misc/generalized_multiple_zeta_values.py b/surface_dynamics/misc/generalized_multiple_zeta_values.py deleted file mode 100644 index 081d6614..00000000 --- a/surface_dynamics/misc/generalized_multiple_zeta_values.py +++ /dev/null @@ -1,985 +0,0 @@ -r""" -Generalized multiple zeta values - -TODO: - -- turns warning for convergence check (ie off/warning/error) - -- introduce linear combinations of *divergent* multiple zeta values - that we possibly simplify at the end. -""" -from __future__ import absolute_import - -import itertools -from collections import defaultdict -import cypari2 - -from sage.all import ZZ, QQ, matrix, bernoulli_polynomial, prod, FreeModule -# Do not use binomial and factorial from sage.all that are slow and broken -from sage.arith.all import binomial, factorial -from sage.sets.disjoint_set import DisjointSet -from sage.misc.cachefunc import cached_method, cached_function - -try: - from sage.modular.multiple_zeta import Multizetas -except (ImportError, cypari2.handle_error.PariError): - def Multizetas(*args, **kwds): - raise ValueError('your sage version does not support multiple zeta values') - -VERBOSE = False - -# pick a new variable and add it in order to the previous ones -# 1 2 3 4 5 6 7 ... -# x y x+y z x+z y+z x+y+z t x+t (y+t) (x+y+t) (z+t) ... -# (1000) (0100) (1100) (0010) (1010) (0110) (1110) (0001) (1001) -def linear_forms(d): - r""" - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import linear_forms - - sage: list(linear_forms(1)) - [(1)] - - sage: list(linear_forms(2)) - [(1, 0), (0, 1), (1, 1)] - - sage: L1 = list(linear_forms(3)) - sage: L2 = L1[:] - sage: L2.sort(key = lambda x: x[::-1]) - sage: assert L1 == L2 - """ - F = FreeModule(ZZ,d) - for n in range(1, 2**d): - v = F(ZZ(n).digits(2,padto=d)) - v.set_immutable() - yield v - -class DivergentZetaError(Exception): - pass - -def handle_term(n, den_tuple): - r""" - Entry point for reduction of generalized multiple zeta values into standard - multiple zeta values. - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import handle_term, is_convergent - - sage: M = Multizetas(QQ) # optional: mzv - - sage: V1 = FreeModule(ZZ, 1) - sage: v = V1((1,)); v.set_immutable() - sage: dt = ((v,3),) - sage: assert is_convergent(1, dt) and handle_term(1, dt) == M((3,)) # optional: mzv - - - sage: V2 = FreeModule(ZZ, 2) - sage: va = V2((1,0)); va.set_immutable() - sage: vb = V2((0,1)); vb.set_immutable() - sage: vc = V2((1,1)); vc.set_immutable() - sage: dt = ((va,2), (vc,3)) - sage: assert is_convergent(2, dt) and handle_term(2, dt) == M((2,3)) # optional: mzv - sage: dt1 = ((va,2),(vb,3)) - sage: dt2 = ((va,3),(vb,2)) - sage: assert is_convergent(2,dt1) and is_convergent(2,dt2) # optional: mzv - sage: assert handle_term(2, ((va,2), (vb,3))) == handle_term(2, ((va,3), (vb,2))) == M((2,)) * M((3,)) # optional: mzv - - sage: V3 = FreeModule(ZZ, 3) - sage: va = V3((1,0,0)); va.set_immutable() - sage: vb = V3((0,1,0)); vb.set_immutable() - sage: vc = V3((0,0,1)); vc.set_immutable() - sage: vd = V3((1,1,0)); vd.set_immutable() - sage: ve = V3((1,0,1)); ve.set_immutable() - sage: vf = V3((0,1,1)); vf.set_immutable() - sage: vg = V3((1,1,1)); vg.set_immutable() - sage: assert handle_term(3, ((va,2), (vd,3), (vg,4))) == M((2,3,4)) # optional: mzv - sage: assert handle_term(3, ((va,2), (vb,3), (vc,4))) == handle_term(3, ((va,3), (vb,2), (vc,4))) # optional mzv - sage: assert handle_term(3, ((va,2), (vb,3), (vc,4))) == M((2,)) * M((3,)) * M((4,)) # optional: mzv - sage: assert handle_term(3, ((va,1), (vc,2), (vd,3))) == handle_term(3, ((va,1), (vb,2), (ve,3))) # optional: mzv - sage: assert handle_term(3, ((va,1), (vc,2), (vd,3))) == handle_term(3, ((va,2), (vb,1), (vf,3))) # optional: mzv - sage: assert handle_term(3, ((va,1), (vc,2), (vd,3))) == M((2,)) * M((1,3)) # optional: mzv - """ - if n == 1: - M = Multizetas(QQ) - ans = M.zero() - for v,p in den_tuple: - # sum (v[0]x)^p - ans += QQ.one() / v[0] **p * M((p,)) - return ans - - elif n == 2: - dat = to_Z2(den_tuple) - if dat is None: - raise NotImplementedError("generalized mzv {}".format(den_tuple)) - return Z2(*dat) - - elif n == 3: - dat = to_Z3(den_tuple) - if dat is None: - raise NotImplementedError("generalized mzv {}".format(den_tuple)) - return Z3(*dat) - - return _handle_term(den_tuple, n) - -@cached_function -def _handle_term(den_tuple, n): - if VERBOSE: - print("handle term({}, {})".format(n, den_tuple)) - assert all(len(v) == n for v,p in den_tuple), (n, den_tuple) - - if any(x < 0 or x > 1 for v,p in den_tuple for x in v): - raise ValueError("unhandled zeta values {}".format(den_tuple)) - - # 0. check convergence - if not is_convergent(n, den_tuple): - raise DivergentZetaError("{} Z({})".format(n, den_tuple)) - - # 1. multizeta - Z = is_multizeta(n, den_tuple) - if Z is not None: - if VERBOSE: - print("standard multizeta") - return Z - - # 2. apply stuffle - P = is_stufflisable(n, den_tuple) - if P is not None: - if VERBOSE: - print("stuffle") - return sum(handle_term(nn, dd) for nn, dd in P) - - # 3. equal rows - # HOW? - - # 4. diminish the number of linear forms without creating convergence problem - data = has_term_sum_of_smaller_terms(n, den_tuple) - if data is not None: - if VERBOSE: - print("relation between linear forms: L_i = sum L_j with L_i > L_j") - return sum(coeff * handle_term(n, new_den_tuple) for coeff, new_den_tuple in kill_relation(n, den_tuple, data[0], data[1])) - - # 5. make "big terms" - data = is_reducible(n, den_tuple) - if data is not None: - if VERBOSE: - print("reduction") - return sum(coeff * handle_term(n, new_den_tuple) for coeff, new_den_tuple in data) - - # 5. 3-variables - if n == 3: - dat = to_Z3(den_tuple) - if dat is None: - raise NotImplementedError("generalized mzv {}".format(den_tuple)) - return Z3(*dat) - - raise NotImplementedError("unhnandled generalized multiple zeta value {}".format(den_tuple)) - -def clean_term(n, den_tuple): - D = {} - for den, p in den_tuple: - if den in D: - D[den] += p - if not D[den]: - del D[den] - else: - D[den] = p - return tuple(sorted(D.items())) - -def to_Z2(den_tuple): - a = b = c = 0 - for v,p in den_tuple: - v = tuple(v) - if len(v) != 2: - raise ValueError - if v == (1,0): - a = p - elif v == (0,1): - b = p - elif v == (1,1): - c = p - else: - return - return a,b,c - -def to_Z3(den_tuple, sort=True): - a = b = c = d = e = f = g = 0 - for v,p in den_tuple: - v = tuple(v) - if len(v) != 3: - raise ValueError - if v == (1,0,0): - a = p - elif v == (0,1,0): - b = p - elif v == (0,0,1): - c = p - elif v == (1,1,0): - d = p - elif v == (1,0,1): - e = p - elif v == (0,1,1): - f = p - elif v == (1,1,1): - g = p - else: - return - # now normalize - if sort: - a,b,c,d,e,f,g = Z3_sort_abc(a,b,c,d,e,f,g) - if b == c: - if a == b: - d,e,f = sorted([d,e,f], reverse=True) - else: - d,e = sorted([d,e], reverse=True) - return a,b,c,d,e,f,g - -def linear_form(R, v): - return sum(R.gen(i) for i,j in enumerate(v) if j) - -def negative_rays(n): - l = [] - v = [0]*n - for i in range(n): - v[i] = -1 - l.append(v[:]) - v[i] = 0 - return l - -def is_convergent(n, den): - r""" - TESTS:: - - sage: import itertools - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import is_convergent - sage: V = FreeModule(ZZ, 3) - sage: va = V((1,0,0)); va.set_immutable() - sage: vb = V((0,1,0)); vb.set_immutable() - sage: vc = V((0,0,1)); vc.set_immutable() - sage: vd = V((1,1,0)); vd.set_immutable() - sage: ve = V((1,0,1)); ve.set_immutable() - sage: vf = V((0,1,1)); vf.set_immutable() - sage: vg = V((1,1,1)); vg.set_immutable() - sage: gens = [va,vb,vc,vd,ve,vf,vg] - sage: N = 0 - sage: for p in itertools.product([0,1,2], repeat=7): # optional: mzv - ....: if sum(map(bool,p)) == 3 and is_convergent(3, list(zip(gens,p))): - ....: print(p) - ....: N += 1 - (0, 0, 0, 0, 1, 1, 2) - (0, 0, 0, 0, 1, 2, 1) - (0, 0, 0, 0, 1, 2, 2) - (0, 0, 0, 0, 2, 1, 1) - (0, 0, 0, 0, 2, 1, 2) - (0, 0, 0, 0, 2, 2, 1) - (0, 0, 0, 0, 2, 2, 2) - (0, 0, 0, 1, 0, 1, 2) - ... - (2, 0, 2, 0, 0, 2, 0) - (2, 0, 2, 2, 0, 0, 0) - (2, 1, 0, 0, 0, 0, 2) - (2, 1, 0, 0, 0, 2, 0) - (2, 2, 0, 0, 0, 0, 2) - (2, 2, 0, 0, 0, 2, 0) - (2, 2, 0, 0, 2, 0, 0) - (2, 2, 2, 0, 0, 0, 0) - sage: print(N) # optional: mzv - 125 - """ - from sage.geometry.polyhedron.constructor import Polyhedron - from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing - - assert all(len(v) == n for v,p in den), (n, den) - - # TODO: fast code path - - R = PolynomialRing(QQ, 'x', n) - x = R.gens() - den_poly = prod(linear_form(R, v)**p for v,p in den) - newton_polytope = Polyhedron(vertices=den_poly.exponents(), rays=negative_rays(n)) - V = newton_polytope.intersection(Polyhedron(rays=[[1]*n])).vertices() - r = max(max(v.vector()) for v in V) - return r > 1 - -def convergent_multizeta(t): - r""" - Multizeta value at a convergent index ``t``. - - TESTS:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import convergent_multizeta - sage: assert all(convergent_multizeta(t) == Multizeta(*t) for t in [(2,),(3,),(1,2),(3,2),(1,1,2)]) # optional: mzv - - sage: convergent_multizeta((0,3)) # optional: mzv - ζ(2) - ζ(3) - sage: convergent_multizeta((0,2,2)) # optional: mzv - ζ(1,2) - ζ(2,2) - sage: convergent_multizeta((1,0,3)) # optional: mzv - ζ(1,2) - ζ(1,3) - ζ(2) + ζ(3) - sage: convergent_multizeta((0,1,3)) # optional: mzv - -ζ(1,3) + ζ(2) - ζ(3) - - sage: convergent_multizeta((0, 4)) # optional: mzv - ζ(3) - ζ(4) - sage: convergent_multizeta((-1, 5)) # optional: mzv - 1/2*ζ(3) - 1/2*ζ(4) - sage: convergent_multizeta((-2, 5)) # optional: mzv - 1/3*ζ(2) - 1/2*ζ(3) + 1/6*ζ(4) - - sage: convergent_multizeta((-1, 3, 4)) # optional: mzv - 1/2*ζ(1,4) - 1/2*ζ(2,4) - sage: convergent_multizeta((-1, -1, 8)) # optional: mzv - 1/8*ζ(4) - 5/12*ζ(5) + 3/8*ζ(6) - 1/12*ζ(7) - sage: convergent_multizeta((-2,-2,10)) # optional: mzv - 1/18*ζ(4) - 4/15*ζ(5) + 31/72*ζ(6) - 1/4*ζ(7) + 1/72*ζ(8) + 1/60*ζ(9) - sage: convergent_multizeta((-1, -2, 10)) # optional: mzv - 1/10*ζ(5) - 3/8*ζ(6) + 5/12*ζ(7) - 1/8*ζ(8) - 1/60*ζ(9) - - sage: convergent_multizeta((4,-4,10)) # optional: mzv - -1/3*ζ(1,10) + 1/30*ζ(3,10) + 1/5*ζ(4,5) - 1/2*ζ(4,6) + 1/3*ζ(4,7) - 1/30*ζ(4,9) - 1/10*ζ(8) - 2/5*ζ(9) + 1/2*ζ(10) - sage: convergent_multizeta((2,-1,8)) # optional: mzv - -1/2*ζ(1,8) + 1/2*ζ(2,6) - 1/2*ζ(2,7) - 1/2*ζ(7) + 1/2*ζ(8) - - sage: convergent_multizeta((0,3,2)) # optional: mzv - ζ(2,2) - ζ(3,2) - - Divergent cases:: - - sage: convergent_multizeta((0,2)) # optional: mzv - Traceback (most recent call last): - ... - DivergentZetaError: divergent multizeta value (0, 2) - sage: convergent_multizeta((0,0,3)) # optional: mzv - Traceback (most recent call last): - ... - DivergentZetaError: divergent multizeta value (0, 0, 3) - sage: convergent_multizeta((1,0,2)) # optional: mzv - Traceback (most recent call last): - ... - DivergentZetaError: divergent multizeta value (1, 0, 2) - sage: convergent_multizeta((0,1,2)) # optional: mzv - Traceback (most recent call last): - ... - DivergentZetaError: divergent multizeta value (0, 1, 2) - """ - if VERBOSE: - print("convergent_multizeta({})".format(t)) - n = len(t) - for i in range(1,n+1): - if sum(t[-i:]) <= i: - raise DivergentZetaError("divergent multizeta value {}".format(t)) - if all(x > 0 for x in t): - return Multizetas(QQ)(t) - else: - # find the first non-positive index and drop the summation on the - # corresponding variable using Faulhaber's formula - i = 0 - while t[i] > 0: - i += 1 - t = list(t) - tt = t[:i] + t[i+1:] - x = QQ['x'].gen() - if t[i] == 0: - # bug: bernoulli_polynomial(x, 0) is an integer - P = x - 1 - else: - P = bernoulli_polynomial(QQ['x'].gen(), -t[i]).integral() - s = Multizetas(QQ).zero() - for c,e in zip(P.coefficients(), P.exponents()): - tt[i] -= e - s += c * convergent_multizeta(tt[:]) - tt[i] += e - if i > 0: - Q = P(x+1) - for c,e in zip(Q.coefficients(), Q.exponents()): - tt[i-1] -= e - s -= c * convergent_multizeta(tt[:]) - tt[i-1] += e - return s - -def is_multizeta(n, den_tuple): - r""" - Return the corresponding zeta values if it is one. - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import linear_forms, is_multizeta - sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) - - sage: is_multizeta(3, ((vb, 2), (vf, 5), (vg, 2))) # optional: mzv - ζ(2,5,2) - - sage: is_multizeta(3, [(vg, 5)]) # optional: mzv - 1/2*ζ(3) - 3/2*ζ(4) + ζ(5) - - sage: assert is_multizeta(3, ((va,3), (vb,3), (vc,3))) is None # optional: mzv - sage: assert is_multizeta(3, ((vb,2), (ve,5), (vg,1))) is None # optional: mzv - """ - assert all(len(v) == n for v,p in den_tuple), (n, den_tuple) - - by_support = [None for _ in range(n+1)] - for vp in den_tuple: - v,p = vp - if any(x > 1 for x in v): - return - m = sum(vp[0]) - if by_support[m] is None: - by_support[m] = vp - else: - return - - k = 0 - linear_order = [None] * n - for step, vp in enumerate(by_support): - if vp is None: - continue - v,p = vp - # discover new variables - for i in range(n): - if v[i] and linear_order[i] is None: - linear_order[i] = k - k += 1 - if k > step: - return - - if all(v is not None for v in linear_order[1:]): - return convergent_multizeta([0 if vp is None else vp[1] for vp in by_support[1:]]) - -def is_stufflisable(n, den_tuple): - m = len(den_tuple) - possibilities = [] - for i in range(n): - for j in range(i): - if all(not (v[i] and v[j]) for v,p in den_tuple): - # see if we can perform a stuffle that will lower rank in next step - for l in range(n): - if all(v[i] + v[j] == v[l] for v,p in den_tuple): - return stuffle(n, den_tuple, i, j) - possibilities.append((i,j)) - - # pick the first candidate ? - # multiple columns shuffling ? - if possibilities: - i, j = possibilities[0] - return stuffle(n, den_tuple, i, j) - -def stuffle(n, den_tuple, i, j): - m = matrix([v for v,p in den_tuple]).transpose() - - # big subsimplex 1 - mm = m.__copy__() - mm[i] += mm[j] - D1 = clean_term(n, tuple((c,p) for c,(v,p) in zip(mm.columns(), den_tuple))) - - # big subsimplex 2 - mm = m.__copy__() - mm[j] += mm[i] - D2 = clean_term(n, tuple((c,p) for c,(v,p) in zip(mm.columns(), den_tuple))) - - # small subsimplex - mm = m.__copy__() - mm[i] += mm[j] - mm = mm.delete_rows([j]) - D3 = clean_term(n-1, tuple((c,p) for c,(v,p) in zip(mm.columns(), den_tuple))) - - return [(n,D1),(n,D2),(n-1,D3)] - -def is_product(n, den_tuple): - r""" - INPUT: - - - ``n`` - number of variables - - - ``den_tuple`` - tuple of pairs ``(vector, power)`` - - TESTS:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import is_product - sage: is_product(3, [((1,0,0),2), ((0,1,0),5), ((1,1,0),1), ((0,0,1),3)]) - [(2, (((0, 1), 5), ((1, 0), 2), ((1, 1), 1))), (1, (((1), 3),))] - sage: is_product(3, [((1,0,0),2), ((0,1,0),3), ((1,0,1),1), ((0,0,1),5)]) - [(2, (((0, 1), 5), ((1, 0), 2), ((1, 1), 1))), (1, (((1), 3),))] - sage: is_product(3, [((1,1,1),3)]) is None - True - """ - D = DisjointSet(n) - assert all(len(v) == n for v,p in den_tuple), (n, den_tuple) - - # 1. product structure - for v,_ in den_tuple: - i0 = 0 - while not v[i0]: - i0 += 1 - i = i0 + 1 - while i < n: - if v[i]: - D.union(i0, i) - i += 1 - if D.number_of_subsets() == 1: - # no way to split variables - return - - # split variables - Rdict = D.root_to_elements_dict() - keys = sorted(Rdict.keys()) - key_indices = {k: i for i,k in enumerate(keys)} - values = [Rdict[k] for k in keys] - values_indices = [{v:i for i,v in enumerate(v)} for v in values] - n_list = [len(J) for J in values] - F = [FreeModule(ZZ, nn) for nn in n_list] - new_terms = [[] for _ in range(len(Rdict))] - for v,p in den_tuple: - i0 = 0 - while not v[i0]: - i0 += 1 - i0 = D.find(i0) - assert all(D.find(i) == i0 for i in range(n) if v[i]), (i0, [D.find(i) for i in range(n) if v[i]]) - k = key_indices[i0] - vv = F[k]() - for i in range(n): - if v[i]: - vv[values_indices[k][i]] = v[i] - vv.set_immutable() - new_terms[k].append((vv,p)) - - return list(zip(n_list, [tuple(sorted(terms)) for terms in new_terms])) - -# TODO: apply a power of a relation at once with multinomial coefficients -def apply_relation(n, den_tuple, i, relation): - # iterator through the pairs (coeff, den_tuples) obtained by applying relation on the i-th term - ans = [] - for j in range(n): - if relation[j] and j != i: - term = list(den_tuple) - v, p = term[i] - term[i] = (v, p+1) - v, p = term[j] - term[j] = (v, p-1) - yield (relation[j], tuple(term)) - -def kill_relation(n, den_tuple, i, relation): - r""" - Make calls to `apply_relation` until we get rid of term - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import kill_relation - - sage: V = FreeModule(ZZ, 2) - sage: va = V((1,0)); va.set_immutable() - sage: vb = V((0,1)); vb.set_immutable() - sage: vc = V((1,1)); vc.set_immutable() - - sage: den_tuple = ((va,2),(vb,3),(vc,4)) - sage: kill_relation(2, den_tuple, 2, [1,1,0]) - [(1, (((0, 1), 3), ((1, 1), 6))), - (2, (((0, 1), 2), ((1, 1), 7))), - (1, (((1, 0), 2), ((1, 1), 7))), - (3, (((0, 1), 1), ((1, 1), 8))), - (3, (((1, 0), 1), ((1, 1), 8)))] - """ - assert len(relation) == len(den_tuple) - assert 0 <= i < len(den_tuple) - assert sum(relation[j] * den_tuple[j][0] for j in range(len(den_tuple))) == den_tuple[i][0] - D = {den_tuple: QQ.one()} - todo = [den_tuple] - while todo: - den_tuple = todo.pop(0) - coeff1 = D.pop(den_tuple) - for coeff, new_den_tuple in apply_relation(n, den_tuple, i, relation): - coeff *= coeff1 - if new_den_tuple in D: - D[new_den_tuple] += coeff - elif any(not new_den_tuple[j][1] for j in range(len(den_tuple))): - new_den_tuple = tuple(x for x in new_den_tuple if x[1]) - if new_den_tuple not in D: - D[new_den_tuple] = coeff - else: - D[new_den_tuple] += coeff - else: - todo.append(new_den_tuple) - D[new_den_tuple] = coeff - - return [(coeff, new_den_tuple) for new_den_tuple, coeff in D.items()] - -def has_term_sum_of_smaller_terms(n, den_tuple): - r""" - Look for a vector v_i and {v_j} with sum(v_j) = v_i - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import linear_forms, has_term_sum_of_smaller_terms - - sage: va,vb,vc = linear_forms(2) - sage: has_term_sum_of_smaller_terms(2, ((va,1),(vb,1),(vc,1))) - [2, [1, 1, 0]] - sage: assert has_term_sum_of_smaller_terms(2, ((va,1),(vc,1))) is None - sage: assert has_term_sum_of_smaller_terms(2, ((vb,1),(vc,1))) is None - sage: assert has_term_sum_of_smaller_terms(2, ((va,1),(vb,1))) is None - - sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) - sage: has_term_sum_of_smaller_terms(3, ((va,1),(vb,1),(vc,1),(vg,1))) - [3, [1, 1, 1, 0]] - sage: has_term_sum_of_smaller_terms(3, ((va,1),(vb,1),(ve,1),(vf,1),(vg,1))) - [4, [0, 1, 1, 0, 0]] - sage: has_term_sum_of_smaller_terms(3, ((vd,1),(ve,1),(vf,1),(vg,1))) - [3, [1/2, 1/2, 1/2, 0]] - sage: has_term_sum_of_smaller_terms(3, ((va,1),(vd,1),(ve,1),(vg,1))) - [3, [-1, 1, 1, 0]] - """ - l = len(den_tuple) - for iv in range(l): - v = den_tuple[iv][0] - if sum(v) == 1: - continue - - candidates = [k for k,(u,_) in enumerate(den_tuple) if sum(u) < sum(v) and all(u[i] <= v[i] for i in range(n))] - M = matrix(QQ, [den_tuple[k][0] for k in candidates]) - try: - coeffs = M.solve_left(v) - except ValueError: - continue - full = [0] * l - for k,c in zip(candidates, coeffs): - full[k] = c - return [iv, full] - -def write_array(level, explanations): - for i,key in enumerate(sorted(level, key=lambda x: 1000 if level[x] is None else level[x])): - print("%s & %s \\\\" % (lin_prod(key),'\\infty' if level[key] is None else level[key])) - -def is_reducible(n, den_tuple): - r""" - If (x1+x2+...+xd) is not present, use a linear relation to create it. Then - try to kill using other forms. Then try to write as P(x1,x2,...,x_{d-1}) (x1+x2+...+xd) and recurse. - - Should solve all Tonrheim - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import linear_forms, is_reducible - sage: va, vb, vd, vc, ve, vf, vg = linear_forms(3) - sage: is_reducible(3, ((va,3),(vb,3),(vc,3))) - [(1, (((0, 0, 1), 3), ((0, 1, 1), 3), ((1, 1, 1), 3))), - (1, (((0, 1, 0), 3), ((0, 1, 1), 3), ((1, 1, 1), 3))), - (3, (((0, 0, 1), 2), ((0, 1, 1), 4), ((1, 1, 1), 3))), - (3, (((0, 1, 0), 2), ((0, 1, 1), 4), ((1, 1, 1), 3))), - ... - (90, (((1, 0, 0), 1), ((1, 0, 1), 1), ((1, 1, 1), 7))), - (90, (((0, 1, 0), 1), ((1, 1, 0), 1), ((1, 1, 1), 7))), - (90, (((1, 0, 0), 1), ((1, 1, 0), 1), ((1, 1, 1), 7)))] - """ - F = FreeModule(ZZ, n) - vmax = F([1] * n) - vmax.set_immutable() - - if len(den_tuple) == 1: - return - - # force the max vector (1,1,...,1) to appear - imax = None - for i,(v,p) in enumerate(den_tuple): - if v == vmax: - imax = i - break - if imax is None: - imax = len(den_tuple) - den_tuple = den_tuple + ((vmax,0),) - if imax != len(den_tuple) - 1: - den_tuple = list(den_tuple) - den_tuple.append(den_tuple.pop(imax)) - den_tuple = tuple(den_tuple) - imax = len(den_tuple) - 1 - - assert den_tuple[imax][0] == vmax - - M = matrix(QQ, [v for v,p in den_tuple if v != vmax]) - try: - relation = M.solve_left(vmax) - except ValueError: - if den_tuple[-1][1] == 0: - return - den_tuple2 = den_tuple[:-1] - variables = set().union(*[[i for i in range(n) if v[i]] for v,p in den_tuple2]) - if len(variables) == n: - return - killed = [(1, den_tuple)] - else: - killed = kill_relation(n, den_tuple, imax, list(relation) + [0]) - - ans = defaultdict(QQ) - for coeff, den_tuple2 in killed: - assert den_tuple2[-1][0] == vmax - pmax = den_tuple2[-1][1] - assert pmax - den_tuple2 = den_tuple2[:-1] - - variables = set().union(*[[i for i in range(n) if v[i]] for v,p in den_tuple2]) - if len(variables) == n: - # removing the maximal vector (1,1,...,1) is not enough to make a variable disappear - data = is_reducible(n, den_tuple2) - if data is None: - data = [(1, den_tuple2)] - else: - # less variables! - nn = len(variables) - variables = sorted(variables) - new_indices = {j:i for i,j in enumerate(variables)} - G = FreeModule(ZZ, nn) - new_den_tuple2 = [] - for v,p in den_tuple2: - vv = G([v[i] for i in variables]) - vv.set_immutable() - new_den_tuple2.append((vv,p)) - new_den_tuple2 = tuple(new_den_tuple2) - data = is_reducible(nn, new_den_tuple2) - if data is None: - data = [(1, new_den_tuple2)] - - # lift to the n variables version - new_data = [] - for coeff3, den_tuple3 in data: - den_tuple3 = [(F([v[new_indices[j]] if j in new_indices else 0 for j in range(n)]), p) for v,p in den_tuple3] - for v,p in den_tuple3: - v.set_immutable() - new_data.append((coeff3, tuple(sorted(den_tuple3)))) - data = new_data - - # update the answer - for coeff3, den_tuple3 in data: - imax = None - for i,(v,p) in enumerate(den_tuple3): - if v == vmax: - imax = i - break - if imax is None: - den_tuple3 = den_tuple3 + ((vmax,pmax),) - else: - den_tuple3 = den_tuple3[:imax] + ((vmax, den_tuple3[imax][1] + pmax),) + den_tuple3[imax+1:] - ans[den_tuple3] += coeff * coeff3 - - if len(ans) > 1: - return [(coeff, den_tuple) for den_tuple, coeff in ans.items()] - - -def Z3_sort_abc(a,b,c,d,e,f,g): - if b > a: - if c > b: - # c > b > a - a,b,c,d,e,f,g = c,b,a,f,e,d,g - elif c > a: - # b >= c > a - a,b,c,d,e,f,g = b,c,a,f,d,e,g - else: - # b > a >= c - a,b,c,d,e,f,g = b,a,c,d,f,e,g - elif c > a: - # c > a >= b - a,b,c,d,e,f,g = c,a,b,e,f,d,g - elif c > b: - # a >= c >= b - a,b,c,d,e,f,g = a,c,b,e,d,f,g - else: - pass - assert a >= b >= c - return a,b,c,d,e,f,g - -def is_Z2_convergent(a,b,c): - r""" - TESTS:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import is_Z2_convergent - - Convergent examples:: - - sage: assert is_Z2_convergent(1,1,1) - sage: assert is_Z2_convergent(2,2,0) - sage: assert is_Z2_convergent(0,0,3) - - Divergent examples:: - - sage: assert not is_Z2_convergent(0,0,2) - sage: assert not is_Z2_convergent(1,2,0) - """ - from sage.geometry.polyhedron.constructor import Polyhedron - x, y = ZZ['x,y'].gens() - poly = x**a * y**b * (x+y)**c - newton_polytope = Polyhedron(vertices=poly.exponents(), rays=[(-1,0),(0,-1)]) - V = newton_polytope.intersection(Polyhedron(rays=[(1,1)])).vertices() - r = max(max(v.vector()) for v in V) - return r > 1 - -def is_Z3_convergent(a,b,c,d,e,f,g): - r""" - TESTS:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import is_Z3_convergent - - Convergent examples:: - - sage: assert is_Z3_convergent(2,0,2,2,0,0,0) - sage: assert is_Z3_convergent(0,0,0,0,0,0,4) - sage: assert is_Z3_convergent(1,0,0,1,0,0,2) - sage: assert is_Z3_convergent(0,1,0,1,0,0,2) - sage: assert is_Z3_convergent(0,1,0,0,0,1,2) - - Divergent examples:: - - sage: assert not is_Z3_convergent(0,0,0,1,1,1,0) - sage: assert not is_Z3_convergent(0,0,0,0,0,0,3) - - """ - from sage.geometry.polyhedron.constructor import Polyhedron - x, y, z = ZZ['x,y,z'].gens() - poly = x**a * y**b * z**c * (x+y)**d * (x+z)**e * (y+z)**f * (x+y+z)**g - newton_polytope = Polyhedron(vertices=poly.exponents(), rays=[(-1,0,0),(0,-1,0),(0,0,-1)]) - V = newton_polytope.intersection(Polyhedron(rays=[(1,1,1)])).vertices() - r = max(max(v.vector()) for v in V) - return r > 1 - -def Z2(a, b, c, check_convergence=True): - M = Multizetas(QQ) - if a == 0 and b == 0: - return M((c-1,)) - M((c,)) - else: - return sum(binomial(a+i-1,i) * M((b-i,c+a+i)) for i in range(b)) + sum(binomial(b+i-1,i) * M((a-i,c+b+i)) for i in range(a)) - -def Z3(a, b, c, d, e, f, g, check_convergence=True): - r""" - The function ``Z3(a,b,c,d,e,f,g)``. - - .. MATH:: - - \sum_{x,y,z \geq 1} x^{-a} y^{-b} z^{-c} (x+y)^{-d} (x+z)^{-e} (y+z)^{-f} (x+y+z)^{-g} - - The reduction algorithm was designed by Bill Allombert. - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import Z3 - - sage: M = Multizetas(QQ) # optional: mzv - - sage: Z3(1,1,1,1,1,1,1) # optional: mzv - 21/2*ζ(1,1,5) + 9/2*ζ(1,2,4) - 3/2*ζ(1,3,3) - 3/2*ζ(1,4,2) + 9/2*ζ(1,6) - sage: Z3(3,0,0,0,0,3,0) # optional: mzv - 6*ζ(1,4) - 12*ζ(1,5) + 3*ζ(2,3) - 6*ζ(2,4) + ζ(3,2) - 2*ζ(3,3) - - sage: assert Z3(2,3,4,0,0,0,0) == M((2,)) * M((3,)) * M((4,)) # optional: mzv - sage: assert Z3(1,0,0,2,0,0,3) == M((1,2,3)) # optional: mzv - - sage: assert Z3(0,0,0,2,0,1,1) == M((4,)) / 2 # optional: mzv - sage: assert Z3(1,0,1,1,0,0,1) == 3 * M((1,1,2)) # optional: mzv - - sage: assert Z3(0,0,0,0,0,0,4) == 1/2 * M((2,)) - 3/2 * M((3,)) + M((4,)) # optional: mzv - - sage: assert Z3(0,0,0,2,0,1,1) == 2 * M((1,1,2)) - M((2,2)) - 3 *M((1,3)) # optional: mzv - """ - M = Multizetas(QQ) - CHECK_CONVERGENCE = False - - # x^a y^b z^c (x+y)^d (x+z)^e (y+z)^f (x+y+z)^g - if VERBOSE: - print("Z3({},{},{},{},{},{},{})".format(a,b,c,d,e,f,g)) - if a < 0 or b < 0 or c < 0 or d < 0 or e < 0 or f < 0 or g < 0: - raise ValueError("invalid exponents for Z3: a={} b={} c={} d={} e={} f={} g={}".format(a,b,c,d,e,f,g)) - - if check_convergence and not is_Z3_convergent(a,b,c,d,e,f,g): - raise DivergentZetaError("divergent Z3({},{},{},{},{},{},{})".format(a,b,c,d,e,f,g)) - - # step 1: try to get rid of the terms (x+y), (x+z), (y+z) - if d and e and f: - if VERBOSE: - print("reduction (x+y+z) = ((x+y) + (x+z) + (y+z)) / 2") - return (Z3(a,b,c,d-1,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c,d,e-1,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c,d,e,f-1,g+1,CHECK_CONVERGENCE)) / 2 - if a and f: - # x^a (y+z)^f -> (x+y+z)^g - # additive version Z3(a-1,b,c,d,e,f,g+1) + Z3(a,b,c,d,e,f-1,g+1) - if VERBOSE: - print("reduction (x+y+z) = (x) + (y+z)") - return sum(binomial(a+k-1, k) * Z3(0,b,c,d,e,f-k,g+a+k,CHECK_CONVERGENCE) for k in range(f)) + \ - sum(binomial(f+k-1, k) * Z3(a-k,b,c,d,e,0,g+f+k,CHECK_CONVERGENCE) for k in range(a)) - if b and e: - # y^b (x+z)^e -> (x+y+z)^g - # additive version Z3(a,b-1,c,d,e,f,g+1) + Z3(a,b,c,d,e-1,f,g+1) - if VERBOSE: - print("reduction (x+y+z) = (y) + (x+z)") - return sum(binomial(b+k-1, k) * Z3(a,0,c,d,e-k,f,g+b+k,CHECK_CONVERGENCE) for k in range(e)) + \ - sum(binomial(e+k-1, k) * Z3(a,b-k,c,d,0,f,g+e+k,CHECK_CONVERGENCE) for k in range(b)) - if c and d: - # z^c (x+y)^d -> (x+y+z)^g - # additive version Z3(a,b,c-1,d,e,f,g+1) + Z3(a,b,c,d-1,e,f,g+1) - if VERBOSE: - print("reduction (x+y+z) = (z) + (x+y)") - return sum(binomial(c+k-1, k) * Z3(a,b,0,d-k,e,f,g+c+k,CHECK_CONVERGENCE) for k in range(d)) + \ - sum(binomial(d+k-1, k) * Z3(a,b,c-k,0,e,f,g+d+k,CHECK_CONVERGENCE) for k in range(c)) - - assert d*e*f == a*f == b*e == c*d == 0 - - a,b,c,d,e,f,g = Z3_sort_abc(a,b,c,d,e,f,g) - - # step 2: kill c - if c: - if VERBOSE: - print("reduction (x+y+z) = (x) + (y) + (z)") - return Z3(a-1,b,c,d,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b-1,c,d,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c-1,d,e,f,g+1,CHECK_CONVERGENCE) - - assert c == 0 - - if b: - # x^a y^b -> (x+y)^d - # additive version Z3(a-1,b,0,d+1,e,f,g) + Z3(a,b-1,0,d+1,e,f,g) - if VERBOSE: - print("reduction (x+y) = (x) + (y)") - return sum(binomial(a+k-1, k) * Z3(0,b-k,c,d+a+k,e,f,g,CHECK_CONVERGENCE) for k in range(b)) + \ - sum(binomial(b+k-1, k) * Z3(a-k,0,c,d+b+k,e,f,g,CHECK_CONVERGENCE) for k in range(a)) - - assert b == c == 0 - - assert b == c == d*e*f == 0 - - if a and f: - raise RuntimeError - - assert b == c == d*e*f == a*f == 0 - - if a == 0: - d,e,f = sorted([d,e,f], reverse=True) - - if a and d and e: - if VERBOSE: - print("reduction (x+y+z) = (x+y) + (x+z) - (x)") - return Z3(a,b,c,d-1,e,f,g+1,CHECK_CONVERGENCE) + Z3(a,b,c,d,e-1,f,g+1,CHECK_CONVERGENCE) - Z3(a-1,b,c,d,e,f,g+1,CHECK_CONVERGENCE) - - assert b == c == d*e*f == a*f == a*d*e == 0 - - # x^a (x+y)^d (x+y+z)^g - if e == f == 0: - return convergent_multizeta((a,d,g)) - # x^a (x+z)^e (x+y+z)^g - if d == f == 0: - return convergent_multizeta((a,e,g)) - if a == f == 0: - if g == 0: - # (x+y)^d (x+z)^e - assert d and e - return M((d-1,e)) - M((d,e)) + M((e-1,d)) - M((e,d)) + M((d+e-1,)) - M((d+e,)) - elif d == e == 1: - if VERBOSE: - print("Bill formulas") - return sum(M((1,g+1-i,i)) for i in range(2,g+1)) - 3*M((1,g+1)) - elif g == 1: - if VERBOSE: - print("[0,0,0,d,0,f,1] = [f,0,d-1,1,0,0,1] - [f,d,1] - zeta([d,1,f]) - zeta([1+d,f])") - f,d = sorted([d,e]) - return Z3(f,0,d-1,1,0,0,1,CHECK_CONVERGENCE) - Z2(f,d,1,CHECK_CONVERGENCE) - M((f,1,d)) - M((f,d+1)) - else: - # THIS LOOKS WRONG - # [0,0,0,d,0,f,g] = [f,0,d,g,0,0,0] - [f,d,0,0,0,0,g] - [f,d,g] - zeta([d,g,f]) - zeta([g+d,f]) - if VERBOSE: - print("[0,0,0,0,e,f,g] = [e,0,f,g,0,0,0] - [e,0,f,0,0,0,g] - Z(e,g,f) - Z(e,f+g) - Z(e+f,g)") - f,d = sorted([d,e]) - return Z3(f,0,d,g,0,0,0,CHECK_CONVERGENCE) - Z3(f,d,0,0,0,0,g,CHECK_CONVERGENCE) - Z2(f,d,g,CHECK_CONVERGENCE) - M((f,g,d)) - M((f,g+d)) - - raise RuntimeError diff --git a/surface_dynamics/misc/groups.py b/surface_dynamics/misc/groups.py new file mode 100644 index 00000000..b81b8c67 --- /dev/null +++ b/surface_dynamics/misc/groups.py @@ -0,0 +1,26 @@ +def relations(x, y): + r""" + Iterate through the relations between ``x`` and ``y`` by length + """ + lx = 'x' + ly = 'y' + lX = 'X' + lY = 'Y' + if x.is_one(): + return lx + if y.is_one(): + return ly + gens = [(lx, x), (ly, y), (lX, ~x), (lY, ~y)] + inv = {lx: lX, lX:lx, ly:lY, lY:ly} + level = gens + while True: + new_level = [] + for w, g in level: + for a, s in gens: + if w[-1] != inv[a]: + gg = g * s + ww = w + a + if gg.is_one(): + yield ww + new_level.append((ww, gg)) + level = new_level diff --git a/surface_dynamics/misc/linalg.py b/surface_dynamics/misc/linalg.py index 9294f834..14c53172 100644 --- a/surface_dynamics/misc/linalg.py +++ b/surface_dynamics/misc/linalg.py @@ -2,7 +2,7 @@ Some linear algebra routines """ #***************************************************************************** -# Copyright (C) 2019 Vincent Delecroix <20100.delecroix@gmail.com> +# Copyright (C) 2019-2023 Vincent Delecroix <20100.delecroix@gmail.com> # # Distributed under the terms of the GNU General Public License (GPL) # as published by the Free Software Foundation; either version 2 of @@ -13,8 +13,10 @@ from sage.arith.all import gcd, lcm from sage.arith.misc import binomial from sage.rings.all import ZZ, QQ +from sage.matrix.constructor import matrix from sage.geometry.polyhedron.constructor import Polyhedron + def relation_space(v): r""" Relation space of the given vector ``v`` @@ -233,3 +235,143 @@ def symbolic_matrix_power(M, n): P *= N p += 1 return result + + +def disjoint_vectors(vectors, min_size=2, max_size=None): + r""" + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import disjoint_vectors + sage: V = FreeModule(ZZ, 4) + sage: v0 = V((1,0,0,0)) + sage: v1 = V((1,1,0,0)) + sage: v2 = V((0,0,1,0)) + sage: v3 = V((0,1,1,1)) + sage: v4 = V((0,1,0,1)) + sage: for subset, s in disjoint_vectors((v0,v1,v2,v3,v4)): + ....: print(subset, s) + [0, 2] (1, 0, 1, 0) + [0, 2, 4] (1, 1, 1, 1) + [0, 3] (1, 1, 1, 1) + [0, 4] (1, 1, 0, 1) + [1, 2] (1, 1, 1, 0) + [2, 4] (0, 1, 1, 1) + sage: for subset, s in disjoint_vectors((v0,v1,v2,v3,v4), min_size=2, max_size=2): + ....: print(subset, s) + [0, 2] (1, 0, 1, 0) + [0, 3] (1, 1, 1, 1) + [0, 4] (1, 1, 0, 1) + [1, 2] (1, 1, 1, 0) + [2, 4] (0, 1, 1, 1) + """ + if max_size is None: + max_size = len(vectors) + F = vectors[0].parent() + n = len(vectors) + d = F.dimension() + stack = [] + i = 0 + c = F.zero() # current sum + while True: + while i < len(vectors) and len(stack) < max_size: + cc = c + vectors[i] + if all(x < 2 for x in cc): + stack.append(i) + c = cc + if len(stack) >= min_size: + c.set_immutable() + yield stack, c + i += 1 + if not stack: + return + i = stack.pop() + c -= vectors[i] + i += 1 + + +def linearly_independent_vectors(vectors, min_size=0, max_size=None): + r""" + Iterate through the subsets of ``vectors`` made of linearly independent vectors. + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import linearly_independent_vectors + sage: V = FreeModule(ZZ, 4) + sage: vecs = [V((1, 0, 0, 0)), V((1, 1, 0, 0)), V((0, 0, 1, 0)), V((0, 1, 1, 1)), V((0, 1, 0, 1))] + sage: for subset, s in linearly_independent_vectors(vecs): + ....: print(subset, s.rows()) + [0] [(1, 0, 0, 0)] + [0, 1] [(1, 0, 0, 0), (1, 1, 0, 0)] + [0, 1, 2] [(1, 0, 0, 0), (1, 1, 0, 0), (0, 0, 1, 0)] + [0, 1, 2, 3] [(1, 0, 0, 0), (1, 1, 0, 0), (0, 0, 1, 0), (0, 1, 1, 1)] + [0, 1, 2, 4] [(1, 0, 0, 0), (1, 1, 0, 0), (0, 0, 1, 0), (0, 1, 0, 1)] + [0, 1, 3] [(1, 0, 0, 0), (1, 1, 0, 0), (0, 1, 1, 1)] + [0, 1, 3, 4] [(1, 0, 0, 0), (1, 1, 0, 0), (0, 1, 1, 1), (0, 1, 0, 1)] + [0, 1, 4] [(1, 0, 0, 0), (1, 1, 0, 0), (0, 1, 0, 1)] + [0, 2] [(1, 0, 0, 0), (0, 0, 1, 0)] + [0, 2, 3] [(1, 0, 0, 0), (0, 0, 1, 0), (0, 1, 1, 1)] + [0, 2, 4] [(1, 0, 0, 0), (0, 0, 1, 0), (0, 1, 0, 1)] + [0, 3] [(1, 0, 0, 0), (0, 1, 1, 1)] + [0, 3, 4] [(1, 0, 0, 0), (0, 1, 1, 1), (0, 1, 0, 1)] + [0, 4] [(1, 0, 0, 0), (0, 1, 0, 1)] + [1] [(1, 1, 0, 0)] + [1, 2] [(1, 1, 0, 0), (0, 0, 1, 0)] + [1, 2, 3] [(1, 1, 0, 0), (0, 0, 1, 0), (0, 1, 1, 1)] + [1, 2, 4] [(1, 1, 0, 0), (0, 0, 1, 0), (0, 1, 0, 1)] + [1, 3] [(1, 1, 0, 0), (0, 1, 1, 1)] + [1, 3, 4] [(1, 1, 0, 0), (0, 1, 1, 1), (0, 1, 0, 1)] + [1, 4] [(1, 1, 0, 0), (0, 1, 0, 1)] + [2] [(0, 0, 1, 0)] + [2, 3] [(0, 0, 1, 0), (0, 1, 1, 1)] + [2, 4] [(0, 0, 1, 0), (0, 1, 0, 1)] + [3] [(0, 1, 1, 1)] + [3, 4] [(0, 1, 1, 1), (0, 1, 0, 1)] + [4] [(0, 1, 0, 1)] + + sage: for subset, s in linearly_independent_vectors(vecs, min_size=2, max_size=2): + ....: print(subset) + [0, 1] + [0, 2] + [0, 3] + [0, 4] + [1, 2] + [1, 3] + [1, 4] + [2, 3] + [2, 4] + [3, 4] + sage: for subset, s in linearly_independent_vectors(vecs, min_size=3, max_size=3): + ....: print(subset) + [0, 1, 2] + [0, 1, 3] + [0, 1, 4] + [0, 2, 3] + [0, 2, 4] + [0, 3, 4] + [1, 2, 3] + [1, 2, 4] + [1, 3, 4] + """ + F = vectors[0].parent() + n = len(vectors) + d = F.dimension() + if max_size is None: + max_size = d + stack = [] + i = 0 + M = matrix(ZZ, max_size, d) # current matrix + while True: + while len(stack) < max_size and i < len(vectors): + M[len(stack)] = vectors[i] + if M.rank() == len(stack) + 1: + stack.append(i) + if len(stack) >= min_size: + yield stack, M[:len(stack)] + else: + M[len(stack)] = F.zero() + i += 1 + if not stack: + return + i = stack.pop() + M[len(stack)] = F.zero() + i += 1 diff --git a/surface_dynamics/misc/linalg.py.orig b/surface_dynamics/misc/linalg.py.orig new file mode 100644 index 00000000..87eb90b6 --- /dev/null +++ b/surface_dynamics/misc/linalg.py.orig @@ -0,0 +1,458 @@ +r""" +Some linear algebra routines +""" +#***************************************************************************** +# Copyright (C) 2019-2021 Vincent Delecroix <20100.delecroix@gmail.com> +# +# Distributed under the terms of the GNU General Public License (GPL) +# as published by the Free Software Foundation; either version 2 of +# the License, or (at your option) any later version. +# https://www.gnu.org/licenses/ +#***************************************************************************** + +<<<<<<< HEAD +from sage.all import ZZ, QQ, vector, Compositions, gcd, lcm, Polyhedron +from sage.numerical.mip import MixedIntegerLinearProgram, MIPSolverException +======= +from sage.arith.all import gcd, lcm +from sage.arith.misc import binomial +from sage.rings.all import ZZ, QQ +from sage.geometry.polyhedron.constructor import Polyhedron +>>>>>>> master + +def relation_space(v): + r""" + Relation space of the given vector ``v`` + + This is the sub vector space of `\QQ^d` given as the kernel of the map + `n \mapsto n \cdot \lambda`. The dimension is `d - rank`. + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import relation_space + sage: K. = QuadraticField(2) + sage: v3 = vector([sqrt2, 1, 1+sqrt2]) + sage: relation_space(v3) + Vector space of degree 3 and dimension 1 over Rational Field + Basis matrix: + [ 1 1 -1] + sage: v4 = vector([sqrt2, 1, 1+sqrt2, 1-sqrt2]) + sage: relation_space(v4) + Vector space of degree 4 and dimension 2 over Rational Field + Basis matrix: + [ 1 0 -1/2 1/2] + [ 0 1 -1/2 -1/2] + + sage: v = vector([1,2,5,3]) + sage: relation_space(v) + Vector space of degree 4 and dimension 3 over Rational Field + Basis matrix: + [ 1 0 0 -1/3] + [ 0 1 0 -2/3] + [ 0 0 1 -5/3] + + The relation space has some covariance relation with respect to matrix + actions:: + + sage: m3 = matrix(3, [1,-1,0,2,-3,4,5,-2,2]) + sage: relation_space(v3 * m3) == relation_space(v3) * ~m3.transpose() + True + sage: relation_space(m3 * v3) == relation_space(v3) * ~m3 + True + + sage: m4 = matrix(4, [1,-1,0,1,2,-3,0,4,5,3,-2,2,1,1,1,1]) + sage: relation_space(v4 * m4) == relation_space(v4) * ~m4.transpose() + True + sage: relation_space(m4 * v4) == relation_space(v4) * ~m4 + True + """ + from sage.matrix.constructor import matrix + try: + m_lengths = matrix([u.vector() for u in v]) + except AttributeError: + v = [QQ.coerce(i) for i in v] + m_lengths = matrix([[i] for i in v]) + return m_lengths.left_kernel() + +def deformation_space(lengths): + r""" + Deformation space of the given ``lengths`` + + This is the smallest vector space defined over `\QQ` that contains + the vector ``lengths``. Its dimension is `rank`. + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import deformation_space + sage: K. = QuadraticField(2) + sage: v3 = vector([sqrt2, 1, 1+sqrt2]) + sage: deformation_space(v3) + Vector space of degree 3 and dimension 2 over Rational Field + Basis matrix: + [1 0 1] + [0 1 1] + sage: v4 = vector([sqrt2, 1, 1+sqrt2, 1-sqrt2]) + sage: deformation_space(v4) + Vector space of degree 4 and dimension 2 over Rational Field + Basis matrix: + [ 1 0 1 -1] + [ 0 1 1 1] + + sage: v = vector([1, 5, 2, 9]) + sage: deformation_space(v) + Vector space of degree 4 and dimension 1 over Rational Field + Basis matrix: + [1 5 2 9] + + The deformation space has some covariance relation with respect to matrix + actions:: + + sage: m3 = matrix(3, [1,-1,0,2,-3,4,5,-2,2]) + sage: deformation_space(v3 * m3) == deformation_space(v3) * m3 + True + sage: deformation_space(m3 * v3) == deformation_space(v3) * m3.transpose() + True + + sage: m4 = matrix(4, [1,-1,0,1,2,-3,0,4,5,3,-2,2,1,1,1,1]) + sage: deformation_space(v4 * m4) == deformation_space(v4) * m4 + True + sage: deformation_space(m4 * v4) == deformation_space(v4) * m4.transpose() + True + """ + from sage.matrix.constructor import matrix + try: + m_lengths = matrix([u.vector() for u in lengths]) + except AttributeError: + lengths = [QQ.coerce(i) for i in lengths] + m_lengths = matrix([[i] for i in lengths]) + + return m_lengths.column_space() + +def pivots(n, length, force_first=False): + r""" + Iterator through increasing sequences of given ``length` in ``{0, 1, ..., n-1}``. + + INPUT: + + - ``force_first`` - whether the first pivot is on the first column + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import pivots + sage: list(pivots(4, 3, force_first=False)) + [[0, 1, 2], [0, 1, 3], [0, 2, 3], [1, 2, 3]] + sage: list(pivots(4, 3, force_first=True)) + [[0, 1, 2], [0, 1, 3], [0, 2, 3]] + """ + if force_first: + for gaps in reversed(Compositions(n, length=length)): + gaps = [1] + list(gaps) + pivots = [gaps[0] - 1] + for i in gaps[1:-1]: + pivots.append(pivots[-1] + i) + yield pivots + else: + for gaps in reversed(Compositions(n + 1, length=length+1)): + pivots = [gaps[0] - 1] + for i in gaps[1:-1]: + pivots.append(pivots[-1] + i) + yield pivots + +def isotropic_subspaces(B, d, bound=1, contains_positive_vector=False): + r""" + Iterator over rational isotropic subspaces of dimension ``d`` of the + anti-symmetric bilinear form ``B``. + + INPUT: + + - ``B`` -- anti-symmetric matrix with integral coefficients + + - ``d`` -- dimension of the isotropic subspace + + - ``bound`` (optional, default 1) -- bound on the coefficients for generators (in echelon + form) + + - ``contains_positive_vector`` (optional, default ``False``) - if ``True``, + then only iterate through subspaces that contain a positive vector + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import isotropic_subspaces + + sage: B = matrix(4, [0, 1, 1, 1, -1, 0, 1, 1, -1, -1, 0, 1, -1, -1, -1, 0]) + sage: for vectors in isotropic_subspaces(B, 2, bound=2, contains_positive_vector=True): + ....: print(matrix(vectors)) + ....: print('*' * 15) + [1 0 0 1] + [0 1 1 0] + *************** + [1 0 0 1] + [0 1 2 0] + *************** + [1 0 0 2] + [0 1 1 2] + *************** + [ 1 0 1 0] + [ 0 1 -2 1] + *************** + [1 0 1 2] + [0 1 0 1] + *************** + [1 0 1 2] + [0 1 2 2] + *************** + [ 1 0 1 2] + [ 0 1 -2 0] + *************** + [ 1 0 2 0] + [ 0 1 -2 1] + *************** + [1 0 2 2] + [0 1 0 1] + *************** + [ 1 0 -2 1] + [ 0 1 1 2] + *************** + [ 1 0 -2 1] + [ 0 1 2 2] + *************** + [ 1 0 -2 2] + [ 0 1 1 0] + *************** + [ 1 0 -2 2] + [ 0 1 2 -1] + *************** + + Some countings:: + + sage: B = matrix(4, [0,1,1,1,-1,0,1,1,-1,-1,0,1,-1,-1,-1,0]) + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=1, contains_positive_vector=True)) + 1 + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=2, contains_positive_vector=True)) + 13 + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=1, contains_positive_vector=False)) + 10 + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=2, contains_positive_vector=False)) + 66 + + sage: B = matrix(5, [0,1,1,1,1,-1,0,1,1,1,-1,-1,0,1,1,-1,-1,-1,0,1,-1,-1,-1,-1,0]) + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=1, contains_positive_vector=True)) + 17 + sage: sum(1 for _ in isotropic_subspaces(B, 3, bound=1, contains_positive_vector=True)) + 0 + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=2, contains_positive_vector=True)) + 204 + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=1, contains_positive_vector=False)) + 175 + sage: sum(1 for _ in isotropic_subspaces(B, 3, bound=1, contains_positive_vector=False)) + 2 + sage: sum(1 for _ in isotropic_subspaces(B, 2, bound=2, contains_positive_vector=False)) + 1649 + """ + assert B.is_square() + n = B.nrows() + + for piv in pivots(n, d, force_first=contains_positive_vector): + positions_to_be_filled = [] + for i in range(d): + positions_to_be_filled.append(list(enumerate(sorted(set(range(piv[i]+1, n)).difference(piv[i+1:]))))) + + # backtracking + # the ordering on entries is 0, 1, -1, 2, -2, 3, -3, .... + vectors = [vector(ZZ, n) for _ in range(d)] + for i, j in enumerate(piv): + vectors[i][j] = 1 + products = [vectors[i] * B for i in range(d)] + + piv.append(n) + + if (all(products[i].dot_product(vectors[j]).is_zero() for i in range(d) for j in range(i)) and + (not contains_positive_vector or has_positive_linear_combination(vectors, n))): + yield vectors + + i = 1 + while True: + # find the next vector at position i + v = vectors[i] + found_next = False + for jpos, j in reversed(positions_to_be_filled[i]): + must_be_positive = False + if contains_positive_vector and piv[i] < j < piv[i+1]: + # v is the last vector with non-zero entry in this column, it must compensate negative coefficients + if all(vectors[ii][j] <= 0 for ii in range(i)): + must_be_positive = True + if v[j] == -bound: + v[j] = 0 + elif v[j] > 0: + if must_be_positive: + if v[j] == bound: + v[j] = 0 + else: + assert v[j] >= 0 + v[j] += 1 + found_next = True + break + else: + v[j] = -v[j] + found_next = True + break + else: + v[j] = -v[j] + 1 + found_next = True + break + + # backtrack ? + if not found_next: + i -= 1 + if i == -1: + break + continue + + # go for the next vector if span(vectors) has positive elements and v[i] is orthogonal to the i-1 first vectors + products[i] = v * B + if all(products[i].dot_product(vectors[j]).is_zero() for j in range(i)): + if i == d - 1: + if not contains_positive_vector or has_positive_linear_combination(vectors, n): + yield vectors + elif not contains_positive_vector or has_positive_linear_combination(vectors[:i+1], piv[i+1]): + i = i + 1 + +def deformation_cone(v): + r""" + Return the deformation cone of the given vector ``v`` + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import deformation_cone + sage: K. = QuadraticField(2) + sage: v3 = vector([sqrt2, 1, 1+sqrt2]) + sage: P = deformation_cone(v3) + sage: P + A 2-dimensional polyhedron in QQ^3 defined as the convex hull of 1 vertex and 2 rays + sage: P.rays_list() + [[1, 0, 1], [0, 1, 1]] + """ + V = deformation_space(v) + P = Polyhedron(lines=deformation_space(v).basis()) + B = Polyhedron(rays=(QQ**V.degree()).basis()) + return P.intersection(B) + +def cone_triangulate(C, hyperplane=None): + r""" + Triangulation of rational cone contained in the positive quadrant. + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import cone_triangulate + sage: P = Polyhedron(rays=[(1,0,0),(0,1,0),(1,0,1),(0,1,1)]) + sage: list(cone_triangulate(P)) # random + [[(0, 1, 1), (0, 1, 0), (1, 0, 0)], [(0, 1, 1), (1, 0, 1), (1, 0, 0)]] + sage: len(_) + 2 + + sage: rays = [(0, 1, 0, -1, 0, 0), + ....: (1, 0, -1, 0, 0, -1), + ....: (0, 1, -1, 0, 0, -1), + ....: (0, 0, 1, 0, 0, 0), + ....: (0, 0, 0, 1, 0, -1), + ....: (1, -1, 0, 0, 1, -1), + ....: (0, 0, 0, 0, 1, -1), + ....: (0, 0, 1, -1, 1, 0), + ....: (0, 0, 1, -1, 0, 0), + ....: (0, 0, 1, 0, -1, 0), + ....: (0, 0, 0, 1, -1, -1), + ....: (1, -1, 0, 0, 0, -1), + ....: (0, 0, 0, 0, 0, -1)] + sage: P = Polyhedron(rays=rays) + sage: list(cone_triangulate(P, hyperplane=(1, 2, 3, -1, 0, -5))) # random + [[(0, 0, 0, 0, 0, -1), + (0, 0, 0, 0, 1, -1), + (0, 0, 0, 1, -1, -1), + (0, 0, 1, 0, 0, 0), + (0, 1, -1, 0, 0, -1), + (1, -1, 0, 0, 1, -1)], + ... + (0, 0, 1, 0, 0, 0), + (0, 1, -1, 0, 0, -1), + (0, 1, 0, -1, 0, 0), + (1, -1, 0, 0, 1, -1), + (1, 0, -1, 0, 0, -1)]] + sage: len(_) + 16 + """ + rays = [r.vector() for r in C.rays()] + dim = len(rays[0]) + if hyperplane is None: + hyperplane = [1] * dim + scalings = [sum(x*h for x,h in zip(r, hyperplane)) for r in rays] + assert all(s > 0 for s in scalings) + normalized_rays = [r / s for r,s in zip(rays, scalings)] + P = Polyhedron(vertices=normalized_rays) + for t in P.triangulate(): + simplex = [P.Vrepresentation(i).vector() for i in t] + yield [(r / gcd(r)).change_ring(ZZ) for r in simplex] + +<<<<<<< HEAD +def has_positive_linear_combination(vectors, d, solver="PPL"): + r""" + Test whether the given ``vectors`` admit a linear combination that is positive. + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import has_positive_linear_combination + sage: has_positive_linear_combination([(1,0), (0,1)], 2) + True + sage: has_positive_linear_combination([(1,-1)], 2) + False + sage: vectors = [(1, 0, 0, 0, 0, 0), (0, 1, 0, 0, 0, -1), (0, 0, 1, 0, -1, 0)] + sage: has_positive_linear_combination(vectors, 3) + True + sage: has_positive_linear_combination(vectors, 4) + False + """ + M = MixedIntegerLinearProgram(solver=solver) + x = M.new_variable() + for i in range(d): + M.add_constraint(M.sum(x[j] * v[i] for j, v in enumerate(vectors)) >= 1) + try: + M.solve(objective_only=True) + except MIPSolverException: + return False + return True +======= +def symbolic_matrix_power(M, n): + r""" + Return the symbolic power ``M^n`` of the unipotent matrix ``M``. + + EXAMPLES:: + + sage: from surface_dynamics.misc.linalg import symbolic_matrix_power + sage: m = matrix(3, [1,1,1,0,1,1,0,0,1]) + sage: n = polygen(QQ, 'n') + sage: symbolic_matrix_power(m, n) + [ 1 n 1/2*n^2 + 1/2*n] + [ 0 1 n] + [ 0 0 1] + + sage: m = matrix(2, [2,1,1,1]) + sage: symbolic_matrix_power(m, n) + Traceback (most recent call last): + ... + NotImplementedError: power only implemented for unipotent matrices + """ + d = M.nrows() + I = M.parent().identity_matrix() + N = M - M.parent().identity_matrix() + char = N.charpoly() + if any(char[i] for i in range(d)): + raise NotImplementedError('power only implemented for unipotent matrices') + + result = I + P = N + p = 1 + while P: + result += binomial(n, p) * P + P *= N + p += 1 + return result +>>>>>>> master diff --git a/tests/test_generalized_multiple_zeta_values.py b/tests/test_generalized_multiple_zeta_values.py index d8458430..39c8ea39 100644 --- a/tests/test_generalized_multiple_zeta_values.py +++ b/tests/test_generalized_multiple_zeta_values.py @@ -279,3 +279,6 @@ def expo_to_den_tuple(e): print("non-convergent relation: {} {}".format(ee,r)) else: assert z1.is_zero() and z2.is_zero(), (z1, z2, ee, r) + +def test_strategies_d3_r3(): + S1 = SolveStrategy(3,3) From 7d3cd8fe0aec76f177fb0084f022cc06a03b749f Mon Sep 17 00:00:00 2001 From: Vincent Delecroix Date: Tue, 24 Jan 2023 22:54:49 +0100 Subject: [PATCH 2/5] fix many doctests related to gmzv --- .../generalized_multiple_zeta_values.py | 53 +++++++++++++++---- .../gmzv_generic_reduction.py | 35 ++---------- .../gmzv_three_variables.py | 4 +- .../gmzv_two_variables.py | 1 + 4 files changed, 50 insertions(+), 43 deletions(-) diff --git a/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py b/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py index d6fe26f4..a89da801 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py +++ b/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py @@ -405,10 +405,9 @@ def convergent_multizeta(t): if VERBOSE: print("convergent_multizeta({})".format(t)) n = len(t) - if not DIVERGENT_MZV: - for i in range(1,n+1): - if sum(t[-i:]) <= i: - raise DivergentZetaError("divergent multizeta value {}".format(t)) + for i in range(1,n+1): + if sum(t[-i:]) <= i: + raise DivergentZetaError("divergent multizeta value {}".format(t)) if all(x > 0 for x in t): M = Multizetas(QQ) W = M.basis().keys() @@ -629,7 +628,7 @@ def kill_relation(n, den_tuple, i, relation): sage: vb = V((0,1)); vb.set_immutable() sage: vc = V((1,1)); vc.set_immutable() - sage: den_tuple = ((va,2),(vb,3),(vc,4)) + sage: den_tuple = ((va,2), (vb,3), (vc,4)) sage: kill_relation(2, den_tuple, 2, [1,1,0]) [(1, (((0, 1), 3), ((1, 1), 6))), (2, (((0, 1), 2), ((1, 1), 7))), @@ -639,7 +638,8 @@ def kill_relation(n, den_tuple, i, relation): """ assert len(relation) == len(den_tuple) assert 0 <= i < len(den_tuple) - assert sum(relation[j] * den_tuple[j][0] for j in range(len(den_tuple))) == den_tuple[i][0] + s = sum(relation[j] * den_tuple[j][0] for j in range(len(den_tuple))) + assert s == den_tuple[i][0], (s, den_tuple[i][0]) D = {den_tuple: QQ.one()} todo = [den_tuple] while todo: @@ -662,6 +662,37 @@ def kill_relation(n, den_tuple, i, relation): return [(coeff, new_den_tuple) for new_den_tuple, coeff in D.items()] +def try_relation(n, den_tuple): + r""" + Assuming that ``den_tuple`` has full rank, make it so that it has only ``n`` columns. + + EXAMPLES:: + + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import try_relation + sage: V = FreeModule(ZZ, 2) + sage: va = V((1,0)); va.set_immutable() + sage: vb = V((0,1)); vb.set_immutable() + sage: vc = V((1,1)); vc.set_immutable() + + sage: den_tuple = ((va,2), (vb,2), (vc,2)) + sage: for r in try_relation(2, den_tuple): + ....: print(r) + [(1, (((0, 1), 2), ((1, 1), 4))), (1, (((1, 0), 2), ((1, 1), 4))), (2, (((0, 1), 1), ((1, 1), 5))), (2, (((1, 0), 1), ((1, 1), 5)))] + """ + # assume it is full rank + if len(den_tuple) <= n: + return + M = matrix([v for v,p in den_tuple]) + for relation in sorted(M.left_kernel().basis(), key=lambda x: sum(bool(cc) for cc in x)): + if sum(x < 0 for x in relation) > sum(x > 0 for x in relation): + relation = -relation + for i in range(len(den_tuple)): + if relation[i] < 0: + relation /= -relation[i] + relation[i] = 0 + yield kill_relation(n, den_tuple, i, relation) + + def has_term_sum_of_smaller_terms(n, den_tuple): r""" Look for a vector v_i and {v_j} with sum(v_j) = v_i @@ -850,7 +881,7 @@ class GeneralizedMultipleZetaFunction: sage: f = GeneralizedMultipleZetaFunction([va, vb, vc, vg]) sage: f GeneralizedMultipleZetaFunction('(0)(1)(2)(0+1+2)') - sage: f((1, 1, 1, 1)) + sage: f((1, 1, 1, 1)) # optional - mzv 6*ζ(1,1,2) """ def __init__(self, *args, as_rows=False, reduced=None): @@ -910,13 +941,15 @@ def __call__(self, e): EXAMPLES:: sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, convergent_multizeta - sage: GeneralizedMultipleZetaFunction([[1,0,0],[1,1,1],[1,0,1]])((2,4,3)) + sage: f = GeneralizedMultipleZetaFunction([[1,0,0],[1,1,1],[1,0,1]]) + sage: f((2,4,3)) # optional - mzv ζ(2,3,4) - sage: GeneralizedMultipleZetaFunction([[1,0,0],[1,1,1]])((2,4)) == convergent_multizeta((2,0,4)) + sage: f = GeneralizedMultipleZetaFunction([[1,0,0],[1,1,1]]) + sage: f((2,4)) == convergent_multizeta((2,0,4)) # optional - mzv True sage: G = GeneralizedMultipleZetaFunction([[1,1,1,1],[1,1,1,0],[1,1,0,0]]) - sage: G((2,2,2)) + sage: G((2,2,2)) # optional - mzv ζ(1,2,2) - ζ(2,2,2) """ M = Multizetas(QQ) diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py index 72e28582..47d49504 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_generic_reduction.py @@ -164,13 +164,12 @@ def stuffle_to_lin_comb(v, rp, cp, stuffle): EXAMPLES:: - sage: from surface_dynamics.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, linear_forms + sage: from surface_dynamics.generalized_multiple_zeta_values.generalized_multiple_zeta_values import GeneralizedMultipleZetaFunction, linear_forms sage: from surface_dynamics.generalized_multiple_zeta_values.gmzv_generic_reduction import stuffle_to_lin_comb sage: va, vb, vc, vd, ve, vf, vg = linear_forms(3) - sage: GeneralizedMultipleZetaFunction([vd, va, vb]) + sage: G = GeneralizedMultipleZetaFunction([vd, va, vb]) sage: S = SymmetricGroup(3) sage: p = S.one() - sage: stuffle_to_lin_comb """ vtmp = v.copy() _, _ = vtmp.canonicalize() @@ -690,34 +689,6 @@ def print_stuffle3(): -def try_relation(n, den_tuple): - r""" - Assuming that ``den_tuple`` has full rank, make it so that it has only ``n`` columns. - - EXAMPLES:: - - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import try_relation - sage: V = FreeModule(ZZ, 2) - sage: va = V((1,0)); va.set_immutable() - sage: vb = V((0,1)); vb.set_immutable() - sage: vc = V((1,1)); vc.set_immutable() - - sage: den_tuple = ((va,2),(vb,2),(vc,2)) - sage: for r in try_relation(2, den_tuple): - ....: print(r) - """ - # assume it is full rank - if len(den_tuple) <= n: - return - M = matrix([v for v,p in den_tuple]) - for relation in sorted(M.left_kernel().basis(), key=lambda x: sum(bool(cc) for cc in x)): - if sum(x < 0 for x in relation) > sum(x > 0 for x in relation): - relation = -relation - for i in range(len(den_tuple)): - if relation[i] < 0: - yield kill_relation(n, den_tuple, i, relation) - - def try_algebraic_forward(S): L = list(S.algebraic(True)) if not L: @@ -775,7 +746,7 @@ def multistuffle_first(d, r, verbose=False): sage: S = gmzv_solvers.multistuffle_first(3, 3) sage: l = GeneralizedMultipleZetaFunction([[1,0,0], [0,1,0], [0,0,1]]) sage: l.set_immutable() - sage: assert S.eval(l, (2,2,2)) == Multizeta(2)**3 + sage: assert S.eval(l, (2,2,2)) == Multizeta(2)**3 # optional - mzv """ S = GMZVSolver(d, r) S.update_multizetas() diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py index 29a03578..7652ffb8 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py @@ -18,6 +18,8 @@ # https://www.gnu.org/licenses/ #***************************************************************************** + +from sage.rings.integer_ring import ZZ from sage.rings.rational_field import QQ from sage.arith.misc import binomial from sage.modular.multiple_zeta import Multizetas @@ -53,7 +55,7 @@ def is_Z3_convergent(a, b, c, d, e, f, g): TESTS:: - sage: from surface_dynamics.misc.generalized_multiple_zeta_values import is_Z3_convergent + sage: from surface_dynamics.generalized_multiple_zeta_values import is_Z3_convergent Convergent examples:: diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py index f6c7b123..89a0147e 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py @@ -17,6 +17,7 @@ #***************************************************************************** +from sage.rings.integer_ring import ZZ from sage.rings.rational_field import QQ from sage.arith.misc import binomial from sage.modular.multiple_zeta import Multizetas From 805c7fd6093d499cb2d4d1340f63c9eeded3645b Mon Sep 17 00:00:00 2001 From: Vincent Delecroix Date: Wed, 25 Jan 2023 08:06:33 +0100 Subject: [PATCH 3/5] adding gmzv to setup.py --- setup.py | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.py b/setup.py index 45238db5..c35b5c74 100755 --- a/setup.py +++ b/setup.py @@ -123,6 +123,7 @@ 'surface_dynamics/topological_recursion', 'surface_dynamics/flat_surfaces', 'surface_dynamics/databases', + 'surface_dynamics/generalized_multiple_zeta_values', 'surface_dynamics/flat_surfaces/origamis', 'surface_dynamics/interval_exchanges'], package_data={ From 11ad3ac1987f61d82d5312b24ec772295a7f4ecf Mon Sep 17 00:00:00 2001 From: Vincent Delecroix Date: Wed, 25 Jan 2023 09:34:33 +0100 Subject: [PATCH 4/5] remove unused file --- .../misc/divergent_multiple_zeta_values.py | 128 ------------------ 1 file changed, 128 deletions(-) delete mode 100644 surface_dynamics/misc/divergent_multiple_zeta_values.py diff --git a/surface_dynamics/misc/divergent_multiple_zeta_values.py b/surface_dynamics/misc/divergent_multiple_zeta_values.py deleted file mode 100644 index 96da9eb6..00000000 --- a/surface_dynamics/misc/divergent_multiple_zeta_values.py +++ /dev/null @@ -1,128 +0,0 @@ -from sage.modular.multiple_zeta import Multizetas - -M = Multizetas(QQ) - -# TODO: without this, it is not possible to create polynomials with Multizeta coefficients -M.is_prime_field = lambda: False - -MT = M['T'] # shuffle regularization -MX = M['X'] # stuffle regularization -T = MT.gen() -X = MX.gen() - -# concatenation product -def prefix(M, p, g): - p = M.basis().keys()(p) - ans = M.zero() - for u, coeff in g.monomial_coefficients().items(): - ans += coeff * M.term(p + u) - return ans - -def suffix(M, f, s): - s = M.basis().keys()(s) - ans = M.zero() - for u , coeff in f.monomial_coefficients().items(): - ans += coeff * M.term(u + s) - return ans - -def concatenate(M, left, right): - r""" - Concatenation product in the free algebra - - EXAMPLES:: - - sage: concatenate(M, M((2,)) + M((3,)), M((5,)) + M((7,))) - ζ(2,5) + ζ(2,7) + ζ(3,5) + ζ(3,7) - """ - ans = M.zero() - for u1, c1 in left.monomial_coefficients(): - for u2, c2 in right.monomial_coefficients(): - ans += c1 * c2 * M.term(u1 + u2) - return ans - -def stuffle_on_basis(M, w1, w2): - r""" - Return the stuffle product of two words - - EXAMPLES:: - - sage: stuffle_on_basis(M, Word([2]), Word([5])) - ζ(2,5) + ζ(5,2) + ζ(7) - sage: stuffle_on_basis(M, Word([2,2]), Word([3,4])) - ζ(2,2,3,4) + ζ(2,3,2,4) + ζ(2,3,4,2) + ζ(2,3,6) + ζ(2,5,4) + ζ(3,2,2,4) + ζ(3,2,4,2) + ζ(3,2,6) + ζ(3,4,2,2) + ζ(3,6,2) + ζ(5,2,4) + ζ(5,4,2) + ζ(5,6) - """ - W = M.basis().keys() - if not w1: - return M.term(W(w2)) - if not w2: - return M.term(W(w1)) - y1 = w1[0] - u1 = w1[1:] - y2 = w2[0] - u2 = w2[1:] - if y1 == y2: - return prefix(M, W([y1]), stuffle_on_basis(M, u1, w2) + stuffle_on_basis(M, w1, u2)) + \ - prefix(M, W([y1+y2]), stuffle_on_basis(M, u1, u2)) - else: - return prefix(M, W([y1]), stuffle_on_basis(M, u1, w2)) + \ - prefix(M, W([y2]), stuffle_on_basis(M, w1, u2)) + \ - prefix(M, W([y1+y2]), stuffle_on_basis(M, u1, u2)) - -def stuffle_product(M, left, right): - ans = M.zero() - for u1, c1 in left.monomial_coefficients().items(): - for u2, c2 in right.monomial_coefficients().items(): - ans += c1 * c2 * stuffle_on_basis(M, u1, u2) - return ans - -def stuffle_power(M, a, n): - if n < 0: - raise ValueError - if n == 0: - return M.one() - apow = a - while not (n & 1): - apow = stuffle_product(M, apow, apow) - n >>= 1 - - res = apow - n >>= 1 - while n: - apow = stuffle_product(M, apow, apow) - if n & 1: - res = stuffle_product(M, apow, res) - n >>= 1 - - return res - - -def stuffle_regularization(s): - r""" - Given a multizeta index s returns it as an expression in MX - - EXAMPLES:: - - sage: divergent_stuffle_regularization((1,)) - sage: divergent_stuffle_regularization((2,)) - sage: divergent_stuffle_regularization((2,1)) - sage: divergent_stuffle_regularization((1,1)) - sage: divergent_stuffle_regularization((2,1,1)) - sage: divergent_stuffle_regularization((1,1,1)) - """ - pass - -def shuffle_regularization(s): - r""" - Given a multizeta index s returns it as an expression in MT - """ - pass - -def shuffle_to_harmonic(n): - r""" - Image of X^n under the morphism rho: M[X] -> M[T] of Ihara-Kaneko-Zagier - """ - -def harmonic_to_shuffle(n): - r""" - Image of T^n under the morphism rho^{-1}: M[T] -> M[X] of Ihara-Kaneko-Zagier - """ From 7998058518cef7abb6438d4767107b7c4ba99902 Mon Sep 17 00:00:00 2001 From: Vincent Delecroix Date: Wed, 25 Jan 2023 09:37:17 +0100 Subject: [PATCH 5/5] catch ImportError when importing Multizetas from sage --- .../generalized_multiple_zeta_values.py | 3 ++- .../gmzv_three_variables.py | 9 ++++++++- .../gmzv_two_variables.py | 9 ++++++++- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py b/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py index a89da801..235fc8e3 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py +++ b/surface_dynamics/generalized_multiple_zeta_values/generalized_multiple_zeta_values.py @@ -26,7 +26,6 @@ import itertools from collections import defaultdict -import cypari2 from sage.misc.cachefunc import cached_method, cached_function from sage.all import ZZ, QQ, matrix, bernoulli_polynomial, prod, FreeModule @@ -38,6 +37,8 @@ from surface_dynamics.misc.linalg import linearly_independent_vectors from .options import VERBOSE, DIVERGENT_MZV +import cypari2.handle_error + try: from sage.modular.multiple_zeta import Multizetas except (ImportError, cypari2.handle_error.PariError): diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py index 7652ffb8..b7e49198 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_three_variables.py @@ -22,7 +22,14 @@ from sage.rings.integer_ring import ZZ from sage.rings.rational_field import QQ from sage.arith.misc import binomial -from sage.modular.multiple_zeta import Multizetas + +import cypari2.handle_error + +try: + from sage.modular.multiple_zeta import Multizetas +except (ImportError, cypari2.handle_error.PariError): + def Multizetas(*args, **kwds): + raise ValueError('your sage version does not support multiple zeta values') from .options import VERBOSE, DIVERGENT_MZV diff --git a/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py b/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py index 89a0147e..eea1ab49 100644 --- a/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py +++ b/surface_dynamics/generalized_multiple_zeta_values/gmzv_two_variables.py @@ -20,7 +20,14 @@ from sage.rings.integer_ring import ZZ from sage.rings.rational_field import QQ from sage.arith.misc import binomial -from sage.modular.multiple_zeta import Multizetas + +import cypari2.handle_error + +try: + from sage.modular.multiple_zeta import Multizetas +except (ImportError, cypari2.handle_error.PariError): + def Multizetas(*args, **kwds): + raise ValueError('your sage version does not support multiple zeta values') from .options import VERBOSE, DIVERGENT_MZV