diff --git a/src/sage/combinat/words/automata.pyx b/src/sage/combinat/words/automata.pyx new file mode 100644 index 00000000000..7b3d853af9a --- /dev/null +++ b/src/sage/combinat/words/automata.pyx @@ -0,0 +1,722 @@ +# coding=utf8 +""" +Finite state machines + +AUTHORS: + +- Paul Mercat +""" + +#***************************************************************************** +# Copyright (C) 2013 Paul Mercat +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# This code is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# The full text of the GPL is available at: +# +# http://www.gnu.org/licenses/ +#***************************************************************************** + +from sage.graphs.digraph import DiGraph +from sage.sets.set import Set + +#test +# +#cdef emonde0_simplify_rec2 (aut, e, a, d, int *n): +# keep = False +# *n = (*n) + 1 +# for u,v,l in aut.outgoing_edges(e): +# if emonde0_simplify_rec(aut, e=v, a=a, d=d, n=n): +# keep = True + + +class Automaton (DiGraph): +# def __init__(self, I, F, **args): +# DiGraph.__init__(self, **args) +# self._I = I +# self._F = F +# return + def __repr__(self): + return "Finite automaton with %d states"%self.num_verts() + + def _latex_(self): + return "Finite automaton with $%d$ states"%self.num_verts() + + def plot (self): + r""" + Plot the automaton. + + EXAMPLES:: + + sage: from sage.combinat.words.morphism import CallableDict + sage: d = CallableDict({1:'one', 2:'zwei', 3:'trois'}) + sage: d(1), d(2), d(3) + ('one', 'zwei', 'trois') + """ + return DiGraph.plot(self, edge_labels=True) + + def copy (self): + a = DiGraph.copy(self) + if hasattr(self, 'I'): + a.I = self.I + if hasattr(self, 'F'): + a.F = self.F + if hasattr(self, 'A'): + a.A = self.A + return a + + def Alphabet (self): + if hasattr(self, 'A'): + return self.A + else: + return set(self.edge_labels()) + + def product2_ (self, A, d=None, verb=False): + L = self.Alphabet() + LA = A.Alphabet() + if d is None: + d = dict([]) + for k in L: + for ka in LA: + d[(k, ka)] = [(k, ka)] #alphabet produit par défaut +# if dA is None: +# dA = dict([]) +# for k in L: +# dA[k] = k +# if d.keys() != dA.keys(): +# print "Error : incompatible keys : %s and %s"%(d,dA) + P = Automaton(loops=True, multiedges=True) + V = self.vertices() + VA = A.vertices() + for v in V: + for va in VA: + P.add_vertex((v,va)) + for v1, v2, l in self.outgoing_edges(v, labels=True): + for va1, va2, la in A.outgoing_edges(va, labels=True): + for e in d[(l, la)]: + if verb: + print "ajout de %s, %s, %s"%((v1, va1), (v2, va2), e) + P.add_edge((v1, va1), (v2, va2), e) + if hasattr(self, 'I') and hasattr(A, 'I'): + P.I = [(u,v) for u in self.I for v in A.I] + if hasattr(self, 'F') and hasattr(A, 'F'): + P.F = [(u,v) for u in self.F for v in A.F] + P.A = set([a for l in L for la in LA if d.has_key((l,la)) for a in d[(l,la)]]) + #[u for u in set(d.values()) if u is not None] + return P + + def product2 (self, A, d=None, verb=False): + L = self.Alphabet() + LA = A.Alphabet() + if d is None: + d = dict([]) + for k in L: + for ka in LA: + d[(k, ka)] = (k, ka) #alphabet produit par défaut + P = Automaton(loops=True, multiedges=True) + S = set([(i, i2) for i in self.I for i2 in A.I]) #ensemble des sommets à explorer + E = set([]) #ensemble des sommets déjà explorés + while len(S) != 0: + (v,va) = S.pop() + P.add_vertex((v, va)) + for v1, v2, l in self.outgoing_edges(v, labels=True): + for va1, va2, la in A.outgoing_edges(va, labels=True): + if d.has_key((l, la)): + #print "has_key (%s,%s)"%(l,la) + for e in d[(l, la)]: + if verb: + print "ajout de %s, %s, %s"%((v1, va1), (v2, va2), e) + P.add_edge((v1, va1), (v2, va2), e) + if len(d[(l, la)]) > 0 and (v2, va2) not in E: + S.add((v2, va2)) + E.add((v,va)) + if hasattr(self, 'I') and hasattr(A, 'I'): + P.I = [(u,v) for u in self.I for v in A.I] + if hasattr(self, 'F') and hasattr(A, 'F'): + P.F = [(u,v) for u in self.F for v in A.F] + #print "A=..." + P.A = set([a for l in L for la in LA if d.has_key((l,la)) for a in d[(l,la)]]) + #print P.A + #[u for u in set(d.values()) if u is not None] + return P + + def product (self, A, d=None, verb=False): + L = self.Alphabet() + LA = A.Alphabet() + if d is None: + d = dict([]) + for k in L: + for ka in LA: + d[(k, ka)] = (k, ka) #alphabet produit par défaut + P = Automaton(loops=True, multiedges=True) + S = set([(i, ia) for i in self.I for ia in A.I]) #ensemble des sommets à explorer + E = set([]) #ensemble des sommets déjà explorés + while len(S) != 0: + if verb: + print S, E + (v,va) = S.pop() + P.add_vertex((v, va)) + for v1, v2, l in self.outgoing_edges(v, labels=True): + for va1, va2, la in A.outgoing_edges(va, labels=True): + if d[(l, la)] is not None: + P.add_edge((v1, va1), (v2, va2), d[(l, la)]) + if (v2, va2) not in E: + S.add((v2, va2)) + E.add((v,va)) + if hasattr(self, 'I') and hasattr(A, 'I'): + P.I = [(u,v) for u in self.I for v in A.I] + if hasattr(self, 'F') and hasattr(A, 'F'): + P.F = [(u,v) for u in self.F for v in A.F] + P.A = [d[(l,la)] for l in L for la in LA if d[(l,la)] is not None] + return P + + def product_ (self, A, d=None): + L = self.Alphabet() + LA = A.Alphabet() + if d is None: + d = dict([]) + for k in L: + for ka in LA: + d[(k, ka)] = (k, ka) #alphabet produit par défaut +# if dA is None: +# dA = dict([]) +# for k in L: +# dA[k] = k +# if d.keys() != dA.keys(): +# print "Error : incompatible keys : %s and %s"%(d,dA) + P = Automaton(loops=True, multiedges=True) + V = self.vertices() + VA = A.vertices() + for v in V: + for va in VA: + P.add_vertex((v,va)) + for v1, v2, l in self.outgoing_edges(v, labels=True): + for va1, va2, la in A.outgoing_edges(va, labels=True): + if d[(l, la)] is not None: + #print "ajout de %s, %s, %s"%((v1, va1), (v2, va2), d[(l, la)]) + P.add_edge((v1, va1), (v2, va2), d[(l, la)]) + if hasattr(self, 'I') and hasattr(A, 'I'): + P.I = [(u,v) for u in self.I for v in A.I] + if hasattr(self, 'F') and hasattr(A, 'F'): + P.F = [(u,v) for u in self.F for v in A.F] + P.A = [d[(l,la)] for l in L for la in LA if d[(l,la)] is not None] + #[u for u in set(d.values()) if u is not None] + return P + + def intersection (self, A, verb=False): + L = self.Alphabet() + LA = A.Alphabet() + d = dict([]) + for k in L: + for ka in LA: + if k == ka: + d[(k, ka)] = k + else: + d[(k, ka)] = None + if verb: + print d + res = self.product(A=A, d=d, verb=verb) + #if self.I is not None and A.I is not None: + #if hasattr(self, 'I') and hasattr(A, 'I'): + # res.I = [(u, v) for u in self.I for v in A.I] + return res + + def emonde0_rec (self, e): + global E + #print "e=%s, E=%s"%(e,E) + #from sage.sets.set import set + E.add(e) + keep = False + O = set([b for a,b,l in self.outgoing_edges(e)]) + #print "O=%s"%O + for b in O: + if self.has_vertex(b): #au cas où le sommet ait été supprimé entre-temps + if b not in E: + self.emonde0_rec (e=b) + if b in E: + keep = True + if not keep: + E.remove(e) + self.delete_vertex(e) + + def emonde0 (self, I=None): #retire les etats puits + if I is None: + if hasattr(self, 'I'): + I = self.I + if I is None: + print "The set I of initial states must be defined !" + return + global E + #from sage.sets.set import set + E = set() + for e in I: + self.emonde0_rec(e=e) + + def emonde0_simplify_rec (self, e, a, d): + keep = False + og = self.outgoing_edges(e) + O = set([v for u,v,l in og]) + cdef int ne = a.add_vertex() + d[e] = ne + for b in O: + if self.has_vertex(b): #au cas où le sommet ait été supprimé entre-temps + if not d.has_key(b): + self.emonde0_simplify_rec (e=b, a=a, d=d) + if d.has_key(b): + keep = True + if not keep: + a.delete_vertex(d.pop(e)) + else: + for u,v,l in og: + if d.has_key(v): + if not a.has_edge(ne, d[v], l): + a.add_edge(ne, d[v], l) + + def emonde0_simplify (self, I=None): #rend l'automate sans états puits et avec des états plus simples + if I is None: + if hasattr(self, 'I'): + I = self.I + if I is None: + raise ValueError("The set I of initial states must be defined !") + a = Automaton(loops=True, multiedges=True) + d = dict() + cdef int n = 0 + for e in I: + #emonde0_simplify_rec(self, e=e, a=a, d=d, n=&n) + self.emonde0_simplify_rec(e=e, a=a, d=d) + a.I = [d[i] for i in I if d.has_key(i)] + if hasattr(self, 'F'): + a.F = [d[f] for f in self.F if d.has_key(f)] + if hasattr(self, 'A'): + a.A = self.A + return a + + def emonde0_simplify_ (self, I=None, verb=False): #rend l'automate sans états puits et avec des états plus simples + if I is None: + if hasattr(self, 'I'): + I = self.I + if I is None: + raise ValueError("The set I of initial states must be defined !") + a = Automaton(loops=True, multiedges=True) + I = list(set(I)) + d = {} + cdef int n = 0 + cdef int ni + if verb: + print "I=%s"%I + while len(I) != 0: + i = I.pop() + if not d.has_key(i): + ni = n + d[i] = ni + a.add_vertex(ni) + n += 1 + else: + ni = d[i] + for u,v,l in self.outgoing_edges(i): + if not d.has_key(v): + d[v] = n + a.add_vertex(n) + I.append(v) + n += 1 + a.add_edge(ni, d[v], l) + a.I = [d[i] for i in self.I] + if hasattr(self, 'F'): + a.F = [d[f] for f in self.F if d.has_key(f)] + if hasattr(self, 'A'): + a.A = self.A + return a + + def emondeI (self, I=None): + #from sage.sets.set import set + if I is None: + if hasattr(self, 'I'): + I = self.I + else: + raise ValueError("The initial states set I must be defined !") + + G = set() + E = set(I) + while E: + #for e in E: + # E.remove(e) + e = E.pop() + if e not in G: + G.add(e) + E = E.union(set([b for (a,b,f) in self.outgoing_edges(e)])) + for e in self.vertices(): + if e not in G: + self.delete_vertex(e) + + def emondeF (self, F=None): + #from sage.sets.set import set + if F is None: + if hasattr(self, 'F'): + F = self.F + else: + raise ValueError("The final states set F must be defined !") + + G = set() + E = set(F) + while E: + #for e in E: + # E.remove(e) + e = E.pop() + if e not in G: + G.add(e) + E = E.union(set([a for (a,b,f) in self.incoming_edges(e)])) + for e in self.vertices(): + if e not in G: + self.delete_vertex(e) + + def emonde (self, I=None, F=None): + if I is None: + if hasattr(self, 'I'): + I = self.I + else: + raise ValueError("The initial states set I must be defined !") + if F is None: + if hasattr(self, 'F'): + F = self.F + else: + raise ValueError("The final states set F must be defined !") + + self.emondeI(I=I) + self.emondeF(F=F) + + def prefix (self, w, i=None, i2=None, verb=False): #retourne un automate dont le langage est w^(-1)L (self est supposé déterministe) + if i is None: + if hasattr(self, 'I'): + i = self.I[0] + else: + raise ValueError("The initial state i must be given !") + a = self.copy() + e = i #etat initial + if i2 is None: + c2 = a.add_vertex() + i2 = c2 + else: + a.add_vertex(i2) + c2 = i2 + for l in w: + found = False + if verb: print "e=%s"%e + for u,v,l2 in self.outgoing_edges(e): + if l2 == l: + found = True + e = v + c = c2 + c2 = a.add_vertex() + a.add_edge(c, c2, l) + break + if not found: + print "Erreur : le mot %s n'est pas reconnu par l'automate !"%w + #connecte + a.add_edge(c, e, l) + a.delete_vertex(c2) + a.I = [i2] + if hasattr(self, 'A'): + a.A = self.A + return a + + def succ (self, i, a): #donne un état atteint en partant de i et en lisant a (rend None si n'existe pas) + if i is None: + return + for u,v,l in self.outgoing_edges(i): + if l == a: + return v + + def transpose (self): + a = Automaton(loops=True, multiedges=True) + for u,v,l in self.edges(): + a.add_edge(v,u,l) + if hasattr(self, 'F'): + a.I = self.F + if hasattr(self, 'I'): + a.F = self.I + if hasattr(self, 'A'): + a.A = self.A + return a + + def save_graphviz (self, file): + a = self.copy() + d=dict() + for u in a.Alphabet(): + d[u]=[str(u)] + a.graphviz_to_file_named(file, edge_labels=True) + + def relabel2 (self, d): + self.allow_multiple_edges(True) + for u,v,l in self.edges(): + self.delete_edge(u, v, l) + for l2 in d[l]: + self.add_edge(u, v, l2) + self.A = [a for A in d.values() for a in A] + + def determinize_rec (self, S, A, a, nof, noempty, verb): + #from sage.sets.set import set + #A2 = set([]) + + if verb: print S + + if not a.has_key(S): + a[S] = dict([]) + o = dict([]) #associe à chaque lettre l'ensemble des états partants de S + for s in S: + for f, d, l in self.outgoing_edges(s): + #A2 = A2.union(set([l])) + #if verb: print "l=%s"%l + if not o.has_key(l): + o[l] = set() + o[l].add(d) + # for l in A: + # if l not in A2: + # if not a[S].has_key(set([])): + # a[S][set([])] = [] + # a[S][set([])] += [l] + + #for l in o.keys(): + for l in A: + if o.has_key(l) or not noempty: + if not o.has_key(l): + o[l] = set() + #if verb: print "S=%s, l=%s, o[l]=%s"%(S, l, o[l]) + if set(nof).isdisjoint(o[l]): + o[l] = Set(o[l]) + if not a[S].has_key(o[l]): + a[S][o[l]] = [] + a[S][o[l]] += [l] + for r in set(a[S].keys()): + a = self.determinize_rec(S=r, A=A, a=a, nof=nof, noempty=noempty, verb=verb) + return a + + def determinize (self, I=None, A=None, nof=set(), noempty=True, verb=False): + if I is None: + if hasattr(self, 'I'): + I = self.I + if I is None: + raise ValueError("The set I of initial states must be defined !") + if A is None: + if hasattr(self, 'A'): + A = self.A + if A is None: + A = set(self.edge_labels()) + #raise ValueError("The alphabet A must be defined !") + #from sage.sets.set import set + + a = dict() + res = Automaton(self.determinize_rec (S=Set(I), A=A, a=a, nof=nof, noempty=noempty, verb=verb), multiedges=True, loops=True) + res.I = [Set(I)] + res.A = A + if hasattr(self, 'F'): + res.F = set([v for v in res.vertices() for f in self.F if f in v]) + #res.F = set() + #for f in self.F: + # for v in res.vertices(): + # if f in v: + # res.F.add(v) + return res + + def determinize2 (self, I=None, A=None, nof=set(), noempty=True, verb=False): + if I is None: + if hasattr(self, 'I'): + I = self.I + if I is None: + raise ValueError("The set I of initial states must be defined !") + if A is None: + A = self.Alphabet() + #raise ValueError("The alphabet A must be defined !") + #from sage.sets.set import set + + nof = set(nof) + + c = 0 + a = Automaton(loops=True, multiedges=True) + SS = set([Set(I)]) + a.I = set([Set(I)]) + a.A = A + a.add_vertex(Set(I)) + while len(SS) != 0: + S = SS.pop() + o = dict([]) + for l in A: + o[l] = set([]) + for s in S: + for f, d, l in self.outgoing_edges(s): + o[l].add(d) + for l in A: + if o[l] != set([]) or not noempty: + if nof.isdisjoint(o[l]): + o[l] = Set(o[l]) + if not a.has_vertex(o[l]): + a.add_vertex(o[l]) + SS.add(o[l]) + c+=1 + if c%1000 == 0 and verb: + print c + a.add_edge(S, o[l], l) + if hasattr(self, 'F'): + a.F = set([v for v in a.vertices() for f in self.F if f in v]) + return a + + def complete (self, name=None, verb=False): + r""" + Complete the given automaton. + If the given automaton is not complete, it will add a hole state named name. + + EXAMPLES:: + + sage: BetaAdicMonoid(3, {0, 1, 3/2}).relations_automaton().complete() + Finite automaton with 4 states + """ + self.allow_multiple_edges(True) + #from sage.sets.set import set + L = self.Alphabet() + V = self.vertices() + ns = True + for v in V: + L2 = [] + for f, d, l in self.outgoing_edges(v): + L2 += [l] + if set(L2) != set(L): + if verb: print "L=%s, L2=%s"%(L,L2) + for l in L.difference(L2): + if ns: + ns = False + if name is None: + name = self.add_vertex() + V.append(name) + elif name in V: + raise ValueError("the name of the hole state must be unused") + self.add_edge(v, name, l) + return name + + def minimize (self, F=None, verb=False): + r""" + Returns the minimal automaton. + The given automaton is supposed to be deterministic. + + EXAMPLES:: + + sage: BetaAdicMonoid(3, {0, 1, 3/2}).relations_automaton().complete().minimize() + Finite automaton with 4 states + """ + puit = None + if F is None: + if hasattr(self, 'F'): + F = self.F + if F is None: + F = self.vertices() + puit = self.complete() + + if verb: print "self=%s"%self + + #from sage.sets.set import set + F = Set(F) + Q = Set(self.vertices()) + P = set([F, Q.difference(F)]) + W = set([F]) + if hasattr(self, 'L'): + L = self.L + else: + L = self.Alphabet() + while W: + if verb: print "P=%s, W=%s"%(P,W) + A = W.pop() + if verb: print "A=%s"%A + X = dict() + for s in A: + for f, d, l in self.incoming_edges(s): + #for l in L2: + if not X.has_key(l): + X[l] = set() + X[l].add(f) + #P2 = set([p for p in P]) + #while P2: + # Y = P2.pop() + stop = False + for l in L: + if verb: print "P=%s, W=%s, l=%s"%(P,W,l) + if X.has_key(l): + X[l] = Set(X[l]) + + cont = True + while cont: + cont=False + for Y in P: + I = Y.intersection(X[l]) + D = Y.difference(X[l]) + if I and D: + if verb: print "I=%s, D=%s"%(I,D) + P.discard(Y) + P = P.union([I,D]) + if Y in W: + W.discard(Y) + W |= set([I,D]) + else: + if len(I) <= len(D): + W = W.union(set([I])) + else: + W = W.union(set([D])) + #stop = True + cont = True + break + #if stop: + # break +# else: +# raise ValueError("the automaton must be complete") + #construct the new automaton + cs = dict() + for S in P: + for s in S: + cs[s] = Set(S) + a = dict() + for S in P: + S = Set(S) + a[S] = dict() + if not S.is_empty(): + for f, d, l in self.outgoing_edges(S.an_element()): + if not a[S].has_key(cs[d]): + a[S][cs[d]] = [] + a[S][cs[d]] += [l] + res = Automaton(a) #, loops=True, multiedges=True) + res.F = [Set(S) for S in P if not set(S).isdisjoint(set(F))] + res.A = L + if hasattr(self, 'I'): + if verb: print "I=%s"%I + res.I = [Set(S) for S in P for i in self.I if i in S] + if verb: print "res.I=%s"%res.I + #for i in self.I: + # for v in res.vertices(): + # if i in v: + # res.I += [v] + if verb: print "res=%s"%res + #delete the puit if it was not + if puit is not None: + if verb: print "effacement de l'etat puit %s"%puit + self.delete_vertex(puit) + for S in P: + if puit in S: + if verb: print "delete %s"%S + res.delete_vertex(Set(S)) + res = res.emonde0_simplify() + return res + + #compute the automaton recognizing the complementary + def complementary (self): + self.complete() + if hasattr(self, 'F'): + self.F = [f for f in self.vertices() if f not in self.F] + else: + self.F = self.vertices() + +# def plot (self, edge_labels=True, **options): +# DiGraph.plot(self, edge_labels=edge_labels, **options) + +# def add_state (self, *args): #ajoute l'état s à l'automate (le réinitialise si déjà présent) +# DiGraph.add_vertex(self, *args) diff --git a/src/sage/monoids/beta_adic_monoid.pyx b/src/sage/monoids/beta_adic_monoid.pyx new file mode 100644 index 00000000000..5e51f9120f8 --- /dev/null +++ b/src/sage/monoids/beta_adic_monoid.pyx @@ -0,0 +1,1521 @@ +# coding=utf8 +r""" +Beta-adic Monoids + +AUTHORS: + +- Paul Mercat (2013) + +Beta-adic monoids are finitely generated monoids with generators of the form + x -> beta*x + c +where beta is a element of a field (for example a complex number), +and c is varying in a finite set of numerals. +It permits to describe beta-adic expansions, that is writing of numbers of the form + x = c_0 + c_1*beta + c_2*beta^2 + ... +for c_i's in a finite set of numerals. +""" + +#***************************************************************************** +# Copyright (C) 2013 Paul Mercat +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# This code is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty +# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +# +# See the GNU General Public License for more details; the full text +# is available at: +# +# http://www.gnu.org/licenses/ +#***************************************************************************** + +from sage.rings.integer import Integer +from sage.rings.number_field.all import * +#from sage.structure.parent_gens import normalize_names +#from free_monoid_element import FreeMonoidElement + +from sage.sets.set import Set +from monoid import Monoid_class +from sage.rings.qqbar import QQbar +from sage.rings.padics.all import * +from sage.combinat.words.automata import Automaton + +#from sage.structure.factory import UniqueFactory +#from sage.misc.cachefunc import cached_method + +#calcul de la valeur absolue p-adique (car non encore implémenté autrement) +def absp (c, p, d): + return ((c.polynomial())(p).norm().abs())**(1/d) + +#garde la composante fortement connexe de 0 +def emonde (a, K): + for s in a.strongly_connected_components_subgraphs(): + if K.zero() in s: + return s + +class BetaAdicMonoid(Monoid_class): + r""" + ``b``-adic monoid with numerals set ``C``. + It is the beta-adic monoid generated by the set of affine transformations ``{x -> b*x + c | c in C}``. + + EXAMPLES:: + + sage: m1 = BetaAdicMonoid(3, {0,1,3}) + sage: m2 = BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + sage: b = (x^3-x-1).roots(ring=QQbar)[0][0] + sage: m3 = BetaAdicMonoid(b, {0,1}) + """ + def __init__ (self, b, C): + r""" + Construction of the b-adic monoid generated by the set of affine transformations ``{x -> b*x + c | c in C}``. + + EXAMPLES:: + sage: m1 = BetaAdicMonoid(3, {0,1,3}) + sage: m2 = BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + sage: b = (x^3-x-1).roots(ring=QQbar)[0][0] + sage: m3 = BetaAdicMonoid(b, {0,1}) + """ + #print "init BAM with (%s,%s)"%(b,C) + if b in QQbar: + # print b + pi = QQbar(b).minpoly() + K = NumberField(pi, 'b', embedding = QQbar(b)) + else: + K = b.parent() + try: + K.places() + except: + print "b=%s must be a algebraic number, ring %s not accepted."%(b,K) + +# print K + self.b = K.gen() #beta (element of an NumberField) +# print "b="; print self.b + self.C = Set([K(c) for c in C]) #set of numerals +# print "C="; print self.C + + def gen (self, i): + r""" + Return the element of C of index i. + """ +# g(x) = self.b*x+self.C[i] + return self.C[i] + + def ngens (self): + r""" + Return the number of elements of C. + """ + return len(self.C) + + def _repr_ (self): + r""" + Returns the string representation of the beta-adic monoid. + + EXAMPLES:: + + sage: BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + Monoid of b-adic expansion with b root of x^2 - x - 1 and numerals set {0, 1} + sage: BetaAdicMonoid(3, {0,1,3}) + Monoid of 3-adic expansion with numerals set {0, 1, 3} + + TESTS:: + + sage: m=BetaAdicMonoid(3/sqrt(2), {0,1}) + sage: repr(m) + 'Monoid of b-adic expansion with b root of x^2 - 9/2 and numerals set {0, 1}' + """ + + str = "" + if hasattr(self, 'ss'): + if self.ss is not None: + str=" with subshift of %s states"%self.ss.num_verts() + + from sage.rings.rational_field import QQ + if self.b in QQ: + return "Monoid of %s-adic expansion with numerals set %s"%(self.b,self.C) + str + else: + K = self.b.parent() + if K.base_field() == QQ: + return "Monoid of b-adic expansion with b root of %s and numerals set %s"%(self.b.minpoly(),self.C) + str + else: + if K.characteristic() != 0: + return "Monoid of b-adic expansion with b root of %s and numerals set %s, in characteristic %s"%(self.b.minpoly(), self.C, K.characteristic()) + str + else: + return "Monoid of b-adic expansion with b root of %s and numerals set %s"%(K.modulus(),self.C) + str + + def default_ss (self): + r""" + Returns the full subshift (given by an Automaton) corresponding to the beta-adic monoid. + + EXAMPLES:: + + sage: m=BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + sage: m.default_ss() + Finite automaton with 1 states + """ + C = self.C + ss = Automaton() + ss.allow_multiple_edges(True) + ss.allow_loops(True) + ss.add_vertex(0) + for c in C: + ss.add_edge(0, 0, c) + ss.I = [0] + ss.F = [0] + ss.A = C + return ss + + def points_exact (self, n=None, ss=None, iss=None): + r""" + Returns a set of exacts values (in the number field of beta) corresponding to the drawing of the limit set of the beta-adic monoid. + + INPUT: + + - ``n`` - integer (default: ``None``) + The number of iterations used to plot the fractal. + Default values: between ``5`` and ``16`` depending on the number of generators. + + - ``ss`` - Automaton (default: ``None``) + The subshift to associate to the beta-adic monoid for this drawing. + + - ``iss`` - set of initial states of the automaton ss (default: ``None``) + + OUTPUT: + + list of exact values + + EXAMPLES: + + #. The dragon fractal:: + + sage: m=BetaAdicMonoid(1/(1+i), {0,1}) + sage: P = m.points_exact() + sage: len(P) + 65536 + """ + #global co + + C = self.C + K = self.b.parent() + b = self.b + ng = C.cardinality() + + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + else: + ss = self.default_ss() + + if iss is None: + if hasattr(ss, 'I'): + iss = [i for i in ss.I][0] + if iss is None: + print "Donner l'état de départ iss de l'automate ss !" + return + + if n is None: + if ng == 2: + n = 16 + elif ng == 3: + n = 9 + else: + n = 5 + if n == 0: + #donne un point au hasard dans l'ensemble limite + #co = co+1 + return [0] + else: + orbit_points = set() + V = set([v for c in C for v in [ss.succ(iss, c)] if v is not None]) + orbit_points0 = dict() + for v in V: + orbit_points0[v] = self.points_exact(n=n-1, ss=ss, iss=v) + for c in C: + v = ss.succ(iss, c) + if v is not None: + orbit_points.update([b*p+c for p in orbit_points0[v]]) + #orbit_points = orbit_points.union([b*p+c for p in self.points_exact(n=n-1, ss=ss, iss=v)]) + #orbit_points0 = self.points_exact(n-1) + #orbit_points = Set([]) + #for c in C: + # v = self.succ(i, c) + # if v is not None: + # orbit_points = orbit_points.union(Set([b*p+c for p in orbit_points0])) + #print "no=%s"%orbit_points.cardinality() + return orbit_points + + def points (self, n=None, place=None, ss=None, iss=None, prec=53): + r""" + Returns a set of values (real or complex) corresponding to the drawing of the limit set of the beta-adic monoid. + + INPUT: + + - ``n`` - integer (default: ``None``) + The number of iterations used to plot the fractal. + Default values: between ``5`` and ``16`` depending on the number of generators. + + - ``place`` - place of the number field of beta (default: ``None``) + The place we should use to evaluate elements of the number field given by points_exact() + + - ``ss`` - Automaton (default: ``None``) + The subshift to associate to the beta-adic monoid for this drawing. + + - ``iss`` - set of initial states of the automaton ss (default: ``None``) + + - ``prec`` - precision of returned values (default: ``53``) + + OUTPUT: + + list of real or complex numbers + + EXAMPLES: + + #. The dragon fractal:: + + sage: m=BetaAdicMonoid(1/(1+I), {0,1}) + sage: P = m.points() + sage: len(P) + 32768 + """ + + C = self.C + K = self.b.parent() + b = self.b + ng = C.cardinality() + + if n is None: + if ng == 2: + n = 18 + elif ng == 3: + n = 14 + elif ng == 4: + n = 10 + elif ng == 5: + n = 7 + else: + n = 5 + from sage.functions.log import log + n = int(5.2/-log(abs(self.b.N(prec=prec)))) + + from sage.rings.complex_field import ComplexField + CC = ComplexField(prec) + if place is None: + if abs(b) < 1: + #garde la place courante + #place = lambda x: CC(x.n()) + return [CC(c).N(prec) for c in self.points_exact(n=n, ss=ss, iss=iss)] + else: + #choisis une place + places = K.places() + place = places[0] + for p in places: + if abs(p(b)) < 1: + place = p + #break + + #from sage.rings.qqbar import QQbar + #from sage.rings.qqbar import QQbar, AA + #if QQbar(self.b) not in AA: + # #print "not in AA !" + # return [(place(c).conjugate().N().real(), place(c).conjugate().N().imag()) for c in self.points_exact(n=n, ss=ss, iss=iss)] + #else: + # #print "in AA !" + # return [place(c).conjugate().N() for c in self.points_exact(n=n, ss=ss, iss=iss)] + return [place(c).N(prec) for c in self.points_exact(n=n, ss=ss, iss=iss)] + +# if n == 0: +# #donne un point au hasard dans l'ensemble limite +# return [0] +# else: +# orbit_points0 = self.points(n-1) +# orbit_points = Set([]) +# for c in C: +# orbit_points = orbit_points.union(Set([place(b)*p+place(c) for p in orbit_points0])) +# return orbit_points + + def plot (self, n=None, place=None, ss=None, iss=None, prec=53, point_size=None, color='blue', verb=False): + r""" + Draw the limit set of the beta-adic monoid (with or without subshift). + + INPUT: + + - ``n`` - integer (default: ``None``) + The number of iterations used to plot the fractal. + Default values: between ``5`` and ``16`` depending on the number of generators. + + - ``place`` - place of the number field of beta (default: ``None``) + The place we should use to evaluate elements of the number field given by points_exact() + + - ``ss`` - Automaton (default: ``None``) + The subshift to associate to the beta-adic monoid for this drawing. + + - ``iss`` - set of initial states of the automaton ss (default: ``None``) + + - ``prec`` - precision of returned values (default: ``53``) + + - ``point_size`` - real (default: ``None``) + Size of the plotted points. + + - ``verb`` - bool (default: ``False``) + Print informations for debugging. + + OUTPUT: + + A Graphics object. + + EXAMPLES: + + #. The dragon fractal:: + + sage: m=BetaAdicMonoid(1/(1+I), {0,1}) + sage: m.plot() # long time + + #. The Rauzy fractal of the Tribonacci substitution:: + + sage: s = WordMorphism('1->12,2->13,3->1') + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: m.plot() # long time + + #. A non-Pisot Rauzy fractal:: + + sage: s = WordMorphism({1:[3,2], 2:[3,3], 3:[4], 4:[1]}) + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: m.b = 1/m.b + sage: m.ss = m.ss.transpose().determinize().minimize() + sage: m.plot() # long time + + #. The dragon fractal and its boundary:: + + sage: m = BetaAdicMonoid(1/(1+I), {0,1}) + sage: p1 = m.plot() # long time + sage: ssi = m.intersection_words(w1=[0], w2=[1]) # long time + sage: p2 = m.plot(ss = ssi, n=18) # long time + sage: p1+p2 # long time + + #. The "Hokkaido" fractal and its boundary:: + + sage: s = WordMorphism('a->ab,b->c,c->d,d->e,e->a') + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: p1 = m.plot() # long time + sage: ssi = m.intersection_words(w1=[0], w2=[1]) # long time + sage: p2 = m.plot(ss=ssi, n=40) # long time + sage: p1+p2 # long time + + #. A limit set that look like a tiling:: + + sage: P=x^4 + x^3 - x + 1 + sage: b = P.roots(ring=QQbar)[2][0] + sage: m = BetaAdicMonoid(b, {0,1}) + sage: m.plot(18) # long time + + """ + + global co + + co = 0 + orbit_points = self.points(n=n, place=place, ss=ss, iss=iss, prec=prec) + if verb: print "co=%s"%co + + # Plot points size + if point_size is None: + point_size = 1 + + # Make graphics + from sage.plot.plot import Graphics + G = Graphics() + + #dim = self.b.minpoly().degree() + + from sage.rings.qqbar import QQbar, AA + if QQbar(self.b) not in AA: #2D plots + from sage.all import points + G = points(orbit_points, size=point_size, color=color) + else: # 1D plots + from sage.all import plot + G += plot(orbit_points, thickness=point_size, color=color) +# if plotbasis: +# from matplotlib import cm +# from sage.plot.arrow import arrow +# canonicalbasis_proj = self.rauzy_fractal_projection(eig=eig, prec=prec) +# for i,a in enumerate(alphabet): +# x = canonicalbasis_proj[a] +# G += arrow((-1.1,0), (-1.1,x[0]), color=cm.__dict__["gist_gray"](0.75*float(i)/float(size_alphabet))[:3]) +# else: +# print "dimension too large !" + G.set_aspect_ratio(1) + + return G + + def _relations_automaton_rec (self, current_state, di, parch, pultra, m, Cd, ext, verb=False, niter=-1): + r""" + Used by relations_automaton() + """ + + if niter == 0: + return di + + global count + if verb and count%10000 == 0: print count + if count == 0: + return di + count -= 1 + + b = self.b + if not di.has_key(current_state): + di[current_state] = dict([]) + for c in Cd: #parcours les transitions partant de current_state + #e = b*current_state + c #calcule l'état obtenu en suivant la transition c + e = (current_state + c)/b #calcule l'état obtenu en suivant la transition c + #if verb: print "b=%s, e=%s, cur=%s, c=%s, di=%s"%(b, e, current_state, c, di) + if not di.has_key(e): #détermine si l'état est déjà dans le dictionnaire + ok = True + #calcule les valeurs abolues pour déterminer si l'état n'est pas trop grand + for p in parch: + if not ext: + if p(e).abs() >= m[p]: + ok = False + break + else: + if p(e).abs() > m[p]+.000000001: + ok = False + break + if not ok: + continue #cesse de considérer cette transition + for p, d in pultra: + if absp(e, p, d) > m[p]: + #if verb: print "abs(%s)=%s trop grand !"%(e, absp(e, p, d)) + ok = False + break + if ok: + #on ajoute l'état et la transition à l'automate + di[current_state][e] = c + di = self._relations_automaton_rec (current_state=e, di=di, parch=parch, pultra=pultra, m=m, Cd=Cd, ext=ext, verb=verb, niter=niter-1) + else: + #ajoute la transition + di[current_state][e] = c + return di + + def relations_automaton (self, ext=False, ss=None, noss=False, verb=False, step=100, limit=None, niter=None): + r""" + Compute the relations automaton of the beta-adic monoid (with or without subshift). + See http://www.latp.univ-mrs.fr/~paul.mercat/Publis/Semi-groupes%20fortement%20automatiques.pdf for a definition of such automaton (without subshift). + + INPUT: + + - ``ext`` - bool (default: ``False``) + If True, compute the extended relations automaton (which permit to describe infinite words in the monoid). + + - ``ss`` - Automaton (default: ``None``) + The subshift to associate to the beta-adic monoid for this operation. + + - ``noss`` - bool (default: ``False``) + + + - ``verb`` - bool (default: ``False``) + If True, print informations for debugging. + + - ``step`` - int (default: ``100``) + Stop to an intermediate state of the computing to verify that all is right. + + - ``limit``- int (default: None) + Stop the computing after a number of states limit. + + OUTPUT: + + A Automaton. + + EXAMPLES:: + + sage: m = BetaAdicMonoid(3, {0,1,3}) + sage: m.relations_automaton() + Finite automaton with 3 states + + sage: b = (x^3-x-1).roots(ring=QQbar)[0][0] + sage: m = BetaAdicMonoid(b, {0,1}) + sage: m.relations_automaton() + Finite automaton with 179 states + + REFERENCES: + + .. [Me13] Mercat P. + Bull. SMF 141, fascicule 3, 2013. + + """ + + if not noss: + a = self.relations_automaton(ext=ext, ss=None, noss=True, verb=verb, step=step, limit=limit) + if not step: + return a + step = step-1 + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + else: + return a #pas de sous-shift + if not step: + return ss + step = step-1 + d=dict() + for u in self.C: + for v in self.C: + if not d.has_key(u-v): + d[u-v] = [] + d[u - v] += [(u,v)] + if not step: + return d + step = step-1 + ss = ss.emonde0_simplify() + P = ss.product(A=ss) + #P = P.emonde0_simplify() + if not step: + return P + step = step-1 + a.relabel2(d) + if not step: + return a + step = step-1 + a = a.intersection(A=P) + if not step: + return a + step = step-1 + a = a.emonde0_simplify() + if not step: + return a + step = step-1 + if not ext: + a.emondeF() + if not step: + return a + step = step-1 + #a = a.determinize(A=a.A, noempty=True) + #if not step: + # return a + #step = step-1 + #return a + return a.minimize() + + K = self.C[0].parent() + b = self.b + + if verb: print K + + #détermine les places qu'il faut considérer + parch = [] + for p in K.places(): #places archimédiennes + if p(b).abs() < 1: + parch += [p] + pi = K.defining_polynomial() + from sage.rings.arith import gcd + pi = pi/gcd(pi.list()) #rend le polynôme à coefficients entiers et de contenu 1 +# den = pi.constant_coefficient().denominator() +# lp = (pi.list()[pi.degree()].numerator()*den).prime_divisors() #liste des nombres premiers concernés + lp = (Integer(pi.list()[0])).prime_divisors() #liste des nombres premiers concernés + pultra = [] #liste des places ultramétriques considérées + for p in lp: + #détermine toutes les places au dessus de p dans le corps de nombres K + k = Qp(p) + Kp = k['a'] + a = Kp.gen() + for f in pi(a).factor(): + kp = f[0].root_field('e') +# c = kp.gen() + if kp == k: + c = f[0].roots(kp)[0][0] + else: + c = kp.gen() + if verb: print "c=%s (abs=%s)"%(c, (c.norm().abs())**(1/f[0].degree())) + if (c.norm().abs())**(1/f[0].degree()) < 1: #absp(c, c, f[0].degree()) > 1: + pultra += [(c, f[0].degree())] + + if verb: print "places: "; print parch; print pultra + + #calcule les bornes max pour chaque valeur absolue + Cd = Set([c-c2 for c in self.C for c2 in self.C]) + if verb: print "Cd = %s"%Cd + m = dict([]) + for p in parch: + m[p] = max([p(c).abs() for c in Cd])/abs(1-p(p.domain().gen()).abs()) + for p, d in pultra: + m[p] = max([absp(c, p, d) for c in Cd]) + + if verb: print "m = %s"%m + + if verb: print K.zero().parent() + + global count + #print limit + if limit is None: + count = -1 + else: + count = limit + if niter is None: + niter = -1 + #print count + if verb: print "Parcours..." + di = self._relations_automaton_rec (current_state=K.zero(), di=dict([]), parch=parch, pultra=pultra, m=m, Cd=Cd, ext=ext, verb=verb, niter=niter) + + if count == 0: + print "Nombre max d'états atteint." + else: + if verb: + if limit is None: + print "%s états parcourus."%(-1-count) + else: + print "%s états parcourus."%(limit-count) + + #a = Automaton([K.zero()], [K.zero()], di) + + #if verb: print "di = %s"%di + + res = Automaton(di, loops=True) #, multiedges=True) + + if verb: print "Avant emondation : %s"%res + + res.I = [K.zero()] + res.A = Set([c-c2 for c in self.C for c2 in self.C]) + if verb: print "Emondation..." + if not ext: + res.F = [K.zero()] + res.emonde() + else: + #res = res.emonde0_simplify() #pour retirer les états puits + res.emonde0() + res.F = res.vertices() + return res + + def relations_automaton2 (self, verb=False, step=100, limit=None, niter=None): + r""" + + Do the same as relations_automaton, but avoid recursivity in order to avoid the crash of sage. + + INPUT: + + - ``verb`` - bool (default: ``False``) + If True, print informations for debugging. + + - ``step`` - int (default: ``100``) + Stop to an intermediate state of the computing to verify that all is right. + + - ``limit``- int (default: None) + Stop the computing after a number of states limit. + + OUTPUT: + + A Automaton. + """ + + K = self.C[0].parent() + b = self.b + + if verb: print K + + #détermine les places qu'il faut considérer + parch = [] + for p in K.places(): #places archimédiennes + if p(b).abs() < 1: + parch += [p] + pi = K.defining_polynomial() + from sage.rings.arith import gcd + pi = pi/gcd(pi.list()) #rend le polynôme à coefficients entiers et de contenu 1 + lp = (Integer(pi.list()[0])).prime_divisors() #liste des nombres premiers concernés + pultra = [] #liste des places ultramétriques considérées + for p in lp: + #détermine toutes les places au dessus de p dans le corps de nombres K + k = Qp(p) + Kp = k['a'] + a = Kp.gen() + for f in pi(a).factor(): + kp = f[0].root_field('e') + if kp == k: + c = f[0].roots(kp)[0][0] + else: + c = kp.gen() + if verb: print "c=%s (abs=%s)"%(c, (c.norm().abs())**(1/f[0].degree())) + if (c.norm().abs())**(1/f[0].degree()) < 1: #absp(c, c, f[0].degree()) > 1: + pultra += [(c, f[0].degree())] + + if verb: print "places: "; print parch; print pultra + + #calcule les bornes max pour chaque valeur absolue + Cd = Set([c-c2 for c in self.C for c2 in self.C]) + if verb: print "Cd = %s"%Cd + m = dict([]) + for p in parch: + m[p] = max([p(c).abs() for c in Cd])/abs(1-p(p.domain().gen()).abs()) + for p, d in pultra: + m[p] = max([absp(c, p, d) for c in Cd]) + + if verb: print "m = %s"%m + + if verb: print K.zero().parent() + + if limit is None: + count = -1 + else: + count = limit + if niter is None: + niter = -1 + + if verb: print "Parcours..." + + di=dict([]) + S = [K.zero()] #set of states to look at + iter = 0 + while len(S) != 0: + if iter == niter: + break + for s in S: + S.remove(s) + if not di.has_key(s): + di[s] = dict([]) + if count%10000 == 0: + print count + count-=1 + if count == 0: + iter = niter-1 #to break the main loop + break + for c in Cd: #parcours les transitions partant de current_state + e = (s + c)/b #calcule l'état obtenu en suivant la transition c + #if verb: print "b=%s, e=%s, cur=%s, c=%s, di=%s"%(b, e, current_state, c, di) + if not di.has_key(e): #détermine si l'état est déjà dans le dictionnaire + ok = True + #calcule les valeurs abolues pour déterminer si l'état n'est pas trop grand + for p in parch: + if p(e).abs() >= m[p]: + ok = False + break + if not ok: + continue #cesse de considérer cette transition + for p, d in pultra: + if absp(e, p, d) > m[p]: + #if verb: print "abs(%s)=%s trop grand !"%(e, absp(e, p, d)) + ok = False + break + if ok: + #on ajoute l'état et la transition à l'automate + di[s][e] = c + S.append(e) + else: + #ajoute la transition + di[s][e] = c + iter+=1 + + if count == 0: + print "Nombre max d'états atteint." + return + else: + if verb: + if limit is None: + print "%s états parcourus."%(-1-count) + else: + print "%s états parcourus."%(limit-count) + + res = Automaton(di, loops=True) #, multiedges=True) + + if verb: print "Avant emondation : %s"%res + + res.I = [K.zero()] + res.A = Set([c-c2 for c in self.C for c2 in self.C]) + res.F = [K.zero()] + if verb: print "Emondation..." + res.emonde() + return res + + def critical_exponent_aprox (self, niter=10, verb=False): + b = self.b + K = b.parent() + C = self.C + S = set([K.zero()]) + for i in range(niter): + S2 = set([]) + for s in S: + for c in C: + S2.add((s+c)/b) + #intervertit S et S2 + S3 = S2 + S2 = S + S = S3 + if verb: print len(S) + from sage.functions.log import log + print "%s"%(log(len(S)).n()/(niter*abs(log(abs(b.n()))))) + + def complexity (self, verb = False): + r""" + Return a estimation of an upper bound of the number of states of the relations automaton. + + INPUT: + + - ``verb`` - Boolean (default: False) Display informations for debug. + + OUTPUT: + + A positive real number. + + EXAMPLES:: + + sage: m = BetaAdicMonoid(3, {0,1,3}) + sage: m.complexity() + 7.06858347... + """ + K = self.C[0].parent() + b = self.b + + if verb: print K + + #détermine les places qu'il faut considérer + parch = K.places() + r = len(parch) + pi = K.defining_polynomial() + from sage.rings.arith import gcd + pi = pi/gcd(pi.list()) #rend le polynôme à coefficients entiers et de contenu 1 +# den = pi.constant_coefficient().denominator() +# lp = (pi.list()[pi.degree()].numerator()*den).prime_divisors() #liste des nombres premiers concernés + lp = (Integer(pi.list()[0])).prime_divisors() #liste des nombres premiers concernés + pultra = [] #liste des places ultramétriques considérées + for p in lp: + #détermine toutes les places au dessus de p dans le corps de nombres K + k = Qp(p) + Kp = k['a'] + a = Kp.gen() + for f in pi(a).factor(): + kp = f[0].root_field('e') +# c = kp.gen() + if kp == k: + c = f[0].roots(kp)[0][0] + else: + c = kp.gen() + if verb: print "c=%s (abs=%s)"%(c, (c.norm().abs())**(1/f[0].degree())) + if c.norm().abs() != 1: #absp(c, c, f[0].degree()) > 1: + pultra += [(c, f[0].degree())] + + if verb: print "places: "; print parch; print pultra + + #calcule les bornes max pour chaque valeur absolue + Cd = Set([c-c2 for c in self.C for c2 in self.C]) + if verb: print "Cd = %s"%Cd + vol = 1. + #from sage.rings.real_mpfr import RR + for p in parch: + if (p(b)).imag() == 0: + vol *= 2*max([p(c).abs() for c in Cd])/abs(1-p(b).abs()) + if verb: print "place réelle %s"%p + else: + vol *= 3.1415926535*(max([p(c).abs() for c in Cd])/abs(1-p(b).abs()))**2 + if verb: print "place complexe %s"%p + #vol *= max([p(c).abs() for c in Cd])/abs(1-p(p.domain().gen()).abs()) + #vol *= max(1, max([p(c).abs() for c in Cd])/abs(1-p(p.domain().gen()).abs())) + for p, d in pultra: + vol *= max([(c.polynomial())(p).norm().abs() for c in Cd]) + #vol *= max([absp(c, p, d) for c in Cd]) + #vol *= max(1, max([absp(c, p, d) for c in Cd])) + #from sage.functions.other import sqrt + #return vol/(K.regulator()*(sqrt(r+1).n())) + return vol + + #def infinite_relations_automaton (self, verb=False): + # a = self.relations_automaton (ext=True, verb=verb) + # #retire la composante fortement connexe de l'etat initial + # K = self.b.parent() + # for s in a.strongly_connected_components_subgraphs(): + # if K.zero() in s: + # a.delete_vertices(a.strongly_connected_components_subgraphs()[0].vertices()) + # return a + + def intersection (self, ss, ss2=None, Iss=None, Iss2=None, ext=True, verb = False): #calcule le sous-shift correspondant à l'intersection des deux monoïdes avec sous-shifts + r""" + Compute the intersection of two beta-adic monoid with subshifts + + INPUT: + + - ``ss``- Automaton (default: ``None``) + The first subshift to associate to the beta-adic monoid for this operation. + + - ``ss2``- Automaton (default: ``None``) + The second subshift to associate to the beta-adic monoid for this operation. + + - ``Iss``- set of states of ss (default: ``None``) + + - ``Iss2``- set of states of ss2 (default: ``None``) + + - ``ext`` - bool (default: ``True``) + If True, compute the extended relations automaton (which permit to describe infinite words in the monoid). + + - ``verb``- bool (default: ``False``) + If True, print informations for debugging. + + OUTPUT: + + A Automaton. + + EXAMPLES: + + #. Compute the boundary of the dragon fractal (see intersection_words for a easier way) :: + + sage: m = BetaAdicMonoid(1/(1+I), {0,1}) + sage: ss=m.default_ss() + sage: iss=ss.I[0] + sage: ss0 = ss.prefix(w=[0], i=iss) + sage: ss1 = ss.prefix(w=[1], i=iss) + sage: ssi = m.intersection(ss=ss0, ss2=ss1) + sage: ssd = ssi.determinize(A=m.C, noempty=True) + sage: ssd = ssd.emonde0_simplify() + sage: m.plot(ss = ssd, n=19) # long time + """ + + m = None + + if ss2 is None: + if hasattr(self, 'ss'): + ss2 = self.ss + else: + raise ValueError("Only one sub-shift given, I need two !") + if type(ss) == type(BetaAdicMonoid(2,{0,1})): + m = ss + ss = m.ss + if m.b != self.b: + raise ValueError("Cannot compute the intersection of two beta-adic monoids with differents beta.") + m.C = m.C.union(self.C) + self.C = self.C.union(m.C) + if hasattr(m, 'ss'): + m.ss.A = m.C + else: + raise ValueError("Only one sub-shift given, I need two !") + self.ss.A = self.C + + if Iss is None: + if hasattr(ss, 'I'): + Iss = ss.I + if Iss is None: + Iss = [ss.vertices()[0]] + if Iss2 is None: + if hasattr(ss2, 'I'): + Iss2 = ss2.I + if Iss2 is None: + Iss2 = [ss2.vertices()[0]] + if verb: print "Iss = %s, Iss2 = %s"%(Iss, Iss2) + + a = ss.product(ss2) + if verb: print "Product = %s"%a + + ar = self.relations_automaton(ext=ext, noss=True) + if verb: print "Arel = %s"%ar + + #maps actual edges to the list of corresponding couple + m = dict([]) + for c in self.C: + for c2 in self.C: + if m.has_key(c-c2): + m[c-c2] += [(c,c2)] + else: + m[c-c2] = [(c,c2)] + if verb: print "m = %s"%m + + L = a.Alphabet() #a.edge_labels() + LA = ar.Alphabet() #ar.edge_labels() + d = dict([]) + for u, v in L: + for ka in LA: + for u2, v2 in m[ka]: + if u == u2 and v == v2: + d[((u,v), ka)] = u + break + else: + d[((u,v), ka)] = None + if verb: print "d = %s"%d + p = a.product(A=ar, d=d) + #I = [((i,i2),self.b.parent().zero()) for i in Iss for i2 in Iss2] + #if verb: print "I = %s"%I + #p.emondeI(I=I) + #if verb: print "%s"%p + #p.emonde0(I=I) + p.I = [((i,i2), self.b.parent().zero()) for i,i2 in zip(Iss, Iss2)] + p = p.emonde0_simplify() + if m is not None: + ssd = p.determinize2(A=self.C, noempty=True) + ssd = ssd.emonde0_simplify() + return ssd + return p + + def intersection_words (self, w1, w2, ss=None, iss=None): + r""" + Compute the intersection of two beta-adic monoid with subshifts corresponding to two prefix + + INPUT: + + - ``w1``- word + The first prefix. + + - ``w2``- word + The second prefix. + + - ``ss``- Automaton (default: ``None``) + The subshift to associate to the beta-adic monoid for this operation. + + - ``iss``- initial state of ss (default: ``None``) + + OUTPUT: + + A Automaton. + + EXAMPLES: + + #. Compute the boundary of the dragon fractal:: + + sage: m = BetaAdicMonoid(1/(1+I), {0,1}) + sage: m.intersection_words(w1=[0], w2=[1]) + Finite automaton with 21 states + + #. Draw the intersection of two sub-sets of a limit set:: + + sage: m = BetaAdicMonoid(1/(1+I), {0,1}) + sage: ssi = m.intersection_words(w1=[0], w2=[1]) + sage: m.plot(n=10, ss=ssi) # long time + """ + + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + else: + ss = self.default_ss() + if iss is None: + if hasattr(ss, 'I'): + iss = ss.I[0] + if iss is None: + iss = ss.vertices()[0] + ss1 = ss.prefix(w=w1, i=iss) + ss2 = ss.prefix(w=w2, i=iss) + ssi = self.intersection(ss=ss1, ss2=ss2) + ssd = ssi.determinize2(A=self.C, noempty=True) + ssd = ssd.emonde0_simplify() + return ssd + + def reduced_words_automaton (self, ss=None, Iss=None, ext=False, verb=False, step=None, arel=None): + r""" + Compute the reduced words automaton of the beta-adic monoid (with or without subshift). + See http://www.latp.univ-mrs.fr/~paul.mercat/Publis/Semi-groupes%20fortement%20automatiques.pdf for a definition of such automaton (without subshift). + + WARNING: It seems there is a bug : result may be incorrect if ss is not None. + + INPUT: + + - ``ss``- Automaton (default: ``None``) + The first subshift to associate to the beta-adic monoid for this operation. + + - ``Iss``- set of states of ss (default: ``None``) + + - ``ext`` - bool (default: ``True``) + If True, compute the extended relations automaton (which permit to describe infinite words in the monoid). + + - ``verb`` - bool (default: ``False``) + If True, print informations for debugging. + + - ``step`` - int (default: ``None``) + Stop to a intermediate state of the computing to make verifications. + + - ``arel`` - Automaton (default: ``None``) + Automaton of relations. + + OUTPUT: + + A Automaton. + + EXAMPLES: + + #. 3-adic expansion with numerals set {0,1,3}:: + + sage: m = BetaAdicMonoid(3, {0,1,3}) + sage: m.reduced_words_automaton() + Finite automaton with 2 states + + #. phi-adic expansion with numerals set {0,1}:: + + sage: m = BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + sage: m.reduced_words_automaton() + Finite automaton with 3 states + + #. beta-adic expansion with numerals set {0,1} where beta is the plastic number:: + sage: b = (x^3-x-1).roots(ring=QQbar)[0][0] + sage: m = BetaAdicMonoid(b, {0,1}) + sage: m.reduced_words_automaton() # long time + Finite automaton with 5321 states + """ + + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + if hasattr(self.ss, 'I'): + Iss = self.ss.I + + if step is None: + step = 1000 + + K = self.C[0].parent() + + if verb: print "Calcul de l'automate des relations..." + + if arel is None: + a = self.relations_automaton(noss=True) + else: + a = arel + + if verb: print " -> %s"%a + + if step == 1: + return ("automate des relations", a) + + #add a state copy of K.0 (it will be the new initial state) + a.add_vertex('O') + +# #add transitions to K.0 to 'O' +# for f, d, l in a.incoming_edges(K.zero(), labels=True): +# if f == K.zero(): +# a.add_edge('O', 'O', l) +# else: +# a.add_edge(f, 'O', l) + + #subset of positives labels + Cdp = [] + for i in range(self.C.cardinality()): + for j in range(i): + Cdp += [self.C[i] - self.C[j]] + + #redirect positives transitions from K.0 + for f, d, l in a.outgoing_edges(K.zero(), labels=True): + if l in Cdp: +# a.delete_edge(K.zero(), d, l) + #add the edge + a.add_edge('O', d, l) + + a.add_edge('O', 'O', a.edge_label(K.zero(), K.zero())) + + #if verb: print a.incoming_edges(K.zero(), labels=True) + + #remove outgoing edges from K.0 (except from K.0 to K.0) + for f, d, l in a.outgoing_edges(K.zero(), labels=True): + if f != d: + a.delete_edge(f, d, l) + + if step == 2: + return ("automate des relations ordonnées", a) + + a.emondeI(I=['O']) + + if step == 3: + return ("automate des relations ordonnées, émondé", a) + + if ss is not None: #not full sub-shift + if Iss is None: + Iss = [ss.vertices()[0]] + #maps actual edges to the list of corresponding couple + m = dict([]) + for c in self.C: + for c2 in self.C: + if m.has_key(c-c2): + m[c-c2] += [(c,c2)] + else: + m[c-c2] = [(c,c2)] + #if verb: print "m=%s"%m + + #calculate the 'product to the right' of a with ss + d = dict([]) + La = a.edge_labels() + Lss = ss.edge_labels() + for ka in La: + for kss in Lss: + d[(ka,kss)] = None + for k in m[ka]: + if k[1] == kss: + d[(ka,kss)] = k[0] + break + + #if verb: print "d=%s"%d + if verb: print "avant produit : a=%s (%s etats)"%(a, a.num_verts()) + a = a.product(A=ss, d=d) + if verb: print "après produit : a=%s"%a + if step == 4: + return ("automate des mots généraux non réduits", a) + + I = [('O',i) for i in Iss] + nof = Set([(K.zero(),i) for i in ss.vertices()]) + + #if verb: print "I=%s, F=%s"%(I, nof) + + if ext: + #a.emondeI(I=I) + #a.emonde0(I=I) #pour retirer les états puits + a = a.emonde0_simplify(I=I) + else: + a = a.emonde0_simplify(I=I) + a.emonde(I=I, F=nof) + #a.emondeI(I=I) + #a.emondeF(F=nof) + #if step == 4: + # return ("automate des mots généraux non réduits, émondé", a) + #a.emondeF(F=nof) + + if verb: print "après émondation : a=%s"%a + if step == 5: + return ("automate des mots généraux non réduits, émondé", a) + + #return a + else: + #maps actual edges to element of self.C (the left part when writted c-c2) + m = dict([]) + for c in self.C: + for c2 in self.C: + if m.has_key(c-c2): + m[c-c2] += [c] + else: + m[c-c2] = [c] + #if verb: print "m=%s"%m + + a.allow_multiple_edges(True) + #replace each label by its mapping + for f, d, l in a.edges(): + a.delete_edge(f, d, l) + for l2 in m[l]: + a.add_edge(f, d, l2) + + I = ['O'] + nof=Set([K.zero()]) + + a.I = I + a.F = nof + a.C = self.C + + if verb: print "avant determinisation : a=%s"%a + if step == 6: + return ("automate des mots généraux non réduits, émondé", a) + + #rend l'automate plus simple + a = a.emonde0_simplify() + + if verb: print "simplification : a=%s"%a + + if verb: print "Determinization..." + #determinize + ad = a.determinize2(nof=a.F) + #ad = a.determinize(nof=a.F, verb=False) + #ad = a.determinize(I, self.C, nof, verb=verb) + + if verb: print " -> %s"%ad + if step == 7: + return ("automate des mots généraux réduits", ad) + + if ss is not None: #not full sub-shift + #calculate the intersection with ss + ad = ad.emonde0_simplify() + ad = ad.intersection(ss) + if verb: print "apres intersection : a=%s"%ad + + if step == 8: + return ("automate des mots réduits", ad) + + #F2=[e for e in a.vertices() nof in e[0]] + #if verb: print "I2=%s"%I2 #, F2=%s"%(I2,F2) + ad.A = self.C + #ad.emondeI(I=I2) #, F=F2) + ad = ad.emonde0_simplify() + ad.F = ad.vertices() + + if verb: print "apres émondation : a=%s"%ad + if step == 9: + return ("automate des mots réduits, émondé", ad) + + return ad + + def critical_exponent_free (self, ss=None, prec = None, verb=False): + r""" + Compute the critical exponent of the beta-adic monoid (with or without subshift), assuming it is free. + See http://www.latp.univ-mrs.fr/~paul.mercat/Publis/Semi-groupes%20fortement%20automatiques.pdf for a definition (without subshift). + + INPUT: + + - ``ss``- Automaton (default: ``None``) + The first subshift to associate to the beta-adic monoid for this operation. + + - ``prec``- precision (default: ``None``) + + - ``verb``- bool (default: ``False``) + If True, print informations for debugging. + + OUTPUT: + + A real number. + + EXAMPLES: + + #. Hausdorff dimension of limit set of 3-adic expansion with numerals set {0,1,3}:: + + sage: m = BetaAdicMonoid(3, {0,1,3}) + sage: m.critical_exponent_free(m.reduced_words_automaton()) + log(y)/log(|3|) where y is the max root of x^2 - 3*x + 1 + 0.8760357589... + + #. Hausdorff dimension of limit set of phi-adic expansion with numerals set {0,1}:: + + sage: m = BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + sage: m.critical_exponent_free(m.reduced_words_automaton()) + log(y)/log(|b|) where y is the max root of x^2 - x - 1 + 1.0000000000... + + #. Hausdorff dimension of the boundary of the dragon fractal:: + + sage: m = BetaAdicMonoid(1/(1+I), {0,1}) + sage: ssi = m.intersection_words(w1=[0], w2=[1]) + sage: m.critical_exponent_free(ss=ssi) + log(y)/log(|b|) where y is the max root of x^3 - x^2 - 2 + 1.5236270862... + + #. Hausdorff dimension of the boundary of a Rauzy fractal:: + + sage: s = WordMorphism('1->12,2->13,3->1') + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: ssi = m.intersection_words(w1=[0], w2=[1]) + sage: m.critical_exponent_free(ss=ssi) + log(y)/log(|b|) where y is the max root of x^4 - 2*x - 1 + 1.0933641642... + + #. Hausdorff dimension of a non-Pisot Rauzy fractal:: + + sage: s = WordMorphism({1:[3,2], 2:[3,3], 3:[4], 4:[1]}) + sage: m = s.rauzy_fractal_beta_adic_monoid() + sage: m.b = 1/m.b + sage: m.ss = m.ss.transpose().determinize().minimize() + sage: m.critical_exponent_free() + log(y)/log(|1/2*b^2 - 1/2*b + 1/2|) where y is the max root of x^3 - x^2 + x - 2 + 1.5485260383... + """ + + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + if ss is None: + y = self.C.cardinality() + print "log(%s)/log(|%s|)"%(y, self.b) + else: + if verb: print "" + M = ss.adjacency_matrix() + if verb: print "Valeurs propres..." + e = M.eigenvalues() + if verb: print "max..." + y = max(e, key=abs) + if verb: print "" + print "log(y)/log(|%s|) where y is the max root of %s"%(self.b, QQbar(y).minpoly()) + y = y.N(prec) + from sage.functions.log import log + b = self.b.N(prec) + if verb: print "y=%s, b=%s"%(y, b) + return abs(log(y)/log(abs(b))).N() + + def critical_exponent (self, ss=None, prec = None, verb=False): + r""" + Compute the critical exponent of the beta-adic monoid (with or without subshift). + See http://www.latp.univ-mrs.fr/~paul.mercat/Publis/Semi-groupes%20fortement%20automatiques.pdf for a definition (without subshift). + + INPUT: + + - ``ss``- Automaton (default: ``None``) + The first subshift to associate to the beta-adic monoid for this operation. + + - ``prec``- precision (default: ``None``) + + - ``verb``- bool (default: ``False``) + If True, print informations for debugging. + + OUTPUT: + + A real number. + + EXAMPLES: + + #. Hausdorff dimension of limit set of 3-adic expansion with numerals set {0,1,3}:: + + sage: m = BetaAdicMonoid(3, {0,1,3}) + sage: m.critical_exponent() + log(y)/log(|3|) where y is the max root of x^2 - 3*x + 1 + 0.8760357589... + + #. Hausdorff dimension of limit set of phi-adic expansion with numerals set {0,1}:: + + sage: m = BetaAdicMonoid((1+sqrt(5))/2, {0,1}) + sage: m.critical_exponent() + log(y)/log(|b|) where y is the max root of x^2 - x - 1 + 1.0000000000... + + #. Critical exponent that is not the Hausdorff dimension of the limit set:: + + sage: P = x^7 - 2*x^6 + x^3 - 2*x^2 + 2*x - 1 + sage: b = P.roots(ring=QQbar)[3][0] + sage: m = BetaAdicMonoid(b, {0,1}) + sage: m.critical_exponent() # long time + log(y)/log(|b|) where y is the max root of x^11 - 2*x^10 - 4*x^2 + 8*x + 2 + 3.3994454205... + + #. See more examples with doc critical_exponent_free() + + """ + # #. Critical exponent that is not the Hausdorff dimension of the limit set:: + # + # sage: + + if ss is None: + if hasattr(self, 'ss'): + ss = self.ss + if verb: + print "Calcul de l'automate des mots réduits...\n" + a = self.reduced_words_automaton(ss=ss, verb=verb) + return self.critical_exponent_free (prec=prec, ss=a, verb=verb) + + #test if 0 is an inner point of the limit set + def ZeroInner (self, verb=False): + + if not hasattr(self, 'ss'): + self.ss = self.default_ss() + + if verb: print "relations automaton..." + + ar = self.relations_automaton(ext=True) + ar.complementary() + + if verb: print "complementaire : %s"%ar + #return ar + + a = self.default_ss().product(self.ss) + + if verb: print "a = %s"%a + #return a + + #a = ar.intersection(a) + + if verb: print "product..." + + L = a.edge_labels() + Lr = ar.edge_labels() + d = dict([]) + for k in L: + for kr in Lr: + if k == kr and k[0] == 0: + d[(k, kr)] = 0 + else: + d[(k, kr)] = None + a2 = a.product(A=ar, d=d) + + if verb: print "product = %s"%a2 + #return a2 + + #test if there is a cycle in the graph + if a2.is_directed_acyclic(): + print "0 is not an inner point." + else: + ss = self.ss + self.ss = None + print "Zero is an inner point iff the %s has non-empty interior."%self + self.ss = ss + + + + + + + \ No newline at end of file