# load the package badic (you can install it with "sage -pip install badic [--user]")
from badic import *

# define a beta-adic set, with beta the Tribonnacci number
m = BetaAdicSet(x^3-x^2-x-1, [0,1])
m = m.reduced()  # describe each point of the set by an unique word
print(m)
m.plot()

b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 4 states and 2 letters

# define an automaton describing a Sierpinsky gasket
a = DetAutomaton([(0,2,0),(0,6,1),(2,3,1),(2,12,0),(6,7,1),(6,9,0),(3,4,1),(3,5,0),(12,13,1),(12,14,0),(7,8,0),(7,15,1),(9,10,0),(9,11,1),(4,0,0),(5,0,0),(5,0,1),(13,0,0),(13,0,1),(14,0,0),(8,0,0),(8,0,1),(15,0,1),(10,0,1),(11,0,1),(11,0,0)],i=0, final_states=[0])
a.plot()

# automaton describing words starting by 00
a2 = dag.Word([0,0]).concat(dag.AnyWord([0,1]))
a2.plot()

# concatenate word 10 to words described by a
a3 = dag.Word([1,0]).concat(a)
a3.plot()

# do the union
a = a2.union(a3)
a

DetAutomaton with 13 states and an alphabet of 2 letters

# take the beta-adic set for this union
ms = BetaAdicSet(m.b, a)
ms = ms.proj(m) # describe the same set with a sub-language of the one of m
ms.plot()

# compute the domain exchange describing this beta-adic set
lm = ms.domain_exchange()
lm

[(b - 1,
  b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 97 states and 2 letters),
 (1,
  b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 91 states and 2 letters),
 (b,
  b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 195 states and 2 letters),
 (b^2 - 1,
  b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 98 states and 2 letters),
 (b + 1,
  b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 98 states and 2 letters),
 (b^2,
  b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 106 states and 2 letters)]

# the domain exchange before exchange
ms.plot_list([a for t,a in lm])

# domain exchange after exchange
ms.plot_list([a.proj(ms, t) for t,a in lm])

# compute a substitution describing this beta-adic set
%time d = ms.substitution()
d

CPU times: user 33.4 s, sys: 310 ms, total: 33.7 s
Wall time: 33.7 s

{1: [59, 6],
 2: [19],
 3: [19, 54],
 4: [50, 42],
 5: [57, 9, 58, 3],
 6: [59, 6, 40, 48],
 7: [59, 6, 53],
 8: [21, 35, 48, 59, 1],
 9: [19, 55, 5],
 10: [21, 65, 49, 59, 1],
 11: [63, 6, 15, 5],
 12: [59, 6, 60, 49, 59, 1],
 13: [53, 63, 7, 25, 4],
 14: [54, 20, 33, 4],
 15: [59, 18, 38, 3, 37, 46, 58, 2],
 16: [36, 17, 45, 41, 46, 58, 2],
 17: [63, 6, 53, 5],
 18: [59, 6, 53, 63, 1],
 19: [57, 9, 58, 3, 37, 46],
 20: [57, 9, 58, 3, 52],
 21: [34, 11, 58, 3, 37, 46],
 22: [34, 11, 58, 3, 52],
 23: [52, 41, 3, 52, 4],
 24: [63, 18, 43, 41, 46, 58, 2],
 25: [63, 18, 43, 50, 4],
 26: [57, 9, 58, 3, 37, 46, 58, 2],
 27: [57, 9, 58, 3, 52, 41, 2],
 28: [40, 48, 59, 7, 64, 47, 58, 2],
 29: [35, 48, 22, 62, 47, 58, 2],
 30: [34, 11, 58, 3, 37, 46, 58, 2],
 31: [34, 11, 58, 3, 52, 41, 2],
 32: [41, 46, 45, 41, 46, 58, 2],
 33: [41, 46, 45, 50, 4],
 34: [15],
 35: [16],
 36: [24],
 37: [26],
 38: [28],
 39: [29],
 40: [30],
 41: [32],
 42: [50, 42, 50],
 43: [13, 42, 50],
 44: [14, 42, 50],
 45: [23, 42, 50],
 46: [19, 54, 5],
 47: [21, 54, 5],
 48: [63, 6, 40, 48, 59, 1],
 49: [59, 6, 40, 48, 59, 1],
 50: [50, 42, 50, 4],
 51: [23, 42, 50, 4],
 52: [27, 44, 50, 4],
 53: [31, 39, 3, 37, 46, 58, 2],
 54: [51, 41, 3, 37, 46, 58, 2],
 55: [58, 46, 58, 3, 37, 46, 58, 2],
 56: [37, 46, 58, 3, 37, 46, 58, 2],
 57: [55],
 58: [56],
 59: [61, 12],
 60: [59, 49],
 61: [60, 49],
 62: [62, 8],
 63: [67, 10],
 64: [63, 49],
 65: [68, 8],
 66: [64, 8],
 67: [65, 49],
 68: [66, 49]}

# convert it to a WordMorphism
s = WordMorphism(d)
s

WordMorphism: 1->59,6, 10->21,65,49,59,1, 11->63,6,15,5, 12->59,6,60,49,59,1, 13->53,63,7,25,4, 14->54,20,33,4, 15->59,18,38,3,37,46,58,2, 16->36,17,45,41,46,58,2, 17->63,6,53,5, 18->59,6,53,63,1, 19->57,9,58,3,37,46, 2->19, 20->57,9,58,3,52, 21->34,11,58,3,37,46, 22->34,11,58,3,52, 23->52,41,3,52,4, 24->63,18,43,41,46,58,2, 25->63,18,43,50,4, 26->57,9,58,3,37,46,58,2, 27->57,9,58,3,52,41,2, 28->40,48,59,7,64,47,58,2, 29->35,48,22,62,47,58,2, 3->19,54, 30->34,11,58,3,37,46,58,2, 31->34,11,58,3,52,41,2, 32->41,46,45,41,46,58,2, 33->41,46,45,50,4, 34->15, 35->16, 36->24, 37->26, 38->28, 39->29, 4->50,42, 40->30, 41->32, 42->50,42,50, 43->13,42,50, 44->14,42,50, 45->23,42,50, 46->19,54,5, 47->21,54,5, 48->63,6,40,48,59,1, 49->59,6,40,48,59,1, 5->57,9,58,3, 50->50,42,50,4, 51->23,42,50,4, 52->27,44,50,4, 53->31,39,3,37,46,58,2, 54->51,41,3,37,46,58,2, 55->58,46,58,3,37,46,58,2, 56->37,46,58,3,37,46,58,2, 57->55, 58->56, 59->61,12, 6->59,6,40,48, 60->59,49, 61->60,49, 62->62,8, 63->67,10, 64->63,49, 65->68,8, 66->64,8, 67->65,49, 68->66,49, 7->59,6,53, 8->21,35,48,59,1, 9->19,55,5

# check that the Rauzy fractal correspond to the wanted beta-adic set
s.rauzy_fractal_plot()