Add files using upload-large-folder tool

lib/python3.12/site-packages/sympy/benchmarks/bench_discrete_log.py

@@ -0,0 +1,83 @@
+import sys
+from time import time
+from sympy.ntheory.residue_ntheory import (discrete_log,
+        _discrete_log_trial_mul, _discrete_log_shanks_steps,
+        _discrete_log_pollard_rho, _discrete_log_pohlig_hellman)
+
+
+# Cyclic group (Z/pZ)* with p prime, order p - 1 and generator g
+data_set_1 = [
+        # p, p - 1, g
+        [191, 190, 19],
+        [46639, 46638, 6],
+        [14789363, 14789362, 2],
+        [4254225211, 4254225210, 2],
+        [432751500361, 432751500360, 7],
+        [158505390797053, 158505390797052, 2],
+        [6575202655312007, 6575202655312006, 5],
+        [8430573471995353769, 8430573471995353768, 3],
+        [3938471339744997827267, 3938471339744997827266, 2],
+        [875260951364705563393093, 875260951364705563393092, 5],
+    ]
+
+
+# Cyclic sub-groups of (Z/nZ)* with prime order p and generator g
+# (n, p are primes and n = 2 * p + 1)
+data_set_2 = [
+        # n, p, g
+        [227, 113, 3],
+        [2447, 1223, 2],
+        [24527, 12263, 2],
+        [245639, 122819, 2],
+        [2456747, 1228373, 3],
+        [24567899, 12283949, 3],
+        [245679023, 122839511, 2],
+        [2456791307, 1228395653, 3],
+        [24567913439, 12283956719, 2],
+        [245679135407, 122839567703, 2],
+        [2456791354763, 1228395677381, 3],
+        [24567913550903, 12283956775451, 2],
+        [245679135509519, 122839567754759, 2],
+    ]
+
+
+# Cyclic sub-groups of (Z/nZ)* with smooth order o and generator g
+data_set_3 = [
+        # n, o, g
+        [2**118, 2**116, 3],
+    ]
+
+
+def bench_discrete_log(data_set, algo=None):
+    if algo is None:
+        f = discrete_log
+    elif algo == 'trial':
+        f = _discrete_log_trial_mul
+    elif algo == 'shanks':
+        f = _discrete_log_shanks_steps
+    elif algo == 'rho':
+        f = _discrete_log_pollard_rho
+    elif algo == 'ph':
+        f = _discrete_log_pohlig_hellman
+    else:
+        raise ValueError("Argument 'algo' should be one"
+                " of ('trial', 'shanks', 'rho' or 'ph')")
+
+    for i, data in enumerate(data_set):
+        for j, (n, p, g) in enumerate(data):
+            t = time()
+            l = f(n, pow(g, p - 1, n), g, p)
+            t = time() - t
+            print('[%02d-%03d] %15.10f' % (i, j, t))
+            assert l == p - 1
+
+
+if __name__ == '__main__':
+    algo = sys.argv[1] \
+            if len(sys.argv) > 1 else None
+    data_set = [
+            data_set_1,
+            data_set_2,
+            data_set_3,
+        ]
+    bench_discrete_log(data_set, algo)

lib/python3.12/site-packages/sympy/benchmarks/bench_meijerint.py

