Sage and Sage-Combinat demo, Affine Schubert Calculus Workshop, Toronto -- Sage Reference Manual v4.4.4 system:sage

Sage and Sage-Combinat demo, Affine Schubert Calculus Workshop, Toronto

Basic calculations

{{{id=0| %hide pretty_print_default() /// }}}

Arithmetic:

{{{id=1| 1 + 1 /// }}} {{{id=2| factor(x^100 - 1) /// }}}

Symbolic calculations:

{{{id=3| var('y') f = sin(x) - cos(x*y) + 1 / (x^3+1) f /// }}} {{{id=4| f.integrate(x) /// }}}

Statistics:

{{{id=5| print r.summary(r.c(1,2,3,111,2,3,2,3,2,5,4)) /// }}}

Documentation

Type plot( and then press the TAB key:

{{{id=6| plot( /// }}}

Click on Help at the top of this page

There is also Live Documentation in which you can evaluate the examples.

Introspection:

{{{id=7| m = matrix([[1/2,1],[2,1]]); m /// }}} {{{id=8| m. /// }}}

Looking at the sources:

{{{id=9| m.determinant /// }}}

Algebraic constructions and categories

{{{id=10| Px.<x> = QQ[] Fx = Px.fraction_field() z = polygen(Fx,'z') Fz = Qx.extension(z^2 - x, 'z') Fz in Fields() # todo: not implemented /// True }}} {{{id=11| Pu.<u> = Fz[] # todo: not implemented (u*z) * z factor(u^2 - x^3) # todo: not implemented /// (u + x z) (u - x z) }}}
{{{id=12| for category in Fx.categories(): print category /// }}}
{{{id=13| g = sage.categories.category.category_graph() g.set_latex_options(format = "dot2tex") view(g, tightpage = True, viewer = "pdf") /// }}}

A demonstration of Sage + GAP4 + GAP3 + Chevie + Semigroupe

Let us create the Coxeter group W:

{{{id=14| W = CoxeterGroup(["H",4]); W /// }}}

It is constructed as a group of permutations, from root data given by GAP3+Chevie (thanks to Franco’s interface):

{{{id=15| W._gap_group /// CoxeterGroup("H",4) }}} {{{id=16| (W._gap_group).parent() /// Gap3 }}}

with operations on permutations implemented in Sage:

{{{id=17| W.an_element()^3 /// (1,5)(2,62)(3,7)(6,9)(8,12)(11,15)(13,17)(16,20)(18,22)(21,25)(26,29)(28,31)(30,33)(32,35)(34,37)(36,39)(38,41)(42,45)(46,48)(47,49)(50,52)(55,56)(57,58)(61,65)(63,67)(66,69)(68,72)(71,75)(73,77)(76,80)(78,82)(81,85)(86,89)(88,91)(90,93)(92,95)(94,97)(96,99)(98,101)(102,105)(106,108)(107,109)(110,112)(115,116)(117,118) }}}

and group operations implemented in GAP:

{{{id=18| len(W.conjugacy_classes_representatives()) /// 34 }}} {{{id=19| W.cardinality() /// 14400 }}}

Now, assume we want to do intensive computations on this group, requiring heavy access to the left and right Cayley graphs (e.g. Bruhat interval calculations, representation theory, ...). Then we can use Jean-Eric Pin’s Semigroupe, a software written in C:

{{{id=20| S = semigroupe.AutomaticSemigroup(W.semigroup_generators(), W.one(), category = FiniteCoxeterGroups()) /// }}}

The following triggers the full expansion of the group and its Cayley graph in memory:

{{{id=21| S.cardinality() /// 14400 }}}

And we can now iterate through the elements, in length-lexicographic order w.r.t. their reduced word:

{{{id=22| sum( x^p.length() for p in S) /// x^60 + 4*x^59 + 9*x^58 + 16*x^57 + 25*x^56 + 36*x^55 + 49*x^54 + 64*x^53 + 81*x^52 + 100*x^51 + 121*x^50 + 144*x^49 + 168*x^48 + 192*x^47 + 216*x^46 + 240*x^45 + 264*x^44 + 288*x^43 + 312*x^42 + 336*x^41 + 359*x^40 + 380*x^39 + 399*x^38 + 416*x^37 + 431*x^36 + 444*x^35 + 455*x^34 + 464*x^33 + 471*x^32 + 476*x^31 + 478*x^30 + 476*x^29 + 471*x^28 + 464*x^27 + 455*x^26 + 444*x^25 + 431*x^24 + 416*x^23 + 399*x^22 + 380*x^21 + 359*x^20 + 336*x^19 + 312*x^18 + 288*x^17 + 264*x^16 + 240*x^15 + 216*x^14 + 192*x^13 + 168*x^12 + 144*x^11 + 121*x^10 + 100*x^9 + 81*x^8 + 64*x^7 + 49*x^6 + 36*x^5 + 25*x^4 + 16*x^3 + 9*x^2 + 4*x + 1 }}} {{{id=23| S[0:10] /// [[], [0], [1], [2], [3], [0, 1], [0, 2], [0, 3], [1, 0], [1, 2]] }}} {{{id=24| S[-1] /// [0, 1, 0, 1, 0, 2, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 2, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 2, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 2, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3, 2, 0, 1, 0, 1, 2, 0, 1, 0, 2, 3] }}}

The elements of S are handles to C objects from Semigroupe:

{{{id=25| x = S.an_element() x /// [0, 1, 2, 3] }}}

Products are calculated by Semigroupe:

{{{id=26| x * x /// [0, 1, 0, 2, 0, 1, 3, 2] }}}

Powering operations are handled by Sage:

{{{id=27| x^3 /// [0, 1, 0, 2, 0, 1, 0, 2, 3, 2, 0, 1] }}} {{{id=28| x^(10^10000) /// }}}

Altogether, S is a full fledged Sage Coxeter group, which passes all the generic tests:

{{{id=29| TestSuite(S).run(verbose = True, skip = "_test_associativity") /// }}}

And of course it works for general semigroups too, and can further compute much more information about those, like the (Knuth-Bendix completion of the) relations between the generators:

{{{id=30| S.print_relations() /// aa = 1 bb = 1 cb = bc cc = 1 da = ad db = bd dd = 1 cac = aca dcd = cdc dcbabacbabcdcbabacbabcdcbabacbabcdcbabacbabcdc = cdcbabacbabcdcbabacbabcdcbabacbabcdcbabacbabcd }}}

which contains the usual commutation + braid relations.

Let’s try now the 0-Hecke monoid:

{{{id=31| from sage.combinat.j_trivial_monoids import * S = semigroupe.AutomaticSemigroup(W.simple_projections(), W.one(), by_action = True) S.cardinality() /// 14400 }}} {{{id=32| S.print_relations() /// aa = a bb = b ca = ac cc = c da = ad db = bd dd = d cbc = bcb dcd = cdc ababacbabacbabcdcbabacbabcdcbabacbabcdcbabacbabcdcbabacbabcd = 0 }}}

Let us throw in more mathematical information:

{{{id=33| W = CoxeterGroup(["A",3]) S = semigroupe.AutomaticSemigroup(W.simple_projections(), W.one(), by_action = True, category = FiniteJTrivialMonoids()) /// }}} {{{id=34| S.cardinality() /// }}} {{{id=35| H = S.algebra(QQ) H.orthogonal_idempotents() /// }}}

Cython: Python \longmapsto C

Here is a function that computes \sum_{k=0}^N k in pure Python:

{{{id=36| def mysum(N): s = int(0) for k in range(1,N): s += k return s time mysum(10^7) /// }}}

Let us compare this with the Cython version:

{{{id=37| %cython def mysum_cython(N): cdef long long s = 0 cdef int k for k in range(1,N): s += k return s time mysum_cython(10^7) /// }}}

A function to count the number of integer partitions with parts in a given set:

{{{id=38| def buying(coins, total): vlist = [ [0] * len(coins) for _ in range(total + 1) ] for i in range(total + 1): for j, coin in enumerate(coins): if j == 0: if i % coin == 0: vlist[i][j] = 1 else: k = 0 while i - k >= 0: vlist[i][j] += vlist[i - k][j - 1] k += coin return vlist[total][len(coins) - 1] [buying([1,2,5,10], i) for i in [1..20]] /// [1, 2, 2, 3, 4, 5, 6, 7, 8, 11, 12, 15, 16, 19, 22, 25, 28, 31, 34, 40] }}} {{{id=39| [1,3..10] /// [1, 3, 5, 7, 9] }}}

Let’s see how long it takes to find the number of partitions of 500 into odd parts:

{{{id=40| time buying([1,3..500], 500) /// 732986521245024 Time: CPU 3.05 s, Wall: 3.05 s }}}

Make two changes:

{{{id=41| %cython def cybuying(coins, total): cdef int i, j, k, coin vlist = [ [0] * len(coins) for _ in range(total + 1) ] for i in range(total + 1): for j, coin in enumerate(coins): if j == 0: if i % coin == 0: vlist[i][j] = 1 else: k = 0 while i - k >= 0: vlist[i][j] += vlist[i - k][j - 1] k += coin return vlist[total][len(coins) - 1] /// }}}

Surely two tiny changes in some Python code can’t make it much faster:

{{{id=42| time cybuying([1,3..500], 500) /// 732986521245024L Time: CPU 0.08 s, Wall: 0.09 s }}} {{{id=43| 3.05/0.08 /// 38.1250000000000 }}}

Sage-Combinat demonstration

{{{id=44| %hide pretty_print_default() /// }}}

Elementary combinatorics

Combinatorial objects

{{{id=45| p = Partition([3,3,2,1]) p /// }}} {{{id=46| p.pp() /// }}} {{{id=47| p.conjugate().pp() /// }}} {{{id=48| p.conjugate /// }}}
{{{id=49| s = Permutation([5,3,2,6,4,8,9,7,1]) s /// }}} {{{id=50| (p,q) = s.robinson_schensted() /// }}} {{{id=51| p.pp() /// 1 4 7 9 2 6 8 3 5 }}} {{{id=52| q.pp() /// 1 4 6 7 2 5 8 3 9 }}} {{{id=53| p.row_stabilizer() /// Permutation Group with generators [(), (7,9), (6,8), (4,7), (2,6), (1,4)] }}}

Enumerated sets (combinatorial classes)

{{{id=54| P5 = Partitions(5) P5 /// Partitions of the integer 5 }}} {{{id=55| P5.list() /// [[5], [4, 1], [3, 2], [3, 1, 1], [2, 2, 1], [2, 1, 1, 1], [1, 1, 1, 1, 1]] }}} {{{id=56| P5.cardinality() /// 7 }}} {{{id=57| Partitions(100000).cardinality() /// 27493510569775696512677516320986352688173429315980054758203125984302147328114964173055050741660736621590157844774296248940493063070200461792764493033510116079342457190155718943509725312466108452006369558934464248716828789832182345009262853831404597021307130674510624419227311238999702284408609370935531629697851569569892196108480158600569421098519 }}} {{{id=58| Permutations(20).random_element() /// [15, 6, 8, 14, 17, 16, 4, 7, 11, 3, 10, 5, 19, 9, 12, 2, 20, 18, 1, 13] }}} {{{id=59| Compositions(10).unrank(100) # TODO: non stupid algorithm /// [1, 1, 3, 1, 2, 1, 1] }}} {{{id=60| for p in StandardTableaux([3,2]): print "-----------------------------" p.pp() /// ----------------------------- 1 3 5 2 4 ----------------------------- 1 2 5 3 4 ----------------------------- 1 3 4 2 5 ----------------------------- 1 2 4 3 5 ----------------------------- 1 2 3 4 5 }}}

Trees

ToDo

Summary:

  • Every mathematical object (element, set, category, ...) is modeled by a Python object</li>
  • All combinatorial classes share a uniform interface</li>

Constructions

{{{id=61| C = DisjointUnionEnumeratedSets( [ Compositions(4), Permutations(3)] ) C /// Union of Family (Compositions of 4, Standard permutations of 3) }}} {{{id=62| C.cardinality() /// 14 }}} {{{id=63| C.list() /// [[1, 1, 1, 1], [1, 1, 2], [1, 2, 1], [1, 3], [2, 1, 1], [2, 2], [3, 1], [4], [1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]] }}}
{{{id=64| C = CartesianProduct(Compositions(8), Permutations(20)) C /// Cartesian product of Compositions of 8, Standard permutations of 20 }}} {{{id=65| C.cardinality() /// 311411457046609920000 }}} {{{id=66| C.random_element() # todo: broken /// }}}
{{{id=67| F = Family(NonNegativeIntegers(), Permutations) F /// Lazy family (Permutations(i))_{i in Set of non negative integers} }}} {{{id=68| F[1000] /// Standard permutations of 1000 }}} {{{id=69| U = DisjointUnionEnumeratedSets(F) U.cardinality() /// +Infinity }}} {{{id=70| p = iter(U) for i in range(12): print p.next() /// [] [1] [1, 2] [2, 1] [1, 2, 3] [1, 3, 2] [2, 1, 3] [2, 3, 1] [3, 1, 2] [3, 2, 1] [1, 2, 3, 4] [1, 2, 4, 3] }}} {{{id=71| for p in U: print p /// [] [1] [1, 2] [2, 1] [1, 2, 3] [1, 3, 2] [2, 1, 3] [2, 3, 1] [3, 1, 2] }}}

Summary:

  • Basic combinatorial classes + constructions give a flexible toolbox
  • This is made possible by uniform interfaces
  • Lazy algorithms and data structures for large / infinite sets (iterators, ...)

Enumeration kernels

Integer lists:

{{{id=72| IntegerVectors(10, 3, min_part = 2, max_part = 5, inner = [2, 4, 2]).list() /// [[4, 4, 2], [3, 5, 2], [3, 4, 3], [2, 5, 3], [2, 4, 4]] }}} {{{id=73| Compositions(5, max_part = 3, min_length = 2, max_length = 3).list() /// [[1, 1, 3], [1, 2, 2], [1, 3, 1], [2, 1, 2], [2, 2, 1], [2, 3], [3, 1, 1], [3, 2]] }}} {{{id=74| Partitions(5, max_slope = -1).list() /// [[5], [4, 1], [3, 2]] }}} {{{id=75| IntegerListsLex(10, length=3, min_part = 2, max_part = 5, floor = [2, 4, 2]).list() /// [[4, 4, 2], [3, 5, 2], [3, 4, 3], [2, 5, 3], [2, 4, 4]] }}} {{{id=76| IntegerListsLex(5, min_part = 1, max_part = 3, min_length = 2, max_length = 3).list() /// [[3, 2], [3, 1, 1], [2, 3], [2, 2, 1], [2, 1, 2], [1, 3, 1], [1, 2, 2], [1, 1, 3]] }}} {{{id=77| IntegerListsLex(5, min_part = 1, max_slope = -1).list() /// [[5], [4, 1], [3, 2]] }}} {{{id=78| c = Compositions(5)[1] c /// [1, 1, 1, 2] }}} {{{id=79| c = IntegerListsLex(5, min_part = 1)[1] /// }}}

Species / decomposable classes

{{{id=80| from sage.combinat.species.library import * o = var("o") /// }}}

Fibonacci words:

{{{id=81| Eps = EmptySetSpecies() Z0 = SingletonSpecies() Z1 = Eps*SingletonSpecies() FW = CombinatorialSpecies() FW.define(Eps + Z0*FW + Z1*Eps + Z1*Z0*FW) FW /// }}} {{{id=82| L = FW.isotype_generating_series().coefficients(15) L /// }}} {{{id=83| sloane_find(L) /// Searching Sloane's online database... [[45, 'Fibonacci numbers: F(n) = F(n-1) + F(n-2), F(0) = 0, F(1) = 1, F(2) = 1, ...', [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352, 24157817, 39088169]], [24595, 'a(n) = s(1)t(n) + s(2)t(n-1) + ... + s(k)t(n+1-k), where k = [ (n+1)/2 ], s = (F(2), F(3), ...), t = A023533.', [1, 0, 0, 1, 2, 3, 5, 0, 0, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1598, 2586, 4184, 6770, 10954, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28658, 46370, 75028, 121398, 196426]], [25109, 'a(n) = s(1)t(n) + s(2)t(n-1) + ... + s(k)t(n-k+1), where k = [ n/2 ], s = (F(2), F(3), F(4), ...), t = A023533.', [0, 0, 1, 2, 3, 0, 0, 0, 1, 2, 3, 5, 8, 13, 21, 34, 55, 0, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1598, 2586, 4181, 6770, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28658, 46370, 75028, 121398, 196426, 317824, 514250]], [132636, 'Fib(n) mod n^3.', [0, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 1685, 7063, 4323, 4896, 12525, 15937, 19271, 10483, 2060, 22040, 5674, 15621, 2752, 3807, 9340, 432, 46989, 19305, 11932, 62155, 31899, 12088, 22273, 3677, 32420]], [132916, 'a(0)=0; a(1)=1; a(n) = Sum a(n-k), k= 1 ... [n^(1/3)] for n>=2.', [0, 1, 1, 1, 1, 1, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 21892, 39603, 72441, 133936, 245980, 452357, 832273, 1530610, 2815240, 5178123, 9523973, 17517336, 32219432, 59260741, 108997509, 200477682]], [147316, 'A000045 Fibonacci mirror sequence Binet: f(n)=(1/5)*2^(-n) ((5 - 2 *Sqrt[5]) (1 + Sqrt[5])^n + (1 - Sqrt[5])^n(5 + 2 * Sqrt[5])).', [1597, -987, 610, -377, 233, -144, 89, -55, 34, -21, 13, -8, 5, -3, 2, -1, 1, 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597]], [39834, 'a(n+2)=-a(n+1)+a(n) (signed Fibonacci numbers); or Fibonacci numbers (A000045) extended to negative indices.', [1, 1, 0, 1, -1, 2, -3, 5, -8, 13, -21, 34, -55, 89, -144, 233, -377, 610, -987, 1597, -2584, 4181, -6765, 10946, -17711, 28657, -46368, 75025, -121393, 196418, -317811, 514229, -832040, 1346269, -2178309, 3524578, -5702887, 9227465, -14930352, 24157817]], [152163, 'a(n)=a(n-1)+a(n-2), n>1 ; a(0)=1, a(1)=-1 .', [1, -1, 0, -1, -1, -2, -3, -5, -8, -13, -21, -34, -55, -89, -144, -233, -377, -610, -987, -1597, -2584, -4181, -6765, -10946, -17711, -28657, -46368, -75025, -121393, -196418, -317811, -514229, -832040, -1346269, -2178309, -3524578, -5702887]]] }}} {{{id=84| FW3 = FW.isotypes([o]*4); FW3 /// Isomorphism types for Combinatorial species with labels [o, o, o, o] }}} {{{id=85| FW3.list() /// [o*(o*(o*(o*{}))), o*(o*(o*(({}*o)*{}))), o*(o*((({}*o)*o)*{})), o*((({}*o)*o)*(o*{})), o*((({}*o)*o)*(({}*o)*{})), (({}*o)*o)*(o*(o*{})), (({}*o)*o)*(o*(({}*o)*{})), (({}*o)*o)*((({}*o)*o)*{})] }}}

Binary trees:

{{{id=86| BT = CombinatorialSpecies() Leaf = SingletonSpecies() BT.define(Leaf+(BT*BT)) BT5 = BT.isotypes([o]*5) /// }}} {{{id=87| BT5.list() /// [o*(o*(o*(o*o))), o*(o*((o*o)*o)), o*((o*o)*(o*o)), o*((o*(o*o))*o), o*(((o*o)*o)*o), (o*o)*(o*(o*o)), (o*o)*((o*o)*o), (o*(o*o))*(o*o), ((o*o)*o)*(o*o), (o*(o*(o*o)))*o, (o*((o*o)*o))*o, ((o*o)*(o*o))*o, ((o*(o*o))*o)*o, (((o*o)*o)*o)*o] }}} {{{id=88| %hide def pbt_to_coordinates(t): e = {} queue = [t] while queue: z = queue.pop() if not isinstance(z[0], int): e[z[1]._labels[0]-1] = z queue.extend(z) coord = [(len(e[i][0]._labels) * len(e[i][1]._labels)) for i in range(len(e))] return sage.geometry.polyhedra.Polytopes.project_1(coord) K4 = Polyhedron(vertices=[pbt_to_coordinates(t) for t in BT.isotypes(range(5))]) K4.show(fill=True).show(frame=False) /// }}}

Juggling automaton:

{{{id=89| F = SingletonSpecies() state_labels = map(tuple, Permutations([0,0,1,1,1]).list()) states = dict((i, CombinatorialSpecies()) for i in state_labels) def successors(state): newstate = state[1:]+(0,) if state[0] == 0: return [(0, newstate)] return [(i+1, newstate[0:i] + (1,) + newstate[i+1:]) for i in range(0, len(state)) if newstate[i] == 0] for state in state_labels: states[state].define( sum( [states[target]*F for (height, target) in successors(state)], Eps )) states /// {(1, 1, 0, 1, 0): Combinatorial species, (0, 1, 1, 0, 1): Combinatorial species, (1, 1, 1, 0, 0): Combinatorial species, (1, 0, 1, 0, 1): Combinatorial species, (0, 1, 0, 1, 1): Combinatorial species, (1, 0, 0, 1, 1): Combinatorial species, (0, 1, 1, 1, 0): Combinatorial species, (1, 0, 1, 1, 0): Combinatorial species, (0, 0, 1, 1, 1): Combinatorial species, (1, 1, 0, 0, 1): Combinatorial species} }}} {{{id=90| states[(1,1,1,0,0)].isotypes([o]*5).cardinality() /// 165 }}}

Lattice points of polytopes

{{{id=91| A=random_matrix(ZZ,3,6,x=7) L=LatticePolytope(A) L.plot3d() /// }}} {{{id=92| L.npoints() # should be cardinality! /// 28 }}}

This example used PALP and J-mol

Graphs up to an isomorphism

{{{id=93| show(graphs(5, lambda G: G.size() <= 4)) /// }}}

Words

An infinite periodic word:

{{{id=94| p = Word([0,1,1,0,1,0,1]) ^ Infinity p /// word: 0110101011010101101010110101011010101101... }}}

The Fibonacci word:

{{{id=95| f = words.FibonacciWord() f /// word: 0100101001001010010100100101001001010010... }}}

The Thue-Morse word:

{{{id=96| t = words.ThueMorseWord() t /// word: 0110100110010110100101100110100110010110... }}}

A random word over the alphabet [0, 1] of length 1000:

{{{id=97| r = words.RandomWord(1000,2) r /// word: 0010101011101100110000000110000111011100... }}}

The fixed point of a morphism:

{{{id=98| m = WordMorphism('a->acabb,b->bcacacbb,c->baba') w = m.fixed_point('a') /// }}} {{{id=99| w /// word: acabbbabaacabbbcacacbbbcacacbbbcacacbbac... }}}

Their prefixes of length 1000:

{{{id=100| pp = p[:1000] ff = f[:1000] tt = t[:1000] ww = w[:1000] /// }}}

A comparison of their complexity function:

{{{id=101| A = list_plot([pp.number_of_factors(i) for i in range(100)], color='red') B = list_plot([ff.number_of_factors(i) for i in range(100)], color='blue') C = list_plot([tt.number_of_factors(i) for i in range(100)], color='green') D = list_plot([r.number_of_factors(i) for i in range(100)], color='black') E = list_plot([ww.number_of_factors(i) for i in range(100)], color='orange') A + B + C + D + E /// }}}

Construction of a permutation and builds its associated Rauzy diagram:

{{{id=102| p = iet.Permutation('a b c d', 'd c b a') p /// a b c d d c b a }}} {{{id=103| r = p.rauzy_diagram() r /// Rauzy diagram with 7 permutations }}} {{{id=104| r.graph().plot() /// }}}

Let us now construct a self-similar interval exchange transformation associated to a loop in the Rauzy diagram:

{{{id=105| g0 = r.path(p, 't', 't', 'b', 't') g1 = r.path(p, 'b', 'b', 't', 'b') g = g0 + g1 m = g.matrix() v = m.eigenvectors_right()[-1][1][0] T = iet.IntervalExchangeTransformation(p, v) /// }}}

We can plot it and all its power:

{{{id=106| T.plot() (T*T).plot() (T*T*T).plot() /// }}}

Check the self similarity of T:

{{{id=107| T.rauzy_diagram(iterations=8).normalize(T.length()) == T /// True }}}

And get the symbolic coding of 0 using the substitution associated to the path:

{{{id=108| s = g.orbit_substitution() s.fixed_point('a') /// word: adbdadcdadbdbdadcdadbdadcdadccdadcdadbda... }}}

Predefined algebraic structures

Root systems, Coxeter groups, ...

{{{id=109| L = RootSystem(['A',2,1]).weight_space() L.plot(size=[[-1..1],[-1..1]],alcovewalks=[[0,2,0,1,2,1,2,0,2,1]]) /// }}} {{{id=110| CartanType.samples() /// [['A', 1], ['A', 5], ['B', 1], ['B', 5], ['C', 1], ['C', 5], ['D', 2], ['D', 3], ['D', 5], ['E', 6], ['E', 7], ['E', 8], ['F', 4], ['G', 2], ['I', 5], ['H', 3], ['H', 4], ['A', 1, 1], ['A', 5, 1], ['B', 1, 1], ['B', 5, 1], ['C', 1, 1], ['C', 5, 1], ['D', 3, 1], ['D', 5, 1], ['E', 6, 1], ['E', 7, 1], ['E', 8, 1], ['F', 4, 1], ['G', 2, 1], ['B', 5, 1]^*, ['C', 4, 1]^*, ['F', 4, 1]^*, ['G', 2, 1]^*, ['BC', 1, 2], ['BC', 5, 2]] }}} {{{id=111| T = CartanType(["E", 8, 1]) print T.dynkin_diagram() /// O 2 | | O---O---O---O---O---O---O---O 1 3 4 5 6 7 8 0 E8~ }}} {{{id=112| T.is_simply_laced(), T.is_finite(), T.is_crystalographic() /// (True, False, True) }}} {{{id=113| W = WeylGroup(["B", 3]) W /// Weyl Group of type ['B', 3] (as a matrix group acting on the ambient space) }}} {{{id=114| W.cayley_graph(side = "left").plot3d(color_by_label = True) /// }}} {{{id=115| print W.character_table() # Thanks GAP! /// CT1 2 4 4 3 3 4 3 1 1 3 4 3 1 . . . . . 1 1 . 1 1a 2a 2b 4a 2c 2d 6a 3a 4b 2e X.1 1 1 1 1 1 1 1 1 1 1 X.2 1 1 1 -1 -1 -1 -1 1 1 -1 X.3 1 1 -1 -1 1 -1 1 1 -1 1 X.4 1 1 -1 1 -1 1 -1 1 -1 -1 X.5 2 2 . . -2 . 1 -1 . -2 X.6 2 2 . . 2 . -1 -1 . 2 X.7 3 -1 1 1 1 -1 . . -1 -3 X.8 3 -1 -1 -1 1 1 . . 1 -3 X.9 3 -1 -1 1 -1 -1 . . 1 3 X.10 3 -1 1 -1 -1 1 . . -1 3 }}} {{{id=116| rho = SymmetricGroupRepresentation([3, 2], "orthogonal"); rho /// Orthogonal representation of the symmetric group corresponding to [3, 2] }}} {{{id=117| rho([1, 3, 2, 4, 5]) /// 1 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{2} & \frac{1}{2} \, \sqrt{3} & 0 & 0 \\ 0 & \frac{1}{2} \, \sqrt{3} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{2} & \frac{1}{2} \, \sqrt{3} \\ 0 & 0 & 0 & \frac{1}{2} \, \sqrt{3} & \frac{1}{2} }}}

Affine Weyl groups:

{{{id=118| W = WeylGroup(["C", 3, 1]) W /// Weyl Group of type ['C', 3, 1] (as a matrix group acting on the root space) }}} {{{id=119| W.category() /// Category of affine weyl groups }}} {{{id=120| W.an_element() /// [-1 1 0 0] [ 0 1 0 0] [ 0 0 1 0] [ 0 0 0 1] }}} {{{id=121| W.from_reduced_word([1,2,3,0,3,0,3,2,1,3,3,2]).stanley_symmetric_function() /// 256*m[1, 1, 1, 1, 1, 1] + 128*m[2, 1, 1, 1, 1] + 64*m[2, 2, 1, 1] + 32*m[2, 2, 2] + 48*m[3, 1, 1, 1] + 24*m[3, 2, 1] + 8*m[3, 3] + 16*m[4, 1, 1] + 8*m[4, 2] + 4*m[5, 1] }}}

Symmetric functions

Classical basis:

{{{id=122| Sym = SymmetricFunctions(QQ) Sym /// Symmetric Functions over Rational Field }}} {{{id=123| s = Sym.schur() h = Sym.complete() e = Sym.elementary() m = Sym.monomial() p = Sym.powersum() /// }}} {{{id=124| m(( ( h[2,1] * ( 1 + 3 * p[2,1]) ) + s[2](s[3]))) /// }}}

Jack polynomials:

{{{id=125| Sym = SymmetricFunctions(QQ['t']) Jack = Sym.jack_polynomials() # todo: not implemented P = Jack.P(); J = Jack.J(); Q = Jack.Q() # todo: not implemented J(P[2,1]) # todo: not implemented /// AttributeError: 'SymmetricFunctions_with_category' object has no attribute 'jack_polynomials' }}}

Macdonald polynomials:

{{{id=126| J = MacdonaldPolynomialsJ(QQ) P = MacdonaldPolynomialsP(QQ) Q = MacdonaldPolynomialsQ(QQ) J /// Macdonald polynomials in the J basis over Fraction Field of Multivariate Polynomial Ring in q, t over Rational Field }}} {{{id=127| f = P(J[2,2] + 3 * Q[3,1]) f /// (q^2*t^6-q^2*t^5-q^2*t^4-q*t^5+q^2*t^3+2*q*t^3+t^3-q*t-t^2-t+1)*McdP[2, 2] + ((3*q^3*t^5-6*q^3*t^4+3*q^3*t^3-3*q^2*t^4+6*q^2*t^3-3*q^2*t^2-3*q*t^3+6*q*t^2-3*q*t+3*t^2-6*t+3)/(q^7*t-2*q^6*t+2*q^4*t-q^4-q^3*t+2*q^3-2*q+1))*McdP[3, 1] }}} {{{id=128| f /// }}} {{{id=129| Sym = SymmetricFunctions(J.base_ring()) s = Sym.s() s(f) /// }}}

Semigroups

{{{id=130| Cat = FiniteSemigroups() Cat /// Category of finite semigroups }}} {{{id=131| Cat.category_graph().plot(layout = "acyclic") /// }}} {{{id=132| S = Cat.example(alphabet = ('a','b','c')) S /// An example of finite semi-group: the left regular band generated by ('a', 'b', 'c') }}} {{{id=133| S.cayley_graph(side = "left", simple = True).plot() /// }}} {{{id=134| S.j_transversal_of_idempotents() /// ['cab', 'ca', 'ab', 'cb', 'a', 'c', 'b'] }}} {{{id=135| S /// }}}

Hopf algebras

{{{id=136| Cat = HopfAlgebrasWithBasis(QQ); Cat /// Category of hopf algebras with basis over Rational Field }}} {{{id=137| g = Cat.category_graph() g.set_latex_options(format = "dot2tex") view(g, tightpage = True, viewer = "pdf") /// }}} {{{id=138| Cat H = Cat.example() H /// The Hopf algebra of the Dihedral group of order 6 as a permutation group over Rational Field }}} {{{id=139| H /// }}}

A real life example

{{{id=140| def path_to_line(path, grid=True): vert = lambda x, y: circle((x, y), .05, rgbcolor=(0, 0, 1), fill=True) pline = [(0,0)] vertices = [vert(0,0)] h = 0 maxh = h ln = len(path) for x, y in enumerate(path): h += y if h > maxh: maxh = h pline += [(x+1, h)] vertices += [vert(x+1, h)] plotted_path = line(pline) + sum(vertices) if grid: gridline = lambda a, b, c, d: line([(a, b), (c,d)], rgbcolor=(0,) * 3, linestyle='--', alpha=.25) # vertical gridlines grid = [gridline(x, 0, x, maxh) for x in [1..ln]] # horizontal gridlines for y in [1..maxh]: grid += [gridline(0, y, ln, y)] plotted_path += sum(grid) plotted_path.set_aspect_ratio(1) return plotted_path from sage.combinat.backtrack import GenericBacktracker class LukPaths(GenericBacktracker): def __init__(self, n, k=1): GenericBacktracker.__init__(self, [], (0, 0)) self._n = n self._k = k if n < 0 or k < 1 or n % (k+1) != 0: def jane_stop_this_crazy_thing(*args): raise StopIteration self._rec = jane_stop_this_crazy_thing def _rec(self, path, state): ln, ht = state if ln < self._n: # if we're high enough, we can yield a path with a # k-downstep at the end if ht >= self._k: yield path + [-self._k], (ln + 1, ht - self._k), False # if the path isn't too high, it can also take an upstep if ht / self._k < (self._n - ln): yield path + [1], (ln + 1, ht + 1), False else: # if length is n, set state to None so we stop trying to # make new paths, and yield what we've got yield path, None, True plots = [path_to_line(p) for p in LukPaths(12, 2)] ga = graphics_array(plots, ceil(len(plots)/6), 6) ga.show(figsize=[9,7]) /// }}}