BetaAdicSet in Sage¶

Author: Paul Mercat, and Dominique Benielli Labex Archimede - I2M - AMU Aix-Marseille Université

The present tool BetaAdicSet is a convenient tool to describe sets like Rauzy fractals and quasicrystals associated to substitutions, and more generally sets described by $\beta$-expansion with list of digits in a regular language.

Introduction and Definitions¶

Subtitution¶

Given a substitution, let’s say for example the Fibonnacci substitution:

\[\begin{split}\left\{ \begin{array}{rcl} a & \mapsto & ab \\ b & \mapsto & a \end{array} \right.\end{split}\]

Up to replace the substitution by a power, a fixed point always exists. Here we have the following unique right-infinite fixed point:

$abaababaabaababaababaabaababaabaababaaba...$

This infinite sequence has a lot of interesting properties in general. If we replace letters by intervals of convenient lengths, we get a self-similar tiling of $\mathbf R_+$. Convenient lengths are given by the non-negative coefficients of a Perron eigenvector of the incidence matrix $M$ of the substitution. Here this matrix is

\[\begin{split}M = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}\end{split}\]

and $\begin{pmatrix} 1 \\ \beta-1 \end{pmatrix}$ is an eigenvector for the Perron eigenvalue $\beta$, which is the golden number.

Self-similar tiling of $\mathbf R_+$¶

The set of points that appears at the boundaries of intervals in this self-similar tiling are elements of $\mathbf Q(\beta)$, and this set of points have very strong properties, because it is a non-periodic self-similar Meyer set. Here, we get the points

$0, 1,\ \beta,\ \beta + 1,\ \beta + 2,\ 2\beta + 1,\ 2\beta + 2,\ 3\beta + 1,\ 3\beta + 2,\ 3\beta + 3, \ 4\beta + 2,\ 4\beta + 3, ...$

To get this set of points, we start from $0$ and read the fixed point: we add $1$ each time we read a letter $a$ and we add $\beta-1$ each time we read a letter $b$.

Automaton¶

This set of points of $\mathbf Q(\beta)$, associated to the Fibonnacci substitution, is described by the automaton $\mathcal{A}$ of the following figure (and such automaton always exists for sets of points coming from substitutions).

The initial state of this automaton $\mathcal{A}$ is the one labelled by letter ‘a’, and the two states ‘a’ and ‘b’ are final states. The language $L_{\mathcal{A}}$ recognized by an automaton $\mathcal{A}$ is defined as the set of sequences of labels of paths from the initial state to a final state. Here, the language of $\mathcal{A}$ is exactly the set of words over the alphabet $\{0,1\}$ that does not contain the subword $11$.

The set of points is obtained from $\beta$ and from the automaton $\mathcal{A}$ by:

\[\{\sum_{i=0}^n a_i \beta^i\ |\ a_n a_{n-1} ... a_1 a_0 \in L_{\mathcal{A}}\}.\]

The self-similar tiling of $\mathbf R_+$ is obtained by taking the embedding of $\mathbf Q(\beta)$ in $\mathbf R$ corresponding to the Perron eigenvalue of the incidence matrix $M$. If we look at the others embeddings, corresponding to eigenvalues less than $1$ we get a bounded set whose adherence is called Rauzy fractal of the substitution. Here there is a unic such embedding, which is a real one, corresponding to the root of $x^2-x-1$ between $-1$ and $1$. The Rauzy fractal is here the interval $[-1, \varphi]$ of $\mathbf R$ where $\varphi$ is the golden number (i.e. greatest root of $x^2-x-1$).

A BetaAdicSet is by definition a set of this form, for a number $\beta$ and an automaton $\mathcal{A}$. In Sage, a BetaAdicSet is the data of a number $b$ (for $\beta$) and a DetAutomaton $\mathcal{a}$ (but we can initialize it from an Automaton or a CAutomaton).

Main Definitions¶

Broken line¶