@@ -0,0 +1,261 @@
+# conceal the implicit import from the code quality tester
+from sympy.core.numbers import (oo, pi)
+from sympy.core.symbol import (Symbol, symbols)
+from sympy.functions.elementary.exponential import exp
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.special.bessel import besseli
+from sympy.functions.special.gamma_functions import gamma
+from sympy.integrals.integrals import integrate
+from sympy.integrals.transforms import (mellin_transform,
+    inverse_fourier_transform, inverse_mellin_transform,
+    laplace_transform, inverse_laplace_transform, fourier_transform)
+
+LT = laplace_transform
+FT = fourier_transform
+MT = mellin_transform
+IFT = inverse_fourier_transform
+ILT = inverse_laplace_transform
+IMT = inverse_mellin_transform
+
+from sympy.abc import x, y
+nu, beta, rho = symbols('nu beta rho')
+
+apos, bpos, cpos, dpos, posk, p = symbols('a b c d k p', positive=True)
+k = Symbol('k', real=True)
+negk = Symbol('k', negative=True)
+
+mu1, mu2 = symbols('mu1 mu2', real=True, nonzero=True, finite=True)
+sigma1, sigma2 = symbols('sigma1 sigma2', real=True, nonzero=True,
+                         finite=True, positive=True)
+rate = Symbol('lambda', positive=True)
+
+
+def normal(x, mu, sigma):
+    return 1/sqrt(2*pi*sigma**2)*exp(-(x - mu)**2/2/sigma**2)
+
+
+def exponential(x, rate):
+    return rate*exp(-rate*x)
+alpha, beta = symbols('alpha beta', positive=True)
+betadist = x**(alpha - 1)*(1 + x)**(-alpha - beta)*gamma(alpha + beta) \
+    /gamma(alpha)/gamma(beta)
+kint = Symbol('k', integer=True, positive=True)
+chi = 2**(1 - kint/2)*x**(kint - 1)*exp(-x**2/2)/gamma(kint/2)
+chisquared = 2**(-k/2)/gamma(k/2)*x**(k/2 - 1)*exp(-x/2)
+dagum = apos*p/x*(x/bpos)**(apos*p)/(1 + x**apos/bpos**apos)**(p + 1)
+d1, d2 = symbols('d1 d2', positive=True)
+f = sqrt(((d1*x)**d1 * d2**d2)/(d1*x + d2)**(d1 + d2))/x \
+    /gamma(d1/2)/gamma(d2/2)*gamma((d1 + d2)/2)
+nupos, sigmapos = symbols('nu sigma', positive=True)
+rice = x/sigmapos**2*exp(-(x**2 + nupos**2)/2/sigmapos**2)*besseli(0, x*
+                         nupos/sigmapos**2)
+mu = Symbol('mu', real=True)
+laplace = exp(-abs(x - mu)/bpos)/2/bpos
+
+u = Symbol('u', polar=True)
+tpos = Symbol('t', positive=True)
+
+
+def E(expr):
+    integrate(expr*exponential(x, rate)*normal(y, mu1, sigma1),
+                     (x, 0, oo), (y, -oo, oo), meijerg=True)
+    integrate(expr*exponential(x, rate)*normal(y, mu1, sigma1),
+                     (y, -oo, oo), (x, 0, oo), meijerg=True)
+
+bench = [
+    'MT(x**nu*Heaviside(x - 1), x, s)',
+    'MT(x**nu*Heaviside(1 - x), x, s)',
+    'MT((1-x)**(beta - 1)*Heaviside(1-x), x, s)',
+    'MT((x-1)**(beta - 1)*Heaviside(x-1), x, s)',
+    'MT((1+x)**(-rho), x, s)',
+    'MT(abs(1-x)**(-rho), x, s)',
+    'MT((1-x)**(beta-1)*Heaviside(1-x) + a*(x-1)**(beta-1)*Heaviside(x-1), x, s)',
+    'MT((x**a-b**a)/(x-b), x, s)',
+    'MT((x**a-bpos**a)/(x-bpos), x, s)',
+    'MT(exp(-x), x, s)',
+    'MT(exp(-1/x), x, s)',
+    'MT(log(x)**4*Heaviside(1-x), x, s)',
+    'MT(log(x)**3*Heaviside(x-1), x, s)',
+    'MT(log(x + 1), x, s)',
+    'MT(log(1/x + 1), x, s)',
+    'MT(log(abs(1 - x)), x, s)',
+    'MT(log(abs(1 - 1/x)), x, s)',
+    'MT(log(x)/(x+1), x, s)',
+    'MT(log(x)**2/(x+1), x, s)',
+    'MT(log(x)/(x+1)**2, x, s)',
+    'MT(erf(sqrt(x)), x, s)',
+
+    'MT(besselj(a, 2*sqrt(x)), x, s)',
+    'MT(sin(sqrt(x))*besselj(a, sqrt(x)), x, s)',
+    'MT(cos(sqrt(x))*besselj(a, sqrt(x)), x, s)',
+    'MT(besselj(a, sqrt(x))**2, x, s)',
+    'MT(besselj(a, sqrt(x))*besselj(-a, sqrt(x)), x, s)',
+    'MT(besselj(a - 1, sqrt(x))*besselj(a, sqrt(x)), x, s)',
+    'MT(besselj(a, sqrt(x))*besselj(b, sqrt(x)), x, s)',
+    'MT(besselj(a, sqrt(x))**2 + besselj(-a, sqrt(x))**2, x, s)',
+    'MT(bessely(a, 2*sqrt(x)), x, s)',
+    'MT(sin(sqrt(x))*bessely(a, sqrt(x)), x, s)',
+    'MT(cos(sqrt(x))*bessely(a, sqrt(x)), x, s)',
+    'MT(besselj(a, sqrt(x))*bessely(a, sqrt(x)), x, s)',
+    'MT(besselj(a, sqrt(x))*bessely(b, sqrt(x)), x, s)',
+    'MT(bessely(a, sqrt(x))**2, x, s)',
+
+    'MT(besselk(a, 2*sqrt(x)), x, s)',
+    'MT(besselj(a, 2*sqrt(2*sqrt(x)))*besselk(a, 2*sqrt(2*sqrt(x))), x, s)',
+    'MT(besseli(a, sqrt(x))*besselk(a, sqrt(x)), x, s)',
+    'MT(besseli(b, sqrt(x))*besselk(a, sqrt(x)), x, s)',
+    'MT(exp(-x/2)*besselk(a, x/2), x, s)',
+
+    # later: ILT, IMT
+
+    'LT((t-apos)**bpos*exp(-cpos*(t-apos))*Heaviside(t-apos), t, s)',
+    'LT(t**apos, t, s)',
+    'LT(Heaviside(t), t, s)',
+    'LT(Heaviside(t - apos), t, s)',
+    'LT(1 - exp(-apos*t), t, s)',
+    'LT((exp(2*t)-1)*exp(-bpos - t)*Heaviside(t)/2, t, s, noconds=True)',
+    'LT(exp(t), t, s)',
+    'LT(exp(2*t), t, s)',
+    'LT(exp(apos*t), t, s)',
+    'LT(log(t/apos), t, s)',
+    'LT(erf(t), t, s)',
+    'LT(sin(apos*t), t, s)',
+    'LT(cos(apos*t), t, s)',
+    'LT(exp(-apos*t)*sin(bpos*t), t, s)',
+    'LT(exp(-apos*t)*cos(bpos*t), t, s)',
+    'LT(besselj(0, t), t, s, noconds=True)',
+    'LT(besselj(1, t), t, s, noconds=True)',
+
+    'FT(Heaviside(1 - abs(2*apos*x)), x, k)',
+    'FT(Heaviside(1-abs(apos*x))*(1-abs(apos*x)), x, k)',
+    'FT(exp(-apos*x)*Heaviside(x), x, k)',
+    'IFT(1/(apos + 2*pi*I*x), x, posk, noconds=False)',
+    'IFT(1/(apos + 2*pi*I*x), x, -posk, noconds=False)',
+    'IFT(1/(apos + 2*pi*I*x), x, negk)',
+    'FT(x*exp(-apos*x)*Heaviside(x), x, k)',
+    'FT(exp(-apos*x)*sin(bpos*x)*Heaviside(x), x, k)',
+    'FT(exp(-apos*x**2), x, k)',
+    'IFT(sqrt(pi/apos)*exp(-(pi*k)**2/apos), k, x)',
+    'FT(exp(-apos*abs(x)), x, k)',
+
+    'integrate(normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
+    'integrate(x*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
+    'integrate(x**2*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
+    'integrate(x**3*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
+    'integrate(normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '          (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate(x*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '          (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate(y*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '          (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate(x*y*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '          (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate((x+y+1)*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '          (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate((x+y-1)*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '                   (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate(x**2*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '                (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate(y**2*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
+    '          (x, -oo, oo), (y, -oo, oo), meijerg=True)',
+    'integrate(exponential(x, rate), (x, 0, oo), meijerg=True)',
+    'integrate(x*exponential(x, rate), (x, 0, oo), meijerg=True)',
+    'integrate(x**2*exponential(x, rate), (x, 0, oo), meijerg=True)',
+    'E(1)',
+    'E(x*y)',
+    'E(x*y**2)',
+    'E((x+y+1)**2)',
+    'E(x+y+1)',
+    'E((x+y-1)**2)',
+    'integrate(betadist, (x, 0, oo), meijerg=True)',
+    'integrate(x*betadist, (x, 0, oo), meijerg=True)',
+    'integrate(x**2*betadist, (x, 0, oo), meijerg=True)',
+    'integrate(chi, (x, 0, oo), meijerg=True)',
+    'integrate(x*chi, (x, 0, oo), meijerg=True)',
+    'integrate(x**2*chi, (x, 0, oo), meijerg=True)',
+    'integrate(chisquared, (x, 0, oo), meijerg=True)',
+    'integrate(x*chisquared, (x, 0, oo), meijerg=True)',
+    'integrate(x**2*chisquared, (x, 0, oo), meijerg=True)',
+    'integrate(((x-k)/sqrt(2*k))**3*chisquared, (x, 0, oo), meijerg=True)',
+    'integrate(dagum, (x, 0, oo), meijerg=True)',
+    'integrate(x*dagum, (x, 0, oo), meijerg=True)',
+    'integrate(x**2*dagum, (x, 0, oo), meijerg=True)',
+    'integrate(f, (x, 0, oo), meijerg=True)',
+    'integrate(x*f, (x, 0, oo), meijerg=True)',
+    'integrate(x**2*f, (x, 0, oo), meijerg=True)',
+    'integrate(rice, (x, 0, oo), meijerg=True)',
+    'integrate(laplace, (x, -oo, oo), meijerg=True)',
+    'integrate(x*laplace, (x, -oo, oo), meijerg=True)',
+    'integrate(x**2*laplace, (x, -oo, oo), meijerg=True)',
+    'integrate(log(x) * x**(k-1) * exp(-x) / gamma(k), (x, 0, oo))',
+
+    'integrate(sin(z*x)*(x**2-1)**(-(y+S(1)/2)), (x, 1, oo), meijerg=True)',
+    'integrate(besselj(0,x)*besselj(1,x)*exp(-x**2), (x, 0, oo), meijerg=True)',
+    'integrate(besselj(0,x)*besselj(1,x)*besselk(0,x), (x, 0, oo), meijerg=True)',
+    'integrate(besselj(0,x)*besselj(1,x)*exp(-x**2), (x, 0, oo), meijerg=True)',
+    'integrate(besselj(a,x)*besselj(b,x)/x, (x,0,oo), meijerg=True)',
+
+    'hyperexpand(meijerg((-s - a/2 + 1, -s + a/2 + 1), (-a/2 - S(1)/2, -s + a/2 + S(3)/2), (a/2, -a/2), (-a/2 - S(1)/2, -s + a/2 + S(3)/2), 1))',
+    "gammasimp(S('2**(2*s)*(-pi*gamma(-a + 1)*gamma(a + 1)*gamma(-a - s + 1)*gamma(-a + s - 1/2)*gamma(a - s + 3/2)*gamma(a + s + 1)/(a*(a + s)) - gamma(-a - 1/2)*gamma(-a + 1)*gamma(a + 1)*gamma(a + 3/2)*gamma(-s + 3/2)*gamma(s - 1/2)*gamma(-a + s + 1)*gamma(a - s + 1)/(a*(-a + s)))*gamma(-2*s + 1)*gamma(s + 1)/(pi*s*gamma(-a - 1/2)*gamma(a + 3/2)*gamma(-s + 1)*gamma(-s + 3/2)*gamma(s - 1/2)*gamma(-a - s + 1)*gamma(-a + s - 1/2)*gamma(a - s + 1)*gamma(a - s + 3/2))'))",
+
+    'mellin_transform(E1(x), x, s)',
+    'inverse_mellin_transform(gamma(s)/s, s, x, (0, oo))',
+    'mellin_transform(expint(a, x), x, s)',
+    'mellin_transform(Si(x), x, s)',
+    'inverse_mellin_transform(-2**s*sqrt(pi)*gamma((s + 1)/2)/(2*s*gamma(-s/2 + 1)), s, x, (-1, 0))',
+    'mellin_transform(Ci(sqrt(x)), x, s)',
+    'inverse_mellin_transform(-4**s*sqrt(pi)*gamma(s)/(2*s*gamma(-s + S(1)/2)),s, u, (0, 1))',
+    'laplace_transform(Ci(x), x, s)',
+    'laplace_transform(expint(a, x), x, s)',
+    'laplace_transform(expint(1, x), x, s)',
+    'laplace_transform(expint(2, x), x, s)',
+    'inverse_laplace_transform(-log(1 + s**2)/2/s, s, u)',
+    'inverse_laplace_transform(log(s + 1)/s, s, x)',
+    'inverse_laplace_transform((s - log(s + 1))/s**2, s, x)',
+    'laplace_transform(Chi(x), x, s)',
+    'laplace_transform(Shi(x), x, s)',
+
+    'integrate(exp(-z*x)/x, (x, 1, oo), meijerg=True, conds="none")',
+    'integrate(exp(-z*x)/x**2, (x, 1, oo), meijerg=True, conds="none")',
+    'integrate(exp(-z*x)/x**3, (x, 1, oo), meijerg=True,conds="none")',
+    'integrate(-cos(x)/x, (x, tpos, oo), meijerg=True)',
+    'integrate(-sin(x)/x, (x, tpos, oo), meijerg=True)',
+    'integrate(sin(x)/x, (x, 0, z), meijerg=True)',
+    'integrate(sinh(x)/x, (x, 0, z), meijerg=True)',
+    'integrate(exp(-x)/x, x, meijerg=True)',
+    'integrate(exp(-x)/x**2, x, meijerg=True)',
+    'integrate(cos(u)/u, u, meijerg=True)',
+    'integrate(cosh(u)/u, u, meijerg=True)',
+    'integrate(expint(1, x), x, meijerg=True)',
+    'integrate(expint(2, x), x, meijerg=True)',
+    'integrate(Si(x), x, meijerg=True)',
+    'integrate(Ci(u), u, meijerg=True)',
+    'integrate(Shi(x), x, meijerg=True)',
+    'integrate(Chi(u), u, meijerg=True)',
+    'integrate(Si(x)*exp(-x), (x, 0, oo), meijerg=True)',
+    'integrate(expint(1, x)*sin(x), (x, 0, oo), meijerg=True)'
+]
+
+from time import time
+from sympy.core.cache import clear_cache
+import sys
+
+timings = []
+
+if __name__ == '__main__':
+    for n, string in enumerate(bench):
+        clear_cache()
+        _t = time()
+        exec(string)
+        _t = time() - _t
+        timings += [(_t, string)]
+        sys.stdout.write('.')
+        sys.stdout.flush()
+        if n % (len(bench) // 10) == 0:
+            sys.stdout.write('%s' % (10*n // len(bench)))
+    print()
+
+    timings.sort(key=lambda x: -x[0])
+
+    for ti, string in timings:
+        print('%.2fs %s' % (ti, string))

lib/python3.12/site-packages/sympy/benchmarks/bench_symbench.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python
+from sympy.core.random import random
+from sympy.core.numbers import (I, Integer, pi)
+from sympy.core.symbol import Symbol
+from sympy.core.sympify import sympify
+from sympy.functions.elementary.miscellaneous import sqrt
+from sympy.functions.elementary.trigonometric import sin
+from sympy.polys.polytools import factor
+from sympy.simplify.simplify import simplify
+from sympy.abc import x, y, z
+from timeit import default_timer as clock
+
+
+def bench_R1():
+    "real(f(f(f(f(f(f(f(f(f(f(i/2)))))))))))"
+    def f(z):
+        return sqrt(Integer(1)/3)*z**2 + I/3
+    f(f(f(f(f(f(f(f(f(f(I/2)))))))))).as_real_imag()[0]
+
+
+def bench_R2():
+    "Hermite polynomial hermite(15, y)"
+    def hermite(n, y):
+        if n == 1:
+            return 2*y
+        if n == 0:
+            return 1
+        return (2*y*hermite(n - 1, y) - 2*(n - 1)*hermite(n - 2, y)).expand()
+
+    hermite(15, y)
+
+
+def bench_R3():
+    "a = [bool(f==f) for _ in range(10)]"
+    f = x + y + z
+    [bool(f == f) for _ in range(10)]
+
+
+def bench_R4():
+    # we don't have Tuples
+    pass
+
+
+def bench_R5():
+    "blowup(L, 8); L=uniq(L)"
+    def blowup(L, n):
+        for i in range(n):
+            L.append( (L[i] + L[i + 1]) * L[i + 2] )
+
+    def uniq(x):
+        v = set(x)
+        return v
+    L = [x, y, z]
+    blowup(L, 8)
+    L = uniq(L)
+
+
+def bench_R6():
+    "sum(simplify((x+sin(i))/x+(x-sin(i))/x) for i in range(100))"
+    sum(simplify((x + sin(i))/x + (x - sin(i))/x) for i in range(100))
+
+
+def bench_R7():
+    "[f.subs(x, random()) for _ in range(10**4)]"
+    f = x**24 + 34*x**12 + 45*x**3 + 9*x**18 + 34*x**10 + 32*x**21
+    [f.subs(x, random()) for _ in range(10**4)]
+
+
+def bench_R8():
+    "right(x^2,0,5,10^4)"
+    def right(f, a, b, n):
+        a = sympify(a)
+        b = sympify(b)
+        n = sympify(n)
+        x = f.atoms(Symbol).pop()
+        Deltax = (b - a)/n
+        c = a
+        est = 0
+        for i in range(n):
+            c += Deltax
+            est += f.subs(x, c)
+        return est*Deltax
+
+    right(x**2, 0, 5, 10**4)
+
+
+def _bench_R9():
+    "factor(x^20 - pi^5*y^20)"
+    factor(x**20 - pi**5*y**20)
+
+
+def bench_R10():
+    "v = [-pi,-pi+1/10..,pi]"
+    def srange(min, max, step):
+        v = [min]
+        while (max - v[-1]).evalf() > 0:
+            v.append(v[-1] + step)
+        return v[:-1]
+    srange(-pi, pi, sympify(1)/10)
+
+
+def bench_R11():
+    "a = [random() + random()*I for w in [0..1000]]"
+    [random() + random()*I for w in range(1000)]
+
+
+def bench_S1():
+    "e=(x+y+z+1)**7;f=e*(e+1);f.expand()"
+    e = (x + y + z + 1)**7
+    f = e*(e + 1)
+    f.expand()
+
+
+if __name__ == '__main__':
+    benchmarks = [
+        bench_R1,
+        bench_R2,
+        bench_R3,
+        bench_R5,
+        bench_R6,
+        bench_R7,
+        bench_R8,
+        #_bench_R9,
+        bench_R10,
+        bench_R11,
+        #bench_S1,
+    ]
+
+    report = []
+    for b in benchmarks:
+        t = clock()
+        b()
+        t = clock() - t
+        print("%s%65s: %f" % (b.__name__, b.__doc__, t))

lib/python3.12/site-packages/sympy/categories/__init__.py

@@ -0,0 +1,33 @@
+"""
+Category Theory module.
+
+Provides some of the fundamental category-theory-related classes,
+including categories, morphisms, diagrams.  Functors are not
+implemented yet.
+
+The general reference work this module tries to follow is
+
+  [JoyOfCats] J. Adamek, H. Herrlich. G. E. Strecker: Abstract and
+              Concrete Categories. The Joy of Cats.
+
+The latest version of this book should be available for free download
+from
+
+   katmat.math.uni-bremen.de/acc/acc.pdf
+
+"""
+
+from .baseclasses import (Object, Morphism, IdentityMorphism,
+                         NamedMorphism, CompositeMorphism, Category,
+                         Diagram)
+
+from .diagram_drawing import (DiagramGrid, XypicDiagramDrawer,
+                             xypic_draw_diagram, preview_diagram)
+
+__all__ = [
+    'Object', 'Morphism', 'IdentityMorphism', 'NamedMorphism',
+    'CompositeMorphism', 'Category', 'Diagram',
+
+    'DiagramGrid', 'XypicDiagramDrawer', 'xypic_draw_diagram',
+    'preview_diagram',
+]

lib/python3.12/site-packages/sympy/categories/baseclasses.py

@@ -0,0 +1,978 @@
+from sympy.core import S, Basic, Dict, Symbol, Tuple, sympify
+from sympy.core.symbol import Str
+from sympy.sets import Set, FiniteSet, EmptySet
+from sympy.utilities.iterables import iterable
+
+
+class Class(Set):
+    r"""
+    The base class for any kind of class in the set-theoretic sense.
+
+    Explanation
+    ===========
+
+    In axiomatic set theories, everything is a class.  A class which
+    can be a member of another class is a set.  A class which is not a
+    member of another class is a proper class.  The class `\{1, 2\}`
+    is a set; the class of all sets is a proper class.
+
+    This class is essentially a synonym for :class:`sympy.core.Set`.
+    The goal of this class is to assure easier migration to the
+    eventual proper implementation of set theory.
+    """
+    is_proper = False
+
+
+class Object(Symbol):
+    """
+    The base class for any kind of object in an abstract category.
+
+    Explanation
+    ===========
+
+    While technically any instance of :class:`~.Basic` will do, this
+    class is the recommended way to create abstract objects in
+    abstract categories.
+    """
+
+
+class Morphism(Basic):
+    """
+    The base class for any morphism in an abstract category.
+
+    Explanation
+    ===========
+
+    In abstract categories, a morphism is an arrow between two
+    category objects.  The object where the arrow starts is called the
+    domain, while the object where the arrow ends is called the
+    codomain.
+
+    Two morphisms between the same pair of objects are considered to
+    be the same morphisms.  To distinguish between morphisms between
+    the same objects use :class:`NamedMorphism`.
+
+    It is prohibited to instantiate this class.  Use one of the
+    derived classes instead.
+
+    See Also
+    ========
+
+    IdentityMorphism, NamedMorphism, CompositeMorphism
+    """
+    def __new__(cls, domain, codomain):
+        raise(NotImplementedError(
+            "Cannot instantiate Morphism.  Use derived classes instead."))
+
+    @property
+    def domain(self):
+        """
+        Returns the domain of the morphism.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> f.domain
+        Object("A")
+
+        """
+        return self.args[0]
+
+    @property
+    def codomain(self):
+        """
+        Returns the codomain of the morphism.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> f.codomain
+        Object("B")
+
+        """
+        return self.args[1]
+
+    def compose(self, other):
+        r"""
+        Composes self with the supplied morphism.
+
+        The order of elements in the composition is the usual order,
+        i.e., to construct `g\circ f` use ``g.compose(f)``.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> g * f
+        CompositeMorphism((NamedMorphism(Object("A"), Object("B"), "f"),
+        NamedMorphism(Object("B"), Object("C"), "g")))
+        >>> (g * f).domain
+        Object("A")
+        >>> (g * f).codomain
+        Object("C")
+
+        """
+        return CompositeMorphism(other, self)
+
+    def __mul__(self, other):
+        r"""
+        Composes self with the supplied morphism.
+
+        The semantics of this operation is given by the following
+        equation: ``g * f == g.compose(f)`` for composable morphisms
+        ``g`` and ``f``.
+
+        See Also
+        ========
+
+        compose
+        """
+        return self.compose(other)
+
+
+class IdentityMorphism(Morphism):
+    """
+    Represents an identity morphism.
+
+    Explanation
+    ===========
+
+    An identity morphism is a morphism with equal domain and codomain,
+    which acts as an identity with respect to composition.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism, IdentityMorphism
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> id_A = IdentityMorphism(A)
+    >>> id_B = IdentityMorphism(B)
+    >>> f * id_A == f
+    True
+    >>> id_B * f == f
+    True
+
+    See Also
+    ========
+
+    Morphism
+    """
+    def __new__(cls, domain):
+        return Basic.__new__(cls, domain)
+
+    @property
+    def codomain(self):
+        return self.domain
+
+
+class NamedMorphism(Morphism):
+    """
+    Represents a morphism which has a name.
+
+    Explanation
+    ===========
+
+    Names are used to distinguish between morphisms which have the
+    same domain and codomain: two named morphisms are equal if they
+    have the same domains, codomains, and names.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> f
+    NamedMorphism(Object("A"), Object("B"), "f")
+    >>> f.name
+    'f'
+
+    See Also
+    ========
+
+    Morphism
+    """
+    def __new__(cls, domain, codomain, name):
+        if not name:
+            raise ValueError("Empty morphism names not allowed.")
+
+        if not isinstance(name, Str):
+            name = Str(name)
+
+        return Basic.__new__(cls, domain, codomain, name)
+
+    @property
+    def name(self):
+        """
+        Returns the name of the morphism.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> f.name
+        'f'
+
+        """
+        return self.args[2].name
+
+
+class CompositeMorphism(Morphism):
+    r"""
+    Represents a morphism which is a composition of other morphisms.
+
+    Explanation
+    ===========
+
+    Two composite morphisms are equal if the morphisms they were
+    obtained from (components) are the same and were listed in the
+    same order.
+
+    The arguments to the constructor for this class should be listed
+    in diagram order: to obtain the composition `g\circ f` from the
+    instances of :class:`Morphism` ``g`` and ``f`` use
+    ``CompositeMorphism(f, g)``.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism, CompositeMorphism
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> g * f
+    CompositeMorphism((NamedMorphism(Object("A"), Object("B"), "f"),
+    NamedMorphism(Object("B"), Object("C"), "g")))
+    >>> CompositeMorphism(f, g) == g * f
+    True
+
+    """
+    @staticmethod
+    def _add_morphism(t, morphism):
+        """
+        Intelligently adds ``morphism`` to tuple ``t``.
+
+        Explanation
+        ===========
+
+        If ``morphism`` is a composite morphism, its components are
+        added to the tuple.  If ``morphism`` is an identity, nothing
+        is added to the tuple.
+
+        No composability checks are performed.
+        """
+        if isinstance(morphism, CompositeMorphism):
+            # ``morphism`` is a composite morphism; we have to
+            # denest its components.
+            return t + morphism.components
+        elif isinstance(morphism, IdentityMorphism):
+            # ``morphism`` is an identity.  Nothing happens.
+            return t
+        else:
+            return t + Tuple(morphism)
+
+    def __new__(cls, *components):
+        if components and not isinstance(components[0], Morphism):
+            # Maybe the user has explicitly supplied a list of
+            # morphisms.
+            return CompositeMorphism.__new__(cls, *components[0])
+
+        normalised_components = Tuple()
+
+        for current, following in zip(components, components[1:]):
+            if not isinstance(current, Morphism) or \
+                    not isinstance(following, Morphism):
+                raise TypeError("All components must be morphisms.")
+
+            if current.codomain != following.domain:
+                raise ValueError("Uncomposable morphisms.")
+
+            normalised_components = CompositeMorphism._add_morphism(
+                normalised_components, current)
+
+        # We haven't added the last morphism to the list of normalised
+        # components.  Add it now.
+        normalised_components = CompositeMorphism._add_morphism(
+            normalised_components, components[-1])
+
+        if not normalised_components:
+            # If ``normalised_components`` is empty, only identities
+            # were supplied.  Since they all were composable, they are
+            # all the same identities.
+            return components[0]
+        elif len(normalised_components) == 1:
+            # No sense to construct a whole CompositeMorphism.
+            return normalised_components[0]
+
+        return Basic.__new__(cls, normalised_components)
+
+    @property
+    def components(self):
+        """
+        Returns the components of this composite morphism.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> (g * f).components
+        (NamedMorphism(Object("A"), Object("B"), "f"),
+        NamedMorphism(Object("B"), Object("C"), "g"))
+
+        """
+        return self.args[0]
+
+    @property
+    def domain(self):
+        """
+        Returns the domain of this composite morphism.
+
+        The domain of the composite morphism is the domain of its
+        first component.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> (g * f).domain
+        Object("A")
+
+        """
+        return self.components[0].domain
+
+    @property
+    def codomain(self):
+        """
+        Returns the codomain of this composite morphism.
+
+        The codomain of the composite morphism is the codomain of its
+        last component.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> (g * f).codomain
+        Object("C")
+
+        """
+        return self.components[-1].codomain
+
+    def flatten(self, new_name):
+        """
+        Forgets the composite structure of this morphism.
+
+        Explanation
+        ===========
+
+        If ``new_name`` is not empty, returns a :class:`NamedMorphism`
+        with the supplied name, otherwise returns a :class:`Morphism`.
+        In both cases the domain of the new morphism is the domain of
+        this composite morphism and the codomain of the new morphism
+        is the codomain of this composite morphism.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> (g * f).flatten("h")
+        NamedMorphism(Object("A"), Object("C"), "h")
+
+        """
+        return NamedMorphism(self.domain, self.codomain, new_name)
+
+
+class Category(Basic):
+    r"""
+    An (abstract) category.
+
+    Explanation
+    ===========
+
+    A category [JoyOfCats] is a quadruple `\mbox{K} = (O, \hom, id,
+    \circ)` consisting of
+
+    * a (set-theoretical) class `O`, whose members are called
+      `K`-objects,
+
+    * for each pair `(A, B)` of `K`-objects, a set `\hom(A, B)` whose
+      members are called `K`-morphisms from `A` to `B`,
+
+    * for a each `K`-object `A`, a morphism `id:A\rightarrow A`,
+      called the `K`-identity of `A`,
+
+    * a composition law `\circ` associating with every `K`-morphisms
+      `f:A\rightarrow B` and `g:B\rightarrow C` a `K`-morphism `g\circ
+      f:A\rightarrow C`, called the composite of `f` and `g`.
+
+    Composition is associative, `K`-identities are identities with
+    respect to composition, and the sets `\hom(A, B)` are pairwise
+    disjoint.
+
+    This class knows nothing about its objects and morphisms.
+    Concrete cases of (abstract) categories should be implemented as
+    classes derived from this one.
+
+    Certain instances of :class:`Diagram` can be asserted to be
+    commutative in a :class:`Category` by supplying the argument
+    ``commutative_diagrams`` in the constructor.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism, Diagram, Category
+    >>> from sympy import FiniteSet
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> d = Diagram([f, g])
+    >>> K = Category("K", commutative_diagrams=[d])
+    >>> K.commutative_diagrams == FiniteSet(d)
+    True
+
+    See Also
+    ========
+
+    Diagram
+    """
+    def __new__(cls, name, objects=EmptySet, commutative_diagrams=EmptySet):
+        if not name:
+            raise ValueError("A Category cannot have an empty name.")
+
+        if not isinstance(name, Str):
+            name = Str(name)
+
+        if not isinstance(objects, Class):
+            objects = Class(objects)
+
+        new_category = Basic.__new__(cls, name, objects,
+                                     FiniteSet(*commutative_diagrams))
+        return new_category
+
+    @property
+    def name(self):
+        """
+        Returns the name of this category.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Category
+        >>> K = Category("K")
+        >>> K.name
+        'K'
+
+        """
+        return self.args[0].name
+
+    @property
+    def objects(self):
+        """
+        Returns the class of objects of this category.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, Category
+        >>> from sympy import FiniteSet
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> K = Category("K", FiniteSet(A, B))
+        >>> K.objects
+        Class({Object("A"), Object("B")})
+
+        """
+        return self.args[1]
+
+    @property
+    def commutative_diagrams(self):
+        """
+        Returns the :class:`~.FiniteSet` of diagrams which are known to
+        be commutative in this category.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism, Diagram, Category
+        >>> from sympy import FiniteSet
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g])
+        >>> K = Category("K", commutative_diagrams=[d])
+        >>> K.commutative_diagrams == FiniteSet(d)
+        True
+
+        """
+        return self.args[2]
+
+    def hom(self, A, B):
+        raise NotImplementedError(
+            "hom-sets are not implemented in Category.")
+
+    def all_morphisms(self):
+        raise NotImplementedError(
+            "Obtaining the class of morphisms is not implemented in Category.")
+
+
+class Diagram(Basic):
+    r"""
+    Represents a diagram in a certain category.
+
+    Explanation
+    ===========
+
+    Informally, a diagram is a collection of objects of a category and
+    certain morphisms between them.  A diagram is still a monoid with
+    respect to morphism composition; i.e., identity morphisms, as well
+    as all composites of morphisms included in the diagram belong to
+    the diagram.  For a more formal approach to this notion see
+    [Pare1970].
+
+    The components of composite morphisms are also added to the
+    diagram.  No properties are assigned to such morphisms by default.
+
+    A commutative diagram is often accompanied by a statement of the
+    following kind: "if such morphisms with such properties exist,
+    then such morphisms which such properties exist and the diagram is
+    commutative".  To represent this, an instance of :class:`Diagram`
+    includes a collection of morphisms which are the premises and
+    another collection of conclusions.  ``premises`` and
+    ``conclusions`` associate morphisms belonging to the corresponding
+    categories with the :class:`~.FiniteSet`'s of their properties.
+
+    The set of properties of a composite morphism is the intersection
+    of the sets of properties of its components.  The domain and
+    codomain of a conclusion morphism should be among the domains and
+    codomains of the morphisms listed as the premises of a diagram.
+
+    No checks are carried out of whether the supplied object and
+    morphisms do belong to one and the same category.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism, Diagram
+    >>> from sympy import pprint, default_sort_key
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> d = Diagram([f, g])
+    >>> premises_keys = sorted(d.premises.keys(), key=default_sort_key)
+    >>> pprint(premises_keys, use_unicode=False)
+    [g*f:A-->C, id:A-->A, id:B-->B, id:C-->C, f:A-->B, g:B-->C]
+    >>> pprint(d.premises, use_unicode=False)
+    {g*f:A-->C: EmptySet, id:A-->A: EmptySet, id:B-->B: EmptySet,
+     id:C-->C: EmptySet, f:A-->B: EmptySet, g:B-->C: EmptySet}
+    >>> d = Diagram([f, g], {g * f: "unique"})
+    >>> pprint(d.conclusions,use_unicode=False)
+    {g*f:A-->C: {unique}}
+
+    References
+    ==========
+
+    [Pare1970] B. Pareigis: Categories and functors.  Academic Press, 1970.
+
+    """
+    @staticmethod
+    def _set_dict_union(dictionary, key, value):
+        """
+        If ``key`` is in ``dictionary``, set the new value of ``key``
+        to be the union between the old value and ``value``.
+        Otherwise, set the value of ``key`` to ``value.
+
+        Returns ``True`` if the key already was in the dictionary and
+        ``False`` otherwise.
+        """
+        if key in dictionary:
+            dictionary[key] = dictionary[key] | value
+            return True
+        else:
+            dictionary[key] = value
+            return False
+
+    @staticmethod
+    def _add_morphism_closure(morphisms, morphism, props, add_identities=True,
+                              recurse_composites=True):
+        """
+        Adds a morphism and its attributes to the supplied dictionary
+        ``morphisms``.  If ``add_identities`` is True, also adds the
+        identity morphisms for the domain and the codomain of
+        ``morphism``.
+        """
+        if not Diagram._set_dict_union(morphisms, morphism, props):
+            # We have just added a new morphism.
+
+            if isinstance(morphism, IdentityMorphism):
+                if props:
+                    # Properties for identity morphisms don't really
+                    # make sense, because very much is known about
+                    # identity morphisms already, so much that they
+                    # are trivial.  Having properties for identity
+                    # morphisms would only be confusing.
+                    raise ValueError(
+                        "Instances of IdentityMorphism cannot have properties.")
+                return
+
+            if add_identities:
+                empty = EmptySet
+
+                id_dom = IdentityMorphism(morphism.domain)
+                id_cod = IdentityMorphism(morphism.codomain)
+
+                Diagram._set_dict_union(morphisms, id_dom, empty)
+                Diagram._set_dict_union(morphisms, id_cod, empty)
+
+            for existing_morphism, existing_props in list(morphisms.items()):
+                new_props = existing_props & props
+                if morphism.domain == existing_morphism.codomain:
+                    left = morphism * existing_morphism
+                    Diagram._set_dict_union(morphisms, left, new_props)
+                if morphism.codomain == existing_morphism.domain:
+                    right = existing_morphism * morphism
+                    Diagram._set_dict_union(morphisms, right, new_props)
+
+            if isinstance(morphism, CompositeMorphism) and recurse_composites:
+                # This is a composite morphism, add its components as
+                # well.
+                empty = EmptySet
+                for component in morphism.components:
+                    Diagram._add_morphism_closure(morphisms, component, empty,
+                                                  add_identities)
+
+    def __new__(cls, *args):
+        """
+        Construct a new instance of Diagram.
+
+        Explanation
+        ===========
+
+        If no arguments are supplied, an empty diagram is created.
+
+        If at least an argument is supplied, ``args[0]`` is
+        interpreted as the premises of the diagram.  If ``args[0]`` is
+        a list, it is interpreted as a list of :class:`Morphism`'s, in
+        which each :class:`Morphism` has an empty set of properties.
+        If ``args[0]`` is a Python dictionary or a :class:`Dict`, it
+        is interpreted as a dictionary associating to some
+        :class:`Morphism`'s some properties.
+
+        If at least two arguments are supplied ``args[1]`` is
+        interpreted as the conclusions of the diagram.  The type of
+        ``args[1]`` is interpreted in exactly the same way as the type
+        of ``args[0]``.  If only one argument is supplied, the diagram
+        has no conclusions.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import IdentityMorphism, Diagram
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g])
+        >>> IdentityMorphism(A) in d.premises.keys()
+        True
+        >>> g * f in d.premises.keys()
+        True
+        >>> d = Diagram([f, g], {g * f: "unique"})
+        >>> d.conclusions[g * f]
+        {unique}
+
+        """
+        premises = {}
+        conclusions = {}
+
+        # Here we will keep track of the objects which appear in the
+        # premises.
+        objects = EmptySet
+
+        if len(args) >= 1:
+            # We've got some premises in the arguments.
+            premises_arg = args[0]
+
+            if isinstance(premises_arg, list):
+                # The user has supplied a list of morphisms, none of
+                # which have any attributes.
+                empty = EmptySet
+
+                for morphism in premises_arg:
+                    objects |= FiniteSet(morphism.domain, morphism.codomain)
+                    Diagram._add_morphism_closure(premises, morphism, empty)
+            elif isinstance(premises_arg, (dict, Dict)):
+                # The user has supplied a dictionary of morphisms and
+                # their properties.
+                for morphism, props in premises_arg.items():
+                    objects |= FiniteSet(morphism.domain, morphism.codomain)
+                    Diagram._add_morphism_closure(
+                        premises, morphism, FiniteSet(*props) if iterable(props) else FiniteSet(props))
+
+        if len(args) >= 2:
+            # We also have some conclusions.
+            conclusions_arg = args[1]
+
+            if isinstance(conclusions_arg, list):
+                # The user has supplied a list of morphisms, none of
+                # which have any attributes.
+                empty = EmptySet
+
+                for morphism in conclusions_arg:
+                    # Check that no new objects appear in conclusions.
+                    if ((sympify(objects.contains(morphism.domain)) is S.true) and
+                        (sympify(objects.contains(morphism.codomain)) is S.true)):
+                        # No need to add identities and recurse
+                        # composites this time.
+                        Diagram._add_morphism_closure(
+                            conclusions, morphism, empty, add_identities=False,
+                            recurse_composites=False)
+            elif isinstance(conclusions_arg, (dict, Dict)):
+                # The user has supplied a dictionary of morphisms and
+                # their properties.
+                for morphism, props in conclusions_arg.items():
+                    # Check that no new objects appear in conclusions.
+                    if (morphism.domain in objects) and \
+                       (morphism.codomain in objects):
+                        # No need to add identities and recurse
+                        # composites this time.
+                        Diagram._add_morphism_closure(
+                            conclusions, morphism, FiniteSet(*props) if iterable(props) else FiniteSet(props),
+                            add_identities=False, recurse_composites=False)
+
+        return Basic.__new__(cls, Dict(premises), Dict(conclusions), objects)
+
+    @property
+    def premises(self):
+        """
+        Returns the premises of this diagram.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import IdentityMorphism, Diagram
+        >>> from sympy import pretty
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> id_A = IdentityMorphism(A)
+        >>> id_B = IdentityMorphism(B)
+        >>> d = Diagram([f])
+        >>> print(pretty(d.premises, use_unicode=False))
+        {id:A-->A: EmptySet, id:B-->B: EmptySet, f:A-->B: EmptySet}
+
+        """
+        return self.args[0]
+
+    @property
+    def conclusions(self):
+        """
+        Returns the conclusions of this diagram.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import IdentityMorphism, Diagram
+        >>> from sympy import FiniteSet
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g])
+        >>> IdentityMorphism(A) in d.premises.keys()
+        True
+        >>> g * f in d.premises.keys()
+        True
+        >>> d = Diagram([f, g], {g * f: "unique"})
+        >>> d.conclusions[g * f] == FiniteSet("unique")
+        True
+
+        """
+        return self.args[1]
+
+    @property
+    def objects(self):
+        """
+        Returns the :class:`~.FiniteSet` of objects that appear in this
+        diagram.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism, Diagram
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g])
+        >>> d.objects
+        {Object("A"), Object("B"), Object("C")}
+
+        """
+        return self.args[2]
+
+    def hom(self, A, B):
+        """
+        Returns a 2-tuple of sets of morphisms between objects ``A`` and
+        ``B``: one set of morphisms listed as premises, and the other set
+        of morphisms listed as conclusions.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism, Diagram
+        >>> from sympy import pretty
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g], {g * f: "unique"})
+        >>> print(pretty(d.hom(A, C), use_unicode=False))
+        ({g*f:A-->C}, {g*f:A-->C})
+
+        See Also
+        ========
+        Object, Morphism
+        """
+        premises = EmptySet
+        conclusions = EmptySet
+
+        for morphism in self.premises.keys():
+            if (morphism.domain == A) and (morphism.codomain == B):
+                premises |= FiniteSet(morphism)
+        for morphism in self.conclusions.keys():
+            if (morphism.domain == A) and (morphism.codomain == B):
+                conclusions |= FiniteSet(morphism)
+
+        return (premises, conclusions)
+
+    def is_subdiagram(self, diagram):
+        """
+        Checks whether ``diagram`` is a subdiagram of ``self``.
+        Diagram `D'` is a subdiagram of `D` if all premises
+        (conclusions) of `D'` are contained in the premises
+        (conclusions) of `D`.  The morphisms contained
+        both in `D'` and `D` should have the same properties for `D'`
+        to be a subdiagram of `D`.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism, Diagram
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g], {g * f: "unique"})
+        >>> d1 = Diagram([f])
+        >>> d.is_subdiagram(d1)
+        True
+        >>> d1.is_subdiagram(d)
+        False
+        """
+        premises = all((m in self.premises) and
+                       (diagram.premises[m] == self.premises[m])
+                       for m in diagram.premises)
+        if not premises:
+            return False
+
+        conclusions = all((m in self.conclusions) and
+                          (diagram.conclusions[m] == self.conclusions[m])
+                          for m in diagram.conclusions)
+
+        # Premises is surely ``True`` here.
+        return conclusions
+
+    def subdiagram_from_objects(self, objects):
+        """
+        If ``objects`` is a subset of the objects of ``self``, returns
+        a diagram which has as premises all those premises of ``self``
+        which have a domains and codomains in ``objects``, likewise
+        for conclusions.  Properties are preserved.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism, Diagram
+        >>> from sympy import FiniteSet
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> d = Diagram([f, g], {f: "unique", g*f: "veryunique"})
+        >>> d1 = d.subdiagram_from_objects(FiniteSet(A, B))
+        >>> d1 == Diagram([f], {f: "unique"})
+        True
+        """
+        if not objects.is_subset(self.objects):
+            raise ValueError(
+                "Supplied objects should all belong to the diagram.")
+
+        new_premises = {}
+        for morphism, props in self.premises.items():
+            if ((sympify(objects.contains(morphism.domain)) is S.true) and
+                (sympify(objects.contains(morphism.codomain)) is S.true)):
+                new_premises[morphism] = props
+
+        new_conclusions = {}
+        for morphism, props in self.conclusions.items():
+            if ((sympify(objects.contains(morphism.domain)) is S.true) and
+                (sympify(objects.contains(morphism.codomain)) is S.true)):
+                new_conclusions[morphism] = props
+
+        return Diagram(new_premises, new_conclusions)

lib/python3.12/site-packages/sympy/categories/diagram_drawing.py

@@ -0,0 +1,2582 @@
+r"""
+This module contains the functionality to arrange the nodes of a
+diagram on an abstract grid, and then to produce a graphical
+representation of the grid.
+
+The currently supported back-ends are Xy-pic [Xypic].
+
+Layout Algorithm
+================
+
+This section provides an overview of the algorithms implemented in
+:class:`DiagramGrid` to lay out diagrams.
+
+The first step of the algorithm is the removal composite and identity
+morphisms which do not have properties in the supplied diagram.  The
+premises and conclusions of the diagram are then merged.
+
+The generic layout algorithm begins with the construction of the
+"skeleton" of the diagram.  The skeleton is an undirected graph which
+has the objects of the diagram as vertices and has an (undirected)
+edge between each pair of objects between which there exist morphisms.
+The direction of the morphisms does not matter at this stage.  The
+skeleton also includes an edge between each pair of vertices `A` and
+`C` such that there exists an object `B` which is connected via
+a morphism to `A`, and via a morphism to `C`.
+
+The skeleton constructed in this way has the property that every
+object is a vertex of a triangle formed by three edges of the
+skeleton.  This property lies at the base of the generic layout
+algorithm.
+
+After the skeleton has been constructed, the algorithm lists all
+triangles which can be formed.  Note that some triangles will not have
+all edges corresponding to morphisms which will actually be drawn.
+Triangles which have only one edge or less which will actually be
+drawn are immediately discarded.
+
+The list of triangles is sorted according to the number of edges which
+correspond to morphisms, then the triangle with the least number of such
+edges is selected.  One of such edges is picked and the corresponding
+objects are placed horizontally, on a grid.  This edge is recorded to
+be in the fringe.  The algorithm then finds a "welding" of a triangle
+to the fringe.  A welding is an edge in the fringe where a triangle
+could be attached.  If the algorithm succeeds in finding such a
+welding, it adds to the grid that vertex of the triangle which was not
+yet included in any edge in the fringe and records the two new edges in
+the fringe.  This process continues iteratively until all objects of
+the diagram has been placed or until no more weldings can be found.
+
+An edge is only removed from the fringe when a welding to this edge
+has been found, and there is no room around this edge to place
+another vertex.
+
+When no more weldings can be found, but there are still triangles
+left, the algorithm searches for a possibility of attaching one of the
+remaining triangles to the existing structure by a vertex.  If such a
+possibility is found, the corresponding edge of the found triangle is
+placed in the found space and the iterative process of welding
+triangles restarts.
+
+When logical groups are supplied, each of these groups is laid out
+independently.  Then a diagram is constructed in which groups are
+objects and any two logical groups between which there exist morphisms
+are connected via a morphism.  This diagram is laid out.  Finally,
+the grid which includes all objects of the initial diagram is
+constructed by replacing the cells which contain logical groups with
+the corresponding laid out grids, and by correspondingly expanding the
+rows and columns.
+
+The sequential layout algorithm begins by constructing the
+underlying undirected graph defined by the morphisms obtained after
+simplifying premises and conclusions and merging them (see above).
+The vertex with the minimal degree is then picked up and depth-first
+search is started from it.  All objects which are located at distance
+`n` from the root in the depth-first search tree, are positioned in
+the `n`-th column of the resulting grid.  The sequential layout will
+therefore attempt to lay the objects out along a line.
+
+References
+==========
+
+.. [Xypic] https://xy-pic.sourceforge.net/
+
+"""
+from sympy.categories import (CompositeMorphism, IdentityMorphism,
+                              NamedMorphism, Diagram)
+from sympy.core import Dict, Symbol, default_sort_key
+from sympy.printing.latex import latex
+from sympy.sets import FiniteSet
+from sympy.utilities.iterables import iterable
+from sympy.utilities.decorator import doctest_depends_on
+
+from itertools import chain
+
+
+__doctest_requires__ = {('preview_diagram',): 'pyglet'}
+
+
+class _GrowableGrid:
+    """
+    Holds a growable grid of objects.
+
+    Explanation
+    ===========
+
+    It is possible to append or prepend a row or a column to the grid
+    using the corresponding methods.  Prepending rows or columns has
+    the effect of changing the coordinates of the already existing
+    elements.
+
+    This class currently represents a naive implementation of the
+    functionality with little attempt at optimisation.
+    """
+    def __init__(self, width, height):
+        self._width = width
+        self._height = height
+
+        self._array = [[None for j in range(width)] for i in range(height)]
+
+    @property
+    def width(self):
+        return self._width
+
+    @property
+    def height(self):
+        return self._height
+
+    def __getitem__(self, i_j):
+        """
+        Returns the element located at in the i-th line and j-th
+        column.
+        """
+        i, j = i_j
+        return self._array[i][j]
+
+    def __setitem__(self, i_j, newvalue):
+        """
+        Sets the element located at in the i-th line and j-th
+        column.
+        """
+        i, j = i_j
+        self._array[i][j] = newvalue
+
+    def append_row(self):
+        """
+        Appends an empty row to the grid.
+        """
+        self._height += 1
+        self._array.append([None for j in range(self._width)])
+
+    def append_column(self):
+        """
+        Appends an empty column to the grid.
+        """
+        self._width += 1
+        for i in range(self._height):
+            self._array[i].append(None)
+
+    def prepend_row(self):
+        """
+        Prepends the grid with an empty row.
+        """
+        self._height += 1
+        self._array.insert(0, [None for j in range(self._width)])
+
+    def prepend_column(self):
+        """
+        Prepends the grid with an empty column.
+        """
+        self._width += 1
+        for i in range(self._height):
+            self._array[i].insert(0, None)
+
+
+class DiagramGrid:
+    r"""
+    Constructs and holds the fitting of the diagram into a grid.
+
+    Explanation
+    ===========
+
+    The mission of this class is to analyse the structure of the
+    supplied diagram and to place its objects on a grid such that,
+    when the objects and the morphisms are actually drawn, the diagram
+    would be "readable", in the sense that there will not be many
+    intersections of moprhisms.  This class does not perform any
+    actual drawing.  It does strive nevertheless to offer sufficient
+    metadata to draw a diagram.
+
+    Consider the following simple diagram.
+
+    >>> from sympy.categories import Object, NamedMorphism
+    >>> from sympy.categories import Diagram, DiagramGrid
+    >>> from sympy import pprint
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> diagram = Diagram([f, g])
+
+    The simplest way to have a diagram laid out is the following:
+
+    >>> grid = DiagramGrid(diagram)
+    >>> (grid.width, grid.height)
+    (2, 2)
+    >>> pprint(grid)
+    A  B
+    <BLANKLINE>
+       C
+
+    Sometimes one sees the diagram as consisting of logical groups.
+    One can advise ``DiagramGrid`` as to such groups by employing the
+    ``groups`` keyword argument.
+
+    Consider the following diagram:
+
+    >>> D = Object("D")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> h = NamedMorphism(D, A, "h")
+    >>> k = NamedMorphism(D, B, "k")
+    >>> diagram = Diagram([f, g, h, k])
+
+    Lay it out with generic layout:
+
+    >>> grid = DiagramGrid(diagram)
+    >>> pprint(grid)
+    A  B  D
+    <BLANKLINE>
+       C
+
+    Now, we can group the objects `A` and `D` to have them near one
+    another:
+
+    >>> grid = DiagramGrid(diagram, groups=[[A, D], B, C])
+    >>> pprint(grid)
+    B     C
+    <BLANKLINE>
+    A  D
+
+    Note how the positioning of the other objects changes.
+
+    Further indications can be supplied to the constructor of
+    :class:`DiagramGrid` using keyword arguments.  The currently
+    supported hints are explained in the following paragraphs.
+
+    :class:`DiagramGrid` does not automatically guess which layout
+    would suit the supplied diagram better.  Consider, for example,
+    the following linear diagram:
+
+    >>> E = Object("E")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> h = NamedMorphism(C, D, "h")
+    >>> i = NamedMorphism(D, E, "i")
+    >>> diagram = Diagram([f, g, h, i])
+
+    When laid out with the generic layout, it does not get to look
+    linear:
+
+    >>> grid = DiagramGrid(diagram)
+    >>> pprint(grid)
+    A  B
+    <BLANKLINE>
+       C  D
+    <BLANKLINE>
+          E
+
+    To get it laid out in a line, use ``layout="sequential"``:
+
+    >>> grid = DiagramGrid(diagram, layout="sequential")
+    >>> pprint(grid)
+    A  B  C  D  E
+
+    One may sometimes need to transpose the resulting layout.  While
+    this can always be done by hand, :class:`DiagramGrid` provides a
+    hint for that purpose:
+
+    >>> grid = DiagramGrid(diagram, layout="sequential", transpose=True)
+    >>> pprint(grid)
+    A
+    <BLANKLINE>
+    B
+    <BLANKLINE>
+    C
+    <BLANKLINE>
+    D
+    <BLANKLINE>
+    E
+
+    Separate hints can also be provided for each group.  For an
+    example, refer to ``tests/test_drawing.py``, and see the different
+    ways in which the five lemma [FiveLemma] can be laid out.
+
+    See Also
+    ========
+
+    Diagram
+
+    References
+    ==========
+
+    .. [FiveLemma] https://en.wikipedia.org/wiki/Five_lemma
+    """
+    @staticmethod
+    def _simplify_morphisms(morphisms):
+        """
+        Given a dictionary mapping morphisms to their properties,
+        returns a new dictionary in which there are no morphisms which
+        do not have properties, and which are compositions of other
+        morphisms included in the dictionary.  Identities are dropped
+        as well.
+        """
+        newmorphisms = {}
+        for morphism, props in morphisms.items():
+            if isinstance(morphism, CompositeMorphism) and not props:
+                continue
+            elif isinstance(morphism, IdentityMorphism):
+                continue
+            else:
+                newmorphisms[morphism] = props
+        return newmorphisms
+
+    @staticmethod
+    def _merge_premises_conclusions(premises, conclusions):
+        """
+        Given two dictionaries of morphisms and their properties,
+        produces a single dictionary which includes elements from both
+        dictionaries.  If a morphism has some properties in premises
+        and also in conclusions, the properties in conclusions take
+        priority.
+        """
+        return dict(chain(premises.items(), conclusions.items()))
+
+    @staticmethod
+    def _juxtapose_edges(edge1, edge2):
+        """
+        If ``edge1`` and ``edge2`` have precisely one common endpoint,
+        returns an edge which would form a triangle with ``edge1`` and
+        ``edge2``.
+
+        If ``edge1`` and ``edge2`` do not have a common endpoint,
+        returns ``None``.
+
+        If ``edge1`` and ``edge`` are the same edge, returns ``None``.
+        """
+        intersection = edge1 & edge2
+        if len(intersection) != 1:
+            # The edges either have no common points or are equal.
+            return None
+
+        # The edges have a common endpoint.  Extract the different
+        # endpoints and set up the new edge.
+        return (edge1 - intersection) | (edge2 - intersection)
+
+    @staticmethod
+    def _add_edge_append(dictionary, edge, elem):
+        """
+        If ``edge`` is not in ``dictionary``, adds ``edge`` to the
+        dictionary and sets its value to ``[elem]``.  Otherwise
+        appends ``elem`` to the value of existing entry.
+
+        Note that edges are undirected, thus `(A, B) = (B, A)`.
+        """
+        if edge in dictionary:
+            dictionary[edge].append(elem)
+        else:
+            dictionary[edge] = [elem]
+
+    @staticmethod
+    def _build_skeleton(morphisms):
+        """
+        Creates a dictionary which maps edges to corresponding
+        morphisms.  Thus for a morphism `f:A\rightarrow B`, the edge
+        `(A, B)` will be associated with `f`.  This function also adds
+        to the list those edges which are formed by juxtaposition of
+        two edges already in the list.  These new edges are not
+        associated with any morphism and are only added to assure that
+        the diagram can be decomposed into triangles.
+        """
+        edges = {}
+        # Create edges for morphisms.
+        for morphism in morphisms:
+            DiagramGrid._add_edge_append(
+                edges, frozenset([morphism.domain, morphism.codomain]), morphism)
+
+        # Create new edges by juxtaposing existing edges.
+        edges1 = dict(edges)
+        for w in edges1:
+            for v in edges1:
+                wv = DiagramGrid._juxtapose_edges(w, v)
+                if wv and wv not in edges:
+                    edges[wv] = []
+
+        return edges
+
+    @staticmethod
+    def _list_triangles(edges):
+        """
+        Builds the set of triangles formed by the supplied edges.  The
+        triangles are arbitrary and need not be commutative.  A
+        triangle is a set that contains all three of its sides.
+        """
+        triangles = set()
+
+        for w in edges:
+            for v in edges:
+                wv = DiagramGrid._juxtapose_edges(w, v)
+                if wv and wv in edges:
+                    triangles.add(frozenset([w, v, wv]))
+
+        return triangles
+
+    @staticmethod
+    def _drop_redundant_triangles(triangles, skeleton):
+        """
+        Returns a list which contains only those triangles who have
+        morphisms associated with at least two edges.
+        """
+        return [tri for tri in triangles
+                if len([e for e in tri if skeleton[e]]) >= 2]
+
+    @staticmethod
+    def _morphism_length(morphism):
+        """
+        Returns the length of a morphism.  The length of a morphism is
+        the number of components it consists of.  A non-composite
+        morphism is of length 1.
+        """
+        if isinstance(morphism, CompositeMorphism):
+            return len(morphism.components)
+        else:
+            return 1
+
+    @staticmethod
+    def _compute_triangle_min_sizes(triangles, edges):
+        r"""
+        Returns a dictionary mapping triangles to their minimal sizes.
+        The minimal size of a triangle is the sum of maximal lengths
+        of morphisms associated to the sides of the triangle.  The
+        length of a morphism is the number of components it consists
+        of.  A non-composite morphism is of length 1.
+
+        Sorting triangles by this metric attempts to address two
+        aspects of layout.  For triangles with only simple morphisms
+        in the edge, this assures that triangles with all three edges
+        visible will get typeset after triangles with less visible
+        edges, which sometimes minimizes the necessity in diagonal
+        arrows.  For triangles with composite morphisms in the edges,
+        this assures that objects connected with shorter morphisms
+        will be laid out first, resulting the visual proximity of
+        those objects which are connected by shorter morphisms.
+        """
+        triangle_sizes = {}
+        for triangle in triangles:
+            size = 0
+            for e in triangle:
+                morphisms = edges[e]
+                if morphisms:
+                    size += max(DiagramGrid._morphism_length(m)
+                                for m in morphisms)
+            triangle_sizes[triangle] = size
+        return triangle_sizes
+
+    @staticmethod
+    def _triangle_objects(triangle):
+        """
+        Given a triangle, returns the objects included in it.
+        """
+        # A triangle is a frozenset of three two-element frozensets
+        # (the edges).  This chains the three edges together and
+        # creates a frozenset from the iterator, thus producing a
+        # frozenset of objects of the triangle.
+        return frozenset(chain(*tuple(triangle)))
+
+    @staticmethod
+    def _other_vertex(triangle, edge):
+        """
+        Given a triangle and an edge of it, returns the vertex which
+        opposes the edge.
+        """
+        # This gets the set of objects of the triangle and then
+        # subtracts the set of objects employed in ``edge`` to get the
+        # vertex opposite to ``edge``.
+        return list(DiagramGrid._triangle_objects(triangle) - set(edge))[0]
+
+    @staticmethod
+    def _empty_point(pt, grid):
+        """
+        Checks if the cell at coordinates ``pt`` is either empty or
+        out of the bounds of the grid.
+        """
+        if (pt[0] < 0) or (pt[1] < 0) or \
+           (pt[0] >= grid.height) or (pt[1] >= grid.width):
+            return True
+        return grid[pt] is None
+
+    @staticmethod
+    def _put_object(coords, obj, grid, fringe):
+        """
+        Places an object at the coordinate ``cords`` in ``grid``,
+        growing the grid and updating ``fringe``, if necessary.
+        Returns (0, 0) if no row or column has been prepended, (1, 0)
+        if a row was prepended, (0, 1) if a column was prepended and
+        (1, 1) if both a column and a row were prepended.
+        """
+        (i, j) = coords
+        offset = (0, 0)
+        if i == -1:
+            grid.prepend_row()
+            i = 0
+            offset = (1, 0)
+            for k in range(len(fringe)):
+                ((i1, j1), (i2, j2)) = fringe[k]
+                fringe[k] = ((i1 + 1, j1), (i2 + 1, j2))
+        elif i == grid.height:
+            grid.append_row()
+
+        if j == -1:
+            j = 0
+            offset = (offset[0], 1)
+            grid.prepend_column()
+            for k in range(len(fringe)):
+                ((i1, j1), (i2, j2)) = fringe[k]
+                fringe[k] = ((i1, j1 + 1), (i2, j2 + 1))
+        elif j == grid.width:
+            grid.append_column()
+
+        grid[i, j] = obj
+        return offset
+
+    @staticmethod
+    def _choose_target_cell(pt1, pt2, edge, obj, skeleton, grid):
+        """
+        Given two points, ``pt1`` and ``pt2``, and the welding edge
+        ``edge``, chooses one of the two points to place the opposing
+        vertex ``obj`` of the triangle.  If neither of this points
+        fits, returns ``None``.
+        """
+        pt1_empty = DiagramGrid._empty_point(pt1, grid)
+        pt2_empty = DiagramGrid._empty_point(pt2, grid)
+
+        if pt1_empty and pt2_empty:
+            # Both cells are empty.  Of these two, choose that cell
+            # which will assure that a visible edge of the triangle
+            # will be drawn perpendicularly to the current welding
+            # edge.
+
+            A = grid[edge[0]]
+
+            if skeleton.get(frozenset([A, obj])):
+                return pt1
+            else:
+                return pt2
+        if pt1_empty:
+            return pt1
+        elif pt2_empty:
+            return pt2
+        else:
+            return None
+
+    @staticmethod
+    def _find_triangle_to_weld(triangles, fringe, grid):
+        """
+        Finds, if possible, a triangle and an edge in the ``fringe`` to
+        which the triangle could be attached.  Returns the tuple
+        containing the triangle and the index of the corresponding
+        edge in the ``fringe``.
+
+        This function relies on the fact that objects are unique in
+        the diagram.
+        """
+        for triangle in triangles:
+            for (a, b) in fringe:
+                if frozenset([grid[a], grid[b]]) in triangle:
+                    return (triangle, (a, b))
+        return None
+
+    @staticmethod
+    def _weld_triangle(tri, welding_edge, fringe, grid, skeleton):
+        """
+        If possible, welds the triangle ``tri`` to ``fringe`` and
+        returns ``False``.  If this method encounters a degenerate
+        situation in the fringe and corrects it such that a restart of
+        the search is required, it returns ``True`` (which means that
+        a restart in finding triangle weldings is required).
+
+        A degenerate situation is a situation when an edge listed in
+        the fringe does not belong to the visual boundary of the
+        diagram.
+        """
+        a, b = welding_edge
+        target_cell = None
+
+        obj = DiagramGrid._other_vertex(tri, (grid[a], grid[b]))
+
+        # We now have a triangle and an edge where it can be welded to
+        # the fringe.  Decide where to place the other vertex of the
+        # triangle and check for degenerate situations en route.
+
+        if (abs(a[0] - b[0]) == 1) and (abs(a[1] - b[1]) == 1):
+            # A diagonal edge.
+            target_cell = (a[0], b[1])
+            if grid[target_cell]:
+                # That cell is already occupied.
+                target_cell = (b[0], a[1])
+
+                if grid[target_cell]:
+                    # Degenerate situation, this edge is not
+                    # on the actual fringe.  Correct the
+                    # fringe and go on.
+                    fringe.remove((a, b))
+                    return True
+        elif a[0] == b[0]:
+            # A horizontal edge.  We first attempt to build the
+            # triangle in the downward direction.
+
+            down_left = a[0] + 1, a[1]
+            down_right = a[0] + 1, b[1]
+
+            target_cell = DiagramGrid._choose_target_cell(
+                down_left, down_right, (a, b), obj, skeleton, grid)
+
+            if not target_cell:
+                # No room below this edge.  Check above.
+                up_left = a[0] - 1, a[1]
+                up_right = a[0] - 1, b[1]
+
+                target_cell = DiagramGrid._choose_target_cell(
+                    up_left, up_right, (a, b), obj, skeleton, grid)
+
+                if not target_cell:
+                    # This edge is not in the fringe, remove it
+                    # and restart.
+                    fringe.remove((a, b))
+                    return True
+        elif a[1] == b[1]:
+            # A vertical edge.  We will attempt to place the other
+            # vertex of the triangle to the right of this edge.
+            right_up = a[0], a[1] + 1
+            right_down = b[0], a[1] + 1
+
+            target_cell = DiagramGrid._choose_target_cell(
+                right_up, right_down, (a, b), obj, skeleton, grid)
+
+            if not target_cell:
+                # No room to the left.  See what's to the right.
+                left_up = a[0], a[1] - 1
+                left_down = b[0], a[1] - 1
+
+                target_cell = DiagramGrid._choose_target_cell(
+                    left_up, left_down, (a, b), obj, skeleton, grid)
+
+                if not target_cell:
+                    # This edge is not in the fringe, remove it
+                    # and restart.
+                    fringe.remove((a, b))
+                    return True
+
+        # We now know where to place the other vertex of the
+        # triangle.
+        offset = DiagramGrid._put_object(target_cell, obj, grid, fringe)
+
+        # Take care of the displacement of coordinates if a row or
+        # a column was prepended.
+        target_cell = (target_cell[0] + offset[0],
+                       target_cell[1] + offset[1])
+        a = (a[0] + offset[0], a[1] + offset[1])
+        b = (b[0] + offset[0], b[1] + offset[1])
+
+        fringe.extend([(a, target_cell), (b, target_cell)])
+
+        # No restart is required.
+        return False
+
+    @staticmethod
+    def _triangle_key(tri, triangle_sizes):
+        """
+        Returns a key for the supplied triangle.  It should be the
+        same independently of the hash randomisation.
+        """
+        objects = sorted(
+            DiagramGrid._triangle_objects(tri), key=default_sort_key)
+        return (triangle_sizes[tri], default_sort_key(objects))
+
+    @staticmethod
+    def _pick_root_edge(tri, skeleton):
+        """
+        For a given triangle always picks the same root edge.  The
+        root edge is the edge that will be placed first on the grid.
+        """
+        candidates = [sorted(e, key=default_sort_key)
+                      for e in tri if skeleton[e]]
+        sorted_candidates = sorted(candidates, key=default_sort_key)
+        # Don't forget to assure the proper ordering of the vertices
+        # in this edge.
+        return tuple(sorted(sorted_candidates[0], key=default_sort_key))
+
+    @staticmethod
+    def _drop_irrelevant_triangles(triangles, placed_objects):
+        """
+        Returns only those triangles whose set of objects is not
+        completely included in ``placed_objects``.
+        """
+        return [tri for tri in triangles if not placed_objects.issuperset(
+            DiagramGrid._triangle_objects(tri))]
+
+    @staticmethod
+    def _grow_pseudopod(triangles, fringe, grid, skeleton, placed_objects):
+        """
+        Starting from an object in the existing structure on the ``grid``,
+        adds an edge to which a triangle from ``triangles`` could be
+        welded.  If this method has found a way to do so, it returns
+        the object it has just added.
+
+        This method should be applied when ``_weld_triangle`` cannot
+        find weldings any more.
+        """
+        for i in range(grid.height):
+            for j in range(grid.width):
+                obj = grid[i, j]
+                if not obj:
+                    continue
+
+                # Here we need to choose a triangle which has only
+                # ``obj`` in common with the existing structure.  The
+                # situations when this is not possible should be
+                # handled elsewhere.
+
+                def good_triangle(tri):
+                    objs = DiagramGrid._triangle_objects(tri)
+                    return obj in objs and \
+                        placed_objects & (objs - {obj}) == set()
+
+                tris = [tri for tri in triangles if good_triangle(tri)]
+                if not tris:
+                    # This object is not interesting.
+                    continue
+
+                # Pick the "simplest" of the triangles which could be
+                # attached.  Remember that the list of triangles is
+                # sorted according to their "simplicity" (see
+                # _compute_triangle_min_sizes for the metric).
+                #
+                # Note that ``tris`` are sequentially built from
+                # ``triangles``, so we don't have to worry about hash
+                # randomisation.
+                tri = tris[0]
+
+                # We have found a triangle which could be attached to
+                # the existing structure by a vertex.
+
+                candidates = sorted([e for e in tri if skeleton[e]],
+                                    key=lambda e: FiniteSet(*e).sort_key())
+                edges = [e for e in candidates if obj in e]
+
+                # Note that a meaningful edge (i.e., and edge that is
+                # associated with a morphism) containing ``obj``
+                # always exists.  That's because all triangles are
+                # guaranteed to have at least two meaningful edges.
+                # See _drop_redundant_triangles.
+
+                # Get the object at the other end of the edge.
+                edge = edges[0]
+                other_obj = tuple(edge - frozenset([obj]))[0]
+
+                # Now check for free directions.  When checking for
+                # free directions, prefer the horizontal and vertical
+                # directions.
+                neighbours = [(i - 1, j), (i, j + 1), (i + 1, j), (i, j - 1),
+                              (i - 1, j - 1), (i - 1, j + 1), (i + 1, j - 1), (i + 1, j + 1)]
+
+                for pt in neighbours:
+                    if DiagramGrid._empty_point(pt, grid):
+                        # We have a found a place to grow the
+                        # pseudopod into.
+                        offset = DiagramGrid._put_object(
+                            pt, other_obj, grid, fringe)
+
+                        i += offset[0]
+                        j += offset[1]
+                        pt = (pt[0] + offset[0], pt[1] + offset[1])
+                        fringe.append(((i, j), pt))
+
+                        return other_obj
+
+        # This diagram is actually cooler that I can handle.  Fail cowardly.
+        return None
+
+    @staticmethod
+    def _handle_groups(diagram, groups, merged_morphisms, hints):
+        """
+        Given the slightly preprocessed morphisms of the diagram,
+        produces a grid laid out according to ``groups``.
+
+        If a group has hints, it is laid out with those hints only,
+        without any influence from ``hints``.  Otherwise, it is laid
+        out with ``hints``.
+        """
+        def lay_out_group(group, local_hints):
+            """
+            If ``group`` is a set of objects, uses a ``DiagramGrid``
+            to lay it out and returns the grid.  Otherwise returns the
+            object (i.e., ``group``).  If ``local_hints`` is not
+            empty, it is supplied to ``DiagramGrid`` as the dictionary
+            of hints.  Otherwise, the ``hints`` argument of
+            ``_handle_groups`` is used.
+            """
+            if isinstance(group, FiniteSet):
+                # Set up the corresponding object-to-group
+                # mappings.
+                for obj in group:
+                    obj_groups[obj] = group
+
+                # Lay out the current group.
+                if local_hints:
+                    groups_grids[group] = DiagramGrid(
+                        diagram.subdiagram_from_objects(group), **local_hints)
+                else:
+                    groups_grids[group] = DiagramGrid(
+                        diagram.subdiagram_from_objects(group), **hints)
+            else:
+                obj_groups[group] = group
+
+        def group_to_finiteset(group):
+            """
+            Converts ``group`` to a :class:``FiniteSet`` if it is an
+            iterable.
+            """
+            if iterable(group):
+                return FiniteSet(*group)
+            else:
+                return group
+
+        obj_groups = {}
+        groups_grids = {}
+
+        # We would like to support various containers to represent
+        # groups.  To achieve that, before laying each group out, it
+        # should be converted to a FiniteSet, because that is what the
+        # following code expects.
+
+        if isinstance(groups, (dict, Dict)):
+            finiteset_groups = {}
+            for group, local_hints in groups.items():
+                finiteset_group = group_to_finiteset(group)
+                finiteset_groups[finiteset_group] = local_hints
+                lay_out_group(group, local_hints)
+            groups = finiteset_groups
+        else:
+            finiteset_groups = []
+            for group in groups:
+                finiteset_group = group_to_finiteset(group)
+                finiteset_groups.append(finiteset_group)
+                lay_out_group(finiteset_group, None)
+            groups = finiteset_groups
+
+        new_morphisms = []
+        for morphism in merged_morphisms:
+            dom = obj_groups[morphism.domain]
+            cod = obj_groups[morphism.codomain]
+            # Note that we are not really interested in morphisms
+            # which do not employ two different groups, because
+            # these do not influence the layout.
+            if dom != cod:
+                # These are essentially unnamed morphisms; they are
+                # not going to mess in the final layout.  By giving
+                # them the same names, we avoid unnecessary
+                # duplicates.
+                new_morphisms.append(NamedMorphism(dom, cod, "dummy"))
+
+        # Lay out the new diagram.  Since these are dummy morphisms,
+        # properties and conclusions are irrelevant.
+        top_grid = DiagramGrid(Diagram(new_morphisms))
+
+        # We now have to substitute the groups with the corresponding
+        # grids, laid out at the beginning of this function.  Compute
+        # the size of each row and column in the grid, so that all
+        # nested grids fit.
+
+        def group_size(group):
+            """
+            For the supplied group (or object, eventually), returns
+            the size of the cell that will hold this group (object).
+            """
+            if group in groups_grids:
+                grid = groups_grids[group]
+                return (grid.height, grid.width)
+            else:
+                return (1, 1)
+
+        row_heights = [max(group_size(top_grid[i, j])[0]
+                           for j in range(top_grid.width))
+                       for i in range(top_grid.height)]
+
+        column_widths = [max(group_size(top_grid[i, j])[1]
+                             for i in range(top_grid.height))
+                         for j in range(top_grid.width)]
+
+        grid = _GrowableGrid(sum(column_widths), sum(row_heights))
+
+        real_row = 0
+        real_column = 0
+        for logical_row in range(top_grid.height):
+            for logical_column in range(top_grid.width):
+                obj = top_grid[logical_row, logical_column]
+
+                if obj in groups_grids:
+                    # This is a group.  Copy the corresponding grid in
+                    # place.
+                    local_grid = groups_grids[obj]
+                    for i in range(local_grid.height):
+                        for j in range(local_grid.width):
+                            grid[real_row + i,
+                                real_column + j] = local_grid[i, j]
+                else:
+                    # This is an object.  Just put it there.
+                    grid[real_row, real_column] = obj
+
+                real_column += column_widths[logical_column]
+            real_column = 0
+            real_row += row_heights[logical_row]
+
+        return grid
+
+    @staticmethod
+    def _generic_layout(diagram, merged_morphisms):
+        """
+        Produces the generic layout for the supplied diagram.
+        """
+        all_objects = set(diagram.objects)
+        if len(all_objects) == 1:
+            # There only one object in the diagram, just put in on 1x1
+            # grid.
+            grid = _GrowableGrid(1, 1)
+            grid[0, 0] = tuple(all_objects)[0]
+            return grid
+
+        skeleton = DiagramGrid._build_skeleton(merged_morphisms)
+
+        grid = _GrowableGrid(2, 1)
+
+        if len(skeleton) == 1:
+            # This diagram contains only one morphism.  Draw it
+            # horizontally.
+            objects = sorted(all_objects, key=default_sort_key)
+            grid[0, 0] = objects[0]
+            grid[0, 1] = objects[1]
+
+            return grid
+
+        triangles = DiagramGrid._list_triangles(skeleton)
+        triangles = DiagramGrid._drop_redundant_triangles(triangles, skeleton)
+        triangle_sizes = DiagramGrid._compute_triangle_min_sizes(
+            triangles, skeleton)
+
+        triangles = sorted(triangles, key=lambda tri:
+                           DiagramGrid._triangle_key(tri, triangle_sizes))
+
+        # Place the first edge on the grid.
+        root_edge = DiagramGrid._pick_root_edge(triangles[0], skeleton)
+        grid[0, 0], grid[0, 1] = root_edge
+        fringe = [((0, 0), (0, 1))]
+
+        # Record which objects we now have on the grid.
+        placed_objects = set(root_edge)
+
+        while placed_objects != all_objects:
+            welding = DiagramGrid._find_triangle_to_weld(
+                triangles, fringe, grid)
+
+            if welding:
+                (triangle, welding_edge) = welding
+
+                restart_required = DiagramGrid._weld_triangle(
+                    triangle, welding_edge, fringe, grid, skeleton)
+                if restart_required:
+                    continue
+
+                placed_objects.update(
+                    DiagramGrid._triangle_objects(triangle))
+            else:
+                # No more weldings found.  Try to attach triangles by
+                # vertices.
+                new_obj = DiagramGrid._grow_pseudopod(
+                    triangles, fringe, grid, skeleton, placed_objects)
+
+                if not new_obj:
+                    # No more triangles can be attached, not even by
+                    # the edge.  We will set up a new diagram out of
+                    # what has been left, laid it out independently,
+                    # and then attach it to this one.
+
+                    remaining_objects = all_objects - placed_objects
+
+                    remaining_diagram = diagram.subdiagram_from_objects(
+                        FiniteSet(*remaining_objects))
+                    remaining_grid = DiagramGrid(remaining_diagram)
+
+                    # Now, let's glue ``remaining_grid`` to ``grid``.
+                    final_width = grid.width + remaining_grid.width
+                    final_height = max(grid.height, remaining_grid.height)
+                    final_grid = _GrowableGrid(final_width, final_height)
+
+                    for i in range(grid.width):
+                        for j in range(grid.height):
+                            final_grid[i, j] = grid[i, j]
+
+                    start_j = grid.width
+                    for i in range(remaining_grid.height):
+                        for j in range(remaining_grid.width):
+                            final_grid[i, start_j + j] = remaining_grid[i, j]
+
+                    return final_grid
+
+                placed_objects.add(new_obj)
+
+            triangles = DiagramGrid._drop_irrelevant_triangles(
+                triangles, placed_objects)
+
+        return grid
+
+    @staticmethod
+    def _get_undirected_graph(objects, merged_morphisms):
+        """
+        Given the objects and the relevant morphisms of a diagram,
+        returns the adjacency lists of the underlying undirected
+        graph.
+        """
+        adjlists = {}
+        for obj in objects:
+            adjlists[obj] = []
+
+        for morphism in merged_morphisms:
+            adjlists[morphism.domain].append(morphism.codomain)
+            adjlists[morphism.codomain].append(morphism.domain)
+
+        # Assure that the objects in the adjacency list are always in
+        # the same order.
+        for obj in adjlists.keys():
+            adjlists[obj].sort(key=default_sort_key)
+
+        return adjlists
+
+    @staticmethod
+    def _sequential_layout(diagram, merged_morphisms):
+        r"""
+        Lays out the diagram in "sequential" layout.  This method
+        will attempt to produce a result as close to a line as
+        possible.  For linear diagrams, the result will actually be a
+        line.
+        """
+        objects = diagram.objects
+        sorted_objects = sorted(objects, key=default_sort_key)
+
+        # Set up the adjacency lists of the underlying undirected
+        # graph of ``merged_morphisms``.
+        adjlists = DiagramGrid._get_undirected_graph(objects, merged_morphisms)
+
+        root = min(sorted_objects, key=lambda x: len(adjlists[x]))
+        grid = _GrowableGrid(1, 1)
+        grid[0, 0] = root
+
+        placed_objects = {root}
+
+        def place_objects(pt, placed_objects):
+            """
+            Does depth-first search in the underlying graph of the
+            diagram and places the objects en route.
+            """
+            # We will start placing new objects from here.
+            new_pt = (pt[0], pt[1] + 1)
+
+            for adjacent_obj in adjlists[grid[pt]]:
+                if adjacent_obj in placed_objects:
+                    # This object has already been placed.
+                    continue
+
+                DiagramGrid._put_object(new_pt, adjacent_obj, grid, [])
+                placed_objects.add(adjacent_obj)
+                placed_objects.update(place_objects(new_pt, placed_objects))
+
+                new_pt = (new_pt[0] + 1, new_pt[1])
+
+            return placed_objects
+
+        place_objects((0, 0), placed_objects)
+
+        return grid
+
+    @staticmethod
+    def _drop_inessential_morphisms(merged_morphisms):
+        r"""
+        Removes those morphisms which should appear in the diagram,
+        but which have no relevance to object layout.
+
+        Currently this removes "loop" morphisms: the non-identity
+        morphisms with the same domains and codomains.
+        """
+        morphisms = [m for m in merged_morphisms if m.domain != m.codomain]
+        return morphisms
+
+    @staticmethod
+    def _get_connected_components(objects, merged_morphisms):
+        """
+        Given a container of morphisms, returns a list of connected
+        components formed by these morphisms.  A connected component
+        is represented by a diagram consisting of the corresponding
+        morphisms.
+        """
+        component_index = {}
+        for o in objects:
+            component_index[o] = None
+
+        # Get the underlying undirected graph of the diagram.
+        adjlist = DiagramGrid._get_undirected_graph(objects, merged_morphisms)
+
+        def traverse_component(object, current_index):
+            """
+            Does a depth-first search traversal of the component
+            containing ``object``.
+            """
+            component_index[object] = current_index
+            for o in adjlist[object]:
+                if component_index[o] is None:
+                    traverse_component(o, current_index)
+
+        # Traverse all components.
+        current_index = 0
+        for o in adjlist:
+            if component_index[o] is None:
+                traverse_component(o, current_index)
+                current_index += 1
+
+        # List the objects of the components.
+        component_objects = [[] for i in range(current_index)]
+        for o, idx in component_index.items():
+            component_objects[idx].append(o)
+
+        # Finally, list the morphisms belonging to each component.
+        #
+        # Note: If some objects are isolated, they will not get any
+        # morphisms at this stage, and since the layout algorithm
+        # relies, we are essentially going to lose this object.
+        # Therefore, check if there are isolated objects and, for each
+        # of them, provide the trivial identity morphism.  It will get
+        # discarded later, but the object will be there.
+
+        component_morphisms = []
+        for component in component_objects:
+            current_morphisms = {}
+            for m in merged_morphisms:
+                if (m.domain in component) and (m.codomain in component):
+                    current_morphisms[m] = merged_morphisms[m]
+
+            if len(component) == 1:
+                # Let's add an identity morphism, for the sake of
+                # surely having morphisms in this component.
+                current_morphisms[IdentityMorphism(component[0])] = FiniteSet()
+
+            component_morphisms.append(Diagram(current_morphisms))
+
+        return component_morphisms
+
+    def __init__(self, diagram, groups=None, **hints):
+        premises = DiagramGrid._simplify_morphisms(diagram.premises)
+        conclusions = DiagramGrid._simplify_morphisms(diagram.conclusions)
+        all_merged_morphisms = DiagramGrid._merge_premises_conclusions(
+            premises, conclusions)
+        merged_morphisms = DiagramGrid._drop_inessential_morphisms(
+            all_merged_morphisms)
+
+        # Store the merged morphisms for later use.
+        self._morphisms = all_merged_morphisms
+
+        components = DiagramGrid._get_connected_components(
+            diagram.objects, all_merged_morphisms)
+
+        if groups and (groups != diagram.objects):
+            # Lay out the diagram according to the groups.
+            self._grid = DiagramGrid._handle_groups(
+                diagram, groups, merged_morphisms, hints)
+        elif len(components) > 1:
+            # Note that we check for connectedness _before_ checking
+            # the layout hints because the layout strategies don't
+            # know how to deal with disconnected diagrams.
+
+            # The diagram is disconnected.  Lay out the components
+            # independently.
+            grids = []
+
+            # Sort the components to eventually get the grids arranged
+            # in a fixed, hash-independent order.
+            components = sorted(components, key=default_sort_key)
+
+            for component in components:
+                grid = DiagramGrid(component, **hints)
+                grids.append(grid)
+
+            # Throw the grids together, in a line.
+            total_width = sum(g.width for g in grids)
+            total_height = max(g.height for g in grids)
+
+            grid = _GrowableGrid(total_width, total_height)
+            start_j = 0
+            for g in grids:
+                for i in range(g.height):
+                    for j in range(g.width):
+                        grid[i, start_j + j] = g[i, j]
+
+                start_j += g.width
+
+            self._grid = grid
+        elif "layout" in hints:
+            if hints["layout"] == "sequential":
+                self._grid = DiagramGrid._sequential_layout(
+                    diagram, merged_morphisms)
+        else:
+            self._grid = DiagramGrid._generic_layout(diagram, merged_morphisms)
+
+        if hints.get("transpose"):
+            # Transpose the resulting grid.
+            grid = _GrowableGrid(self._grid.height, self._grid.width)
+            for i in range(self._grid.height):
+                for j in range(self._grid.width):
+                    grid[j, i] = self._grid[i, j]
+            self._grid = grid
+
+    @property
+    def width(self):
+        """
+        Returns the number of columns in this diagram layout.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import Diagram, DiagramGrid
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> diagram = Diagram([f, g])
+        >>> grid = DiagramGrid(diagram)
+        >>> grid.width
+        2
+
+        """
+        return self._grid.width
+
+    @property
+    def height(self):
+        """
+        Returns the number of rows in this diagram layout.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import Diagram, DiagramGrid
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> diagram = Diagram([f, g])
+        >>> grid = DiagramGrid(diagram)
+        >>> grid.height
+        2
+
+        """
+        return self._grid.height
+
+    def __getitem__(self, i_j):
+        """
+        Returns the object placed in the row ``i`` and column ``j``.
+        The indices are 0-based.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import Diagram, DiagramGrid
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> diagram = Diagram([f, g])
+        >>> grid = DiagramGrid(diagram)
+        >>> (grid[0, 0], grid[0, 1])
+        (Object("A"), Object("B"))
+        >>> (grid[1, 0], grid[1, 1])
+        (None, Object("C"))
+
+        """
+        i, j = i_j
+        return self._grid[i, j]
+
+    @property
+    def morphisms(self):
+        """
+        Returns those morphisms (and their properties) which are
+        sufficiently meaningful to be drawn.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import Diagram, DiagramGrid
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> diagram = Diagram([f, g])
+        >>> grid = DiagramGrid(diagram)
+        >>> grid.morphisms
+        {NamedMorphism(Object("A"), Object("B"), "f"): EmptySet,
+        NamedMorphism(Object("B"), Object("C"), "g"): EmptySet}
+
+        """
+        return self._morphisms
+
+    def __str__(self):
+        """
+        Produces a string representation of this class.
+
+        This method returns a string representation of the underlying
+        list of lists of objects.
+
+        Examples
+        ========
+
+        >>> from sympy.categories import Object, NamedMorphism
+        >>> from sympy.categories import Diagram, DiagramGrid
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> diagram = Diagram([f, g])
+        >>> grid = DiagramGrid(diagram)
+        >>> print(grid)
+        [[Object("A"), Object("B")],
+        [None, Object("C")]]
+
+        """
+        return repr(self._grid._array)
+
+
+class ArrowStringDescription:
+    r"""
+    Stores the information necessary for producing an Xy-pic
+    description of an arrow.
+
+    The principal goal of this class is to abstract away the string
+    representation of an arrow and to also provide the functionality
+    to produce the actual Xy-pic string.
+
+    ``unit`` sets the unit which will be used to specify the amount of
+    curving and other distances.  ``horizontal_direction`` should be a
+    string of ``"r"`` or ``"l"`` specifying the horizontal offset of the
+    target cell of the arrow relatively to the current one.
+    ``vertical_direction`` should  specify the vertical offset using a
+    series of either ``"d"`` or ``"u"``.  ``label_position`` should be
+    either ``"^"``, ``"_"``,  or ``"|"`` to specify that the label should
+    be positioned above the arrow, below the arrow or just over the arrow,
+    in a break.  Note that the notions "above" and "below" are relative
+    to arrow direction.  ``label`` stores the morphism label.
+
+    This works as follows (disregard the yet unexplained arguments):
+
+    >>> from sympy.categories.diagram_drawing import ArrowStringDescription
+    >>> astr = ArrowStringDescription(
+    ... unit="mm", curving=None, curving_amount=None,
+    ... looping_start=None, looping_end=None, horizontal_direction="d",
+    ... vertical_direction="r", label_position="_", label="f")
+    >>> print(str(astr))
+    \ar[dr]_{f}
+
+    ``curving`` should be one of ``"^"``, ``"_"`` to specify in which
+    direction the arrow is going to curve. ``curving_amount`` is a number
+    describing how many ``unit``'s the morphism is going to curve:
+
+    >>> astr = ArrowStringDescription(
+    ... unit="mm", curving="^", curving_amount=12,
+    ... looping_start=None, looping_end=None, horizontal_direction="d",
+    ... vertical_direction="r", label_position="_", label="f")
+    >>> print(str(astr))
+    \ar@/^12mm/[dr]_{f}
+
+    ``looping_start`` and ``looping_end`` are currently only used for
+    loop morphisms, those which have the same domain and codomain.
+    These two attributes should store a valid Xy-pic direction and
+    specify, correspondingly, the direction the arrow gets out into
+    and the direction the arrow gets back from:
+
+    >>> astr = ArrowStringDescription(
+    ... unit="mm", curving=None, curving_amount=None,
+    ... looping_start="u", looping_end="l", horizontal_direction="",
+    ... vertical_direction="", label_position="_", label="f")
+    >>> print(str(astr))
+    \ar@(u,l)[]_{f}
+
+    ``label_displacement`` controls how far the arrow label is from
+    the ends of the arrow.  For example, to position the arrow label
+    near the arrow head, use ">":
+
+    >>> astr = ArrowStringDescription(
+    ... unit="mm", curving="^", curving_amount=12,
+    ... looping_start=None, looping_end=None, horizontal_direction="d",
+    ... vertical_direction="r", label_position="_", label="f")
+    >>> astr.label_displacement = ">"
+    >>> print(str(astr))
+    \ar@/^12mm/[dr]_>{f}
+
+    Finally, ``arrow_style`` is used to specify the arrow style.  To
+    get a dashed arrow, for example, use "{-->}" as arrow style:
+
+    >>> astr = ArrowStringDescription(
+    ... unit="mm", curving="^", curving_amount=12,
+    ... looping_start=None, looping_end=None, horizontal_direction="d",
+    ... vertical_direction="r", label_position="_", label="f")
+    >>> astr.arrow_style = "{-->}"
+    >>> print(str(astr))
+    \ar@/^12mm/@{-->}[dr]_{f}
+
+    Notes
+    =====
+
+    Instances of :class:`ArrowStringDescription` will be constructed
+    by :class:`XypicDiagramDrawer` and provided for further use in
+    formatters.  The user is not expected to construct instances of
+    :class:`ArrowStringDescription` themselves.
+
+    To be able to properly utilise this class, the reader is encouraged
+    to checkout the Xy-pic user guide, available at [Xypic].
+
+    See Also
+    ========
+
+    XypicDiagramDrawer
+
+    References
+    ==========
+
+    .. [Xypic] https://xy-pic.sourceforge.net/
+    """
+    def __init__(self, unit, curving, curving_amount, looping_start,
+                 looping_end, horizontal_direction, vertical_direction,
+                 label_position, label):
+        self.unit = unit
+        self.curving = curving
+        self.curving_amount = curving_amount
+        self.looping_start = looping_start
+        self.looping_end = looping_end
+        self.horizontal_direction = horizontal_direction
+        self.vertical_direction = vertical_direction
+        self.label_position = label_position
+        self.label = label
+
+        self.label_displacement = ""
+        self.arrow_style = ""
+
+        # This flag shows that the position of the label of this
+        # morphism was set while typesetting a curved morphism and
+        # should not be modified later.
+        self.forced_label_position = False
+
+    def __str__(self):
+        if self.curving:
+            curving_str = "@/%s%d%s/" % (self.curving, self.curving_amount,
+                                         self.unit)
+        else:
+            curving_str = ""
+
+        if self.looping_start and self.looping_end:
+            looping_str = "@(%s,%s)" % (self.looping_start, self.looping_end)
+        else:
+            looping_str = ""
+
+        if self.arrow_style:
+
+            style_str = "@" + self.arrow_style
+        else:
+            style_str = ""
+
+        return "\\ar%s%s%s[%s%s]%s%s{%s}" % \
+               (curving_str, looping_str, style_str, self.horizontal_direction,
+                self.vertical_direction, self.label_position,
+                self.label_displacement, self.label)
+
+
+class XypicDiagramDrawer:
+    r"""
+    Given a :class:`~.Diagram` and the corresponding
+    :class:`DiagramGrid`, produces the Xy-pic representation of the
+    diagram.
+
+    The most important method in this class is ``draw``.  Consider the
+    following triangle diagram:
+
+    >>> from sympy.categories import Object, NamedMorphism, Diagram
+    >>> from sympy.categories import DiagramGrid, XypicDiagramDrawer
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> diagram = Diagram([f, g], {g * f: "unique"})
+
+    To draw this diagram, its objects need to be laid out with a
+    :class:`DiagramGrid`::
+
+    >>> grid = DiagramGrid(diagram)
+
+    Finally, the drawing:
+
+    >>> drawer = XypicDiagramDrawer()
+    >>> print(drawer.draw(diagram, grid))
+    \xymatrix{
+    A \ar[d]_{g\circ f} \ar[r]^{f} & B \ar[ld]^{g} \\
+    C &
+    }
+
+    For further details see the docstring of this method.
+
+    To control the appearance of the arrows, formatters are used.  The
+    dictionary ``arrow_formatters`` maps morphisms to formatter
+    functions.  A formatter is accepts an
+    :class:`ArrowStringDescription` and is allowed to modify any of
+    the arrow properties exposed thereby.  For example, to have all
+    morphisms with the property ``unique`` appear as dashed arrows,
+    and to have their names prepended with `\exists !`, the following
+    should be done:
+
+    >>> def formatter(astr):
+    ...   astr.label = r"\exists !" + astr.label
+    ...   astr.arrow_style = "{-->}"
+    >>> drawer.arrow_formatters["unique"] = formatter
+    >>> print(drawer.draw(diagram, grid))
+    \xymatrix{
+    A \ar@{-->}[d]_{\exists !g\circ f} \ar[r]^{f} & B \ar[ld]^{g} \\
+    C &
+    }
+
+    To modify the appearance of all arrows in the diagram, set
+    ``default_arrow_formatter``.  For example, to place all morphism
+    labels a little bit farther from the arrow head so that they look
+    more centred, do as follows:
+
+    >>> def default_formatter(astr):
+    ...   astr.label_displacement = "(0.45)"
+    >>> drawer.default_arrow_formatter = default_formatter
+    >>> print(drawer.draw(diagram, grid))
+    \xymatrix{
+    A \ar@{-->}[d]_(0.45){\exists !g\circ f} \ar[r]^(0.45){f} & B \ar[ld]^(0.45){g} \\
+    C &
+    }
+
+    In some diagrams some morphisms are drawn as curved arrows.
+    Consider the following diagram:
+
+    >>> D = Object("D")
+    >>> E = Object("E")
+    >>> h = NamedMorphism(D, A, "h")
+    >>> k = NamedMorphism(D, B, "k")
+    >>> diagram = Diagram([f, g, h, k])
+    >>> grid = DiagramGrid(diagram)
+    >>> drawer = XypicDiagramDrawer()
+    >>> print(drawer.draw(diagram, grid))
+    \xymatrix{
+    A \ar[r]_{f} & B \ar[d]^{g} & D \ar[l]^{k} \ar@/_3mm/[ll]_{h} \\
+    & C &
+    }
+
+    To control how far the morphisms are curved by default, one can
+    use the ``unit`` and ``default_curving_amount`` attributes:
+
+    >>> drawer.unit = "cm"
+    >>> drawer.default_curving_amount = 1
+    >>> print(drawer.draw(diagram, grid))
+    \xymatrix{
+    A \ar[r]_{f} & B \ar[d]^{g} & D \ar[l]^{k} \ar@/_1cm/[ll]_{h} \\
+    & C &
+    }
+
+    In some diagrams, there are multiple curved morphisms between the
+    same two objects.  To control by how much the curving changes
+    between two such successive morphisms, use
+    ``default_curving_step``:
+
+    >>> drawer.default_curving_step = 1
+    >>> h1 = NamedMorphism(A, D, "h1")
+    >>> diagram = Diagram([f, g, h, k, h1])
+    >>> grid = DiagramGrid(diagram)
+    >>> print(drawer.draw(diagram, grid))
+    \xymatrix{
+    A \ar[r]_{f} \ar@/^1cm/[rr]^{h_{1}} & B \ar[d]^{g} & D \ar[l]^{k} \ar@/_2cm/[ll]_{h} \\
+    & C &
+    }
+
+    The default value of ``default_curving_step`` is 4 units.
+
+    See Also
+    ========
+
+    draw, ArrowStringDescription
+    """
+    def __init__(self):
+        self.unit = "mm"
+        self.default_curving_amount = 3
+        self.default_curving_step = 4
+
+        # This dictionary maps properties to the corresponding arrow
+        # formatters.
+        self.arrow_formatters = {}
+
+        # This is the default arrow formatter which will be applied to
+        # each arrow independently of its properties.
+        self.default_arrow_formatter = None
+
+    @staticmethod
+    def _process_loop_morphism(i, j, grid, morphisms_str_info, object_coords):
+        """
+        Produces the information required for constructing the string
+        representation of a loop morphism.  This function is invoked
+        from ``_process_morphism``.
+
+        See Also
+        ========
+
+        _process_morphism
+        """
+        curving = ""
+        label_pos = "^"
+        looping_start = ""
+        looping_end = ""
+
+        # This is a loop morphism.  Count how many morphisms stick
+        # in each of the four quadrants.  Note that straight
+        # vertical and horizontal morphisms count in two quadrants
+        # at the same time (i.e., a morphism going up counts both
+        # in the first and the second quadrants).
+
+        # The usual numbering (counterclockwise) of quadrants
+        # applies.
+        quadrant = [0, 0, 0, 0]
+
+        obj = grid[i, j]
+
+        for m, m_str_info in morphisms_str_info.items():
+            if (m.domain == obj) and (m.codomain == obj):
+                # That's another loop morphism.  Check how it
+                # loops and mark the corresponding quadrants as
+                # busy.
+                (l_s, l_e) = (m_str_info.looping_start, m_str_info.looping_end)
+
+                if (l_s, l_e) == ("r", "u"):
+                    quadrant[0] += 1
+                elif (l_s, l_e) == ("u", "l"):
+                    quadrant[1] += 1
+                elif (l_s, l_e) == ("l", "d"):
+                    quadrant[2] += 1
+                elif (l_s, l_e) == ("d", "r"):
+                    quadrant[3] += 1
+
+                continue
+            if m.domain == obj:
+                (end_i, end_j) = object_coords[m.codomain]
+                goes_out = True
+            elif m.codomain == obj:
+                (end_i, end_j) = object_coords[m.domain]
+                goes_out = False
+            else:
+                continue
+
+            d_i = end_i - i
+            d_j = end_j - j
+            m_curving = m_str_info.curving
+
+            if (d_i != 0) and (d_j != 0):
+                # This is really a diagonal morphism.  Detect the
+                # quadrant.
+                if (d_i > 0) and (d_j > 0):
+                    quadrant[0] += 1
+                elif (d_i > 0) and (d_j < 0):
+                    quadrant[1] += 1
+                elif (d_i < 0) and (d_j < 0):
+                    quadrant[2] += 1
+                elif (d_i < 0) and (d_j > 0):
+                    quadrant[3] += 1
+            elif d_i == 0:
+                # Knowing where the other end of the morphism is
+                # and which way it goes, we now have to decide
+                # which quadrant is now the upper one and which is
+                # the lower one.
+                if d_j > 0:
+                    if goes_out:
+                        upper_quadrant = 0
+                        lower_quadrant = 3
+                    else:
+                        upper_quadrant = 3
+                        lower_quadrant = 0
+                else:
+                    if goes_out:
+                        upper_quadrant = 2
+                        lower_quadrant = 1
+                    else:
+                        upper_quadrant = 1
+                        lower_quadrant = 2
+
+                if m_curving:
+                    if m_curving == "^":
+                        quadrant[upper_quadrant] += 1
+                    elif m_curving == "_":
+                        quadrant[lower_quadrant] += 1
+                else:
+                    # This morphism counts in both upper and lower
+                    # quadrants.
+                    quadrant[upper_quadrant] += 1
+                    quadrant[lower_quadrant] += 1
+            elif d_j == 0:
+                # Knowing where the other end of the morphism is
+                # and which way it goes, we now have to decide
+                # which quadrant is now the left one and which is
+                # the right one.
+                if d_i < 0:
+                    if goes_out:
+                        left_quadrant = 1
+                        right_quadrant = 0
+                    else:
+                        left_quadrant = 0
+                        right_quadrant = 1
+                else:
+                    if goes_out:
+                        left_quadrant = 3
+                        right_quadrant = 2
+                    else:
+                        left_quadrant = 2
+                        right_quadrant = 3
+
+                if m_curving:
+                    if m_curving == "^":
+                        quadrant[left_quadrant] += 1
+                    elif m_curving == "_":
+                        quadrant[right_quadrant] += 1
+                else:
+                    # This morphism counts in both upper and lower
+                    # quadrants.
+                    quadrant[left_quadrant] += 1
+                    quadrant[right_quadrant] += 1
+
+        # Pick the freest quadrant to curve our morphism into.
+        freest_quadrant = 0
+        for i in range(4):
+            if quadrant[i] < quadrant[freest_quadrant]:
+                freest_quadrant = i
+
+        # Now set up proper looping.
+        (looping_start, looping_end) = [("r", "u"), ("u", "l"), ("l", "d"),
+                                        ("d", "r")][freest_quadrant]
+
+        return (curving, label_pos, looping_start, looping_end)
+
+    @staticmethod
+    def _process_horizontal_morphism(i, j, target_j, grid, morphisms_str_info,
+                                     object_coords):
+        """
+        Produces the information required for constructing the string
+        representation of a horizontal morphism.  This function is
+        invoked from ``_process_morphism``.
+
+        See Also
+        ========
+
+        _process_morphism
+        """
+        # The arrow is horizontal.  Check if it goes from left to
+        # right (``backwards == False``) or from right to left
+        # (``backwards == True``).
+        backwards = False
+        start = j
+        end = target_j
+        if end < start:
+            (start, end) = (end, start)
+            backwards = True
+
+        # Let's see which objects are there between ``start`` and
+        # ``end``, and then count how many morphisms stick out
+        # upwards, and how many stick out downwards.
+        #
+        # For example, consider the situation:
+        #
+        #    B1 C1
+        #    |  |
+        # A--B--C--D
+        #    |
+        #    B2
+        #
+        # Between the objects `A` and `D` there are two objects:
+        # `B` and `C`.  Further, there are two morphisms which
+        # stick out upward (the ones between `B1` and `B` and
+        # between `C` and `C1`) and one morphism which sticks out
+        # downward (the one between `B and `B2`).
+        #
+        # We need this information to decide how to curve the
+        # arrow between `A` and `D`.  First of all, since there
+        # are two objects between `A` and `D``, we must curve the
+        # arrow.  Then, we will have it curve downward, because
+        # there is more space (less morphisms stick out downward
+        # than upward).
+        up = []
+        down = []
+        straight_horizontal = []
+        for k in range(start + 1, end):
+            obj = grid[i, k]
+            if not obj:
+                continue
+
+            for m in morphisms_str_info:
+                if m.domain == obj:
+                    (end_i, end_j) = object_coords[m.codomain]
+                elif m.codomain == obj:
+                    (end_i, end_j) = object_coords[m.domain]
+                else:
+                    continue
+
+                if end_i > i:
+                    down.append(m)
+                elif end_i < i:
+                    up.append(m)
+                elif not morphisms_str_info[m].curving:
+                    # This is a straight horizontal morphism,
+                    # because it has no curving.
+                    straight_horizontal.append(m)
+
+        if len(up) < len(down):
+            # More morphisms stick out downward than upward, let's
+            # curve the morphism up.
+            if backwards:
+                curving = "_"
+                label_pos = "_"
+            else:
+                curving = "^"
+                label_pos = "^"
+
+            # Assure that the straight horizontal morphisms have
+            # their labels on the lower side of the arrow.
+            for m in straight_horizontal:
+                (i1, j1) = object_coords[m.domain]
+                (i2, j2) = object_coords[m.codomain]
+
+                m_str_info = morphisms_str_info[m]
+                if j1 < j2:
+                    m_str_info.label_position = "_"
+                else:
+                    m_str_info.label_position = "^"
+
+                # Don't allow any further modifications of the
+                # position of this label.
+                m_str_info.forced_label_position = True
+        else:
+            # More morphisms stick out downward than upward, let's
+            # curve the morphism up.
+            if backwards:
+                curving = "^"
+                label_pos = "^"
+            else:
+                curving = "_"
+                label_pos = "_"
+
+            # Assure that the straight horizontal morphisms have
+            # their labels on the upper side of the arrow.
+            for m in straight_horizontal:
+                (i1, j1) = object_coords[m.domain]
+                (i2, j2) = object_coords[m.codomain]
+
+                m_str_info = morphisms_str_info[m]
+                if j1 < j2:
+                    m_str_info.label_position = "^"
+                else:
+                    m_str_info.label_position = "_"
+
+                # Don't allow any further modifications of the
+                # position of this label.
+                m_str_info.forced_label_position = True
+
+        return (curving, label_pos)
+
+    @staticmethod
+    def _process_vertical_morphism(i, j, target_i, grid, morphisms_str_info,
+                                   object_coords):
+        """
+        Produces the information required for constructing the string
+        representation of a vertical morphism.  This function is
+        invoked from ``_process_morphism``.
+
+        See Also
+        ========
+
+        _process_morphism
+        """
+        # This arrow is vertical.  Check if it goes from top to
+        # bottom (``backwards == False``) or from bottom to top
+        # (``backwards == True``).
+        backwards = False
+        start = i
+        end = target_i
+        if end < start:
+            (start, end) = (end, start)
+            backwards = True
+
+        # Let's see which objects are there between ``start`` and
+        # ``end``, and then count how many morphisms stick out to
+        # the left, and how many stick out to the right.
+        #
+        # See the corresponding comment in the previous branch of
+        # this if-statement for more details.
+        left = []
+        right = []
+        straight_vertical = []
+        for k in range(start + 1, end):
+            obj = grid[k, j]
+            if not obj:
+                continue
+
+            for m in morphisms_str_info:
+                if m.domain == obj:
+                    (end_i, end_j) = object_coords[m.codomain]
+                elif m.codomain == obj:
+                    (end_i, end_j) = object_coords[m.domain]
+                else:
+                    continue
+
+                if end_j > j:
+                    right.append(m)
+                elif end_j < j:
+                    left.append(m)
+                elif not morphisms_str_info[m].curving:
+                    # This is a straight vertical morphism,
+                    # because it has no curving.
+                    straight_vertical.append(m)
+
+        if len(left) < len(right):
+            # More morphisms stick out to the left than to the
+            # right, let's curve the morphism to the right.
+            if backwards:
+                curving = "^"
+                label_pos = "^"
+            else:
+                curving = "_"
+                label_pos = "_"
+
+            # Assure that the straight vertical morphisms have
+            # their labels on the left side of the arrow.
+            for m in straight_vertical:
+                (i1, j1) = object_coords[m.domain]
+                (i2, j2) = object_coords[m.codomain]
+
+                m_str_info = morphisms_str_info[m]
+                if i1 < i2:
+                    m_str_info.label_position = "^"
+                else:
+                    m_str_info.label_position = "_"
+
+                # Don't allow any further modifications of the
+                # position of this label.
+                m_str_info.forced_label_position = True
+        else:
+            # More morphisms stick out to the right than to the
+            # left, let's curve the morphism to the left.
+            if backwards:
+                curving = "_"
+                label_pos = "_"
+            else:
+                curving = "^"
+                label_pos = "^"
+
+            # Assure that the straight vertical morphisms have
+            # their labels on the right side of the arrow.
+            for m in straight_vertical:
+                (i1, j1) = object_coords[m.domain]
+                (i2, j2) = object_coords[m.codomain]
+
+                m_str_info = morphisms_str_info[m]
+                if i1 < i2:
+                    m_str_info.label_position = "_"
+                else:
+                    m_str_info.label_position = "^"
+
+                # Don't allow any further modifications of the
+                # position of this label.
+                m_str_info.forced_label_position = True
+
+        return (curving, label_pos)
+
+    def _process_morphism(self, diagram, grid, morphism, object_coords,
+                          morphisms, morphisms_str_info):
+        """
+        Given the required information, produces the string
+        representation of ``morphism``.
+        """
+        def repeat_string_cond(times, str_gt, str_lt):
+            """
+            If ``times > 0``, repeats ``str_gt`` ``times`` times.
+            Otherwise, repeats ``str_lt`` ``-times`` times.
+            """
+            if times > 0:
+                return str_gt * times
+            else:
+                return str_lt * (-times)
+
+        def count_morphisms_undirected(A, B):
+            """
+            Counts how many processed morphisms there are between the
+            two supplied objects.
+            """
+            return len([m for m in morphisms_str_info
+                        if {m.domain, m.codomain} == {A, B}])
+
+        def count_morphisms_filtered(dom, cod, curving):
+            """
+            Counts the processed morphisms which go out of ``dom``
+            into ``cod`` with curving ``curving``.
+            """
+            return len([m for m, m_str_info in morphisms_str_info.items()
+                        if (m.domain, m.codomain) == (dom, cod) and
+                        (m_str_info.curving == curving)])
+
+        (i, j) = object_coords[morphism.domain]
+        (target_i, target_j) = object_coords[morphism.codomain]
+
+        # We now need to determine the direction of
+        # the arrow.
+        delta_i = target_i - i
+        delta_j = target_j - j
+        vertical_direction = repeat_string_cond(delta_i,
+                                                "d", "u")
+        horizontal_direction = repeat_string_cond(delta_j,
+                                                  "r", "l")
+
+        curving = ""
+        label_pos = "^"
+        looping_start = ""
+        looping_end = ""
+
+        if (delta_i == 0) and (delta_j == 0):
+            # This is a loop morphism.
+            (curving, label_pos, looping_start,
+             looping_end) = XypicDiagramDrawer._process_loop_morphism(
+                 i, j, grid, morphisms_str_info, object_coords)
+        elif (delta_i == 0) and (abs(j - target_j) > 1):
+            # This is a horizontal morphism.
+            (curving, label_pos) = XypicDiagramDrawer._process_horizontal_morphism(
+                i, j, target_j, grid, morphisms_str_info, object_coords)
+        elif (delta_j == 0) and (abs(i - target_i) > 1):
+            # This is a vertical morphism.
+            (curving, label_pos) = XypicDiagramDrawer._process_vertical_morphism(
+                i, j, target_i, grid, morphisms_str_info, object_coords)
+
+        count = count_morphisms_undirected(morphism.domain, morphism.codomain)
+        curving_amount = ""
+        if curving:
+            # This morphisms should be curved anyway.
+            curving_amount = self.default_curving_amount + count * \
+                self.default_curving_step
+        elif count:
+            # There are no objects between the domain and codomain of
+            # the current morphism, but this is not there already are
+            # some morphisms with the same domain and codomain, so we
+            # have to curve this one.
+            curving = "^"
+            filtered_morphisms = count_morphisms_filtered(
+                morphism.domain, morphism.codomain, curving)
+            curving_amount = self.default_curving_amount + \
+                filtered_morphisms * \
+                self.default_curving_step
+
+        # Let's now get the name of the morphism.
+        morphism_name = ""
+        if isinstance(morphism, IdentityMorphism):
+            morphism_name = "id_{%s}" + latex(grid[i, j])
+        elif isinstance(morphism, CompositeMorphism):
+            component_names = [latex(Symbol(component.name)) for
+                               component in morphism.components]
+            component_names.reverse()
+            morphism_name = "\\circ ".join(component_names)
+        elif isinstance(morphism, NamedMorphism):
+            morphism_name = latex(Symbol(morphism.name))
+
+        return ArrowStringDescription(
+            self.unit, curving, curving_amount, looping_start,
+            looping_end, horizontal_direction, vertical_direction,
+            label_pos, morphism_name)
+
+    @staticmethod
+    def _check_free_space_horizontal(dom_i, dom_j, cod_j, grid):
+        """
+        For a horizontal morphism, checks whether there is free space
+        (i.e., space not occupied by any objects) above the morphism
+        or below it.
+        """
+        if dom_j < cod_j:
+            (start, end) = (dom_j, cod_j)
+            backwards = False
+        else:
+            (start, end) = (cod_j, dom_j)
+            backwards = True
+
+        # Check for free space above.
+        if dom_i == 0:
+            free_up = True
+        else:
+            free_up = all(grid[dom_i - 1, j] for j in
+                          range(start, end + 1))
+
+        # Check for free space below.
+        if dom_i == grid.height - 1:
+            free_down = True
+        else:
+            free_down = not any(grid[dom_i + 1, j] for j in
+                                range(start, end + 1))
+
+        return (free_up, free_down, backwards)
+
+    @staticmethod
+    def _check_free_space_vertical(dom_i, cod_i, dom_j, grid):
+        """
+        For a vertical morphism, checks whether there is free space
+        (i.e., space not occupied by any objects) to the left of the
+        morphism or to the right of it.
+        """
+        if dom_i < cod_i:
+            (start, end) = (dom_i, cod_i)
+            backwards = False
+        else:
+            (start, end) = (cod_i, dom_i)
+            backwards = True
+
+        # Check if there's space to the left.
+        if dom_j == 0:
+            free_left = True
+        else:
+            free_left = not any(grid[i, dom_j - 1] for i in
+                                range(start, end + 1))
+
+        if dom_j == grid.width - 1:
+            free_right = True
+        else:
+            free_right = not any(grid[i, dom_j + 1] for i in
+                                 range(start, end + 1))
+
+        return (free_left, free_right, backwards)
+
+    @staticmethod
+    def _check_free_space_diagonal(dom_i, cod_i, dom_j, cod_j, grid):
+        """
+        For a diagonal morphism, checks whether there is free space
+        (i.e., space not occupied by any objects) above the morphism
+        or below it.
+        """
+        def abs_xrange(start, end):
+            if start < end:
+                return range(start, end + 1)
+            else:
+                return range(end, start + 1)
+
+        if dom_i < cod_i and dom_j < cod_j:
+            # This morphism goes from top-left to
+            # bottom-right.
+            (start_i, start_j) = (dom_i, dom_j)
+            (end_i, end_j) = (cod_i, cod_j)
+            backwards = False
+        elif dom_i > cod_i and dom_j > cod_j:
+            # This morphism goes from bottom-right to
+            # top-left.
+            (start_i, start_j) = (cod_i, cod_j)
+            (end_i, end_j) = (dom_i, dom_j)
+            backwards = True
+        if dom_i < cod_i and dom_j > cod_j:
+            # This morphism goes from top-right to
+            # bottom-left.
+            (start_i, start_j) = (dom_i, dom_j)
+            (end_i, end_j) = (cod_i, cod_j)
+            backwards = True
+        elif dom_i > cod_i and dom_j < cod_j:
+            # This morphism goes from bottom-left to
+            # top-right.
+            (start_i, start_j) = (cod_i, cod_j)
+            (end_i, end_j) = (dom_i, dom_j)
+            backwards = False
+
+        # This is an attempt at a fast and furious strategy to
+        # decide where there is free space on the two sides of
+        # a diagonal morphism.  For a diagonal morphism
+        # starting at ``(start_i, start_j)`` and ending at
+        # ``(end_i, end_j)`` the rectangle defined by these
+        # two points is considered.  The slope of the diagonal
+        # ``alpha`` is then computed.  Then, for every cell
+        # ``(i, j)`` within the rectangle, the slope
+        # ``alpha1`` of the line through ``(start_i,
+        # start_j)`` and ``(i, j)`` is considered.  If
+        # ``alpha1`` is between 0 and ``alpha``, the point
+        # ``(i, j)`` is above the diagonal, if ``alpha1`` is
+        # between ``alpha`` and infinity, the point is below
+        # the diagonal.  Also note that, with some beforehand
+        # precautions, this trick works for both the main and
+        # the secondary diagonals of the rectangle.
+
+        # I have considered the possibility to only follow the
+        # shorter diagonals immediately above and below the
+        # main (or secondary) diagonal.  This, however,
+        # wouldn't have resulted in much performance gain or
+        # better detection of outer edges, because of
+        # relatively small sizes of diagram grids, while the
+        # code would have become harder to understand.
+
+        alpha = float(end_i - start_i)/(end_j - start_j)
+        free_up = True
+        free_down = True
+        for i in abs_xrange(start_i, end_i):
+            if not free_up and not free_down:
+                break
+
+            for j in abs_xrange(start_j, end_j):
+                if not free_up and not free_down:
+                    break
+
+                if (i, j) == (start_i, start_j):
+                    continue
+
+                if j == start_j:
+                    alpha1 = "inf"
+                else:
+                    alpha1 = float(i - start_i)/(j - start_j)
+
+                if grid[i, j]:
+                    if (alpha1 == "inf") or (abs(alpha1) > abs(alpha)):
+                        free_down = False
+                    elif abs(alpha1) < abs(alpha):
+                        free_up = False
+
+        return (free_up, free_down, backwards)
+
+    def _push_labels_out(self, morphisms_str_info, grid, object_coords):
+        """
+        For all straight morphisms which form the visual boundary of
+        the laid out diagram, puts their labels on their outer sides.
+        """
+        def set_label_position(free1, free2, pos1, pos2, backwards, m_str_info):
+            """
+            Given the information about room available to one side and
+            to the other side of a morphism (``free1`` and ``free2``),
+            sets the position of the morphism label in such a way that
+            it is on the freer side.  This latter operations involves
+            choice between ``pos1`` and ``pos2``, taking ``backwards``
+            in consideration.
+
+            Thus this function will do nothing if either both ``free1
+            == True`` and ``free2 == True`` or both ``free1 == False``
+            and ``free2 == False``.  In either case, choosing one side
+            over the other presents no advantage.
+            """
+            if backwards:
+                (pos1, pos2) = (pos2, pos1)
+
+            if free1 and not free2:
+                m_str_info.label_position = pos1
+            elif free2 and not free1:
+                m_str_info.label_position = pos2
+
+        for m, m_str_info in morphisms_str_info.items():
+            if m_str_info.curving or m_str_info.forced_label_position:
+                # This is either a curved morphism, and curved
+                # morphisms have other magic, or the position of this
+                # label has already been fixed.
+                continue
+
+            if m.domain == m.codomain:
+                # This is a loop morphism, their labels, again have a
+                # different magic.
+                continue
+
+            (dom_i, dom_j) = object_coords[m.domain]
+            (cod_i, cod_j) = object_coords[m.codomain]
+
+            if dom_i == cod_i:
+                # Horizontal morphism.
+                (free_up, free_down,
+                 backwards) = XypicDiagramDrawer._check_free_space_horizontal(
+                     dom_i, dom_j, cod_j, grid)
+
+                set_label_position(free_up, free_down, "^", "_",
+                                   backwards, m_str_info)
+            elif dom_j == cod_j:
+                # Vertical morphism.
+                (free_left, free_right,
+                 backwards) = XypicDiagramDrawer._check_free_space_vertical(
+                     dom_i, cod_i, dom_j, grid)
+
+                set_label_position(free_left, free_right, "_", "^",
+                                   backwards, m_str_info)
+            else:
+                # A diagonal morphism.
+                (free_up, free_down,
+                 backwards) = XypicDiagramDrawer._check_free_space_diagonal(
+                     dom_i, cod_i, dom_j, cod_j, grid)
+
+                set_label_position(free_up, free_down, "^", "_",
+                                   backwards, m_str_info)
+
+    @staticmethod
+    def _morphism_sort_key(morphism, object_coords):
+        """
+        Provides a morphism sorting key such that horizontal or
+        vertical morphisms between neighbouring objects come
+        first, then horizontal or vertical morphisms between more
+        far away objects, and finally, all other morphisms.
+        """
+        (i, j) = object_coords[morphism.domain]
+        (target_i, target_j) = object_coords[morphism.codomain]
+
+        if morphism.domain == morphism.codomain:
+            # Loop morphisms should get after diagonal morphisms
+            # so that the proper direction in which to curve the
+            # loop can be determined.
+            return (3, 0, default_sort_key(morphism))
+
+        if target_i == i:
+            return (1, abs(target_j - j), default_sort_key(morphism))
+
+        if target_j == j:
+            return (1, abs(target_i - i), default_sort_key(morphism))
+
+        # Diagonal morphism.
+        return (2, 0, default_sort_key(morphism))
+
+    @staticmethod
+    def _build_xypic_string(diagram, grid, morphisms,
+                            morphisms_str_info, diagram_format):
+        """
+        Given a collection of :class:`ArrowStringDescription`
+        describing the morphisms of a diagram and the object layout
+        information of a diagram, produces the final Xy-pic picture.
+        """
+        # Build the mapping between objects and morphisms which have
+        # them as domains.
+        object_morphisms = {}
+        for obj in diagram.objects:
+            object_morphisms[obj] = []
+        for morphism in morphisms:
+            object_morphisms[morphism.domain].append(morphism)
+
+        result = "\\xymatrix%s{\n" % diagram_format
+
+        for i in range(grid.height):
+            for j in range(grid.width):
+                obj = grid[i, j]
+                if obj:
+                    result += latex(obj) + " "
+
+                    morphisms_to_draw = object_morphisms[obj]
+                    for morphism in morphisms_to_draw:
+                        result += str(morphisms_str_info[morphism]) + " "
+
+                # Don't put the & after the last column.
+                if j < grid.width - 1:
+                    result += "& "
+
+            # Don't put the line break after the last row.
+            if i < grid.height - 1:
+                result += "\\\\"
+            result += "\n"
+
+        result += "}\n"
+
+        return result
+
+    def draw(self, diagram, grid, masked=None, diagram_format=""):
+        r"""
+        Returns the Xy-pic representation of ``diagram`` laid out in
+        ``grid``.
+
+        Consider the following simple triangle diagram.
+
+        >>> from sympy.categories import Object, NamedMorphism, Diagram
+        >>> from sympy.categories import DiagramGrid, XypicDiagramDrawer
+        >>> A = Object("A")
+        >>> B = Object("B")
+        >>> C = Object("C")
+        >>> f = NamedMorphism(A, B, "f")
+        >>> g = NamedMorphism(B, C, "g")
+        >>> diagram = Diagram([f, g], {g * f: "unique"})
+
+        To draw this diagram, its objects need to be laid out with a
+        :class:`DiagramGrid`::
+
+        >>> grid = DiagramGrid(diagram)
+
+        Finally, the drawing:
+
+        >>> drawer = XypicDiagramDrawer()
+        >>> print(drawer.draw(diagram, grid))
+        \xymatrix{
+        A \ar[d]_{g\circ f} \ar[r]^{f} & B \ar[ld]^{g} \\
+        C &
+        }
+
+        The argument ``masked`` can be used to skip morphisms in the
+        presentation of the diagram:
+
+        >>> print(drawer.draw(diagram, grid, masked=[g * f]))
+        \xymatrix{
+        A \ar[r]^{f} & B \ar[ld]^{g} \\
+        C &
+        }
+
+        Finally, the ``diagram_format`` argument can be used to
+        specify the format string of the diagram.  For example, to
+        increase the spacing by 1 cm, proceeding as follows:
+
+        >>> print(drawer.draw(diagram, grid, diagram_format="@+1cm"))
+        \xymatrix@+1cm{
+        A \ar[d]_{g\circ f} \ar[r]^{f} & B \ar[ld]^{g} \\
+        C &
+        }
+
+        """
+        # This method works in several steps.  It starts by removing
+        # the masked morphisms, if necessary, and then maps objects to
+        # their positions in the grid (coordinate tuples).  Remember
+        # that objects are unique in ``Diagram`` and in the layout
+        # produced by ``DiagramGrid``, so every object is mapped to a
+        # single coordinate pair.
+        #
+        # The next step is the central step and is concerned with
+        # analysing the morphisms of the diagram and deciding how to
+        # draw them.  For example, how to curve the arrows is decided
+        # at this step.  The bulk of the analysis is implemented in
+        # ``_process_morphism``, to the result of which the
+        # appropriate formatters are applied.
+        #
+        # The result of the previous step is a list of
+        # ``ArrowStringDescription``.  After the analysis and
+        # application of formatters, some extra logic tries to assure
+        # better positioning of morphism labels (for example, an
+        # attempt is made to avoid the situations when arrows cross
+        # labels).  This functionality constitutes the next step and
+        # is implemented in ``_push_labels_out``.  Note that label
+        # positions which have been set via a formatter are not
+        # affected in this step.
+        #
+        # Finally, at the closing step, the array of
+        # ``ArrowStringDescription`` and the layout information
+        # incorporated in ``DiagramGrid`` are combined to produce the
+        # resulting Xy-pic picture.  This part of code lies in
+        # ``_build_xypic_string``.
+
+        if not masked:
+            morphisms_props = grid.morphisms
+        else:
+            morphisms_props = {}
+            for m, props in grid.morphisms.items():
+                if m in masked:
+                    continue
+                morphisms_props[m] = props
+
+        # Build the mapping between objects and their position in the
+        # grid.
+        object_coords = {}
+        for i in range(grid.height):
+            for j in range(grid.width):
+                if grid[i, j]:
+                    object_coords[grid[i, j]] = (i, j)
+
+        morphisms = sorted(morphisms_props,
+                           key=lambda m: XypicDiagramDrawer._morphism_sort_key(
+                               m, object_coords))
+
+        # Build the tuples defining the string representations of
+        # morphisms.
+        morphisms_str_info = {}
+        for morphism in morphisms:
+            string_description = self._process_morphism(
+                diagram, grid, morphism, object_coords, morphisms,
+                morphisms_str_info)
+
+            if self.default_arrow_formatter:
+                self.default_arrow_formatter(string_description)
+
+            for prop in morphisms_props[morphism]:
+                # prop is a Symbol.  TODO: Find out why.
+                if prop.name in self.arrow_formatters:
+                    formatter = self.arrow_formatters[prop.name]
+                    formatter(string_description)
+
+            morphisms_str_info[morphism] = string_description
+
+        # Reposition the labels a bit.
+        self._push_labels_out(morphisms_str_info, grid, object_coords)
+
+        return XypicDiagramDrawer._build_xypic_string(
+            diagram, grid, morphisms, morphisms_str_info, diagram_format)
+
+
+def xypic_draw_diagram(diagram, masked=None, diagram_format="",
+                       groups=None, **hints):
+    r"""
+    Provides a shortcut combining :class:`DiagramGrid` and
+    :class:`XypicDiagramDrawer`.  Returns an Xy-pic presentation of
+    ``diagram``.  The argument ``masked`` is a list of morphisms which
+    will be not be drawn.  The argument ``diagram_format`` is the
+    format string inserted after "\xymatrix".  ``groups`` should be a
+    set of logical groups.  The ``hints`` will be passed directly to
+    the constructor of :class:`DiagramGrid`.
+
+    For more information about the arguments, see the docstrings of
+    :class:`DiagramGrid` and ``XypicDiagramDrawer.draw``.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism, Diagram
+    >>> from sympy.categories import xypic_draw_diagram
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> diagram = Diagram([f, g], {g * f: "unique"})
+    >>> print(xypic_draw_diagram(diagram))
+    \xymatrix{
+    A \ar[d]_{g\circ f} \ar[r]^{f} & B \ar[ld]^{g} \\
+    C &
+    }
+
+    See Also
+    ========
+
+    XypicDiagramDrawer, DiagramGrid
+    """
+    grid = DiagramGrid(diagram, groups, **hints)
+    drawer = XypicDiagramDrawer()
+    return drawer.draw(diagram, grid, masked, diagram_format)
+
+
+@doctest_depends_on(exe=('latex', 'dvipng'), modules=('pyglet',))
+def preview_diagram(diagram, masked=None, diagram_format="", groups=None,
+                    output='png', viewer=None, euler=True, **hints):
+    """
+    Combines the functionality of ``xypic_draw_diagram`` and
+    ``sympy.printing.preview``.  The arguments ``masked``,
+    ``diagram_format``, ``groups``, and ``hints`` are passed to
+    ``xypic_draw_diagram``, while ``output``, ``viewer, and ``euler``
+    are passed to ``preview``.
+
+    Examples
+    ========
+
+    >>> from sympy.categories import Object, NamedMorphism, Diagram
+    >>> from sympy.categories import preview_diagram
+    >>> A = Object("A")
+    >>> B = Object("B")
+    >>> C = Object("C")
+    >>> f = NamedMorphism(A, B, "f")
+    >>> g = NamedMorphism(B, C, "g")
+    >>> d = Diagram([f, g], {g * f: "unique"})
+    >>> preview_diagram(d)
+
+    See Also
+    ========
+
+    XypicDiagramDrawer
+    """
+    from sympy.printing import preview
+    latex_output = xypic_draw_diagram(diagram, masked, diagram_format,
+                                      groups, **hints)
+    preview(latex_output, output, viewer, euler, ("xypic",))

lib/python3.12/site-packages/sympy/categories/tests/test_baseclasses.py

@@ -0,0 +1,209 @@
+from sympy.categories import (Object, Morphism, IdentityMorphism,
+                              NamedMorphism, CompositeMorphism,
+                              Diagram, Category)
+from sympy.categories.baseclasses import Class
+from sympy.testing.pytest import raises
+from sympy.core.containers import (Dict, Tuple)
+from sympy.sets import EmptySet
+from sympy.sets.sets import FiniteSet
+
+
+def test_morphisms():
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+
+    # Test the base morphism.
+    f = NamedMorphism(A, B, "f")
+    assert f.domain == A
+    assert f.codomain == B
+    assert f == NamedMorphism(A, B, "f")
+
+    # Test identities.
+    id_A = IdentityMorphism(A)
+    id_B = IdentityMorphism(B)
+    assert id_A.domain == A
+    assert id_A.codomain == A
+    assert id_A == IdentityMorphism(A)
+    assert id_A != id_B
+
+    # Test named morphisms.
+    g = NamedMorphism(B, C, "g")
+    assert g.name == "g"
+    assert g != f
+    assert g == NamedMorphism(B, C, "g")
+    assert g != NamedMorphism(B, C, "f")
+
+    # Test composite morphisms.
+    assert f == CompositeMorphism(f)
+
+    k = g.compose(f)
+    assert k.domain == A
+    assert k.codomain == C
+    assert k.components == Tuple(f, g)
+    assert g * f == k
+    assert CompositeMorphism(f, g) == k
+
+    assert CompositeMorphism(g * f) == g * f
+
+    # Test the associativity of composition.
+    h = NamedMorphism(C, D, "h")
+
+    p = h * g
+    u = h * g * f
+
+    assert h * k == u
+    assert p * f == u
+    assert CompositeMorphism(f, g, h) == u
+
+    # Test flattening.
+    u2 = u.flatten("u")
+    assert isinstance(u2, NamedMorphism)
+    assert u2.name == "u"
+    assert u2.domain == A
+    assert u2.codomain == D
+
+    # Test identities.
+    assert f * id_A == f
+    assert id_B * f == f
+    assert id_A * id_A == id_A
+    assert CompositeMorphism(id_A) == id_A
+
+    # Test bad compositions.
+    raises(ValueError, lambda: f * g)
+
+    raises(TypeError, lambda: f.compose(None))
+    raises(TypeError, lambda: id_A.compose(None))
+    raises(TypeError, lambda: f * None)
+    raises(TypeError, lambda: id_A * None)
+
+    raises(TypeError, lambda: CompositeMorphism(f, None, 1))
+
+    raises(ValueError, lambda: NamedMorphism(A, B, ""))
+    raises(NotImplementedError, lambda: Morphism(A, B))
+
+
+def test_diagram():
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    id_A = IdentityMorphism(A)
+    id_B = IdentityMorphism(B)
+
+    empty = EmptySet
+
+    # Test the addition of identities.
+    d1 = Diagram([f])
+
+    assert d1.objects == FiniteSet(A, B)
+    assert d1.hom(A, B) == (FiniteSet(f), empty)
+    assert d1.hom(A, A) == (FiniteSet(id_A), empty)
+    assert d1.hom(B, B) == (FiniteSet(id_B), empty)
+
+    assert d1 == Diagram([id_A, f])
+    assert d1 == Diagram([f, f])
+
+    # Test the addition of composites.
+    d2 = Diagram([f, g])
+    homAC = d2.hom(A, C)[0]
+
+    assert d2.objects == FiniteSet(A, B, C)
+    assert g * f in d2.premises.keys()
+    assert homAC == FiniteSet(g * f)
+
+    # Test equality, inequality and hash.
+    d11 = Diagram([f])
+
+    assert d1 == d11
+    assert d1 != d2
+    assert hash(d1) == hash(d11)
+
+    d11 = Diagram({f: "unique"})
+    assert d1 != d11
+
+    # Make sure that (re-)adding composites (with new properties)
+    # works as expected.
+    d = Diagram([f, g], {g * f: "unique"})
+    assert d.conclusions == Dict({g * f: FiniteSet("unique")})
+
+    # Check the hom-sets when there are premises and conclusions.
+    assert d.hom(A, C) == (FiniteSet(g * f), FiniteSet(g * f))
+    d = Diagram([f, g], [g * f])
+    assert d.hom(A, C) == (FiniteSet(g * f), FiniteSet(g * f))
+
+    # Check how the properties of composite morphisms are computed.
+    d = Diagram({f: ["unique", "isomorphism"], g: "unique"})
+    assert d.premises[g * f] == FiniteSet("unique")
+
+    # Check that conclusion morphisms with new objects are not allowed.
+    d = Diagram([f], [g])
+    assert d.conclusions == Dict({})
+
+    # Test an empty diagram.
+    d = Diagram()
+    assert d.premises == Dict({})
+    assert d.conclusions == Dict({})
+    assert d.objects == empty
+
+    # Check a SymPy Dict object.
+    d = Diagram(Dict({f: FiniteSet("unique", "isomorphism"), g: "unique"}))
+    assert d.premises[g * f] == FiniteSet("unique")
+
+    # Check the addition of components of composite morphisms.
+    d = Diagram([g * f])
+    assert f in d.premises
+    assert g in d.premises
+
+    # Check subdiagrams.
+    d = Diagram([f, g], {g * f: "unique"})
+
+    d1 = Diagram([f])
+    assert d.is_subdiagram(d1)
+    assert not d1.is_subdiagram(d)
+
+    d = Diagram([NamedMorphism(B, A, "f'")])
+    assert not d.is_subdiagram(d1)
+    assert not d1.is_subdiagram(d)
+
+    d1 = Diagram([f, g], {g * f: ["unique", "something"]})
+    assert not d.is_subdiagram(d1)
+    assert not d1.is_subdiagram(d)
+
+    d = Diagram({f: "blooh"})
+    d1 = Diagram({f: "bleeh"})
+    assert not d.is_subdiagram(d1)
+    assert not d1.is_subdiagram(d)
+
+    d = Diagram([f, g], {f: "unique", g * f: "veryunique"})
+    d1 = d.subdiagram_from_objects(FiniteSet(A, B))
+    assert d1 == Diagram([f], {f: "unique"})
+    raises(ValueError, lambda: d.subdiagram_from_objects(FiniteSet(A,
+           Object("D"))))
+
+    raises(ValueError, lambda: Diagram({IdentityMorphism(A): "unique"}))
+
+
+def test_category():
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+
+    d1 = Diagram([f, g])
+    d2 = Diagram([f])
+
+    objects = d1.objects | d2.objects
+
+    K = Category("K", objects, commutative_diagrams=[d1, d2])
+
+    assert K.name == "K"
+    assert K.objects == Class(objects)
+    assert K.commutative_diagrams == FiniteSet(d1, d2)
+
+    raises(ValueError, lambda: Category(""))

lib/python3.12/site-packages/sympy/categories/tests/test_drawing.py

@@ -0,0 +1,919 @@
+from sympy.categories.diagram_drawing import _GrowableGrid, ArrowStringDescription
+from sympy.categories import (DiagramGrid, Object, NamedMorphism,
+                              Diagram, XypicDiagramDrawer, xypic_draw_diagram)
+from sympy.sets.sets import FiniteSet
+
+
+def test_GrowableGrid():
+    grid = _GrowableGrid(1, 2)
+
+    # Check dimensions.
+    assert grid.width == 1
+    assert grid.height == 2
+
+    # Check initialization of elements.
+    assert grid[0, 0] is None
+    assert grid[1, 0] is None
+
+    # Check assignment to elements.
+    grid[0, 0] = 1
+    grid[1, 0] = "two"
+
+    assert grid[0, 0] == 1
+    assert grid[1, 0] == "two"
+
+    # Check appending a row.
+    grid.append_row()
+
+    assert grid.width == 1
+    assert grid.height == 3
+
+    assert grid[0, 0] == 1
+    assert grid[1, 0] == "two"
+    assert grid[2, 0] is None
+
+    # Check appending a column.
+    grid.append_column()
+    assert grid.width == 2
+    assert grid.height == 3
+
+    assert grid[0, 0] == 1
+    assert grid[1, 0] == "two"
+    assert grid[2, 0] is None
+
+    assert grid[0, 1] is None
+    assert grid[1, 1] is None
+    assert grid[2, 1] is None
+
+    grid = _GrowableGrid(1, 2)
+    grid[0, 0] = 1
+    grid[1, 0] = "two"
+
+    # Check prepending a row.
+    grid.prepend_row()
+    assert grid.width == 1
+    assert grid.height == 3
+
+    assert grid[0, 0] is None
+    assert grid[1, 0] == 1
+    assert grid[2, 0] == "two"
+
+    # Check prepending a column.
+    grid.prepend_column()
+    assert grid.width == 2
+    assert grid.height == 3
+
+    assert grid[0, 0] is None
+    assert grid[1, 0] is None
+    assert grid[2, 0] is None
+
+    assert grid[0, 1] is None
+    assert grid[1, 1] == 1
+    assert grid[2, 1] == "two"
+
+
+def test_DiagramGrid():
+    # Set up some objects and morphisms.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(D, A, "h")
+    k = NamedMorphism(D, B, "k")
+
+    # A one-morphism diagram.
+    d = Diagram([f])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 2
+    assert grid.height == 1
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid.morphisms == {f: FiniteSet()}
+
+    # A triangle.
+    d = Diagram([f, g], {g * f: "unique"})
+    grid = DiagramGrid(d)
+
+    assert grid.width == 2
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[1, 0] == C
+    assert grid[1, 1] is None
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(),
+                              g * f: FiniteSet("unique")}
+
+    # A triangle with a "loop" morphism.
+    l_A = NamedMorphism(A, A, "l_A")
+    d = Diagram([f, g, l_A])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 2
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[1, 0] is None
+    assert grid[1, 1] == C
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), l_A: FiniteSet()}
+
+    # A simple diagram.
+    d = Diagram([f, g, h, k])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 3
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] == D
+    assert grid[1, 0] is None
+    assert grid[1, 1] == C
+    assert grid[1, 2] is None
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), h: FiniteSet(),
+                              k: FiniteSet()}
+
+    assert str(grid) == '[[Object("A"), Object("B"), Object("D")], ' \
+        '[None, Object("C"), None]]'
+
+    # A chain of morphisms.
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(C, D, "h")
+    k = NamedMorphism(D, E, "k")
+    d = Diagram([f, g, h, k])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 3
+    assert grid.height == 3
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] is None
+    assert grid[1, 0] is None
+    assert grid[1, 1] == C
+    assert grid[1, 2] == D
+    assert grid[2, 0] is None
+    assert grid[2, 1] is None
+    assert grid[2, 2] == E
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), h: FiniteSet(),
+                              k: FiniteSet()}
+
+    # A square.
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, D, "g")
+    h = NamedMorphism(A, C, "h")
+    k = NamedMorphism(C, D, "k")
+    d = Diagram([f, g, h, k])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 2
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[1, 0] == C
+    assert grid[1, 1] == D
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), h: FiniteSet(),
+                              k: FiniteSet()}
+
+    # A strange diagram which resulted from a typo when creating a
+    # test for five lemma, but which allowed to stop one extra problem
+    # in the algorithm.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+    A_ = Object("A'")
+    B_ = Object("B'")
+    C_ = Object("C'")
+    D_ = Object("D'")
+    E_ = Object("E'")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(C, D, "h")
+    i = NamedMorphism(D, E, "i")
+
+    # These 4 morphisms should be between primed objects.
+    j = NamedMorphism(A, B, "j")
+    k = NamedMorphism(B, C, "k")
+    l = NamedMorphism(C, D, "l")
+    m = NamedMorphism(D, E, "m")
+
+    o = NamedMorphism(A, A_, "o")
+    p = NamedMorphism(B, B_, "p")
+    q = NamedMorphism(C, C_, "q")
+    r = NamedMorphism(D, D_, "r")
+    s = NamedMorphism(E, E_, "s")
+
+    d = Diagram([f, g, h, i, j, k, l, m, o, p, q, r, s])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 3
+    assert grid.height == 4
+    assert grid[0, 0] is None
+    assert grid[0, 1] == A
+    assert grid[0, 2] == A_
+    assert grid[1, 0] == C
+    assert grid[1, 1] == B
+    assert grid[1, 2] == B_
+    assert grid[2, 0] == C_
+    assert grid[2, 1] == D
+    assert grid[2, 2] == D_
+    assert grid[3, 0] is None
+    assert grid[3, 1] == E
+    assert grid[3, 2] == E_
+
+    morphisms = {}
+    for m in [f, g, h, i, j, k, l, m, o, p, q, r, s]:
+        morphisms[m] = FiniteSet()
+    assert grid.morphisms == morphisms
+
+    # A cube.
+    A1 = Object("A1")
+    A2 = Object("A2")
+    A3 = Object("A3")
+    A4 = Object("A4")
+    A5 = Object("A5")
+    A6 = Object("A6")
+    A7 = Object("A7")
+    A8 = Object("A8")
+
+    # The top face of the cube.
+    f1 = NamedMorphism(A1, A2, "f1")
+    f2 = NamedMorphism(A1, A3, "f2")
+    f3 = NamedMorphism(A2, A4, "f3")
+    f4 = NamedMorphism(A3, A4, "f3")
+
+    # The bottom face of the cube.
+    f5 = NamedMorphism(A5, A6, "f5")
+    f6 = NamedMorphism(A5, A7, "f6")
+    f7 = NamedMorphism(A6, A8, "f7")
+    f8 = NamedMorphism(A7, A8, "f8")
+
+    # The remaining morphisms.
+    f9 = NamedMorphism(A1, A5, "f9")
+    f10 = NamedMorphism(A2, A6, "f10")
+    f11 = NamedMorphism(A3, A7, "f11")
+    f12 = NamedMorphism(A4, A8, "f11")
+
+    d = Diagram([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 4
+    assert grid.height == 3
+    assert grid[0, 0] is None
+    assert grid[0, 1] == A5
+    assert grid[0, 2] == A6
+    assert grid[0, 3] is None
+    assert grid[1, 0] is None
+    assert grid[1, 1] == A1
+    assert grid[1, 2] == A2
+    assert grid[1, 3] is None
+    assert grid[2, 0] == A7
+    assert grid[2, 1] == A3
+    assert grid[2, 2] == A4
+    assert grid[2, 3] == A8
+
+    morphisms = {}
+    for m in [f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12]:
+        morphisms[m] = FiniteSet()
+    assert grid.morphisms == morphisms
+
+    # A line diagram.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(C, D, "h")
+    i = NamedMorphism(D, E, "i")
+    d = Diagram([f, g, h, i])
+    grid = DiagramGrid(d, layout="sequential")
+
+    assert grid.width == 5
+    assert grid.height == 1
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] == C
+    assert grid[0, 3] == D
+    assert grid[0, 4] == E
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), h: FiniteSet(),
+                              i: FiniteSet()}
+
+    # Test the transposed version.
+    grid = DiagramGrid(d, layout="sequential", transpose=True)
+
+    assert grid.width == 1
+    assert grid.height == 5
+    assert grid[0, 0] == A
+    assert grid[1, 0] == B
+    assert grid[2, 0] == C
+    assert grid[3, 0] == D
+    assert grid[4, 0] == E
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), h: FiniteSet(),
+                              i: FiniteSet()}
+
+    # A pullback.
+    m1 = NamedMorphism(A, B, "m1")
+    m2 = NamedMorphism(A, C, "m2")
+    s1 = NamedMorphism(B, D, "s1")
+    s2 = NamedMorphism(C, D, "s2")
+    f1 = NamedMorphism(E, B, "f1")
+    f2 = NamedMorphism(E, C, "f2")
+    g = NamedMorphism(E, A, "g")
+
+    d = Diagram([m1, m2, s1, s2, f1, f2], {g: "unique"})
+    grid = DiagramGrid(d)
+
+    assert grid.width == 3
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] == E
+    assert grid[1, 0] == C
+    assert grid[1, 1] == D
+    assert grid[1, 2] is None
+
+    morphisms = {g: FiniteSet("unique")}
+    for m in [m1, m2, s1, s2, f1, f2]:
+        morphisms[m] = FiniteSet()
+    assert grid.morphisms == morphisms
+
+    # Test the pullback with sequential layout, just for stress
+    # testing.
+    grid = DiagramGrid(d, layout="sequential")
+
+    assert grid.width == 5
+    assert grid.height == 1
+    assert grid[0, 0] == D
+    assert grid[0, 1] == B
+    assert grid[0, 2] == A
+    assert grid[0, 3] == C
+    assert grid[0, 4] == E
+    assert grid.morphisms == morphisms
+
+    # Test a pullback with object grouping.
+    grid = DiagramGrid(d, groups=FiniteSet(E, FiniteSet(A, B, C, D)))
+
+    assert grid.width == 3
+    assert grid.height == 2
+    assert grid[0, 0] == E
+    assert grid[0, 1] == A
+    assert grid[0, 2] == B
+    assert grid[1, 0] is None
+    assert grid[1, 1] == C
+    assert grid[1, 2] == D
+    assert grid.morphisms == morphisms
+
+    # Five lemma, actually.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+    A_ = Object("A'")
+    B_ = Object("B'")
+    C_ = Object("C'")
+    D_ = Object("D'")
+    E_ = Object("E'")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(C, D, "h")
+    i = NamedMorphism(D, E, "i")
+
+    j = NamedMorphism(A_, B_, "j")
+    k = NamedMorphism(B_, C_, "k")
+    l = NamedMorphism(C_, D_, "l")
+    m = NamedMorphism(D_, E_, "m")
+
+    o = NamedMorphism(A, A_, "o")
+    p = NamedMorphism(B, B_, "p")
+    q = NamedMorphism(C, C_, "q")
+    r = NamedMorphism(D, D_, "r")
+    s = NamedMorphism(E, E_, "s")
+
+    d = Diagram([f, g, h, i, j, k, l, m, o, p, q, r, s])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 5
+    assert grid.height == 3
+    assert grid[0, 0] is None
+    assert grid[0, 1] == A
+    assert grid[0, 2] == A_
+    assert grid[0, 3] is None
+    assert grid[0, 4] is None
+    assert grid[1, 0] == C
+    assert grid[1, 1] == B
+    assert grid[1, 2] == B_
+    assert grid[1, 3] == C_
+    assert grid[1, 4] is None
+    assert grid[2, 0] == D
+    assert grid[2, 1] == E
+    assert grid[2, 2] is None
+    assert grid[2, 3] == D_
+    assert grid[2, 4] == E_
+
+    morphisms = {}
+    for m in [f, g, h, i, j, k, l, m, o, p, q, r, s]:
+        morphisms[m] = FiniteSet()
+    assert grid.morphisms == morphisms
+
+    # Test the five lemma with object grouping.
+    grid = DiagramGrid(d, FiniteSet(
+        FiniteSet(A, B, C, D, E), FiniteSet(A_, B_, C_, D_, E_)))
+
+    assert grid.width == 6
+    assert grid.height == 3
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] is None
+    assert grid[0, 3] == A_
+    assert grid[0, 4] == B_
+    assert grid[0, 5] is None
+    assert grid[1, 0] is None
+    assert grid[1, 1] == C
+    assert grid[1, 2] == D
+    assert grid[1, 3] is None
+    assert grid[1, 4] == C_
+    assert grid[1, 5] == D_
+    assert grid[2, 0] is None
+    assert grid[2, 1] is None
+    assert grid[2, 2] == E
+    assert grid[2, 3] is None
+    assert grid[2, 4] is None
+    assert grid[2, 5] == E_
+    assert grid.morphisms == morphisms
+
+    # Test the five lemma with object grouping, but mixing containers
+    # to represent groups.
+    grid = DiagramGrid(d, [(A, B, C, D, E), {A_, B_, C_, D_, E_}])
+
+    assert grid.width == 6
+    assert grid.height == 3
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] is None
+    assert grid[0, 3] == A_
+    assert grid[0, 4] == B_
+    assert grid[0, 5] is None
+    assert grid[1, 0] is None
+    assert grid[1, 1] == C
+    assert grid[1, 2] == D
+    assert grid[1, 3] is None
+    assert grid[1, 4] == C_
+    assert grid[1, 5] == D_
+    assert grid[2, 0] is None
+    assert grid[2, 1] is None
+    assert grid[2, 2] == E
+    assert grid[2, 3] is None
+    assert grid[2, 4] is None
+    assert grid[2, 5] == E_
+    assert grid.morphisms == morphisms
+
+    # Test the five lemma with object grouping and hints.
+    grid = DiagramGrid(d, {
+        FiniteSet(A, B, C, D, E): {"layout": "sequential",
+                                   "transpose": True},
+        FiniteSet(A_, B_, C_, D_, E_): {"layout": "sequential",
+                                        "transpose": True}},
+        transpose=True)
+
+    assert grid.width == 5
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] == C
+    assert grid[0, 3] == D
+    assert grid[0, 4] == E
+    assert grid[1, 0] == A_
+    assert grid[1, 1] == B_
+    assert grid[1, 2] == C_
+    assert grid[1, 3] == D_
+    assert grid[1, 4] == E_
+    assert grid.morphisms == morphisms
+
+    # A two-triangle disconnected diagram.
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    f_ = NamedMorphism(A_, B_, "f")
+    g_ = NamedMorphism(B_, C_, "g")
+    d = Diagram([f, g, f_, g_], {g * f: "unique", g_ * f_: "unique"})
+    grid = DiagramGrid(d)
+
+    assert grid.width == 4
+    assert grid.height == 2
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] == A_
+    assert grid[0, 3] == B_
+    assert grid[1, 0] == C
+    assert grid[1, 1] is None
+    assert grid[1, 2] == C_
+    assert grid[1, 3] is None
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet(), f_: FiniteSet(),
+                              g_: FiniteSet(), g * f: FiniteSet("unique"),
+                              g_ * f_: FiniteSet("unique")}
+
+    # A two-morphism disconnected diagram.
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(C, D, "g")
+    d = Diagram([f, g])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 4
+    assert grid.height == 1
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+    assert grid[0, 2] == C
+    assert grid[0, 3] == D
+    assert grid.morphisms == {f: FiniteSet(), g: FiniteSet()}
+
+    # Test a one-object diagram.
+    f = NamedMorphism(A, A, "f")
+    d = Diagram([f])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 1
+    assert grid.height == 1
+    assert grid[0, 0] == A
+
+    # Test a two-object disconnected diagram.
+    g = NamedMorphism(B, B, "g")
+    d = Diagram([f, g])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 2
+    assert grid.height == 1
+    assert grid[0, 0] == A
+    assert grid[0, 1] == B
+
+
+def test_DiagramGrid_pseudopod():
+    # Test a diagram in which even growing a pseudopod does not
+    # eventually help.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+    F = Object("F")
+    A_ = Object("A'")
+    B_ = Object("B'")
+    C_ = Object("C'")
+    D_ = Object("D'")
+    E_ = Object("E'")
+
+    f1 = NamedMorphism(A, B, "f1")
+    f2 = NamedMorphism(A, C, "f2")
+    f3 = NamedMorphism(A, D, "f3")
+    f4 = NamedMorphism(A, E, "f4")
+    f5 = NamedMorphism(A, A_, "f5")
+    f6 = NamedMorphism(A, B_, "f6")
+    f7 = NamedMorphism(A, C_, "f7")
+    f8 = NamedMorphism(A, D_, "f8")
+    f9 = NamedMorphism(A, E_, "f9")
+    f10 = NamedMorphism(A, F, "f10")
+    d = Diagram([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10])
+    grid = DiagramGrid(d)
+
+    assert grid.width == 5
+    assert grid.height == 3
+    assert grid[0, 0] == E
+    assert grid[0, 1] == C
+    assert grid[0, 2] == C_
+    assert grid[0, 3] == E_
+    assert grid[0, 4] == F
+    assert grid[1, 0] == D
+    assert grid[1, 1] == A
+    assert grid[1, 2] == A_
+    assert grid[1, 3] is None
+    assert grid[1, 4] is None
+    assert grid[2, 0] == D_
+    assert grid[2, 1] == B
+    assert grid[2, 2] == B_
+    assert grid[2, 3] is None
+    assert grid[2, 4] is None
+
+    morphisms = {}
+    for f in [f1, f2, f3, f4, f5, f6, f7, f8, f9, f10]:
+        morphisms[f] = FiniteSet()
+    assert grid.morphisms == morphisms
+
+
+def test_ArrowStringDescription():
+    astr = ArrowStringDescription("cm", "", None, "", "", "d", "r", "_", "f")
+    assert str(astr) == "\\ar[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "", 12, "", "", "d", "r", "_", "f")
+    assert str(astr) == "\\ar[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "^", 12, "", "", "d", "r", "_", "f")
+    assert str(astr) == "\\ar@/^12cm/[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "", 12, "r", "", "d", "r", "_", "f")
+    assert str(astr) == "\\ar[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "", 12, "r", "u", "d", "r", "_", "f")
+    assert str(astr) == "\\ar@(r,u)[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "", 12, "r", "u", "d", "r", "_", "f")
+    assert str(astr) == "\\ar@(r,u)[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "", 12, "r", "u", "d", "r", "_", "f")
+    astr.arrow_style = "{-->}"
+    assert str(astr) == "\\ar@(r,u)@{-->}[dr]_{f}"
+
+    astr = ArrowStringDescription("cm", "_", 12, "", "", "d", "r", "_", "f")
+    astr.arrow_style = "{-->}"
+    assert str(astr) == "\\ar@/_12cm/@{-->}[dr]_{f}"
+
+
+def test_XypicDiagramDrawer_line():
+    # A linear diagram.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(C, D, "h")
+    i = NamedMorphism(D, E, "i")
+    d = Diagram([f, g, h, i])
+    grid = DiagramGrid(d, layout="sequential")
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[r]^{f} & B \\ar[r]^{g} & C \\ar[r]^{h} & D \\ar[r]^{i} & E \n" \
+        "}\n"
+
+    # The same diagram, transposed.
+    grid = DiagramGrid(d, layout="sequential", transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[d]^{f} \\\\\n" \
+        "B \\ar[d]^{g} \\\\\n" \
+        "C \\ar[d]^{h} \\\\\n" \
+        "D \\ar[d]^{i} \\\\\n" \
+        "E \n" \
+        "}\n"
+
+
+def test_XypicDiagramDrawer_triangle():
+    # A triangle diagram.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+
+    d = Diagram([f, g], {g * f: "unique"})
+    grid = DiagramGrid(d)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[d]_{g\\circ f} \\ar[r]^{f} & B \\ar[ld]^{g} \\\\\n" \
+        "C & \n" \
+        "}\n"
+
+    # The same diagram, transposed.
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[r]^{g\\circ f} \\ar[d]_{f} & C \\\\\n" \
+        "B \\ar[ru]_{g} & \n" \
+        "}\n"
+
+    # The same diagram, with a masked morphism.
+    assert drawer.draw(d, grid, masked=[g]) == "\\xymatrix{\n" \
+        "A \\ar[r]^{g\\circ f} \\ar[d]_{f} & C \\\\\n" \
+        "B & \n" \
+        "}\n"
+
+    # The same diagram with a formatter for "unique".
+    def formatter(astr):
+        astr.label = "\\exists !" + astr.label
+        astr.arrow_style = "{-->}"
+
+    drawer.arrow_formatters["unique"] = formatter
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar@{-->}[r]^{\\exists !g\\circ f} \\ar[d]_{f} & C \\\\\n" \
+        "B \\ar[ru]_{g} & \n" \
+        "}\n"
+
+    # The same diagram with a default formatter.
+    def default_formatter(astr):
+        astr.label_displacement = "(0.45)"
+
+    drawer.default_arrow_formatter = default_formatter
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar@{-->}[r]^(0.45){\\exists !g\\circ f} \\ar[d]_(0.45){f} & C \\\\\n" \
+        "B \\ar[ru]_(0.45){g} & \n" \
+        "}\n"
+
+    # A triangle diagram with a lot of morphisms between the same
+    # objects.
+    f1 = NamedMorphism(B, A, "f1")
+    f2 = NamedMorphism(A, B, "f2")
+    g1 = NamedMorphism(C, B, "g1")
+    g2 = NamedMorphism(B, C, "g2")
+    d = Diagram([f, f1, f2, g, g1, g2], {f1 * g1: "unique", g2 * f2: "unique"})
+
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid, masked=[f1*g1*g2*f2, g2*f2*f1*g1]) == \
+        "\\xymatrix{\n" \
+        "A \\ar[r]^{g_{2}\\circ f_{2}} \\ar[d]_{f} \\ar@/^3mm/[d]^{f_{2}} " \
+        "& C \\ar@/^3mm/[l]^{f_{1}\\circ g_{1}} \\ar@/^3mm/[ld]^{g_{1}} \\\\\n" \
+        "B \\ar@/^3mm/[u]^{f_{1}} \\ar[ru]_{g} \\ar@/^3mm/[ru]^{g_{2}} & \n" \
+        "}\n"
+
+
+def test_XypicDiagramDrawer_cube():
+    # A cube diagram.
+    A1 = Object("A1")
+    A2 = Object("A2")
+    A3 = Object("A3")
+    A4 = Object("A4")
+    A5 = Object("A5")
+    A6 = Object("A6")
+    A7 = Object("A7")
+    A8 = Object("A8")
+
+    # The top face of the cube.
+    f1 = NamedMorphism(A1, A2, "f1")
+    f2 = NamedMorphism(A1, A3, "f2")
+    f3 = NamedMorphism(A2, A4, "f3")
+    f4 = NamedMorphism(A3, A4, "f3")
+
+    # The bottom face of the cube.
+    f5 = NamedMorphism(A5, A6, "f5")
+    f6 = NamedMorphism(A5, A7, "f6")
+    f7 = NamedMorphism(A6, A8, "f7")
+    f8 = NamedMorphism(A7, A8, "f8")
+
+    # The remaining morphisms.
+    f9 = NamedMorphism(A1, A5, "f9")
+    f10 = NamedMorphism(A2, A6, "f10")
+    f11 = NamedMorphism(A3, A7, "f11")
+    f12 = NamedMorphism(A4, A8, "f11")
+
+    d = Diagram([f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12])
+    grid = DiagramGrid(d)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "& A_{5} \\ar[r]^{f_{5}} \\ar[ldd]_{f_{6}} & A_{6} \\ar[rdd]^{f_{7}} " \
+        "& \\\\\n" \
+        "& A_{1} \\ar[r]^{f_{1}} \\ar[d]^{f_{2}} \\ar[u]^{f_{9}} & A_{2} " \
+        "\\ar[d]^{f_{3}} \\ar[u]_{f_{10}} & \\\\\n" \
+        "A_{7} \\ar@/_3mm/[rrr]_{f_{8}} & A_{3} \\ar[r]^{f_{3}} \\ar[l]_{f_{11}} " \
+        "& A_{4} \\ar[r]^{f_{11}} & A_{8} \n" \
+        "}\n"
+
+    # The same diagram, transposed.
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "& & A_{7} \\ar@/^3mm/[ddd]^{f_{8}} \\\\\n" \
+        "A_{5} \\ar[d]_{f_{5}} \\ar[rru]^{f_{6}} & A_{1} \\ar[d]^{f_{1}} " \
+        "\\ar[r]^{f_{2}} \\ar[l]^{f_{9}} & A_{3} \\ar[d]_{f_{3}} " \
+        "\\ar[u]^{f_{11}} \\\\\n" \
+        "A_{6} \\ar[rrd]_{f_{7}} & A_{2} \\ar[r]^{f_{3}} \\ar[l]^{f_{10}} " \
+        "& A_{4} \\ar[d]_{f_{11}} \\\\\n" \
+        "& & A_{8} \n" \
+        "}\n"
+
+
+def test_XypicDiagramDrawer_curved_and_loops():
+    # A simple diagram, with a curved arrow.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(D, A, "h")
+    k = NamedMorphism(D, B, "k")
+    d = Diagram([f, g, h, k])
+    grid = DiagramGrid(d)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[r]_{f} & B \\ar[d]^{g} & D \\ar[l]^{k} \\ar@/_3mm/[ll]_{h} \\\\\n" \
+        "& C & \n" \
+        "}\n"
+
+    # The same diagram, transposed.
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[d]^{f} & \\\\\n" \
+        "B \\ar[r]^{g} & C \\\\\n" \
+        "D \\ar[u]_{k} \\ar@/^3mm/[uu]^{h} & \n" \
+        "}\n"
+
+    # The same diagram, larger and rotated.
+    assert drawer.draw(d, grid, diagram_format="@+1cm@dr") == \
+        "\\xymatrix@+1cm@dr{\n" \
+        "A \\ar[d]^{f} & \\\\\n" \
+        "B \\ar[r]^{g} & C \\\\\n" \
+        "D \\ar[u]_{k} \\ar@/^3mm/[uu]^{h} & \n" \
+        "}\n"
+
+    # A simple diagram with three curved arrows.
+    h1 = NamedMorphism(D, A, "h1")
+    h2 = NamedMorphism(A, D, "h2")
+    k = NamedMorphism(D, B, "k")
+    d = Diagram([f, g, h, k, h1, h2])
+    grid = DiagramGrid(d)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[r]_{f} \\ar@/^3mm/[rr]^{h_{2}} & B \\ar[d]^{g} & D \\ar[l]^{k} " \
+        "\\ar@/_7mm/[ll]_{h} \\ar@/_11mm/[ll]_{h_{1}} \\\\\n" \
+        "& C & \n" \
+        "}\n"
+
+    # The same diagram, transposed.
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[d]^{f} \\ar@/_3mm/[dd]_{h_{2}} & \\\\\n" \
+        "B \\ar[r]^{g} & C \\\\\n" \
+        "D \\ar[u]_{k} \\ar@/^7mm/[uu]^{h} \\ar@/^11mm/[uu]^{h_{1}} & \n" \
+        "}\n"
+
+    # The same diagram, with "loop" morphisms.
+    l_A = NamedMorphism(A, A, "l_A")
+    l_D = NamedMorphism(D, D, "l_D")
+    l_C = NamedMorphism(C, C, "l_C")
+    d = Diagram([f, g, h, k, h1, h2, l_A, l_D, l_C])
+    grid = DiagramGrid(d)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[r]_{f} \\ar@/^3mm/[rr]^{h_{2}} \\ar@(u,l)[]^{l_{A}} " \
+        "& B \\ar[d]^{g} & D \\ar[l]^{k} \\ar@/_7mm/[ll]_{h} " \
+        "\\ar@/_11mm/[ll]_{h_{1}} \\ar@(r,u)[]^{l_{D}} \\\\\n" \
+        "& C \\ar@(l,d)[]^{l_{C}} & \n" \
+        "}\n"
+
+    # The same diagram with "loop" morphisms, transposed.
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[d]^{f} \\ar@/_3mm/[dd]_{h_{2}} \\ar@(r,u)[]^{l_{A}} & \\\\\n" \
+        "B \\ar[r]^{g} & C \\ar@(r,u)[]^{l_{C}} \\\\\n" \
+        "D \\ar[u]_{k} \\ar@/^7mm/[uu]^{h} \\ar@/^11mm/[uu]^{h_{1}} " \
+        "\\ar@(l,d)[]^{l_{D}} & \n" \
+        "}\n"
+
+    # The same diagram with two "loop" morphisms per object.
+    l_A_ = NamedMorphism(A, A, "n_A")
+    l_D_ = NamedMorphism(D, D, "n_D")
+    l_C_ = NamedMorphism(C, C, "n_C")
+    d = Diagram([f, g, h, k, h1, h2, l_A, l_D, l_C, l_A_, l_D_, l_C_])
+    grid = DiagramGrid(d)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[r]_{f} \\ar@/^3mm/[rr]^{h_{2}} \\ar@(u,l)[]^{l_{A}} " \
+        "\\ar@/^3mm/@(l,d)[]^{n_{A}} & B \\ar[d]^{g} & D \\ar[l]^{k} " \
+        "\\ar@/_7mm/[ll]_{h} \\ar@/_11mm/[ll]_{h_{1}} \\ar@(r,u)[]^{l_{D}} " \
+        "\\ar@/^3mm/@(d,r)[]^{n_{D}} \\\\\n" \
+        "& C \\ar@(l,d)[]^{l_{C}} \\ar@/^3mm/@(d,r)[]^{n_{C}} & \n" \
+        "}\n"
+
+    # The same diagram with two "loop" morphisms per object, transposed.
+    grid = DiagramGrid(d, transpose=True)
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == "\\xymatrix{\n" \
+        "A \\ar[d]^{f} \\ar@/_3mm/[dd]_{h_{2}} \\ar@(r,u)[]^{l_{A}} " \
+        "\\ar@/^3mm/@(u,l)[]^{n_{A}} & \\\\\n" \
+        "B \\ar[r]^{g} & C \\ar@(r,u)[]^{l_{C}} \\ar@/^3mm/@(d,r)[]^{n_{C}} \\\\\n" \
+        "D \\ar[u]_{k} \\ar@/^7mm/[uu]^{h} \\ar@/^11mm/[uu]^{h_{1}} " \
+        "\\ar@(l,d)[]^{l_{D}} \\ar@/^3mm/@(d,r)[]^{n_{D}} & \n" \
+        "}\n"
+
+
+def test_xypic_draw_diagram():
+    # A linear diagram.
+    A = Object("A")
+    B = Object("B")
+    C = Object("C")
+    D = Object("D")
+    E = Object("E")
+
+    f = NamedMorphism(A, B, "f")
+    g = NamedMorphism(B, C, "g")
+    h = NamedMorphism(C, D, "h")
+    i = NamedMorphism(D, E, "i")
+    d = Diagram([f, g, h, i])
+
+    grid = DiagramGrid(d, layout="sequential")
+    drawer = XypicDiagramDrawer()
+    assert drawer.draw(d, grid) == xypic_draw_diagram(d, layout="sequential")

lib/python3.12/site-packages/sympy/liealgebras/__init__.py

@@ -0,0 +1,3 @@
+from sympy.liealgebras.cartan_type import CartanType
+
+__all__ = ['CartanType']

lib/python3.12/site-packages/sympy/liealgebras/cartan_matrix.py

@@ -0,0 +1,25 @@
+from .cartan_type import CartanType
+
+def CartanMatrix(ct):
+    """Access the Cartan matrix of a specific Lie algebra
+
+    Examples
+    ========
+
+    >>> from sympy.liealgebras.cartan_matrix import CartanMatrix
+    >>> CartanMatrix("A2")
+    Matrix([
+    [ 2, -1],
+    [-1,  2]])
+
+    >>> CartanMatrix(['C', 3])
+    Matrix([
+    [ 2, -1,  0],
+    [-1,  2, -1],
+    [ 0, -2,  2]])
+
+    This method works by returning the Cartan matrix
+    which corresponds to Cartan type t.
+    """
+
+    return CartanType(ct).cartan_matrix()

lib/python3.12/site-packages/sympy/liealgebras/cartan_type.py

@@ -0,0 +1,73 @@
+from sympy.core import Atom, Basic
+
+
+class CartanType_generator():
+    """
+    Constructor for actually creating things
+    """
+    def __call__(self, *args):
+        c = args[0]
+        if isinstance(c, list):
+            letter, n = c[0], int(c[1])
+        elif isinstance(c, str):
+            letter, n = c[0], int(c[1:])
+        else:
+            raise TypeError("Argument must be a string (e.g. 'A3') or a list (e.g. ['A', 3])")
+
+        if n < 0:
+            raise ValueError("Lie algebra rank cannot be negative")
+        if letter == "A":
+            from . import type_a
+            return type_a.TypeA(n)
+        if letter == "B":
+            from . import type_b
+            return type_b.TypeB(n)
+
+        if letter == "C":
+            from . import type_c
+            return type_c.TypeC(n)
+
+        if letter == "D":
+            from . import type_d
+            return type_d.TypeD(n)
+
+        if letter == "E":
+            if n >= 6 and n <= 8:
+                from . import type_e
+                return type_e.TypeE(n)
+
+        if letter == "F":
+            if n == 4:
+                from . import type_f
+                return type_f.TypeF(n)
+
+        if letter == "G":
+            if n == 2:
+                from . import type_g
+                return type_g.TypeG(n)
+
+CartanType = CartanType_generator()
+
+
+class Standard_Cartan(Atom):
+    """
+    Concrete base class for Cartan types such as A4, etc
+    """
+
+    def __new__(cls, series, n):
+        obj = Basic.__new__(cls)
+        obj.n = n
+        obj.series = series
+        return obj
+
+    def rank(self):
+        """
+        Returns the rank of the Lie algebra
+        """
+        return self.n
+
+    def series(self):
+        """
+        Returns the type of the Lie algebra
+        """
+        return self.series