Consider $s$ a substitution, or in others words, a word morphism over a finite alphabet $A = \{1, 2, ..., n\}$. Up to replace $s$ by a power, we can assume that $s$ has a fixed point $\omega$.

We defined the defi{broken line} associated to $\omega$ as the subset of $\mathbf Z^n$ defined by

\[\begin{split}\{ {\begin{pmatrix} \text{number of occurences of } 1 \text{ in } \omega_k \\ \text{number of occurences of } 2 \text{ in } \omega_k \\ \vdots \\ \text{number of occurences of } n \text{ in } \omega_k \end{pmatrix} \in \mathbf Z^n } {k \in \mathbf N} \}\end{split}\]

where $\omega_k$ is the prefix of length $k$ of the infinite word $\omega$.

This broken line is very interesting since it is a geometrical object which completely encode the substitution and is stable by multiplication by the incidence matrix. And we get the quasicrystal or the Rauzy fractal by projecting this set of point.

Rauzy fractal¶

The Rauzy fractal is the closure of the projection of the broken line to the contracting space along the expanding line.

Expanding line and contracting space¶

The “expanding line” has dimension $1$ for Pisot numbers, but it can have greater dimension for other Perron numbers.

Let $M_s$ be the incidence matrix of the substitution $s$. By definition the coefficient $(i,j)$ of this matrix is the number of occurrences of the letter $j$ in the word $s(i)$. By Perron-Frobenius theorem, there exists an eigenvector $v \in (\mathbf R_+)^n$, unique if the matrix is irreducible, for an eigenvalue $\lambda$ which is the spectral radius of $M_s$, and moreover we can assume that $v \in (\mathbf Q(\lambda))^{n}$.

We can directly define the projection of the broken line in $\mathbf Q(\lambda)$, by the following.

\[\{ Q_\omega = \sum_{k=1}^N v_{a_k}\ |\ N \in \mathbf N,\ a_1 a_2 ... a_N \text{ prefix of } \omega \text{ of length } N \}.\]

This set is invariant by multiplication by the Perron eigenvalue $\lambda$ and gives a self-similar tiling of $\mathbf R_+$. The definition of $Q_\omega$ depends of the choice of an eigenvector. We prefer to choose an eigenvector whose coefficients belongs to the integer ring $\mathcal O_\lambda$, in order to have $Q_\omega \subset \mathcal O_\lambda$.

For $\mathbf Q(\lambda)$, there are natural contracting and expanding spaces for the multiplication by $\lambda$. Indeed, consider the biggest sets $P_+$ and $P_-$ of places (i.e. equivalence class of absolute values) of $\mathbf Q(\lambda)$ such that

\[\forall v \in P_+,\ |{\lambda}|_v > 1 \quad \text{ and } \quad \forall v \in P_-,\ |{\lambda}|_v < 1.\]

If $\lambda$ is an algebraic unit, the set $P_+$ corresponds to roots of the minimal polynomial of $\lambda$ greater than $1$ in absolute value, counting two conjugate complexes only once, and it is the same for $P_-$ with the roots of modulus less than $1$.

For each place $v$, we define a space $E_v$ as the completion of $\mathbf Q(\lambda)$ for the absolute value $v$. If $v$ is a real place (i.e. corresponding to a real root or the minimal polynomial of $\lambda$), then $E_v = \mathbf R$. If $v$ is a complex place (i.e. corresponding to two conjugated complex roots or the minimal polynomial of $\lambda$), then $E_v = \mathbf C$. Otherwise, $E_v$ is a $p$-adic space, which is a finite extension of the $p$-adic field $\mathbf Q_p$ (which is the completion of $\mathbf Q$ for the $p$-adic absolute value).

We can define the expanding space

\[E_\lambda^+ := \prod_{v \in P_+} E_v,\]

and the contracting one

\[E_\lambda^- := \prod_{v \in P_-} E_v.\]

Let’s take $\sigma_+$ and $\sigma_-$ some embeddings of $\mathbf Q(\lambda)$ into the spaces $E_+$ and $E_-$ respectively. We will also denote by $\sigma_\beta$ the maximal real embedding when $\beta$ is a Perron number.