lib/python3.12/site-packages/sympy/liealgebras/dynkin_diagram.py

@@ -0,0 +1,24 @@
+from .cartan_type import CartanType
+
+
+def DynkinDiagram(t):
+    """Display the Dynkin diagram of a given Lie algebra
+
+    Works by generating the CartanType for the input, t, and then returning the
+    Dynkin diagram method from the individual classes.
+
+    Examples
+    ========
+
+    >>> from sympy.liealgebras.dynkin_diagram import DynkinDiagram
+    >>> print(DynkinDiagram("A3"))
+    0---0---0
+    1   2   3
+
+    >>> print(DynkinDiagram("B4"))
+    0---0---0=>=0
+    1   2   3   4
+
+    """
+
+    return CartanType(t).dynkin_diagram()

lib/python3.12/site-packages/sympy/liealgebras/root_system.py

@@ -0,0 +1,199 @@
+from .cartan_type import CartanType
+from sympy.core.basic import Atom
+
+class RootSystem(Atom):
+    """Represent the root system of a simple Lie algebra
+
+    Every simple Lie algebra has a unique root system.  To find the root
+    system, we first consider the Cartan subalgebra of g, which is the maximal
+    abelian subalgebra, and consider the adjoint action of g on this
+    subalgebra.  There is a root system associated with this action. Now, a
+    root system over a vector space V is a set of finite vectors Phi (called
+    roots), which satisfy:
+
+    1.  The roots span V
+    2.  The only scalar multiples of x in Phi are x and -x
+    3.  For every x in Phi, the set Phi is closed under reflection
+        through the hyperplane perpendicular to x.
+    4.  If x and y are roots in Phi, then the projection of y onto
+        the line through x is a half-integral multiple of x.
+
+    Now, there is a subset of Phi, which we will call Delta, such that:
+    1.  Delta is a basis of V
+    2.  Each root x in Phi can be written x = sum k_y y for y in Delta
+
+    The elements of Delta are called the simple roots.
+    Therefore, we see that the simple roots span the root space of a given
+    simple Lie algebra.
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Root_system
+    .. [2] Lie Algebras and Representation Theory - Humphreys
+
+    """
+
+    def __new__(cls, cartantype):
+        """Create a new RootSystem object
+
+        This method assigns an attribute called cartan_type to each instance of
+        a RootSystem object.  When an instance of RootSystem is called, it
+        needs an argument, which should be an instance of a simple Lie algebra.
+        We then take the CartanType of this argument and set it as the
+        cartan_type attribute of the RootSystem instance.
+
+        """
+        obj = Atom.__new__(cls)
+        obj.cartan_type = CartanType(cartantype)
+        return obj
+
+    def simple_roots(self):
+        """Generate the simple roots of the Lie algebra
+
+        The rank of the Lie algebra determines the number of simple roots that
+        it has.  This method obtains the rank of the Lie algebra, and then uses
+        the simple_root method from the Lie algebra classes to generate all the
+        simple roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.root_system import RootSystem
+        >>> c = RootSystem("A3")
+        >>> roots = c.simple_roots()
+        >>> roots
+        {1: [1, -1, 0, 0], 2: [0, 1, -1, 0], 3: [0, 0, 1, -1]}
+
+        """
+        n = self.cartan_type.rank()
+        roots = {}
+        for i in range(1, n+1):
+            root = self.cartan_type.simple_root(i)
+            roots[i] = root
+        return roots
+
+
+    def all_roots(self):
+        """Generate all the roots of a given root system
+
+        The result is a dictionary where the keys are integer numbers.  It
+        generates the roots by getting the dictionary of all positive roots
+        from the bases classes, and then taking each root, and multiplying it
+        by -1 and adding it to the dictionary.  In this way all the negative
+        roots are generated.
+
+        """
+        alpha = self.cartan_type.positive_roots()
+        keys = list(alpha.keys())
+        k = max(keys)
+        for val in keys:
+            k += 1
+            root = alpha[val]
+            newroot = [-x for x in root]
+            alpha[k] = newroot
+        return alpha
+
+    def root_space(self):
+        """Return the span of the simple roots
+
+        The root space is the vector space spanned by the simple roots, i.e. it
+        is a vector space with a distinguished basis, the simple roots.  This
+        method returns a string that represents the root space as the span of
+        the simple roots, alpha[1],...., alpha[n].
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.root_system import RootSystem
+        >>> c = RootSystem("A3")
+        >>> c.root_space()
+        'alpha[1] + alpha[2] + alpha[3]'
+
+        """
+        n = self.cartan_type.rank()
+        rs = " + ".join("alpha["+str(i) +"]" for i in range(1, n+1))
+        return rs
+
+    def add_simple_roots(self, root1, root2):
+        """Add two simple roots together
+
+        The function takes as input two integers, root1 and root2.  It then
+        uses these integers as keys in the dictionary of simple roots, and gets
+        the corresponding simple roots, and then adds them together.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.root_system import RootSystem
+        >>> c = RootSystem("A3")
+        >>> newroot = c.add_simple_roots(1, 2)
+        >>> newroot
+        [1, 0, -1, 0]
+
+        """
+
+        alpha = self.simple_roots()
+        if root1 > len(alpha) or root2 > len(alpha):
+            raise ValueError("You've used a root that doesn't exist!")
+        a1 = alpha[root1]
+        a2 = alpha[root2]
+        newroot = [_a1 + _a2 for _a1, _a2 in zip(a1, a2)]
+        return newroot
+
+    def add_as_roots(self, root1, root2):
+        """Add two roots together if and only if their sum is also a root
+
+        It takes as input two vectors which should be roots.  It then computes
+        their sum and checks if it is in the list of all possible roots.  If it
+        is, it returns the sum.  Otherwise it returns a string saying that the
+        sum is not a root.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.root_system import RootSystem
+        >>> c = RootSystem("A3")
+        >>> c.add_as_roots([1, 0, -1, 0], [0, 0, 1, -1])
+        [1, 0, 0, -1]
+        >>> c.add_as_roots([1, -1, 0, 0], [0, 0, -1, 1])
+        'The sum of these two roots is not a root'
+
+        """
+        alpha = self.all_roots()
+        newroot = [r1 + r2 for r1, r2 in zip(root1, root2)]
+        if newroot in alpha.values():
+            return newroot
+        else:
+            return "The sum of these two roots is not a root"
+
+
+    def cartan_matrix(self):
+        """Cartan matrix of Lie algebra associated with this root system
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.root_system import RootSystem
+        >>> c = RootSystem("A3")
+        >>> c.cartan_matrix()
+        Matrix([
+            [ 2, -1,  0],
+            [-1,  2, -1],
+            [ 0, -1,  2]])
+        """
+        return self.cartan_type.cartan_matrix()
+
+    def dynkin_diagram(self):
+        """Dynkin diagram of the Lie algebra associated with this root system
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.root_system import RootSystem
+        >>> c = RootSystem("A3")
+        >>> print(c.dynkin_diagram())
+        0---0---0
+        1   2   3
+        """
+        return self.cartan_type.dynkin_diagram()

lib/python3.12/site-packages/sympy/liealgebras/tests/test_cartan_matrix.py

@@ -0,0 +1,10 @@
+from sympy.liealgebras.cartan_matrix import CartanMatrix
+from sympy.matrices import Matrix
+
+def test_CartanMatrix():
+    c = CartanMatrix("A3")
+    m = Matrix(3, 3, [2, -1, 0, -1, 2, -1, 0, -1, 2])
+    assert c == m
+    a = CartanMatrix(["G",2])
+    mt = Matrix(2, 2, [2, -1, -3, 2])
+    assert a == mt

lib/python3.12/site-packages/sympy/liealgebras/tests/test_cartan_type.py

@@ -0,0 +1,12 @@
+from sympy.liealgebras.cartan_type import CartanType, Standard_Cartan
+
+def test_Standard_Cartan():
+    c = CartanType("A4")
+    assert c.rank() == 4
+    assert c.series == "A"
+    m = Standard_Cartan("A", 2)
+    assert m.rank() == 2
+    assert m.series == "A"
+    b = CartanType("B12")
+    assert b.rank() == 12
+    assert b.series == "B"

lib/python3.12/site-packages/sympy/liealgebras/tests/test_dynkin_diagram.py

@@ -0,0 +1,9 @@
+from sympy.liealgebras.dynkin_diagram import DynkinDiagram
+
+def test_DynkinDiagram():
+    c = DynkinDiagram("A3")
+    diag = "0---0---0\n1   2   3"
+    assert c == diag
+    ct = DynkinDiagram(["B", 3])
+    diag2 = "0---0=>=0\n1   2   3"
+    assert ct == diag2