So Rauzy fractal of the substitution $s$ can be define as the adherence of $\sigma_-(Q_\omega)$ in $E_{\lambda}^-$.

$\beta$-adic sets¶

A $\beta$-adic set, for an number $\beta$, is a subset of $\mathbf Q(\beta)$ of the form

\[{ \mathbf Q_{\beta,L} := \{ \sum_{i=0}^n a_i \beta^i} { n \in \mathbf N,\ a_0 a_1 ... a_n \in L } \}.\]

where $L$ is a regular language over a finite alphabet $\Sigma \subset \mathbf Q(\beta)$. It is represented in Sage by the class BetaAdicSet which contains a number $b$ for $\beta$, and the data of a DetAutomaton $a$ recognizing the language $L$.

Some Properties¶

For a fixed algebraic number $\beta$ with no conjugate of modulus one, the set of $\beta$-adic sets is stable by

intersection
union
complementary (in another $\beta$-adic set)
Minkowski sum (i.e. the sum of two $\beta$-adic sets is a $\beta$-adic set)
multiplication by an element of $\mathbf Q(\beta)$
translation by an element of $\mathbf Q(\beta)$
adherence, interior, boundary, for the topology of $\mathcal O_\beta$ induced by $E_-$.

The fact that $\beta$-adic sets come naturally to describe quasicrystals arising from substitutions and has a lot of nice properties show that it is an interesting fundamental object.

Remark¶

We see from theses properties that we can construct $\beta$-adic sets with any shape in the contracting space $E^-$. This allows us to construct Rauzy fractals of any shape: a $\beta$-adic set comes from a substitution for a Pisot number $\beta^k$ if and only if it is invariant by multiplication by $\beta^k$, it contains zero, and it is a Meyer set in the expanding direction.

Construction of a domain exchange¶

The first step, to construct a substitution from a quasicrystal, is to construct a domain exchange describing the shift on the quasicrystal.

Let $\beta$ be a Pisot number (eventually non unit), and let $Q \subseteq \mathbf Q(\beta)$ such that $\sigma_+(Q)$ is a quasicrystal of $\mathbf R$ or $\mathbf R^+$. Then there exists a domain exchange with a finite number of pieces such that the union of the pieces is $Q$ and such that it is conjugated to the shift on $\sigma_+(Q)$.

Construction of a domain exchange in the unit disk, for the points of the integer ring $\mathcal O_\beta$ that fall in the disk, where $\beta$ is the Tribonnacci number. $\color{red}{-2\beta^2+2\beta}$, $\color{orange}{\beta^2-\beta-1}$, $\color{lime}{\beta-1}$, $\color{green}{1}$, $\color{cyan}{-\beta^2+2\beta+1}$, $\color{blue}{\beta^2-\beta}$, $\color{magenta}{\beta}$

The domain exchange described in the above figure for the open unit disk, gives exactly the list of Pisot numbers (including non-unit ones) of degree $3$ in $\mathbf Q(\beta)$, where $\beta$ is the Tribonnacci number (i.e. greatest root of $x^3-x^2-x-1$). Indeed if $x$ is a Pisot number of degree three in $\mathbf Q(\beta)$, the next Pisot number is obtained by looking in which piece is the conjugate $\overline{x}$, and adding the corresponding translation to $x$.

For a $\beta$-adic set, for a Pisot number $\beta$, we can compute the domain exchange when it is finite. It is done by the function domain_exchange() of the class BetaAdicSet.

Construction of a substitution¶

If we know that a quasicrystal $\sigma_+(Q)$ of $\mathbf R$ or $\mathbf R_+$ comes from the fixed point of a substitution for a Pisot number $\lambda$, it is not difficult to guess what is the substitution. Indeed, it is enough to take intervals between two consecutive points, multiply it by $\lambda$, and see how the result is covered by others intervals.

But we have to take care of the fact that one interval can have several substitutions rules, corresponding to the fact that several letters of a substitution can give intervals of same lengths.