lib/python3.12/site-packages/sympy/liealgebras/tests/test_root_system.py

@@ -0,0 +1,18 @@
+from sympy.liealgebras.root_system import RootSystem
+from sympy.liealgebras.type_a import TypeA
+from sympy.matrices import Matrix
+
+def test_root_system():
+    c = RootSystem("A3")
+    assert c.cartan_type == TypeA(3)
+    assert c.simple_roots() ==  {1: [1, -1, 0, 0], 2: [0, 1, -1, 0], 3: [0, 0, 1, -1]}
+    assert c.root_space() == "alpha[1] + alpha[2] + alpha[3]"
+    assert c.cartan_matrix() == Matrix([[ 2, -1,  0], [-1,  2, -1], [ 0, -1,  2]])
+    assert c.dynkin_diagram() == "0---0---0\n1   2   3"
+    assert c.add_simple_roots(1, 2) == [1, 0, -1, 0]
+    assert c.all_roots() == {1: [1, -1, 0, 0], 2: [1, 0, -1, 0],
+            3: [1, 0, 0, -1], 4: [0, 1, -1, 0], 5: [0, 1, 0, -1],
+            6: [0, 0, 1, -1], 7: [-1, 1, 0, 0], 8: [-1, 0, 1, 0],
+            9: [-1, 0, 0, 1], 10: [0, -1, 1, 0],
+            11: [0, -1, 0, 1], 12: [0, 0, -1, 1]}
+    assert c.add_as_roots([1, 0, -1, 0], [0, 0, 1, -1]) == [1, 0, 0, -1]

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_A.py