If we look at what happens in the contracting space $E^-$, we have to do a sort of induction on $\lambda Q$ for the domain exchange on $Q$, and we have to iterate it up to stabilization. But it’s not really an induction : we have to distinguish between different possible trajectories for points in $\lambda Q$ before they come back to $\lambda Q$, otherwise the induction only give the same domain exchange on $\lambda Q$ than in $Q$.

For a $\beta$-adic set, for a Pisot number $\beta$, we can compute the substitution when it exists. It is done by the method substitution() of the class BetaAdicSet.

Examples of Usage of BetaAdicSet¶

A Sierpinsky gasket¶

Take the Tribonnacci Pisot number $\beta$, root of $x^3 - x^2 - x - 1$, and take $L$ the regular language defined by the following automaton.

This automaton describe a regular language giving a $\beta$-adic set which is a Sierpiński gasket union a set of non-empty interior for $\beta$ the Tribonnacci number.

Obtained by the code:

# automaton that describe 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)

# automaton recognizing a set of non-empty interior
a2 = DetAutomaton([(0,1,0),(1,2,0),(2,2,0),(2,2,1)],i=0, final_states=[2])
# multiply by b^2
a3 = a.unshift(0, final=True).unshift(1)
a = a2.union(a3)
a.plot()

We define and plot this $\beta$-adic set:

m = BetaAdicSet(x^3-x^2-x-1, a) #choose to work with the alphabet {0,1} and with the Tribonnacci polynomial
m.plot(nprec=6)

The domain exchange with $6$ pieces, describing the shift on $\beta$-adic-set can be computed:

# compute a domain exchange
l = m.domain_exchange()
print("Exchange with %s pieces."%len(l))
Exchange with 6 pieces.
# plot it
m.plot_list([a for t,a in l], nprec=6)

And we plot the domains after exchange

# plot it after exchange
m.plot_list([a.proj(m, t) for t,a in l], nprec=6)

Compute the subtitution

# compute a substitution whose broken line is this BetaAdicSet
%time d = m.substitution()
d
CPU times: user 24 s, sys: 156 ms, total: 24.1 s
Wall time: 24.1 s

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

And we can plot directly the Rauzy fractal from a WordMorphism describing this subtitution:

#plot the Rauzy fractal from the substitution
s = WordMorphism(d)
s.rauzy_fractal_plot()

The Dragon Fractal¶

################################################
# The dragon fractal
################################################
m = BetaAdicSet(1/(1+I), [0,1])
m
b-adic set with b root of x^2 - x + 1/2, and an automaton of 1 states and 2 letters

a = m.relations_automaton(ext=True)
a.plot()

mi = m.intersection_words([0], [1])
m.plot_list([mi])

mi.plot(nprec=6)

Compute the Hausdorff dimension.

sage: # compute the Hausdorff dimension
sage: mi.critical_exponent()
log(y)/log(1.414213562373095?) where y is the max root of x^3 - x^2 - 2, and 1.414213562373095? is root of x^2 - 2.
1.523627086202492

Any Shape¶

We can define $\beta$-adic sets approximating any shape. This can be done with the method approx() of the class BetaAdicSet. And we can also draw a $\beta$-adic set with the mouse, with the method user_draw() of the class BetaAdicSet. We present various examples.

Disk¶

We start by taking a simple $\beta$-adic set.

######################################
# BetaAdicSet approximating a disk
######################################
#. BetaAdicSet from Tribonnacci
m = WordMorphism('a->ab,b->ac,c->a').DumontThomas().mirror()
m
b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 3 states and 2 letters

The relation automaton associated

pm = m.b.parent().places()[1]
pm
Ring morphism:
  From: Number Field in b with defining polynomial x^3 - x^2 - x - 1
  To:   Complex Field with 53 bits of precision
  Defn: b |--> -0.419643377607080 + 0.606290729207199*I

We compute an approximation of a disk included in this $\beta$-adic set:

md = m.approx(14, lambda x: (pm(x).real())^2 + (pm(x).imag())^2 < .4)
print(md)
b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 265 states and 2 letters

m.plot_list([md])

We compute the domain exchange describing this $\beta$-adic set:

md1 = md.proj(m)
md1
b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 253 states and 2 letters

# domain exchange for this set
l = md1.domain_exchange()
print(l)
[(1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 99 states and 2 letters), (b^2 - b, b-adic set with b root of x^3 - x^2
- x - 1, and an automaton of 70 states and 2 letters), (b, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 134 states and 2
letters), (b + 1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 99 states and 2 letters), (b^2, b-adic set with b root of
x^3 - x^2 - x - 1, and an automaton of 164 states and 2 letters), (b^2 + 1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of
61 states and 2 letters), (b^2 + b, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 37 states and 2 letters),
(b^2 + b + 1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 56 states and 2 letters)]

md1.plot_list([a for t,a in l], nprec=6)

And the domain exchange after exchange

# plot the domain exchange after exchange
md1.plot_list([a.proj(md, t) for t,a in l], nprec=6)

Compute a substitution whose Rauzy fractal is this approximation of disk

%time d, lm = md.substitution(get_aut=True)
d
CPU times: user 48.6 s, sys: 783 ms, total: 49.4 s
Wall time: 49.2 s

{1: [248, 318, 288, 324],
[168, 272, 92],
[264, 274],
[407, 2],
[117, 99],
[352, 218],
[226, 235, 372, 323],
[415, 134, 309, 232, 380, 291, 93],
[411, 6, 232, 288, 273, 208],
[374, 310, 123, 168, 292, 92],
[232, 169, 273, 208],
[232, 288, 281],
[411, 329, 232, 288, 273, 208],
[415, 339, 232, 380, 291, 93],
[237],
[4],
[5],
[152],
[8],
...
...
[191, 47],
[191, 308],
[383],
[384],
[373, 48, 125],
[373, 252, 125],
[393, 349],
[394, 349],
[345],
[346],
[352],
[390],
[395, 49, 126],
[396, 49, 126],
[401, 48, 231],
[402, 48, 231],
[395],
[374, 164, 123],
[374, 310, 123],
[406, 348],
[406, 350],
[408],
[409],
[407, 344],
[407, 348],
[412],
[413]}

s = WordMorphism(d)
s.rauzy_fractal_plot()

m.plot_list([a for a,t in lm], nprec=6)

Square¶

Rauzy fractal approximating a square.

sage: #########################################
sage: # choose a Pisot number and an alphabet #
sage: #########################################
sage: pi = x^3-x^2-x-1 #Tribonacci
sage: print(pi.roots(ring=CC))
sage: b = pi.roots(ring=QQbar)[1][0] #we take the conjugate of modulus < 1 in order to plot the result
sage:  m = WordMorphism('a->ab,b->ac,c->a').DumontThomas().mirror()
sage: pm = m.b.parent().places()[1]
sage: ########################################
sage: # Rauzy fractal approximating a square #
sage: ########################################
sage: # compute a g-b-set approximating a square
sage: # the first argument of approx() is the precision,
sage: # and the second one is the characteristic function of the shape to approximate
sage: # the shape must be not too big in order to be inside the set of elements that admit a b-expansion
sage: md  = m.approx(15, lambda x: abs(pm(x).real()) < .5 and abs(pm(x).imag()) < .5 )
sage:  m.plot_list([md])
[(1.83928675521416, 1), (-0.419643377607081 - 0.606290729207199*I, 1), (-0.419643377607081 + 0.606290729207199*I, 1)]
-0.4196433776070806? - 0.6062907292071993?*I

md = md.proj(m)
mdb-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 241 states and 2 letters

domain exchange

# domain exchange for this set
l = md.domain_exchange()
print(l)
md.plot_list([a for t,a in l], nprec=6)
[(1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 84 states and 2 letters),
(b, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 142 states and 2 letters),
(b + 1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 128 states and 2 letters),
(b^2, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 168 states and 2 letters),
(b^2 + 1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 115 states and 2 letters),
(b^2 + b, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 84 states and 2 letters),
(b^2 + b + 1, b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 91 states and 2 letters)]

The projection

# plot the domain exchange after exchange
md.plot_list([a.proj(md, t) for t,a in l], nprec=6)

d , lm = md.substitution(get_aut=True)
d
{1: [244, 131, 200, 174],
[265, 160, 206],
[31, 1, 262],
[314, 72, 132, 265, 86],
[77, 174, 121],
[100, 278, 184, 74, 94],
[298, 308, 203, 277, 124, 217, 105],
[189, 244, 40, 77, 186, 121],
[255, 275, 239, 74],
[56, 69, 1, 262],
[25, 154, 241, 261],
[117, 258, 288],
[301, 66, 157],
[301, 66, 26],
[301, 20, 28],
[201, 140, 22],
[201, 290, 8],
[77, 186, 121],
[255, 304, 152, 239, 74, 93],
[196, 296, 291, 217, 105],
[129, 296, 123, 216, 106],
[244, 43, 77, 174, 121],
...
...
[194, 303, 16],
[194, 303, 39],
[274, 178],
[274, 219],
[201, 290, 293, 18],
[201, 290, 294, 18],
[210, 235, 103],
[210, 236, 103],
[297, 309],
[210, 234, 103],
[210, 234, 169],
[210, 234, 252],
[297, 311],
[297, 312]}

s = WordMorphism(d)
s.rauzy_fractal_plot()

m.plot_list([a for a,t in lm], nprec=6)

Imported Image¶

from sage.arith.beta_adic import ImageIn
im = ImageIn('gabian.png')
im
Image of size 573x628
# compute a beta-adic set approximating the image
w = im.width()
h = im.height()
ma = max(w,h)

# approximate the image using Tribonnacci
m = WordMorphism('a->ab,b->ac,c->a').DumontThomas().mirror()
pm = m.b.parent().places()[1]
md  = m.approx(17, lambda x: (pm(x).conjugate()+.5*(1+I))*ma in im)
m.plot_list([md])

m.relations_automaton().plot()

Domain exchange with 17 pieces for the substitution described in figure

# compute a domain exchange
l = md.domain_exchange()
print("Exchange with %s pieces."%len(l))

Exchange with 17 pieces.

md = md.proj(m)
md.plot_list([a for t,a in l], nprec=6)

# plot the domain exchange after exchange
md.plot_list([a.proj(md, t) for t,a in l], nprec=6)

d , lm = md.substitution(get_aut=True)
d
{1: [122, 137, 113, 64, 104, 142, 43, 107, 100, 47, 102, 10, 105, 131, 145, 128, 33, 78, 144, 39, 86, 63, 66, 153, 134, 69, 123, 110, 124,
98],
...
...
153: [121, 114, 89, 64, 103, 141, 81, 143, 99, 149, 11, 51, 148, 124, 79, 126, 76, 40, 63, 102, 136, 113, 64, 104, 96, 38, 107, 100, 20, 9,
101, 74, 125, 84, 127, 147, 36, 130, 90, 112, 64, 103, 150, 16, 111, 143,  59, 17, 15, 58, 41, 6, 31, 140, 16, 129, 108, 37, 122, 137, 113,
64, 104, 142, 43, 107, 100, 18, 8, 102, 45, 131, 64, 103, 128, 33, 78, 144, 12, 63, 66, 153, 134, 69, 123, 110, 124, 16, 129, 92, 93, 102,
137, 113, 64, 104, 96, 1, 32, 107, 56, 67, 135, 101, 74, 125, 29]}

md.plot_list([a for a,t in lm], nprec=6)

s = WordMorphism(d)
s.rauzy_fractal_plot()

User Draw¶