@@ -0,0 +1,17 @@
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+
+def test_type_A():
+    c = CartanType("A3")
+    m = Matrix(3, 3, [2, -1, 0, -1, 2, -1, 0, -1, 2])
+    assert m == c.cartan_matrix()
+    assert c.basis() == 8
+    assert c.roots() == 12
+    assert c.dimension() == 4
+    assert c.simple_root(1) == [1, -1, 0, 0]
+    assert c.highest_root() == [1, 0, 0, -1]
+    assert c.lie_algebra() == "su(4)"
+    diag = "0---0---0\n1   2   3"
+    assert c.dynkin_diagram() == diag
+    assert c.positive_roots() == {1: [1, -1, 0, 0], 2: [1, 0, -1, 0],
+            3: [1, 0, 0, -1], 4: [0, 1, -1, 0], 5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_B.py

@@ -0,0 +1,17 @@
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+
+def test_type_B():
+    c = CartanType("B3")
+    m = Matrix(3, 3, [2, -1, 0, -1, 2, -2, 0, -1, 2])
+    assert m == c.cartan_matrix()
+    assert c.dimension() == 3
+    assert c.roots() == 18
+    assert c.simple_root(3) == [0, 0, 1]
+    assert c.basis() == 3
+    assert c.lie_algebra() == "so(6)"
+    diag = "0---0=>=0\n1   2   3"
+    assert c.dynkin_diagram() == diag
+    assert c.positive_roots() ==  {1: [1, -1, 0], 2: [1, 1, 0], 3: [1, 0, -1],
+            4: [1, 0, 1], 5: [0, 1, -1], 6: [0, 1, 1], 7: [1, 0, 0],
+            8: [0, 1, 0], 9: [0, 0, 1]}

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_C.py

@@ -0,0 +1,22 @@
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+
+def test_type_C():
+    c = CartanType("C4")
+    m = Matrix(4, 4, [2, -1, 0, 0, -1, 2, -1, 0, 0, -1, 2, -1, 0, 0, -2, 2])
+    assert c.cartan_matrix() == m
+    assert c.dimension() == 4
+    assert c.simple_root(4) == [0, 0, 0, 2]
+    assert c.roots() == 32
+    assert c.basis() == 36
+    assert c.lie_algebra() == "sp(8)"
+    t = CartanType(['C', 3])
+    assert t.dimension() == 3
+    diag = "0---0---0=<=0\n1   2   3   4"
+    assert c.dynkin_diagram() == diag
+    assert c.positive_roots() == {1: [1, -1, 0, 0], 2: [1, 1, 0, 0],
+            3: [1, 0, -1, 0], 4: [1, 0, 1, 0], 5: [1, 0, 0, -1],
+            6: [1, 0, 0, 1], 7: [0, 1, -1, 0], 8: [0, 1, 1, 0],
+            9: [0, 1, 0, -1], 10: [0, 1, 0, 1], 11: [0, 0, 1, -1],
+            12: [0, 0, 1, 1], 13: [2, 0, 0, 0], 14: [0, 2, 0, 0], 15: [0, 0, 2, 0],
+            16: [0, 0, 0, 2]}

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_D.py

@@ -0,0 +1,19 @@
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+
+
+
+def test_type_D():
+    c = CartanType("D4")
+    m = Matrix(4, 4, [2, -1, 0, 0, -1, 2, -1, -1, 0, -1, 2, 0, 0, -1, 0, 2])
+    assert c.cartan_matrix() == m
+    assert c.basis() == 6
+    assert c.lie_algebra() == "so(8)"
+    assert c.roots() == 24
+    assert c.simple_root(3) == [0, 0, 1, -1]
+    diag = "    3\n    0\n    |\n    |\n0---0---0\n1   2   4"
+    assert diag == c.dynkin_diagram()
+    assert c.positive_roots() == {1: [1, -1, 0, 0], 2: [1, 1, 0, 0],
+            3: [1, 0, -1, 0], 4: [1, 0, 1, 0], 5: [1, 0, 0, -1], 6: [1, 0, 0, 1],
+            7: [0, 1, -1, 0], 8: [0, 1, 1, 0], 9: [0, 1, 0, -1], 10: [0, 1, 0, 1],
+            11: [0, 0, 1, -1], 12: [0, 0, 1, 1]}

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_E.py

@@ -0,0 +1,19 @@
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+
+def test_type_E():
+    c = CartanType("E6")
+    m = Matrix(6, 6, [2, 0, -1, 0, 0, 0, 0, 2, 0, -1, 0, 0,
+        -1, 0, 2, -1, 0, 0, 0, -1, -1, 2, -1, 0, 0, 0, 0,
+        -1, 2, -1, 0, 0, 0, 0, -1, 2])
+    assert c.cartan_matrix() == m
+    assert c.dimension() == 8
+    assert c.simple_root(6) == [0, 0, 0, -1, 1, 0, 0, 0]
+    assert c.roots() == 72
+    assert c.basis() == 78
+    diag = " "*8 + "2\n" + " "*8 + "0\n" + " "*8 + "|\n" + " "*8 + "|\n"
+    diag += "---".join("0" for i in range(1, 6))+"\n"
+    diag += "1   " + "   ".join(str(i) for i in range(3, 7))
+    assert c.dynkin_diagram() == diag
+    posroots = c.positive_roots()
+    assert posroots[8] == [1, 0, 0, 0, 1, 0, 0, 0]

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_F.py

@@ -0,0 +1,24 @@
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+from sympy.core.backend import S
+
+def test_type_F():
+    c = CartanType("F4")
+    m = Matrix(4, 4, [2, -1, 0, 0, -1, 2, -2, 0, 0, -1, 2, -1, 0, 0, -1, 2])
+    assert c.cartan_matrix() == m
+    assert c.dimension() == 4
+    assert c.simple_root(1) == [1, -1, 0, 0]
+    assert c.simple_root(2) == [0, 1, -1, 0]
+    assert c.simple_root(3) == [0, 0, 0, 1]
+    assert c.simple_root(4) == [-S.Half, -S.Half, -S.Half, -S.Half]
+    assert c.roots() == 48
+    assert c.basis() == 52
+    diag = "0---0=>=0---0\n" + "   ".join(str(i) for i in range(1, 5))
+    assert c.dynkin_diagram() == diag
+    assert c.positive_roots() == {1: [1, -1, 0, 0], 2: [1, 1, 0, 0], 3: [1, 0, -1, 0],
+            4: [1, 0, 1, 0], 5: [1, 0, 0, -1], 6: [1, 0, 0, 1], 7: [0, 1, -1, 0],
+            8: [0, 1, 1, 0], 9: [0, 1, 0, -1], 10: [0, 1, 0, 1], 11: [0, 0, 1, -1],
+            12: [0, 0, 1, 1], 13: [1, 0, 0, 0], 14: [0, 1, 0, 0], 15: [0, 0, 1, 0],
+            16: [0, 0, 0, 1], 17: [S.Half, S.Half, S.Half, S.Half], 18: [S.Half, -S.Half, S.Half, S.Half],
+            19: [S.Half, S.Half, -S.Half, S.Half], 20: [S.Half, S.Half, S.Half, -S.Half], 21: [S.Half, S.Half, -S.Half, -S.Half],
+            22: [S.Half, -S.Half, S.Half, -S.Half], 23: [S.Half, -S.Half, -S.Half, S.Half], 24: [S.Half, -S.Half, -S.Half, -S.Half]}

lib/python3.12/site-packages/sympy/liealgebras/tests/test_type_G.py

@@ -0,0 +1,16 @@
+# coding=utf-8
+from sympy.liealgebras.cartan_type import CartanType
+from sympy.matrices import Matrix
+
+def test_type_G():
+    c = CartanType("G2")
+    m = Matrix(2, 2, [2, -1, -3, 2])
+    assert c.cartan_matrix() == m
+    assert c.simple_root(2) == [1, -2, 1]
+    assert c.basis() == 14
+    assert c.roots() == 12
+    assert c.dimension() == 3
+    diag = "0≡<≡0\n1   2"
+    assert diag == c.dynkin_diagram()
+    assert c.positive_roots() == {1: [0, 1, -1], 2: [1, -2, 1], 3: [1, -1, 0],
+            4: [1, 0, 1], 5: [1, 1, -2], 6: [2, -1, -1]}

lib/python3.12/site-packages/sympy/liealgebras/tests/test_weyl_group.py

@@ -0,0 +1,35 @@
+from sympy.liealgebras.weyl_group import WeylGroup
+from sympy.matrices import Matrix
+
+def test_weyl_group():
+    c = WeylGroup("A3")
+    assert c.matrix_form('r1*r2') == Matrix([[0, 0, 1, 0], [1, 0, 0, 0],
+        [0, 1, 0, 0], [0, 0, 0, 1]])
+    assert c.generators() == ['r1', 'r2', 'r3']
+    assert c.group_order() == 24.0
+    assert c.group_name() == "S4: the symmetric group acting on 4 elements."
+    assert c.coxeter_diagram() == "0---0---0\n1   2   3"
+    assert c.element_order('r1*r2*r3') == 4
+    assert c.element_order('r1*r3*r2*r3') == 3
+    d = WeylGroup("B5")
+    assert d.group_order() == 3840
+    assert d.element_order('r1*r2*r4*r5') == 12
+    assert d.matrix_form('r2*r3') ==  Matrix([[0, 0, 1, 0, 0], [1, 0, 0, 0, 0],
+        [0, 1, 0, 0, 0], [0, 0, 0, 1, 0], [0, 0, 0, 0, 1]])
+    assert d.element_order('r1*r2*r1*r3*r5') == 6
+    e = WeylGroup("D5")
+    assert e.element_order('r2*r3*r5') == 4
+    assert e.matrix_form('r2*r3*r5') == Matrix([[1, 0, 0, 0, 0], [0, 0, 0, 0, -1],
+        [0, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, -1, 0]])
+    f = WeylGroup("G2")
+    assert f.element_order('r1*r2*r1*r2') == 3
+    assert f.element_order('r2*r1*r1*r2') == 1
+
+    assert f.matrix_form('r1*r2*r1*r2') == Matrix([[0, 1, 0], [0, 0, 1], [1, 0, 0]])
+    g = WeylGroup("F4")
+    assert g.matrix_form('r2*r3') == Matrix([[1, 0, 0, 0], [0, 1, 0, 0],
+        [0, 0, 0, -1], [0, 0, 1, 0]])
+
+    assert g.element_order('r2*r3') == 4
+    h = WeylGroup("E6")
+    assert h.group_order() == 51840

lib/python3.12/site-packages/sympy/liealgebras/type_a.py

@@ -0,0 +1,166 @@
+from sympy.liealgebras.cartan_type import Standard_Cartan
+from sympy.core.backend import eye
+
+
+class TypeA(Standard_Cartan):
+    """
+    This class contains the information about
+    the A series of simple Lie algebras.
+    ====
+    """
+
+    def __new__(cls, n):
+        if n < 1:
+            raise ValueError("n cannot be less than 1")
+        return Standard_Cartan.__new__(cls, "A", n)
+
+
+    def dimension(self):
+        """Dimension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A4")
+        >>> c.dimension()
+        5
+        """
+        return self.n+1
+
+
+    def basic_root(self, i, j):
+        """
+        This is a method just to generate roots
+        with a 1 iin the ith position and a -1
+        in the jth position.
+
+        """
+
+        n = self.n
+        root = [0]*(n+1)
+        root[i] = 1
+        root[j] = -1
+        return root
+
+    def simple_root(self, i):
+        """
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        In A_n the ith simple root is the root which has a 1
+        in the ith position, a -1 in the (i+1)th position,
+        and zeroes elsewhere.
+
+        This method returns the ith simple root for the A series.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A4")
+        >>> c.simple_root(1)
+        [1, -1, 0, 0, 0]
+
+        """
+
+        return self.basic_root(i-1, i)
+
+    def positive_roots(self):
+        """
+        This method generates all the positive roots of
+        A_n.  This is half of all of the roots of A_n;
+        by multiplying all the positive roots by -1 we
+        get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+        """
+
+        n = self.n
+        posroots = {}
+        k = 0
+        for i in range(0, n):
+            for j in range(i+1, n+1):
+               k += 1
+               posroots[k] = self.basic_root(i, j)
+        return posroots
+
+    def highest_root(self):
+        """
+        Returns the highest weight root for A_n
+        """
+
+        return self.basic_root(0, self.n)
+
+    def roots(self):
+        """
+        Returns the total number of roots for A_n
+        """
+        n = self.n
+        return n*(n+1)
+
+    def cartan_matrix(self):
+        """
+        Returns the Cartan matrix for A_n.
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType('A4')
+        >>> c.cartan_matrix()
+        Matrix([
+        [ 2, -1,  0,  0],
+        [-1,  2, -1,  0],
+        [ 0, -1,  2, -1],
+        [ 0,  0, -1,  2]])
+
+        """
+
+        n = self.n
+        m = 2 * eye(n)
+        i = 1
+        while i < n-1:
+            m[i, i+1] = -1
+            m[i, i-1] = -1
+            i += 1
+        m[0,1] = -1
+        m[n-1, n-2] = -1
+        return m
+
+    def basis(self):
+        """
+        Returns the number of independent generators of A_n
+        """
+        n = self.n
+        return n**2 - 1
+
+    def lie_algebra(self):
+        """
+        Returns the Lie algebra associated with A_n
+        """
+        n = self.n
+        return "su(" + str(n + 1) + ")"
+
+    def dynkin_diagram(self):
+        n = self.n
+        diag = "---".join("0" for i in range(1, n+1)) + "\n"
+        diag += "   ".join(str(i) for i in range(1, n+1))
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/type_b.py

@@ -0,0 +1,172 @@
+from .cartan_type import Standard_Cartan
+from sympy.core.backend import eye
+
+class TypeB(Standard_Cartan):
+
+    def __new__(cls, n):
+        if n < 2:
+            raise ValueError("n cannot be less than 2")
+        return Standard_Cartan.__new__(cls, "B", n)
+
+    def dimension(self):
+        """Dimension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("B3")
+        >>> c.dimension()
+        3
+        """
+
+        return self.n
+
+    def basic_root(self, i, j):
+        """
+        This is a method just to generate roots
+        with a 1 iin the ith position and a -1
+        in the jth position.
+
+        """
+        root = [0]*self.n
+        root[i] = 1
+        root[j] = -1
+        return root
+
+    def simple_root(self, i):
+        """
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        In B_n the first n-1 simple roots are the same as the
+        roots in A_(n-1) (a 1 in the ith position, a -1 in
+        the (i+1)th position, and zeroes elsewhere).  The n-th
+        simple root is the root with a 1 in the nth position
+        and zeroes elsewhere.
+
+        This method returns the ith simple root for the B series.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("B3")
+        >>> c.simple_root(2)
+        [0, 1, -1]
+
+        """
+        n = self.n
+        if i < n:
+            return self.basic_root(i-1, i)
+        else:
+            root = [0]*self.n
+            root[n-1] = 1
+            return root
+
+    def positive_roots(self):
+        """
+        This method generates all the positive roots of
+        A_n.  This is half of all of the roots of B_n;
+        by multiplying all the positive roots by -1 we
+        get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+        """
+
+        n = self.n
+        posroots = {}
+        k = 0
+        for i in range(0, n-1):
+            for j in range(i+1, n):
+               k += 1
+               posroots[k] = self.basic_root(i, j)
+               k += 1
+               root = self.basic_root(i, j)
+               root[j] = 1
+               posroots[k] = root
+
+        for i in range(0, n):
+            k += 1
+            root = [0]*n
+            root[i] = 1
+            posroots[k] = root
+
+        return posroots
+
+    def roots(self):
+        """
+        Returns the total number of roots for B_n"
+        """
+
+        n = self.n
+        return 2*(n**2)
+
+    def cartan_matrix(self):
+        """
+        Returns the Cartan matrix for B_n.
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType('B4')
+        >>> c.cartan_matrix()
+        Matrix([
+        [ 2, -1,  0,  0],
+        [-1,  2, -1,  0],
+        [ 0, -1,  2, -2],
+        [ 0,  0, -1,  2]])
+
+        """
+
+        n = self.n
+        m = 2* eye(n)
+        i = 1
+        while i < n-1:
+            m[i, i+1] = -1
+            m[i, i-1] = -1
+            i += 1
+        m[0, 1] = -1
+        m[n-2, n-1] = -2
+        m[n-1, n-2] = -1
+        return m
+
+    def basis(self):
+        """
+        Returns the number of independent generators of B_n
+        """
+
+        n = self.n
+        return (n**2 - n)/2
+
+    def lie_algebra(self):
+        """
+        Returns the Lie algebra associated with B_n
+        """
+
+        n = self.n
+        return "so(" + str(2*n) + ")"
+
+    def dynkin_diagram(self):
+        n = self.n
+        diag = "---".join("0" for i in range(1, n)) + "=>=0\n"
+        diag += "   ".join(str(i) for i in range(1, n+1))
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/type_c.py

@@ -0,0 +1,171 @@
+from .cartan_type import Standard_Cartan
+from sympy.core.backend import eye
+
+class TypeC(Standard_Cartan):
+
+    def __new__(cls, n):
+        if n < 3:
+            raise ValueError("n cannot be less than 3")
+        return Standard_Cartan.__new__(cls, "C", n)
+
+
+    def dimension(self):
+        """Dimension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("C3")
+        >>> c.dimension()
+        3
+        """
+        n = self.n
+        return n
+
+    def basic_root(self, i, j):
+        """Generate roots with 1 in ith position and a -1 in jth position
+        """
+        n = self.n
+        root = [0]*n
+        root[i] = 1
+        root[j] = -1
+        return root
+
+    def simple_root(self, i):
+        """The ith simple root for the C series
+
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        In C_n, the first n-1 simple roots are the same as
+        the roots in A_(n-1) (a 1 in the ith position, a -1
+        in the (i+1)th position, and zeroes elsewhere).  The
+        nth simple root is the root in which there is a 2 in
+        the nth position and zeroes elsewhere.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("C3")
+        >>> c.simple_root(2)
+        [0, 1, -1]
+
+        """
+
+        n = self.n
+        if i < n:
+            return self.basic_root(i-1,i)
+        else:
+            root = [0]*self.n
+            root[n-1] = 2
+            return root
+
+
+    def positive_roots(self):
+        """Generates all the positive roots of A_n
+
+        This is half of all of the roots of C_n; by multiplying all the
+        positive roots by -1 we get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+
+        """
+
+        n = self.n
+        posroots = {}
+        k = 0
+        for i in range(0, n-1):
+            for j in range(i+1, n):
+               k += 1
+               posroots[k] = self.basic_root(i, j)
+               k += 1
+               root = self.basic_root(i, j)
+               root[j] = 1
+               posroots[k] = root
+
+        for i in range(0, n):
+            k += 1
+            root = [0]*n
+            root[i] = 2
+            posroots[k] = root
+
+        return posroots
+
+    def roots(self):
+        """
+        Returns the total number of roots for C_n"
+        """
+
+        n = self.n
+        return 2*(n**2)
+
+    def cartan_matrix(self):
+        """The Cartan matrix for C_n
+
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType('C4')
+        >>> c.cartan_matrix()
+        Matrix([
+        [ 2, -1,  0,  0],
+        [-1,  2, -1,  0],
+        [ 0, -1,  2, -1],
+        [ 0,  0, -2,  2]])
+
+        """
+
+        n = self.n
+        m = 2 * eye(n)
+        i = 1
+        while i < n-1:
+           m[i, i+1] = -1
+           m[i, i-1] = -1
+           i += 1
+        m[0,1] = -1
+        m[n-1, n-2] = -2
+        return m
+
+
+    def basis(self):
+        """
+        Returns the number of independent generators of C_n
+        """
+
+        n = self.n
+        return n*(2*n + 1)
+
+    def lie_algebra(self):
+        """
+        Returns the Lie algebra associated with C_n"
+        """
+
+        n = self.n
+        return "sp(" + str(2*n) + ")"
+
+    def dynkin_diagram(self):
+        n = self.n
+        diag = "---".join("0" for i in range(1, n)) + "=<=0\n"
+        diag += "   ".join(str(i) for i in range(1, n+1))
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/type_d.py

@@ -0,0 +1,175 @@
+from .cartan_type import Standard_Cartan
+from sympy.core.backend import eye
+
+class TypeD(Standard_Cartan):
+
+    def __new__(cls, n):
+        if n < 3:
+            raise ValueError("n cannot be less than 3")
+        return Standard_Cartan.__new__(cls, "D", n)
+
+
+    def dimension(self):
+        """Dmension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("D4")
+        >>> c.dimension()
+        4
+        """
+
+        return self.n
+
+    def basic_root(self, i, j):
+        """
+        This is a method just to generate roots
+        with a 1 iin the ith position and a -1
+        in the jth position.
+
+        """
+
+        n = self.n
+        root = [0]*n
+        root[i] = 1
+        root[j] = -1
+        return root
+
+    def simple_root(self, i):
+        """
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        In D_n, the first n-1 simple roots are the same as
+        the roots in A_(n-1) (a 1 in the ith position, a -1
+        in the (i+1)th position, and zeroes elsewhere).
+        The nth simple root is the root in which there 1s in
+        the nth and (n-1)th positions, and zeroes elsewhere.
+
+        This method returns the ith simple root for the D series.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("D4")
+        >>> c.simple_root(2)
+        [0, 1, -1, 0]
+
+        """
+
+        n = self.n
+        if i < n:
+            return self.basic_root(i-1, i)
+        else:
+            root = [0]*n
+            root[n-2] = 1
+            root[n-1] = 1
+            return root
+
+
+    def positive_roots(self):
+        """
+        This method generates all the positive roots of
+        A_n.  This is half of all of the roots of D_n
+        by multiplying all the positive roots by -1 we
+        get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+        """
+
+        n = self.n
+        posroots = {}
+        k = 0
+        for i in range(0, n-1):
+            for j in range(i+1, n):
+               k += 1
+               posroots[k] = self.basic_root(i, j)
+               k += 1
+               root = self.basic_root(i, j)
+               root[j] = 1
+               posroots[k] = root
+        return posroots
+
+    def roots(self):
+        """
+        Returns the total number of roots for D_n"
+        """
+
+        n = self.n
+        return 2*n*(n-1)
+
+    def cartan_matrix(self):
+        """
+        Returns the Cartan matrix for D_n.
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType('D4')
+        >>> c.cartan_matrix()
+            Matrix([
+            [ 2, -1,  0,  0],
+            [-1,  2, -1, -1],
+            [ 0, -1,  2,  0],
+            [ 0, -1,  0,  2]])
+
+        """
+
+        n = self.n
+        m = 2*eye(n)
+        i = 1
+        while i < n-2:
+           m[i,i+1] = -1
+           m[i,i-1] = -1
+           i += 1
+        m[n-2, n-3] = -1
+        m[n-3, n-1] = -1
+        m[n-1, n-3] = -1
+        m[0, 1] = -1
+        return m
+
+    def basis(self):
+        """
+        Returns the number of independent generators of D_n
+        """
+        n = self.n
+        return n*(n-1)/2
+
+    def lie_algebra(self):
+        """
+        Returns the Lie algebra associated with D_n"
+        """
+
+        n = self.n
+        return "so(" + str(2*n) + ")"
+
+    def dynkin_diagram(self):
+        n = self.n
+        diag = " "*4*(n-3) + str(n-1) + "\n"
+        diag += " "*4*(n-3) + "0\n"
+        diag += " "*4*(n-3) +"|\n"
+        diag += " "*4*(n-3) + "|\n"
+        diag += "---".join("0" for i in range(1,n)) + "\n"
+        diag += "   ".join(str(i) for i in range(1, n-1)) + "   "+str(n)
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/type_e.py

@@ -0,0 +1,287 @@
+from .cartan_type import Standard_Cartan
+from sympy.core.backend import eye, Rational
+
+
+class TypeE(Standard_Cartan):
+
+    def __new__(cls, n):
+        if n < 6 or n > 8:
+            raise ValueError("Invalid value of n")
+        return Standard_Cartan.__new__(cls, "E", n)
+
+    def dimension(self):
+        """Dimension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("E6")
+        >>> c.dimension()
+        8
+        """
+
+        return 8
+
+    def basic_root(self, i, j):
+        """
+        This is a method just to generate roots
+        with a -1 in the ith position and a 1
+        in the jth position.
+
+        """
+
+        root = [0]*8
+        root[i] = -1
+        root[j] = 1
+        return root
+
+    def simple_root(self, i):
+        """
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        This method returns the ith simple root for E_n.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("E6")
+        >>> c.simple_root(2)
+        [1, 1, 0, 0, 0, 0, 0, 0]
+        """
+        n = self.n
+        if i == 1:
+            root = [-0.5]*8
+            root[0] = 0.5
+            root[7] = 0.5
+            return root
+        elif i == 2:
+            root = [0]*8
+            root[1] = 1
+            root[0] = 1
+            return root
+        else:
+            if i in (7, 8) and n == 6:
+                raise ValueError("E6 only has six simple roots!")
+            if i == 8 and n == 7:
+                raise ValueError("E7 has only 7 simple roots!")
+
+            return self.basic_root(i - 3, i - 2)
+
+    def positive_roots(self):
+        """
+        This method generates all the positive roots of
+        A_n.  This is half of all of the roots of E_n;
+        by multiplying all the positive roots by -1 we
+        get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+        """
+        n = self.n
+        if n == 6:
+            posroots = {}
+            k = 0
+            for i in range(n-1):
+                for j in range(i+1, n-1):
+                    k += 1
+                    root = self.basic_root(i, j)
+                    posroots[k] = root
+                    k += 1
+                    root = self.basic_root(i, j)
+                    root[i] = 1
+                    posroots[k] = root
+
+            root = [Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2),
+                    Rational(-1, 2), Rational(-1, 2), Rational(1, 2)]
+            for a in range(0, 2):
+                for b in range(0, 2):
+                    for c in range(0, 2):
+                        for d in range(0, 2):
+                            for e in range(0, 2):
+                                if (a + b + c + d + e)%2 == 0:
+                                    k += 1
+                                    if a == 1:
+                                        root[0] = Rational(-1, 2)
+                                    if b == 1:
+                                        root[1] = Rational(-1, 2)
+                                    if c == 1:
+                                        root[2] = Rational(-1, 2)
+                                    if d == 1:
+                                        root[3] = Rational(-1, 2)
+                                    if e == 1:
+                                        root[4] = Rational(-1, 2)
+                                    posroots[k] = root
+
+            return posroots
+        if n == 7:
+            posroots = {}
+            k = 0
+            for i in range(n-1):
+                for j in range(i+1, n-1):
+                    k += 1
+                    root = self.basic_root(i, j)
+                    posroots[k] = root
+                    k += 1
+                    root = self.basic_root(i, j)
+                    root[i] = 1
+                    posroots[k] = root
+
+            k += 1
+            posroots[k] = [0, 0, 0, 0, 0, 1, 1, 0]
+            root = [Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2),
+                    Rational(-1, 2), Rational(-1, 2), Rational(1, 2)]
+            for a in range(0, 2):
+                for b in range(0, 2):
+                    for c in range(0, 2):
+                        for d in range(0, 2):
+                            for e in range(0, 2):
+                                for f in range(0, 2):
+                                    if (a + b + c + d + e + f)%2 == 0:
+                                        k += 1
+                                        if a == 1:
+                                            root[0] = Rational(-1, 2)
+                                        if b == 1:
+                                            root[1] = Rational(-1, 2)
+                                        if c == 1:
+                                            root[2] = Rational(-1, 2)
+                                        if d == 1:
+                                            root[3] = Rational(-1, 2)
+                                        if e == 1:
+                                            root[4] = Rational(-1, 2)
+                                        if f == 1:
+                                            root[5] = Rational(1, 2)
+                                        posroots[k] = root
+
+            return posroots
+        if n == 8:
+            posroots = {}
+            k = 0
+            for i in range(n):
+                for j in range(i+1, n):
+                    k += 1
+                    root = self.basic_root(i, j)
+                    posroots[k] = root
+                    k += 1
+                    root = self.basic_root(i, j)
+                    root[i] = 1
+                    posroots[k] = root
+
+            root = [Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2),
+                    Rational(-1, 2), Rational(-1, 2), Rational(1, 2)]
+            for a in range(0, 2):
+                for b in range(0, 2):
+                    for c in range(0, 2):
+                        for d in range(0, 2):
+                            for e in range(0, 2):
+                                for f in range(0, 2):
+                                    for g in range(0, 2):
+                                        if (a + b + c + d + e + f + g)%2 == 0:
+                                            k += 1
+                                            if a == 1:
+                                                root[0] = Rational(-1, 2)
+                                            if b == 1:
+                                                root[1] = Rational(-1, 2)
+                                            if c == 1:
+                                                root[2] = Rational(-1, 2)
+                                            if d == 1:
+                                                root[3] = Rational(-1, 2)
+                                            if e == 1:
+                                                root[4] = Rational(-1, 2)
+                                            if f == 1:
+                                                root[5] = Rational(1, 2)
+                                            if g == 1:
+                                                root[6] = Rational(1, 2)
+                                            posroots[k] = root
+
+            return posroots
+
+
+
+    def roots(self):
+        """
+        Returns the total number of roots of E_n
+        """
+
+        n = self.n
+        if n == 6:
+            return 72
+        if n == 7:
+            return 126
+        if n == 8:
+            return 240
+
+
+    def cartan_matrix(self):
+        """
+        Returns the Cartan matrix for G_2
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType('A4')
+        >>> c.cartan_matrix()
+        Matrix([
+        [ 2, -1,  0,  0],
+        [-1,  2, -1,  0],
+        [ 0, -1,  2, -1],
+        [ 0,  0, -1,  2]])
+
+
+        """
+
+        n = self.n
+        m = 2*eye(n)
+        i = 3
+        while i < n-1:
+            m[i, i+1] = -1
+            m[i, i-1] = -1
+            i += 1
+        m[0, 2] = m[2, 0] = -1
+        m[1, 3] = m[3, 1] = -1
+        m[2, 3] = -1
+        m[n-1, n-2] = -1
+        return m
+
+
+    def basis(self):
+        """
+        Returns the number of independent generators of E_n
+        """
+
+        n = self.n
+        if n == 6:
+            return 78
+        if n == 7:
+            return 133
+        if n == 8:
+            return 248
+
+    def dynkin_diagram(self):
+        n = self.n
+        diag = " "*8 + str(2) + "\n"
+        diag += " "*8 + "0\n"
+        diag += " "*8 + "|\n"
+        diag += " "*8 + "|\n"
+        diag += "---".join("0" for i in range(1, n)) + "\n"
+        diag += "1   " + "   ".join(str(i) for i in range(3, n+1))
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/type_f.py