sage: #########################################
sage: # choose a Pisot number and an alphabet #
sage: #########################################
sage: pi = x^3-x^2-x-1 #Tribonacci
sage: b = pi.roots(ring=QQbar)[1][0] #we take the conjugate of modulus < 1 in order to plot the result
sage: m = BetaAdicSet(b, {0,1}) #choose the alphabet
sage: pm = m.b.parent().places()[1]
sage: aoc = m.user_draw()
[(1.83928675521416, 1), (-0.419643377607081 - 0.606290729207199*I, 1), (-0.419643377607081 + 0.606290729207199*I, 1)]
-0.4196433776070806? - 0.6062907292071993?*I
sage: aoc.plot()

Compute the domain exchange

# compute a domain exchange
l = aoc.domain_exchange()
print("Exchange with %s pieces."%len(l))
Exchange with 13 pieces.

# plot it
aoc.plot_list([a for t,a in l], nprec=6)

# plot it after exchange
aoc.plot_list([a.proj(m, t) for t,a in l], nprec=6)

aoc.relations_automaton()
DetAutomaton with 7 states and an alphabet of 3 letters
aoc.relations_automaton().plot()

aoc.reduced_words_automaton()
DetAutomaton with 43 states and an alphabet of 2 letters
aoc.reduced_words_automaton().plot()

d , lm = aoc.substitution(get_aut=True)
d
{'a': ['o', 'm', 'e', 'r', 'g'],
 'b': ['m', 'b', 'j', 'e', 'r', 'g', 'm', 'p', 'm', 'q'],
 'c': ['m', 'b', 'j', 'e', 'r', 'g', 'm', 'o', 'j', 'e', 'r', 'g', 'm'],
 'd': ['o', 'm', 'q', 'i', 'r', 'g', 'm', 'b', 'j', 'e', 'r', 'g', 'm', 'p', 'm', 'q'],
 'e': ['m', 'b', 'j', 'e', 'r', 'g', 'm', 'o', 'j', 'e', 'r', 'g', 'm', 'o', 'm', 'e', 'r', 'g'],
 'f': ['l', 'd', 'm', 'q'],
 'g': ['p', 'm', 'q'],
 'h': ['o', 'j', 'e', 'r', 'g', 'm'],
 'i': ['o', 'j', 'e', 'r', 'g', 'm', 'o', 'm', 'e', 'r', 'g'],
 'j': ['h', 'a', 'r', 'g', 'j', 'e', 'r', 'f', 'm', 'o', 'm', 'e', 'r', 'g'],
 'k': ['k', 'j', 'e', 'r', 'g', 'm', 'o', 'm', 'e', 'r', 'g', 'm', 'b', 'j', 'e', 'r', 'g', 'm', 'p', 'm', 'q'],
 'l': ['m', 'b', 'j', 'c', 'a', 'r', 'g', 'm'],
 'm': ['i', 'r', 'g', 'j', 'e', 'r', 'g', 'm', 'o', 'm', 'e', 'r', 'g'],
 'n': ['i', 'r', 'g', 'j', 'q', 'i', 'r', 'g', 'm', 'o', 'm', 'e', 'r', 'g'],
 'o': ['m', 'b', 'j', 'c', 'a', 'r', 'g', 'm', 'o', 'm', 'e', 'r', 'g', 'm', 'b', 'j', 'e', 'r', 'g', 'm', 'p', 'm', 'q'],
 'p': ['m', 'b', 'j', 'c', 'a', 'r', 'g', 'm', 'o', 'm', 'q', 'i', 'r', 'g', 'm', 'b', 'j', 'e', 'r', 'g', 'm', 'p', 'm', 'q'],
 'q': ['m', 'b', 'j', 'e', 'r', 'g', 'm'],
 'r': ['n', 'b', 'j', 'e', 'r', 'g', 'm']}

aoc.plot_list([a for a,t in lm], nprec=6)

Intervals¶

For any quadratic Pisot number $\beta$, the $\beta$-adic set $Q_{(-1,1)}$ of all elements of the integer rings that are in the intervalle $(-1,1)$ in the contracting space, can be computed, and then we can compute a substitution describing the quasicrystal.