@@ -0,0 +1,162 @@
+from .cartan_type import Standard_Cartan
+from sympy.core.backend import Matrix, Rational
+
+
+class TypeF(Standard_Cartan):
+
+    def __new__(cls, n):
+        if n != 4:
+            raise ValueError("n should be 4")
+        return Standard_Cartan.__new__(cls, "F", 4)
+
+    def dimension(self):
+        """Dimension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("F4")
+        >>> c.dimension()
+        4
+        """
+
+        return 4
+
+
+    def basic_root(self, i, j):
+        """Generate roots with 1 in ith position and -1 in jth position
+
+        """
+
+        n = self.n
+        root = [0]*n
+        root[i] = 1
+        root[j] = -1
+        return root
+
+    def simple_root(self, i):
+        """The ith simple root of F_4
+
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("F4")
+        >>> c.simple_root(3)
+        [0, 0, 0, 1]
+
+        """
+
+        if i < 3:
+            return self.basic_root(i-1, i)
+        if i == 3:
+            root = [0]*4
+            root[3] = 1
+            return root
+        if i == 4:
+            root = [Rational(-1, 2)]*4
+            return root
+
+    def positive_roots(self):
+        """Generate all the positive roots of A_n
+
+        This is half of all of the roots of F_4; by multiplying all the
+        positive roots by -1 we get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+
+        """
+
+        n = self.n
+        posroots = {}
+        k = 0
+        for i in range(0, n-1):
+            for j in range(i+1, n):
+               k += 1
+               posroots[k] = self.basic_root(i, j)
+               k += 1
+               root = self.basic_root(i, j)
+               root[j] = 1
+               posroots[k] = root
+
+        for i in range(0, n):
+            k += 1
+            root = [0]*n
+            root[i] = 1
+            posroots[k] = root
+
+        k += 1
+        root = [Rational(1, 2)]*n
+        posroots[k] = root
+        for i in range(1, 4):
+            k += 1
+            root = [Rational(1, 2)]*n
+            root[i] = Rational(-1, 2)
+            posroots[k] = root
+
+        posroots[k+1] = [Rational(1, 2), Rational(1, 2), Rational(-1, 2), Rational(-1, 2)]
+        posroots[k+2] = [Rational(1, 2), Rational(-1, 2), Rational(1, 2), Rational(-1, 2)]
+        posroots[k+3] = [Rational(1, 2), Rational(-1, 2), Rational(-1, 2), Rational(1, 2)]
+        posroots[k+4] = [Rational(1, 2), Rational(-1, 2), Rational(-1, 2), Rational(-1, 2)]
+
+        return posroots
+
+
+    def roots(self):
+        """
+        Returns the total number of roots for F_4
+        """
+        return 48
+
+    def cartan_matrix(self):
+        """The Cartan matrix for F_4
+
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType('A4')
+        >>> c.cartan_matrix()
+        Matrix([
+        [ 2, -1,  0,  0],
+        [-1,  2, -1,  0],
+        [ 0, -1,  2, -1],
+        [ 0,  0, -1,  2]])
+        """
+
+        m = Matrix( 4, 4, [2, -1, 0, 0, -1, 2, -2, 0, 0,
+            -1, 2, -1, 0, 0, -1, 2])
+        return m
+
+    def basis(self):
+        """
+        Returns the number of independent generators of F_4
+        """
+        return 52
+
+    def dynkin_diagram(self):
+        diag = "0---0=>=0---0\n"
+        diag += "   ".join(str(i) for i in range(1, 5))
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/type_g.py

@@ -0,0 +1,111 @@
+# -*- coding: utf-8 -*-
+
+from .cartan_type import Standard_Cartan
+from sympy.core.backend import Matrix
+
+class TypeG(Standard_Cartan):
+
+    def __new__(cls, n):
+        if n != 2:
+            raise ValueError("n should be 2")
+        return Standard_Cartan.__new__(cls, "G", 2)
+
+
+    def dimension(self):
+        """Dimension of the vector space V underlying the Lie algebra
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("G2")
+        >>> c.dimension()
+        3
+        """
+        return 3
+
+    def simple_root(self, i):
+        """The ith simple root of G_2
+
+        Every lie algebra has a unique root system.
+        Given a root system Q, there is a subset of the
+        roots such that an element of Q is called a
+        simple root if it cannot be written as the sum
+        of two elements in Q.  If we let D denote the
+        set of simple roots, then it is clear that every
+        element of Q can be written as a linear combination
+        of elements of D with all coefficients non-negative.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("G2")
+        >>> c.simple_root(1)
+        [0, 1, -1]
+
+        """
+        if i == 1:
+            return [0, 1, -1]
+        else:
+            return [1, -2, 1]
+
+    def positive_roots(self):
+        """Generate all the positive roots of A_n
+
+        This is half of all of the roots of A_n; by multiplying all the
+        positive roots by -1 we get the negative roots.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("A3")
+        >>> c.positive_roots()
+        {1: [1, -1, 0, 0], 2: [1, 0, -1, 0], 3: [1, 0, 0, -1], 4: [0, 1, -1, 0],
+                5: [0, 1, 0, -1], 6: [0, 0, 1, -1]}
+
+        """
+
+        roots = {1: [0, 1, -1], 2: [1, -2, 1], 3: [1, -1, 0], 4: [1, 0, 1],
+                5: [1, 1, -2], 6: [2, -1, -1]}
+        return roots
+
+    def roots(self):
+        """
+        Returns the total number of roots of G_2"
+        """
+        return 12
+
+    def cartan_matrix(self):
+        """The Cartan matrix for G_2
+
+        The Cartan matrix matrix for a Lie algebra is
+        generated by assigning an ordering to the simple
+        roots, (alpha[1], ...., alpha[l]).  Then the ijth
+        entry of the Cartan matrix is (<alpha[i],alpha[j]>).
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.cartan_type import CartanType
+        >>> c = CartanType("G2")
+        >>> c.cartan_matrix()
+        Matrix([
+            [ 2, -1],
+            [-3,  2]])
+
+        """
+
+        m = Matrix( 2, 2, [2, -1, -3, 2])
+        return m
+
+    def basis(self):
+        """
+        Returns the number of independent generators of G_2
+        """
+        return 14
+
+    def dynkin_diagram(self):
+        diag = "0≡<≡0\n1   2"
+        return diag

lib/python3.12/site-packages/sympy/liealgebras/weyl_group.py

@@ -0,0 +1,403 @@
+# -*- coding: utf-8 -*-
+
+from .cartan_type import CartanType
+from mpmath import fac
+from sympy.core.backend import Matrix, eye, Rational, igcd
+from sympy.core.basic import Atom
+
+class WeylGroup(Atom):
+
+    """
+    For each semisimple Lie group, we have a Weyl group.  It is a subgroup of
+    the isometry group of the root system.  Specifically, it's the subgroup
+    that is generated by reflections through the hyperplanes orthogonal to
+    the roots.  Therefore, Weyl groups are reflection groups, and so a Weyl
+    group is a finite Coxeter group.
+
+    """
+
+    def __new__(cls, cartantype):
+        obj = Atom.__new__(cls)
+        obj.cartan_type = CartanType(cartantype)
+        return obj
+
+    def generators(self):
+        """
+        This method creates the generating reflections of the Weyl group for
+        a given Lie algebra.  For a Lie algebra of rank n, there are n
+        different generating reflections.  This function returns them as
+        a list.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.weyl_group import WeylGroup
+        >>> c = WeylGroup("F4")
+        >>> c.generators()
+        ['r1', 'r2', 'r3', 'r4']
+        """
+        n = self.cartan_type.rank()
+        generators = []
+        for i in range(1, n+1):
+            reflection = "r"+str(i)
+            generators.append(reflection)
+        return generators
+
+    def group_order(self):
+        """
+        This method returns the order of the Weyl group.
+        For types A, B, C, D, and E the order depends on
+        the rank of the Lie algebra.  For types F and G,
+        the order is fixed.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.weyl_group import WeylGroup
+        >>> c = WeylGroup("D4")
+        >>> c.group_order()
+        192.0
+        """
+        n = self.cartan_type.rank()
+        if self.cartan_type.series == "A":
+            return fac(n+1)
+
+        if self.cartan_type.series in ("B", "C"):
+            return fac(n)*(2**n)
+
+        if self.cartan_type.series == "D":
+            return fac(n)*(2**(n-1))
+
+        if self.cartan_type.series == "E":
+            if n == 6:
+                return 51840
+            if n == 7:
+                return 2903040
+            if n == 8:
+                return 696729600
+        if self.cartan_type.series == "F":
+            return 1152
+
+        if self.cartan_type.series == "G":
+            return 12
+
+    def group_name(self):
+        """
+        This method returns some general information about the Weyl group for
+        a given Lie algebra.  It returns the name of the group and the elements
+        it acts on, if relevant.
+        """
+        n = self.cartan_type.rank()
+        if self.cartan_type.series == "A":
+            return "S"+str(n+1) + ": the symmetric group acting on " + str(n+1) + " elements."
+
+        if self.cartan_type.series in ("B", "C"):
+            return "The hyperoctahedral group acting on " + str(2*n) + " elements."
+
+        if self.cartan_type.series == "D":
+            return "The symmetry group of the " + str(n) + "-dimensional demihypercube."
+
+        if self.cartan_type.series == "E":
+            if n == 6:
+                return "The symmetry group of the 6-polytope."
+
+            if n == 7:
+                return "The symmetry group of the 7-polytope."
+
+            if n == 8:
+                return "The symmetry group of the 8-polytope."
+
+        if self.cartan_type.series == "F":
+            return "The symmetry group of the 24-cell, or icositetrachoron."
+
+        if self.cartan_type.series == "G":
+            return "D6, the dihedral group of order 12, and symmetry group of the hexagon."
+
+    def element_order(self, weylelt):
+        """
+        This method returns the order of a given Weyl group element, which should
+        be specified by the user in the form of products of the generating
+        reflections, i.e. of the form r1*r2 etc.
+
+        For types A-F, this method current works by taking the matrix form of
+        the specified element, and then finding what power of the matrix is the
+        identity.  It then returns this power.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.weyl_group import WeylGroup
+        >>> b = WeylGroup("B4")
+        >>> b.element_order('r1*r4*r2')
+        4
+        """
+        n = self.cartan_type.rank()
+        if self.cartan_type.series == "A":
+            a = self.matrix_form(weylelt)
+            order = 1
+            while a != eye(n+1):
+                a *= self.matrix_form(weylelt)
+                order += 1
+            return order
+
+        if self.cartan_type.series == "D":
+            a = self.matrix_form(weylelt)
+            order = 1
+            while a != eye(n):
+                a *= self.matrix_form(weylelt)
+                order += 1
+            return order
+
+        if self.cartan_type.series == "E":
+            a = self.matrix_form(weylelt)
+            order = 1
+            while a != eye(8):
+                a *= self.matrix_form(weylelt)
+                order += 1
+            return order
+
+        if self.cartan_type.series == "G":
+            elts = list(weylelt)
+            reflections = elts[1::3]
+            m = self.delete_doubles(reflections)
+            while self.delete_doubles(m) != m:
+                m = self.delete_doubles(m)
+                reflections = m
+            if len(reflections) % 2 == 1:
+                return 2
+
+            elif len(reflections) == 0:
+                return 1
+
+            else:
+                if len(reflections) == 1:
+                    return 2
+                else:
+                    m = len(reflections) // 2
+                    lcm = (6 * m)/ igcd(m, 6)
+                order = lcm / m
+                return order
+
+
+        if self.cartan_type.series == 'F':
+            a = self.matrix_form(weylelt)
+            order = 1
+            while a != eye(4):
+                a *= self.matrix_form(weylelt)
+                order += 1
+            return order
+
+
+        if self.cartan_type.series in ("B", "C"):
+            a = self.matrix_form(weylelt)
+            order = 1
+            while a != eye(n):
+                a *= self.matrix_form(weylelt)
+                order += 1
+            return order
+
+    def delete_doubles(self, reflections):
+        """
+        This is a helper method for determining the order of an element in the
+        Weyl group of G2.  It takes a Weyl element and if repeated simple reflections
+        in it, it deletes them.
+        """
+        counter = 0
+        copy = list(reflections)
+        for elt in copy:
+            if counter < len(copy)-1:
+                if copy[counter + 1] == elt:
+                    del copy[counter]
+                    del copy[counter]
+            counter += 1
+
+
+        return copy
+
+
+    def matrix_form(self, weylelt):
+        """
+        This method takes input from the user in the form of products of the
+        generating reflections, and returns the matrix corresponding to the
+        element of the Weyl group.  Since each element of the Weyl group is
+        a reflection of some type, there is a corresponding matrix representation.
+        This method uses the standard representation for all the generating
+        reflections.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.weyl_group import WeylGroup
+        >>> f = WeylGroup("F4")
+        >>> f.matrix_form('r2*r3')
+        Matrix([
+        [1, 0, 0,  0],
+        [0, 1, 0,  0],
+        [0, 0, 0, -1],
+        [0, 0, 1,  0]])
+
+        """
+        elts = list(weylelt)
+        reflections = elts[1::3]
+        n = self.cartan_type.rank()
+        if self.cartan_type.series == 'A':
+            matrixform = eye(n+1)
+            for elt in reflections:
+                a = int(elt)
+                mat = eye(n+1)
+                mat[a-1, a-1] = 0
+                mat[a-1, a] = 1
+                mat[a, a-1] = 1
+                mat[a, a] = 0
+                matrixform *= mat
+            return matrixform
+
+        if self.cartan_type.series == 'D':
+            matrixform = eye(n)
+            for elt in reflections:
+                a = int(elt)
+                mat = eye(n)
+                if a < n:
+                    mat[a-1, a-1] = 0
+                    mat[a-1, a] = 1
+                    mat[a, a-1] = 1
+                    mat[a, a] = 0
+                    matrixform *= mat
+                else:
+                    mat[n-2, n-1] = -1
+                    mat[n-2, n-2] = 0
+                    mat[n-1, n-2] = -1
+                    mat[n-1, n-1] = 0
+                    matrixform *= mat
+            return matrixform
+
+        if self.cartan_type.series == 'G':
+            matrixform = eye(3)
+            for elt in reflections:
+                a = int(elt)
+                if a == 1:
+                    gen1 = Matrix([[1, 0, 0], [0, 0, 1], [0, 1, 0]])
+                    matrixform *= gen1
+                else:
+                    gen2 = Matrix([[Rational(2, 3), Rational(2, 3), Rational(-1, 3)],
+                                   [Rational(2, 3), Rational(-1, 3), Rational(2, 3)],
+                                   [Rational(-1, 3), Rational(2, 3), Rational(2, 3)]])
+                    matrixform *= gen2
+            return matrixform
+
+        if self.cartan_type.series == 'F':
+            matrixform = eye(4)
+            for elt in reflections:
+                a = int(elt)
+                if a == 1:
+                    mat = Matrix([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])
+                    matrixform *= mat
+                elif a == 2:
+                    mat = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
+                    matrixform *= mat
+                elif a == 3:
+                    mat = Matrix([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]])
+                    matrixform *= mat
+                else:
+
+                    mat = Matrix([[Rational(1, 2), Rational(1, 2), Rational(1, 2), Rational(1, 2)],
+                        [Rational(1, 2), Rational(1, 2), Rational(-1, 2), Rational(-1, 2)],
+                        [Rational(1, 2), Rational(-1, 2), Rational(1, 2), Rational(-1, 2)],
+                        [Rational(1, 2), Rational(-1, 2), Rational(-1, 2), Rational(1, 2)]])
+                    matrixform *= mat
+            return matrixform
+
+        if self.cartan_type.series == 'E':
+            matrixform = eye(8)
+            for elt in reflections:
+                a = int(elt)
+                if a == 1:
+                    mat = Matrix([[Rational(3, 4), Rational(1, 4), Rational(1, 4), Rational(1, 4),
+                        Rational(1, 4), Rational(1, 4), Rational(1, 4), Rational(-1, 4)],
+                        [Rational(1, 4), Rational(3, 4), Rational(-1, 4), Rational(-1, 4),
+                            Rational(-1, 4), Rational(-1, 4), Rational(1, 4), Rational(-1, 4)],
+                        [Rational(1, 4), Rational(-1, 4), Rational(3, 4), Rational(-1, 4),
+                        Rational(-1, 4), Rational(-1, 4), Rational(-1, 4), Rational(1, 4)],
+                        [Rational(1, 4), Rational(-1, 4), Rational(-1, 4), Rational(3, 4),
+                        Rational(-1, 4), Rational(-1, 4), Rational(-1, 4), Rational(1, 4)],
+                        [Rational(1, 4), Rational(-1, 4), Rational(-1, 4), Rational(-1, 4),
+                        Rational(3, 4), Rational(-1, 4), Rational(-1, 4), Rational(1, 4)],
+                        [Rational(1, 4), Rational(-1, 4), Rational(-1, 4), Rational(-1, 4),
+                        Rational(-1, 4), Rational(3, 4), Rational(-1, 4), Rational(1, 4)],
+                        [Rational(1, 4), Rational(-1, 4), Rational(-1, 4), Rational(-1, 4),
+                        Rational(-1, 4), Rational(-1, 4), Rational(-3, 4), Rational(1, 4)],
+                        [Rational(1, 4), Rational(-1, 4), Rational(-1, 4), Rational(-1, 4),
+                        Rational(-1, 4), Rational(-1, 4), Rational(-1, 4), Rational(3, 4)]])
+                    matrixform *= mat
+                elif a == 2:
+                    mat = eye(8)
+                    mat[0, 0] = 0
+                    mat[0, 1] = -1
+                    mat[1, 0] = -1
+                    mat[1, 1] = 0
+                    matrixform *= mat
+                else:
+                    mat = eye(8)
+                    mat[a-3, a-3] = 0
+                    mat[a-3, a-2] = 1
+                    mat[a-2, a-3] = 1
+                    mat[a-2, a-2] = 0
+                    matrixform *= mat
+            return matrixform
+
+
+        if self.cartan_type.series in ("B", "C"):
+            matrixform = eye(n)
+            for elt in reflections:
+                a = int(elt)
+                mat = eye(n)
+                if a == 1:
+                    mat[0, 0] = -1
+                    matrixform *= mat
+                else:
+                    mat[a - 2, a - 2] = 0
+                    mat[a-2, a-1] = 1
+                    mat[a - 1, a - 2] = 1
+                    mat[a -1, a - 1] = 0
+                    matrixform *= mat
+            return matrixform
+
+
+
+    def coxeter_diagram(self):
+        """
+        This method returns the Coxeter diagram corresponding to a Weyl group.
+        The Coxeter diagram can be obtained from a Lie algebra's Dynkin diagram
+        by deleting all arrows; the Coxeter diagram is the undirected graph.
+        The vertices of the Coxeter diagram represent the generating reflections
+        of the Weyl group, $s_i$.  An edge is drawn between $s_i$ and $s_j$ if the order
+        $m(i, j)$ of $s_is_j$ is greater than two.  If there is one edge, the order
+        $m(i, j)$ is 3.  If there are two edges, the order $m(i, j)$ is 4, and if there
+        are three edges, the order $m(i, j)$ is 6.
+
+        Examples
+        ========
+
+        >>> from sympy.liealgebras.weyl_group import WeylGroup
+        >>> c = WeylGroup("B3")
+        >>> print(c.coxeter_diagram())
+        0---0===0
+        1   2   3
+        """
+        n = self.cartan_type.rank()
+        if self.cartan_type.series in ("A", "D", "E"):
+            return self.cartan_type.dynkin_diagram()
+
+        if self.cartan_type.series in ("B", "C"):
+            diag = "---".join("0" for i in range(1, n)) + "===0\n"
+            diag += "   ".join(str(i) for i in range(1, n+1))
+            return diag
+
+        if self.cartan_type.series == "F":
+            diag = "0---0===0---0\n"
+            diag += "   ".join(str(i) for i in range(1, 5))
+            return diag
+
+        if self.cartan_type.series == "G":
+            diag = "0≡≡≡0\n1   2"
+            return diag