Asymptotic Analysis of Regular Sequences

In this article, q-regular sequences in the sense of Allouche and Shallit are analysed asymptotically. It is shown that the summatory function of a regular sequence can asymptotically be decomposed as a finite sum of periodic fluctuations multiplied by a scaling factor. Each of these terms corresponds to an eigenvalue of the sum of matrices of a linear representation of the sequence; only the eigenvalues of absolute value larger than the joint spectral radius of the matrices contribute terms which grow faster than the error term. The paper has a particular focus on the Fourier coefficients of the periodic fluctuations: they are expressed as residues of the corresponding Dirichlet generating function. This makes it possible to compute them in an efficient way. The asymptotic analysis deals with Mellin–Perron summations and uses two arguments to overcome convergence issues, namely Hölder regularity of the fluctuations together with a pseudo-Tauberian argument. Apart from the very general result, three examples are discussed in more detail: sequences defined as the sum of outputs written by a transducer when reading a q-ary expansion of the input; the amount of esthetic numbers in the first N natural numbers; and the number of odd entries in the rows of Pascal’s rhombus. For these examples, very precise asymptotic formulæ are presented. In the latter two examples, prior to this analysis only rough estimates were known.

C. Heuberger and D. Krenn are supported by the Austrian Science Fund (FWF): P 28466-N35. The extended abstract [30] (with appendices containing proofs available as arXiv:1802.03266) imposes a restriction on the asymptotic growth. The extended abstract [29] (with appendices containing proofs available as arXiv:1808.00842) lifts this restriction by completely getting rid of the corresponding technical condition. This article now contains the full (majorly restructured) proof covering all cases. It is shorter and simpler. We now use a generating functions approach which also gives additional insights. For example, the cancellations in the proof in [30] seem to be a kind of magic at that point, but with the new approach, it is now clear and no surprise anymore that they have to appear. Besides, the examples investigated in [29,30] are now presented with full details. A new part on computational aspects of the computation of Fourier coefficients is added. Reading strategies for various interests are now outlined in Part I so that readers find their ways through this article. The authors thank Helmut Prodinger for his early involvement and Sara Kropf for her comments on an early version of this paper.
Extended author information available on the last page of the article Keywords Regular sequence · Mellin-Perron summation · Summatory function · Tauberian theorem · Transducer · Esthetic numbers · Pascal's rhombus In this paper, we study the asymptotic behaviour of the summatory function of a qregular sequence x(n). At this point, we give a short overview of the notion of q-regular sequences 1 and our main result.

Mathematics Subject Classification
One characterisation of a q-regular sequence is as follows: the sequence x(n) is said to be q-regular if there are square matrices A 0 , . . . , A q−1 and a vector-valued sequence v(n) such that v(qn + r ) = A r v(n) for 0 ≤ r < q and n ≥ 0 and such that x(n) is the first component of v(n). Regular sequences are intimately related to the q-ary expansion of their arguments. They have been introduced by Allouche and Shallit [2]; see also [3,Chapter 16]. Many special cases have been investigated in the literature; this is also due to their relation to divide-and-conquer algorithms. Moreover, every q-automatic sequencethose sequences are defined by finite automata-is q-regular as well. Take also a look at the book [3] for many examples.
Our main result is, roughly speaking, that the summatory function of a q-regular sequence x(n) has the asymptotic form as N → ∞ for a suitable positive integer J , suitable constants λ j ∈ C, suitable non-negative integers k j , a suitable R and 1-periodic continuous functions k j . The λ j will turn out to be eigenvalues of C := A 0 + · · · + A q−1 , the k j be related to the multiplicities of these eigenvalues and the constant R will be a bound for the joint spectral radius of the matrices A 0 , . . . , A q−1 .
While (1.1) gives the shape of the asymptotic form, gathering as much information as possible on the periodic fluctuations k j is required to have a full picture. To this aim, we will give a description of the Fourier coefficients of the k j which allows to compute them algorithmically and therefore to describe these periodic fluctuations with high precision. In particular, this allows to detect non-vanishing fluctuations. Code 2 is provided to compute the Fourier coefficients.
We close this introductory section by noting that the normalized sum 1 N n<N x(n) enlightens us about the expectation of a random element of the sequence x(n) with respect to uniform distribution on the non-negative integers smaller than a certain N .

How to Read This Paper
This is a long (and perhaps sometimes technical) paper and not all readers might find the time to read it from the very beginning to the very end. We therefore outline reading strategies for various interests.
For the reader who wants to apply our results to a particular problem: Read Sect. 3.1 on the definition of q-regular sequences and Sect. 3.2 containing the main result in a condensed version which should cover most applications. These two sections also have a simple, illustrative and well-known running example. If it turns out that the refined versions of the results are needed, follow the upcoming paragraph below.
For the reader who still wants to apply our results to a particular problem but finds the condensed version insufficient, turn to the overview of the results (Sect. 4.1) and then continue with Sect. 6 where the notations and results are stated in full generality. Formulating them will need quite a number of definitions provided in Sect. 6.2. In order to cut straight to the results themselves, we will refrain from motivations and comments on these definitions and postpone those comments to Sect. 7.
For the reader who wants to determine the asympotics of a regular sequence instead of determining the asymptotics of the summatory function of the regular sequence, advice is given in Sect. 3.3. For the reader who wants to read more about showcase applications of our method yielding new asymptotic results, additionally to Sect. 3 read Sect. 5 where an overview of the examples in this paper is given and then Part II where these examples are discussed in detail. For many more examples to which the methods can be applied, read the original papers [2,4] and the book by Allouche and Shallit [3] which contain many examples of q-regular sequences.
For the reader who wants to compute the Fourier coefficients for a particular application, use the provided code. Read Part IV for more details, in particular, see Sect. 19 for some comments on how to decide whether fluctuations are constant or even vanish.
Moreover, for the reader who is interested in the background on the algorithmic aspects and details of the implementation of the actual computation, we also refer to Part IV; this part will also be useful for the reader who wants to review the code written for SageMath.
For the reader who is interested in the history of the problem, we refer to Sect. 4.4. For the reader who wants to see a heuristic argument why everything works out, there is Sect. 4.2 where it is shown that once one does not care about convergence issues, the Mellin-Perron summation formula of order zero explains the result.
For the reader who wants to understand the idea of the proof, there is Sect. 4.3 with a high level overview of the proof how the above mentioned convergence issues with the Mellin-Perron summation formula can be overcome by a pseudo-Tauberian argument.
For the reader who wants to overcome convergence problems with the Mellin-Perron summation formula in other contexts involving periodic fluctuations, we note that the pseudo-Tauberian argument (Proposition 14.1) is completely independent of our application to q-regular sequences; the only prerequisite is the knowledge on the existence of the fluctuation and sufficient knowledge on analyticity and growth of the Dirichlet generating function. As a consequence, Theorem E has been formulated as an independent result and provisions have been made for several applications of the pseudo-Tauberian argument.
Finally, for the reader who wants to fully understand the proof : We have no other advice than reading the whole introduction, the whole Sect. 6 on results and the whole Part III on the proofs starting with a very short Sect. 11 where a few notations used throughout the proofs are fixed.

q-Regular Sequences
We start by giving a definition of q-regular sequences; see Allouche and Shallit [2]. Let q ≥ 2 be a fixed integer and x be a sequence on Z ≥0 . 3 Then x is said to be (C, q)-regular (briefly: q-regular or simply regular) if the C-vector space generated by its q-kernel x • (n → q j n + r ) : integers j ≥ 0, 0 ≤ r < q j has finite dimension. In other words, x is q-regular if there are an integer D and sequences x 1 , . . . , x D such that for every j ≥ 0 and 0 ≤ r < q j there exist complex numbers c 1 , . . . , c D with x(q j n + r ) = c 1 x 1 (n) + · · · + c D x D (n) for all n ≥ 0.
By Allouche and Shallit [2,Theorem 2.2], the sequence x is q-regular if and only if there exists a vector-valued sequence v whose first component coincides with x and there exist square matrices A 0 , . . . , A q−1 ∈ C d×d such that v(qn + r ) = A r v(n) for 0 ≤ r < q and n ≥ 0. (3.1) This is called a q-linear representation of the sequence x. The best-known example for a 2-regular function is the binary sum-of-digits function.

Example 3.1
For n ≥ 0, let x(n) = s(n) be the binary sum-of-digits of n. We clearly have x(2n + 1) = x(n) + 1 (3.2) for n ≥ 0. Indeed, we have x(2 j n + r ) = x(n) + x(r ) · 1 for integers j ≥ 0, 0 ≤ r < 2 j and n ≥ 0; i.e., the complex vector space generated by the 2-kernel is generated by x and the constant sequence n → 1. Alternatively, we set v = (x, n → 1) and have v(2n) = x(n) for n ≥ 0. Thus (3.1) holds with At this point, we note that a linear representation (3.1) immediately leads to an explicit expression for x(n) by induction. where e 1 = 1 0 . . . 0 .

Condensed Main Result
We are interested in the asymptotic behaviour of the summatory function X (N ) = 0≤n<N x(n). At this point, we give a simplified version of our results. We choose any vector norm · on C d and its induced matrix norm. We set C := 0≤r <q A r . We choose R > 0 such that A r 1 · · · A r = O(R ) holds for all ≥ 0 and r 1 , . . . , r ∈ {0, . . . , q − 1}. In other words, R is an upper bound for the joint spectral radius of A 0 , . . . , A q−1 . The spectrum of C, i.e., the set of eigenvalues of C, is denoted by σ (C). For λ ∈ C, let m(λ) denote the size of the largest Jordan block of C associated with λ; in particular, m(λ) = 0 if λ / ∈ σ (C). Finally, we consider the scalar-valued Dirichlet series X and the vector-valued Dirichlet series V defined by 5 and where v(n) is the vector-valued sequence defined in (3.1). Of course, X (s) is the first component of V(s). The principal value of the complex logarithm is denoted by log. The fractional part of a real number z is denoted by {z} := z − z .
Theorem A (User-friendly all-in-one theorem) With the notations above, we have for suitable 1-periodic continuous functions λk . If there are no eigenvalues λ ∈ σ (C) with |λ| ≤ R, the O-term can be omitted. For |λ| > R and 0 ≤ k < m(λ), the function λk is Hölder continuous with any exponent smaller than log q (|λ|/R).
The Dirichlet series V(s) converges absolutely and uniformly on compact subsets of the half plane s > log q R + 1 and can be continued to a meromorphic function on the half plane s > log q R. It satisfies the functional equation for s > log q R. The right-hand side of (3.4) converges absolutely and uniformly on compact subsets of s > log q R. In particular, V(s) can only have poles where q s ∈ σ (C).
For λ ∈ σ (C) with |λ| > R, the Fourier series converges pointwise for u ∈ R where the Fourier coefficients ϕ λk are defined by the singular expansion 6 This theorem is proved in Sect. 15. We note: • We write λk ({log q N }) to optically emphasise the 1-periodicity; technically, we have λk ({log q N }) = λk (log q N ). • The arguments in the proof could be used to meromophically continue the Dirichlet series to the complex plane, but we do not need this result for our purposes. See [1] for the corresponding argument for automatic sequences. • Sometimes, it will be convenient to write (3.5) in the equivalent explicit formulation In particular, this can be used to algorithmically compute the ϕ λk . • Computing the Fourier coefficients ϕ λk via the explicit formulation (3.6) by reliable numerical arithmetic (see Part IV for details) enables us to detect the non-vanishing of a fluctuation; see also the example below and in Sect. 8 (on sequences defined by transducers) for examples where the fluctuation of the leading term is in fact constant. There, additional arguments are required to actually prove this fact; see Sect. 19 for more details.
We come back to the binary sum of digits.
As A 0 is the identity matrix, any product A r 1 · · · A r has the shape A k 1 = 1 k 0 1 where k is the number of factors A 1 in the product. This implies that R with A r 1 · · · A r = O(R ) may be chosen to be any number greater than 1. As C is a Jordan block itself, we simply read off that the only eigenvalue of C is λ = 2 with m(2) = 2.
Thus Theorem A yields for suitable 1-periodic continuous functions 21 and 20 .
In principle, we can now use the functional equation (3.4) to obtain the Dirichlet series X . Due to the fact that one component of v is the constant sequence where everything is known, it is more efficient to use an ad-hoc calculation for X by splitting the sum according to the parity of the index and using the recurrence relation (3.2) for x(n). We obtain where the Hurwitz zeta function ζ (s, α) := n+α>0 (n + α) −s has been used. We get As the sum of digits is bounded by the length of the expansion, we have x(n) = O(log n). By combining this estimate with we see that the sum in (3.7) converges absolutely for s > 0 and is therefore analytic for s > 0. Therefore, the right-hand side of (3.7) is a meromorphic function for s > 0 whose only pole is simple and at s = 1 which originates from ζ s, 1 2 . Thus, X (s) is a meromorphic function for s > 0 with a double pole at s = 1 and simple poles at 1 + 2 πi log 2 for ∈ Z\{0}. This gives us 21 by (3.6) and (3.7).
We conclude that We will explain in Part IV how to compute rigorous numerical values for the Fourier coefficients, in our case those of the fluctuation 20 which can be deduced from (3.7).
In this particular case of the binary sum-of-digits, simpler and even explicit expressions for the Fourier coefficients have been stated and derived by other authors: they can be obtained in our set-up by rewriting the residues of X (s) in terms of shifted residues of n≥1 (x(n) − x(n − 1)) n −s and by computing the latter explicitly; see [31, Proof of Corollary 2.5]. This yields the well-known result by Delange [9]. It will also turn out that (3.8) being a constant function is an immediate consequence of the fact that 0 1 is a left eigenvector of both A 0 and A 1 associated with the eigenvalue 1; see Theorem B.

Asymptotics of Regular Sequences
This article is written with a focus on the sequence of partial sums of a regular sequence. In this section, however, we explain how to use all material for the regular sequence itself.
Let x(N ) be a q-regular sequence. We may rewrite it as a telescoping sum By [2, Theorems 2.5 and 2.6], the sequence of differences x(n + 1) − x(n) is again q-regular. Conversely, it is also well-known that the summatory function of a q-regular sequence is itself q-regular. (This is an immediate consequence of [2, Theorem 3.1].) Therefore, we might also start to analyse a regular sequence by considering it to be the summatory function of its sequence of differences as in (3.9). In this way, we can apply all of the machinery developed in this article.
We end this short section with some remarks on why focusing on the sequence of partial sums can be rewarding. When modelling a quantity by a regular sequences, its asymptotic behaviour is often not smooth, but the asymptotic behaviour of its summatory function is. Moreover, we will see throughout this work that from a technical perspective, considering partial sums is appropriate. Therefore, we adopt this point of view of summatory functions of q-regular sequences throughout this paper.

Overview of the Results
We have already seen the main results collected in a user-friendly simplified version as Theorem A which was written down in a self-contained way in Sect. 3.2.
In Theorem B the assumptions are refined. In particular, this theorem uses the joint spectral radius R of the matrices in a linear representation of the sequence (instead of a suitable bound for this quantity in Theorem A). Theorem B states the contribution of each eigenvalue of the sum C of matrices of the linear representation-split into the three cases of smaller, equal and larger in absolute value than R, respectively. This is formulated in terms of generalised eigenvectors. As a consequence of this precise breakdown of contributions, Theorem C, which collects the different cases into one result, provides a condition on when the error term vanishes.
Theorem D brings up the full formulation of the functional equation of the Dirichlet series associated to our regular sequence. This is accompanied by a meromorphic continuation as well as bounds on the growth of the Dirichlet series along vertical lines (i.e., points with fixed real value). The analytic properties provided by Theorem D will be used to verify the assumptions of Theorem E.
Theorem E is in fact stated and proved very generally: it is not limited to Dirichlet series coming from matrix products and regular sequences, but it works for general Dirichlet series provided that periodicity and continuity properties of the result are known a priori. This theorem handles the Mellin-Perron summation and the theoretical foundations for the computation of the Fourier coefficients of the appearing fluctuations.
We want to point out that Theorem E can be viewed as a "successful" version of the Mellin-Perron summation formula of order zero. In fact, the theorem states sufficient conditions to provide the analytic justification for the zeroth order formula.
Note that there is another result shown in this article, namely a pseudo-Tauberian theorem for summing up periodic functions. This is formulated as Proposition 14.1, and all the details around this topic are collected in Sect. 14.1. This pseudo-Tauberian argument is an essential step in proving Theorem E.

Heuristic Approach: Mellin-Perron Summation
The purpose of this section is to explain why the formula (3.5) for the Fourier coefficients is expected. The approach here is heuristic and non-rigorous because we do not have the required growth estimates. See also [10].
By the Mellin-Perron summation formula of order 0 (see, for example, [18, Theo- By Remark 3.2 and the definition of R, we have Adding the summand x(0) to match our definition of X (N ) amounts to adding O (1). Shifting the line of integration to the left-we have no analytic justification that this is allowed-and using the location of the poles of X claimed in Theorem A yield for some ε > 0. Expanding N s as and assuming that the remainder integral converges absolutely yield where m λ denotes the order of the pole of X (s)/s at log q λ + 2 πi log q and ϕ λk is as in (3.5). (For λ = 1 and k = 0, the contribution of x(0)/s in (3.5) is absorbed by the error term O(1) here.) Summarising, this heuristic approach explains most of the formulae in Theorem A. Some details (exact error term and order of the poles) are not explained by this approach. A result "repairing" the zeroth order Mellin-Perron formula is known as Landau's theorem; see [5,Sect. 9]. It is not applicable to our situation due to multiple poles along vertical lines which then yield the periodic fluctuations. Instead, we present Theorem E which provides the required justification (not by estimating the relevant quantities, but by reducing the problem to higher order Mellin-Perron summation). The essential assumption is that the summatory function can be decomposed into fluctuations multiplied by some growth factors such as in (3.3).

High Level Overview of the Proof
As we want to use Mellin-Perron summation in some form, we derive properties of the Dirichlet series associated to the regular sequence. In particular, we derive a functional equation which allows to compute the Dirichlet series and its residues with arbitrary precision (Theorem D).
We cannot directly use Mellin-Perron summation of order zero for computing the Fourier coefficients of the fluctuations of interest. As demonstrated in Sect. 4.2, however, our theorems coincide with the results which Mellin-Perron summation of order zero would give if the required growth estimates could be provided. Unfortunately, we are unable to prove these required growth estimates. Therefore, we have to circum-vent the problem by applying a generalisation of the pseudo-Tauberian argument by Flajolet, Grabner, Kirschenhofer, Prodinger and Tichy [18].
In order to use this argument, we have to know that the asymptotic formula has the shape (3.3). Note that a successful application (not directly possible!) of Mellin-Perron summation of order zero would give this directly. Therefore, we first prove (3.3) and the existence of the fluctuations (Theorems B, C). To do so, we decompose the problem into contributions of the eigenspaces of the matrix C = A 0 + · · · + A q−1 . The regular sequence is then expressed as a matrix product. Next, we construct the fluctuations by elementary means: We replace finite sums occurring in the summatory functions by infinite sums involving digits using the factorisation as a matrix product.
Then the pseudo-Tauberian argument states that the summatory function of the fluctuation is again a fluctuation and there is a relation between the Fourier coefficients of these fluctuations. The Fourier coefficients of the summatory function of the fluctuation, however, can be computed by Mellin-Perron summation of order one, so the Fourier coefficients of the original fluctuation can be recovered; see Theorem E.

Relation to Previous Work
The asymptotics of the summatory function of specific examples of regular sequences has been studied in [14,23,24]. There, various methods have been used to show that the fluctuations exist; then the original pseudo-Tauberian argument by Flajolet, Grabner, Kirschenhofer, Prodinger and Tichy [18] is used to compute the Fourier coefficients of the fluctuations.
The first version of the pseudo-Tauberian argument in Theorem E was provided in [18]: there, no logarithmic factors were allowed, only values γ with γ > 0 were allowed and the result contained an error term of o(1) whereas we give a more precise error estimate in order to allow repeated application.
Dumas [12,13] proved the first part of Theorem A using dilation equations. We re-prove it here in a self-contained way because we need more explicit results than obtained by Dumas (e.g., we need explicit expressions for the fluctuations) to explicitly get the precise structure depending on the eigenspaces (Theorem B). Notice that the order of factors in Dumas' paper is inconsistent between his versions of (3.1) and Remark 3.2.
A functional equation for the Dirichlet series of an automatic sequence has been proved by Allouche, Mendès France and Peyrière [1].
In Sect. 8 we study transducers. The sequences there are defined as the output sum of transducer automata in the sense of [31]. They are a special case of regular sequences and are a generalisation of many previously studied concepts. In that case, much more is known (variance, limiting distribution, higher dimensional input); see [31] for references and results. A more detailed comparison can be found in Sect. 8. Divide and conquer recurrences (see [11,32]) can also be seen as special cases of regular sequences.
The present article gives a unified approach which covers all cases of regular sequences. As long as the conditions on the joint spectral radius are met, the main asymptotic terms are not absorbed by the error terms. Otherwise, the regular sequence is so irregular that the summatory function is not smooth enough to allow a result of this shape.

Overview of the Examples
We take a closer look at three particular examples. In this section, we provide an overview of these examples; all details can be found in Part II.
At first gance it seems that these examples are straight-forward applications of the results. However, we have to reformulate the relevant questions in terms of a q-regular sequence and will then provide shortcuts for the computation of the Fourier series. We put a special effort on the details which gives additional insights like dependencies on certain residue classes; see Sect. 5.3. Moreover, the study of these examples also encourages us to investigate symmetries in the eigenvalues; see Sect. 5.4 for an overview and Sect. 6.6 for general considerations.
We start with transducer automata. Transducers have been chosen in order to compare the results here with the previously available results [31]. In some sense, the results complement each other: while the results in [31] also contain information on the variance and the limiting distribution, our approach here yields more terms of the asymptotic expansion of the mean, at least in the general case. Also, it is a class of examples.
We then continue with esthetic numbers. These numbers are an example of an automatic sequence, therefore can be treated by a transducer. However, it turns out that the generic results (the results here and in [31]) degenerate: they are too weak to give a meaningful main term. Therefore a different effort is needed for esthetic numbers. No precise asymptotic results were known previously.
The example on Pascal's Rhombus is a choice of a regular sequence where all components of the vector sequence have some combinatorial meaning. Again, no precise asymptotic results were known previously. Section 5.6 contains further examples. Note that there are the two additional Sects. 5.3 and 5.4 pointing out phenomena appearing in the analysis of our examples.

Transducers
The sum T (n) of the output labels of a complete deterministic finite transducer T when reading the q-ary expansion of an integer n has been investigated in [31]. As this can be seen as a q-regular sequence, we reconsider the problem in the light of our general results in this article; see Sect. 8. For the summatory function, the main terms corresponding to the eigenvalue q can be extracted by both results; if there are further eigenvalues larger than the joint spectral radius, our Corollary F allows to describe more asymptotic terms which are absorbed by the error term in [31]. Note, however, that our approach here does not give any readily available information on the variance (this could somehow be repaired for specific examples because regular sequences are known to form a ring) nor on the limiting distribution.

Esthetic Numbers
In this article, we also contribute a precise asymptotic analysis of q-esthetic numbers; see De Koninck and Doyon [8]. These are numbers whose q-ary digit expansion satisfies the condition that neighboring digits differ by exactly one. The sequence of such numbers turns out to be q-automatic, thus are q-regular and can also be seen as an output sum of a transducer; see the first author's joint work with Kropf and Prodinger [31] or Sect. 8. However, the asymptotics obtained by using the main result of [31] is degenerated in the sense that the provided main term and second order term both equal zero; only an error term remains. On the other hand, using a more direct approach via our main theorem brings up the actual main term and the fluctuation in this main term. We also explicitly compute the Fourier coefficients. The full theorem is formulated in Sect. 9. Prior to this precise analysis, the authors of [8] only performed an analysis of esthetic numbers by digit-length (and not by the number itself).
The approach used in the analysis of q-esthetic numbers can easily be adapted to numbers defined by other conditions on the word of digits of their q-ary expansion.

Dependence on Residue Classes
The analysis of q-esthetic numbers also brings another aspect into the light of day, namely a quite interesting dependence of the behaviour with respect to q on different moduli: • The dimensions in the matrix approach of [8] need to be increased for certain residue classes of q modulo 4 in order to get a formulation as a q-automatic and q-regular sequence, respectively. • The main result in [8] already depends on the parity of q (i.e., on q modulo 2). This reflects our Corollary G by having 2-periodic fluctuations (in contrast to 1-periodic fluctuations in the main Theorem A).
• Surprisingly, the error term in the resulting formula of Corollary G depends on the residue class of q modulo 3. This can be seen in the spectrum of the matrix C = 0≤r <q A r : there is an appearance of an eigenvalue 1 in certain cases. • As an interesting side-note: In the spectrum of C, the algebraic multiplicity of the eigenvalue 0 changes again only modulo 2.

Symmetrically Arranged Eigenvalues
Fluctuations with longer periods (like in the second of the four bullet points above) come from a particular configuration in the spectrum of C. Whenever eigenvalues are arranged as vertices of a regular polygon, then their influence can be collected; this results in periodic fluctuations with larger period than 1. We elaborate on the influence of such eigenvalues in Sect. 6.6. This is then used in the particular cases of esthetic numbers and in conjunction with the output sum of transducers. More specifically, in the latter example this yields the second order term in Corollary F; see also [31].

Pascal's Rhombus
Beside esthetic numbers, we perform an asymptotic analysis of the number of ones in the rows of Pascal's rhombus. The rhombus is in some sense a variant of Pascal's triangle-its recurrence is similar to that of Pascal's triangle. It turns out that the number of ones in the rows of Pascal's rhombus can be modelled by a 2-regular sequence. The authors of [21] investigate this number of ones, but only for blocks whose number of rows is a power of 2. In the precise analysis in Sect. 10 we not only obtain the asymptotic formula, we also explicitly compute the Fourier coefficients.

Further Examples
There are many further examples of specific q-regular sequences which await precise asymptotic analysis, for example the Stern-Brocot sequence [39, A002487], the denominators of Farey tree fractions [39, A007306], the number of unbordered factors of length n of the Thue-Morse sequence (see [22]).
The Stern-Brocot sequence is a typical example: it is defined by x(0) = 0, x(1) = 1 and i.e., the right-hand sides are linear combinations of shifted versions of the original sequence. Note that recurrence relations like (5.1) are not proper linear representations of regular sequences in the sense of (3.1). The good news, however, is that in general, such a sequence is q-regular. The following remark formulates this more explicitly.

Remark 5.1
Let x(n) be a sequence such that there are fixed integers ≤ 0 ≤ u and constants c rk for 0 ≤ r < q and ≤ k ≤ u such that holds for 0 ≤ r < q and n ≥ 0. Then the sequence x(n) is q-regular with q-linear representation for v(n) = x(n + ), . . . , x(n), . . . , x(n + u ) where Note that if < 0, then a simple permutation of the components of v(n) brings x(n) to its first component (so that the above is indeed a proper linear representation as defined in Sect. 3.1).

Full Results
In this section, we fully formulate our results. As pointed out in Remark 3.2, regular sequences can essentially be seen as matrix products. Therefore, we will study these matrix products instead of regular sequences. Theorem A can then be proved as a simple corollary of the results for matrix products; see Sect. 15.

Problem Statement
Let q ≥ 2, d ≥ 1 be fixed integers and A 0 , . . . , A q−1 ∈ C d×d . We investigate the sequence f of d × d matrices such that f (qn + r ) = A r f (n) for 0 ≤ r < q, 0 ≤ n with qn + r = 0 (6.1) and f (0) = I . Let n be an integer with q-ary expansion r −1 . . . r 0 . Then it is easily seen that (6.1) implies that We are interested in the asymptotic behaviour of F(N ) := 0≤n<N f (n).

Definitions and Notations
In this section, we give all definitions and notations which are required in order to state the results. For the sake of conciseness, we do not give any motivations for our definitions here; those are deferred to Sect. 7.
The following notations are essential: • Let · denote a fixed norm on C d and its induced matrix norm on C d×d .
• We set B r := 0≤r <r A r for 0 ≤ r < q and C := 0≤r <q A r .
• The joint spectral radius of A 0 , . . . , A q−1 is denoted by If the set of matrices A 0 , . . . , A q−1 has the finiteness property, i.e., there is an > 0 such that then we set R = ρ. Otherwise, we choose R > ρ in such a way that there is no eigenvalue λ of C with ρ < |λ| ≤ R. • The spectrum of C, i.e., the set of eigenvalues of C, is denoted by σ (C).
• For a positive integer n 0 , let F n 0 be the matrix-valued Dirichlet series defined by log q for k ∈ Z. In the formulation of Theorems B and C, the following constants are needed additionally: • Choose a regular matrix T such that T CT −1 =: J is in Jordan form.
• Let D be the diagonal matrix whose jth diagonal element is 1 if the jth diagonal element of J is not equal to 1; otherwise the jth diagonal element of D is 0.
• For a λ ∈ C, let m(λ) be the size of the largest Jordan block associated with λ. In particular, All implicit O-constants depend on q, d, the matrices A 0 , . . . , A q−1 (and therefore on ρ), as well as on R.

Decomposition into Periodic Fluctuations
Instead of considering F(N ), it is certainly enough to consider wF(N ) for all generalised left eigenvectors w of C, e.g., the rows of T . The result for F(N ) then follows by taking appropriate linear combinations.
Theorem B Let w be a generalised left eigenvector of rank m of C corresponding to the eigenvalue λ.

If |λ| < R, then
for N ≥ q m−1 . The function k is Hölder continuous with any exponent smaller than log q |λ|/R. If, additionally, the left eigenvector w(C − λI ) m−1 of C happens to be a left eigenvector to each matrix A 0 , . . . , A q−1 associated with the eigenvalue 1, then Here, wK = 0 for λ = 1 and wϑ m = 0 for λ = 1.
This theorem is proved in Sect. 12. Note that in general, the three summands in the theorem have different growths: a constant, a logarithmic term and a term whose growth depends essentially on the joint spectral radius and the eigenvalues larger than the joint spectral radius, respectively. The vector w is not directly visible in front of the third summand; instead, the vectors of its Jordan chain are part of the function k . Expressing the identity matrix as linear combinations of generalised left eigenvalues and summing up the contributions of Theorem B essentially yields the following corollary.

Theorem C With the notations above, we have
for suitable 1-periodic continuous functions λk . If 1 is not an eigenvalue of C, then ϑ = 0. If there are no eigenvalues λ ∈ σ (C) with |λ| ≤ ρ, then the O-term can be omitted.
This theorem is proved in Sect. 12.4. Remark 6. 1 We want to point out that the condition |λ| > R is inherent in the problem: Single summands f (n) might be as large as n log q R and must therefore be absorbed by the error term in any smooth asymptotic formula for the summatory function.

Dirichlet Series
This section gives the required result on the Dirichlet series F n 0 . For theoretical purposes, it is enough to study F := F 1 ; for numerical purposes, however, convergence improves for larger values of n 0 . This is because for large n 0 and large s, the value of F n 0 (s) is roughly n −s 0 f (n 0 ); see also Part IV. Theorem D Let n 0 be a positive integer. Then the Dirichlet series F n 0 (s) converges absolutely and uniformly on compact subsets of the half plane s > log q ρ + 1, thus is analytic there.
Here, the implicit O-constant also depends on δ.
Note that by the introductory remark on F n 0 (s), the infinite sum over k in (6.4) can be well approximated by a finite sum. Detailed error bounds are discussed in Part IV. Therefore the theorem allows to transfer the information on F n 0 (s) for large s where convergence is unproblematical to values of s where the convergence of the Dirichlet series F n 0 itself is bad. Remark 6.2 By the identity theorem for analytic functions, the meromorphic continuation of F n 0 is unique on the domain given in the theorem. Therefore, the bound (6.5) does not depend on the particular expression for the meromorphic continuation given in (6.3) and (6.4).
Theorem D is proved in Sect. 13. In the proof we translate the linear representation of f into a system of equations involving F n 0 (s) and shifted versions like n≥n 0 f (n)(n + β) −s . We will have to bound the difference between the shifted and unshifted versions of the Dirichlet series. These bounds are provided by the following lemma. It will turn out to be useful to have it as a result listed in this section and not buried in the proofs sections.

Lemma 6.3 Let D(s) = n≥n 0 d(n)/n s be a Dirichlet series with coefficients d(n)
where the series converges absolutely and uniformly on compact sets for s > log q ρ, thus (s, β, D) is analytic there. Moreover, with μ δ as in Theorem D, as | s| → ∞ holds uniformly for log q ρ + δ ≤ s ≤ log q ρ + δ + 1.

Fourier Coefficients
As discussed in Sect. 4.2, we would like to apply the zeroth order Mellin-Perron summation formula but need analytic justification. In the following theorem we prove that whenever it is known that the result is a periodic fluctuation, the use of zeroth order Mellin-Perron summation can be justified. In contrast to the remaining parts of the paper, this theorem does not assume that f (n) is a matrix product.
Theorem E Let f be a sequence on Z >0 , let γ 0 ∈ R\Z ≤0 and γ ∈ C with γ > γ 0 , δ > 0, q > 1 be real numbers with δ ≤ π/(log q) and δ < γ − γ 0 , and let m be a positive integer. Moreover, let j be Hölder continuous (with exponent α with • there is some real number σ abs ≥ γ such that F(s) converges absolutely for s > σ abs ; • the function F(s)/s can be continued to a meromorphic function for s > γ 0 − δ such that poles can only occur at γ + χ for ∈ Z and such that these poles have order at most m and a possible pole at 0; the local expansions are written as with suitable constants ϕ j for j, ∈ Z; • there is some real number η > 0 such that for γ 0 ≤ s ≤ σ abs and |s −γ −χ | ≥ δ for all ∈ Z, we have All implicit O-constants may depend on f , q, m, γ , γ 0 , α, δ, σ abs and η. Then for u ∈ R, ∈ Z and 0 ≤ j < m. If γ 0 < 0 and γ / ∈ 2πi log q Z, then F(0) = 0. This theorem is proved in Sect. 14. The theorem is more general than necessary for q-regular sequences because Theorem D shows that we could use some 0 < η < 1. However, it might be applicable in other cases, so we prefer to state it in this more general form.

Fluctuations of Symmetrically Arranged Eigenvalues
In our main results, the occurring fluctuations are always 1-periodic functions. However, if eigenvalues of the sum of matrices of the linear representation are arranged in a symmetric way, then we can combine summands and get fluctuations with longer periods. This is in particular true if all vertices of a regular polygon (with center 0) are eigenvalues.

Proposition 6.4
Let λ ∈ C, and let k ≥ 0 and p > 0 be integers. Denote by U p the set of pth roots of unity. Suppose for each ζ ∈ U p we have a continuous 1-periodic function whose Fourier coefficients are with a continuous p-periodic function Note that we again write ( p{log q p N }) to optically emphasise the p-periodicity. Moreover, the factor (log q N ) k in (6.10) could be cancelled, however it is there to optically highlight the similarities to the main results (e.g. Theorem A). The proof of Proposition 6.4 can be found in Sect. 16.
The above proposition will be used for proving Corollary F which deals with transducer automata; there, the second order term exhibits a fluctuation with possible period larger than 1. We will also use the proposition for the analysis of esthetic numbers in Sect. 9. Remark 6. 5 We can view Proposition 6.4 from a different perspective: A q-regular sequence is q p -regular as well (by [2, Theorem 2.9]). Then, all eigenvalues ζ λ of the original sequence become eigenvalues λ p whose algebraic multiplicity is the sum of the individual multiplicities but the sizes of the corresponding Jordan blocks do not change. Moreover, the joint spectral radius is also taken to the pth power. We apply, for example, Theorem A in our q p -world and get again 1-period fluctuations. Note that for actually computing the Fourier coefficients, the approach presented in the proposition seems to be more suitable.

Remarks on the Definitions
In this section, we give some motivation for and comments on the definitions listed in Sect. 6.2.

q-Regular Sequences Versus Matrix Products
We note one significant difference between the study of q-regular sequences as in (3.1) and the study of matrix products (6.2). The recurrence (3.1) is supposed to hold for qn + r = 0, too; i.e. v(0) = A 0 v(0). This implies that v(0) is either the zero vector (which is not interesting at all) or that v(0) is a right eigenvector of A 0 associated with the eigenvalue 1.
We do not want to impose this condition in the study of the matrix product (6.2). Therefore, we exclude the case qn + r = 0 in (6.1). This comes at the price of the terms K , ϑ m , ϑ in Theorem B which vanish if multiplied by a right eigenvector to the eigenvalue 1 of A 0 from the right. This is the reason why Theorem A has simpler expressions than those encountered in Theorem B.

Joint Spectral Radius
Let Then the submultiplicativity of the norm and Fekete's subadditivity lemma [15] imply that lim →∞ ρ = inf >0 ρ = ρ; cf. [37]. In view of equivalence of norms, this shows that the joint spectral radius does not depend on the chosen norm. For our purposes, the important point is that the choice of R ensures that there is an 0 > 0 such that For any > 0, we use long division to write = s 0 + r , and by submultiplicativity of the norm, we get A r 1 . . . A r ≤ R s 0 ρ r r and thus for all r j ∈ {0, . . . , q −1} and → ∞. We will only use (7.1) and no further properties of the joint spectral radius. Note that (6.2) and (7.1) imply that for n → ∞.
As mentioned, we say that the set of matrices A 0 , . . . , A q−1 , has the finiteness property if there is an > 0 with ρ = ρ; see [34,35].

Constants for Theorem B
In contrast to usual conventions, we write matrix representations of endomorphisms as multiplications x → x M where x is a (row) vector in C d and M is a matrix. Note that we usually denote this endomorphism by the corresponding calligraphic letter, for example, the endomorphism represented by the matrix M is denoted by M.
Consider the endomorphism C which maps a row vector x ∈ C d to xC and its generalised eigenspaces W λ for λ ∈ C. (These are the generalised left eigenspaces of Let T be the basis formed by the rows of T . Then the matrix representation of C with respect to T is J .
Let now D be the endomorphism of C d which acts as identity on W λ for λ = 1 and as zero on W 1 . Its matrix representation with respect to the basis T is D; its matrix representation with respect to the standard basis is T −1 DT .
Finally, let C be the endomorphism C = C • D. As C and D decompose along C d = λ∈σ (C) W λ and D commutes with every other endomorphism on W λ for all λ, we clearly also have C = D • C. Thus the matrix representation of C with respect to T is D J = J D; its matrix representation with respect to the standard basis is Now consider a generalised left eigenvector w of C. If it is associated to the eigenvalue 1, then and wϑ m = 0. Also note that 1 is not an eigenvalue of C , thus I − C is indeed regular. If 1 is not an eigenvalue of C, then everything is simpler: D is the identity matrix,

Part II: Examples
In this part we investigate three examples in-depth. For an overview, we refer to Sect. 5 where some of the appearing phenomena are discussed as well. Further examples are also mentioned there.

Sequences Defined by Transducer Automata
We discuss the asymptotic analysis related to transducers; see also Sect. 5.1 for an overview.

Transducer and Automata
Let us start with two paragraphs recalling some notions around transducer automata. A transducer automaton has a finite set of states together with transitions (directed edges) between these states. Each transition has an input label and an output label out of the input alphabet and the output alphabet, respectively. A transducer is said to be deterministic and complete if for every state and every letter of the input alphabet, there is exactly one transition starting in this state with this input label. A deterministic and complete transducer processes a word (over the input alphabet) in the following way: • It starts at its unique initial state.
• Then the transducer reads the word letter by letter and for each letter -takes the transition with matching input label, -the output label is written, and -we proceed to the next state (according to the end of the transition).
• Each state has a final output label that is written when we halt in this final state; we call a transducer with this property a subsequential transducer.
We refer to [6, Chapter 1] for a more detailed introduction to transducers and automata. Now we are ready to start with the set-up for our example.

Sums of Output Labels
Let q ≥ 2 be a positive integer. We consider a complete deterministic subsequential transducer T with input alphabet {0, . . . , q − 1} and output alphabet C; see [31]. For a non-negative integer n, let T (n) be the sum of the output labels (including the final output label) encountered when the transducer reads the q-ary expansion of n. Therefore, letters of the input alphabet will from now on be called digits. This concept has been thoroughly studied in [31]: there, T (n) is considered as a random variable defined on the probability space {0, . . . , N − 1} equipped with uniform distribution. The expectation in this model corresponds (up to a factor of N ) to our summatory function 0≤n<N T (n). We remark that in [31], the variance and limiting distribution of the random variable T (n) have also been investigated. Most of the results there are also valid for higher dimensional input.
The purpose of this section is to show that T (n) is a q-regular sequence and to see that the corresponding results in [31] also follow from our more general framework here. We note that the binary sum of digits considered in Example 3.1 is the special case of q = 2 and the transducer consisting of a single state which implements the identity map. For additional special cases of this concept; see [31]. Note that our result here for the summatory function contains (fluctuating) terms for all eigenvalues λ of the adjacency matrix of the underlying digraph with |λ| > 1 whereas in [31] only contributions of those eigenvalues λ with |λ| = q are available, all other contributions are absorbed by the error term there.

Some Perron-Frobenius Theory
We will need the following consequence of Perron-Frobenius theory. By a component of a digraph we always mean a strongly connected component. We call a component final if there are no arcs leaving the component. The period of a component is the greatest common divisor of its cycle lengths. The final period of a digraph is the least common multiple of the periods of its final components. Lemma 8.1 Let D be a directed graph where each vertex has outdegree q. Let M be its adjacency matrix and p be its final period. Then M has spectral radius q, q is an eigenvalue of M and for all eigenvalues λ of M of modulus q, the algebraic and geometric multiplicities coincide and λ = qζ for some pth root of unity ζ .
This lemma follows from setting t = 0 in [31, Lemma 2.3]. As [31,Lemma 2.3] proves more than we need here and depends on the notions of that article, we extract the relevant parts of [31] to provide a self-contained (apart from Perron-Frobenius theorem) proof of Lemma 8.1.
Proof As usual, the condensation of D is the graph resulting from contracting each component of the original digraph to a single new vertex. By construction, the condensation is acyclic.
We choose a refinement of the partial order of the components given by the successor relation in the condensation to a linear order in such a way that the final components come last. Note that this implies that if there is an arc from one component to another, the former component comes before the latter component in our linear order. We then denote the components by C 1 , . . . , C k , C k+1 , . . . , C k+ where the first k components are non-final and the last are final. Without loss of generality, we assume that the vertices of the original digraph D are labeled such that vertices within a component get successive labels and such that the linear order of the components established above is respected.
Therefore, the adjacency matrix M is an upper block triagonal matrix of the shape where M j is the adjacency matrix of the component C j .
Each row of the non-negative square matrix M has sum q by construction. Thus M ∞ = q and therefore the spectral radius of M is bounded from above by q. As the all ones vector is obviously a right eigenvector associated with the eigenvalue q of M, the spectral radius of M equals q. The same argument applies to M k+1 , . . . , M k+ .
By construction, the matrices M k+1 , . . . , M k+ are irreducible. For 1 ≤ j ≤ all eigenvalues λ of M k+ j of modulus q have algebraic and geometric multiplicities 1 by Perron-Frobenius theory and λ = qζ for some p k+ j th root of unity ζ where p k+ j is the period of C k+ j .
By construction, the vertices of the components C j for 1 ≤ j ≤ k have out-degree at most q. We add loops to these vertices to increase their out-degree to q, resulting in C j . The corresponding adjacency matrices are denoted by M j . By the above argument, M j has spectral radius q for 1 ≤ j ≤ k. As M j ≤ M j (component-wise) and M j = M j by construction, the spectral radius of M j is strictly less than q by [20,Theorem 8.8.1].
A left eigenvector v j of M k+ j for 1 ≤ j ≤ can easily be extended to a left eigenvector (0, . . . , 0, v j , 0, . . . , 0) of M. This observation shows that the geometric multiplicity of any eigenvalue of M of modulus q is at least its algebraic multiplicity. This concludes the proof.

Analysis of Output Sums of Transducers
We consider the states of T to be numbered by {1, . . . , d} for some positive integer d ≥ 1 such that the initial state is state 1. We set T j (n) to be the sum of the output labels (including the final output label) encountered when the transducer reads the q-ary expansion of n when starting in state j. By construction, we have 1}-matrix P r in such a way that there is a one in row j, column k if and only if there is a transition from state j to state k with input label r . The vector o r is defined by setting its jth coordinate to be the output label of the transition from state j with input label r .
For n 0 ≥ 1, we set The last Dirichlet series is a truncated version of the Hurwitz zeta function.
Corollary F Let T be a transducer as described at the beginning of this section. Let M be the adjacency matrix and p be the final period of the underlying digraph. For λ ∈ C let m(λ) be the size of the largest Jordan block associated with the eigenvalue λ of M. Then the sequence n → T (n) is a q-regular sequence and for some continuous p-periodic function , some continuous 1-periodic functions λk for λ ∈ σ (M) with 1 < |λ| < q and 0 ≤ k < m(λ) and some constant e T . Furthermore, The Dirichlet series Y n 0 satisfies the functional equation Note that the functional equation (8.2) is preferrable over the functional equation given in Theorem D for the generic case of a regular sequence: the generic functional equation suggests a double pole at s = 1 + χ for all ∈ Z whereas the occurrence of the Hurwitz zeta function in (8.2) shows that there is a double pole s = 1 but single poles at s = 1+χ for all ∈ Z\{0}. Numerically, the same occurrence of the Hurwitz zeta function is also advantageous because it allows to decouple the problem.

Proof of Corollary
and o( j, r ) to be the target state and output label of the unique transition from state j with input label r , respectively. Therefore, for n ≥ 0, 0 ≤ r < q with qn + r > 0. Defining f (n) as in (6.1) for these A r , we see q-Regular Sequence. If we insist on a proper formulation as a regular sequence, we rewrite (8.3) to for n ≥ 0, 0 ≤ r < q.
Eigenvalue 1. By construction, the matrices A r have the shape It is clear that (0, . . . , 0, 1) is a left eigenvector of A r associated with the eigenvalue 1.
Joint Spectral Radius. We claim that A 0 , . . . , A q−1 have joint spectral radius 1. Let · ∞ denote the maximum norm of complex vectors as well as the induced matrix norm, i.e., the maximum row sum norm. Let j 1 , . . . , j ∈ {0, . . . , q − 1}. It is easily shown by induction on that A j 1 · · · A j = P b P 0 1 for some P ∈ C d×d and b P ∈ C d with P ∞ ≤ 1 and b P ∞ ≤ max 0≤r <q o r ∞ . Thus, we obtain As 1 is an eigenvalue of each matrix A r for 0 ≤ r < q, the joint spectral radius equals 1, which proves the claim.
Eigenvectors and Asymptotics. We now consider C = 0≤r <q A r . It has the shape where b M is some complex vector. Let w 1 , . . . , w be a linearly independent system of left eigenvectors of M associated with the eigenvector q. If w j b M = 0 for 1 ≤ j ≤ , then (w 1 , 0), . . . , (w , 0), (0, 1) is a linearly independent system of left eigenvectors of C associated with the eigenvalue q. In that case and because of Lemma 8.1, algebraic and geometric multiplicities of q as an eigenvalue of C are both equal to + 1.
Otherwise, assume without loss of generality that w 1 b M = 1. Then is a linearly independent system of left eigenvectors of C associated with the eigenvalue q. Additionally, (w 1 , 0) is a generalised left eigenvector of rank 2 of C associated with the eigenvalue q with (w 1 , 0)(C − q I ) = (0, 1). As noted above, the vector (0, 1) is a left eigenvector to each matrix A 0 , . . . , A q−1 . Similarly, it is easily seen that any left eigenvector of M associated with some eigenvalue λ = q can be extended uniquely to a left eigenvector of C associated with the same eigenvalue. The same is true for chains of generalised left eigenvectors associated with λ = q.
Therefore, in both of the above cases, Theorem B yields 0≤n<N for some constant e T (which vanishes in the first case) and some 1-periodic continuous functions (qζ ) and λk where ζ runs through the pth roots of unity U p and λ through the eigenvalues of M with 1 < |λ| < q and 0 ≤ k < m(λ). Proposition 6.4 leads to (8.1).
Fourier Coefficients. By Theorem A, we have for a pth root of unity ζ ∈ U p and ∈ Z. Therefore and by noting that T (0) does not contribute to the residue, Proposition 6.4 leads to the Fourier series given in the corollary.

Esthetic Numbers
We discuss the asymptotic analysis of esthetic numbers; see also Sect. 5.2 for an overview. Let again be q ≥ 2 a fixed integer. We call a non-negative integer n a q-esthetic number (or simply an esthetic number) if its q-ary digit expansion r −1 . . . r 0 satisfies |r j − r j−1 | = 1 for all j ∈ {1, . . . , − 1}; see De Koninck and Doyon [8].
In [8] the authors count q-esthetic numbers with a given length of their q-ary digit expansion. They provide an explicit (in form of a sum of q summands) as well as an asymptotic formula for these counts. We aim for a more precise analysis and head for an asymptotic description of the amount of q-esthetic numbers up the an arbitrary value N (in contrast to only powers of q in [8]).

A q-Linear Representation
The language consisting of the q-ary digit expansions (seen as words of digits) which are q-esthetic is a regular language, because it is recognized by the automaton A in Fig. 1. Therefore, the indicator sequence of this language, i.e., the nth entry is 1 if n is q-esthetic and 0 otherwise, is a q-automatic sequence and therefore also q-regular. Let us name this sequence x(n).
Let A 0 , . . . , A q−1 be the transition matrices of the automaton A, i.e., A r is the adjacency matrix of the directed graph induced by a transition with digit r . To make this more explicit, we have the following (q + 1)-dimensional square matrices: Each row and column corresponds to the states 0, 1, . . . , q − 1, I. In matrix A r , the only non-zero entries are in column r ∈ {0, 1, . . . , q − 1}, namely 1 in the rows r − 1 and r + 1 (if available) and in row I as there are transitions from these states to state r in the automaton A.
Let us make this more concrete by considering q = 4. We obtain the matrices We are almost at a q-linear representation of our sequence; we still need vectors on both sides of the matrix products. We have for r −1 . . . r 0 being the q-ary expansion of n and vectors e q+1 = 0 . . . 0 1 and v(0) = 0 1 . . . 1 . As A 0 v(0) = 0 = v(0), this is not a linear representation of a regular sequence. Thus we cannot use Theorem A, but need to use Theorem B. However, the difference is slight: we simply cannot omit the contributions of the constant vector K v(0). However, it will turn out that the joint spectral radius is 1, so the contribution will be absorbed by the error term anyway.
To see that the above holds, we have two different interpretations: the first is that the row vector w(n) = e q+1 A r 0 · · · A r −1 is the unit vector corresponding to the most significant digit of the q-ary expansion of n or, in view of the automaton A, corresponding to the final state. Note that we read the digit expansion from the least significant digit to the most significant one (although it would be possible the other way round as well). We have w(0) = e q+1 which corresponds to the empty word and being in the initial state I in the automaton. The vector v(0) corresponds to the fact that all states of A except 0 are accepting. The other interpretation is: the r th component of the column vector v(n) = A r 0 · · · A r −1 v(0) has the following two meanings: • In the automaton A, we start in state r and then read the digit expansion of n. The r th component is then the indicator function whether we remain esthetic, i.e., end in an accepting state. • To a word ending with r we append the digit expansion of n. The r th component is then the indicator function whether the result is an esthetic word.
At first glance, our problem here seems to be a special case of the transducers studied in Sect. 8. However, the automaton A is not complete. Adding a sink to have a formally complete automaton, however, adds an eigenvalue q and thus a much larger dominant asymptotic term, which would then be multiplied by 0. Therefore, the results of [31] do not apply to this case here.

Full Asymptotics
We now formulate our main result for the amount of esthetic numbers smaller than a given integer N . We abbreviate this amount by and have the following corollary.

Corollary G Fix an integer q ≥ 2. Then the number X (N ) of q-esthetic numbers smaller than N is
with 2-periodic continuous functions j . Moreover, we can effectively compute the Fourier coefficients of each j (as explained in Part IV). If q is even, then the functions j are actually 1-periodic. If q is odd, then the functions j for even j vanish.
If q = 2, then the corollary results in X (N ) = O(log N ). However, for each length, the only word of digits satisfying the esthetic number condition has alternating digits 0 and 1, starting with 1 at its most significant digit. The corresponding numbers n form the so-called Lichtenberg sequence [39, A000975].
Back to a general q: For the asymptotics, the main quantities influencing the growth are the eigenvalues of the matrix C = A 0 +· · ·+ A q−1 . Continuing our example q = 4 above, this matrix is Therefore it turns out that the growth of the main term is N log 4 ( √ 5+1)− 1 2 = N 0.347... , see Fig. 2. The first few Fourier coefficients are shown in Table 1.

Eigenvectors
Before proving Corollary G, we collect information on the eigenvalues of C.
The matrix C = A 0 + · · · + A q−1 has a block decomposition into for vectors 0 (vector of zeros) and 1 (vector of ones) of suitable dimension. Therefore, one eigenvalue of C is 0 and the others are the eigenvalues of M.
In contrast to [8,Sects. 4 and 5], we use the Chebyshev polynomials 8,9 U n of the second kind defined by for n ≥ 1. It is well-known that U n (cos ϕ) = sin((n + 1)ϕ) sin(ϕ) (9.2) and, as a consequence, the roots of U n are given by for n ≥ 1.
The following lemma is similar to [8, Proposition 3]. Proof We write ϕ := kπ/(q + 1). By Lemma 9.1 and (9.2) and a summation similar to the Dirichlet kernel, we have 8 Chebyshev polynomials are frequently occurring phenomena in lattice path analysis, see for instance [7,16]. We have such a lattice path here, so their appearance is not surprising. 9 Up to replacing 2X by X , the polynomials U n used here correspond to the polynomials p n used in [8]. .
As a consequence, 1, v = 0 if and only if k/2 is an integer.

Lemma 9.3 The characteristic polynomial of C is
In particular, all eigenvalues of M apart from 0 are eigenvalues of C with algebraic multiplicity 1. If q is even, then 0 has algebraic multiplicity 1 as an eigenvalue of C; if q is odd, then 0 has algebraic multiplicity 2 as an eigenvalue of C.
Proof The matrix C is a block lower triangular matrix, so the characteristic polynomial is the product of the characteristic polynomials of the matrices M and 0. The statement on the algebraic multiplicities follows from Lemma 9.1.
We can summarise our findings on the eigenvectors and eigenvalues of C as follows.

Proposition 9.4
Let v ∈ C q , w ∈ C, not both 0, and let λ ∈ C.
Then v w = 0 is an eigenvector of C to the eigenvalue λ if and only if one of the following conditions hold: 1. 0 = λ = 2 cos kπ q+1 for some 1 ≤ k ≤ q and k = q+1 2 , v is an eigenvector of M to λ, and w = 0 if k is even and λw = 1, v = 0 if k is odd; , v is an eigenvector of M and w = 0.

Proof The vector v w is an eigenvector if and only if
First assume that λ = 0. Then v = 0 leads to w = 0, contradiction. Therefore, v is an eigenvector of M to the eigenvalue λ and λ = 2 cos kπ q+1 for some 1 ≤ k ≤ q by Lemma 9.1. Then w = 0 if and only if k is even by Lemma 9.2. Now assume that λ = 0 and q is even. Then 0 is not an eigenvalue of M by Lemma 9.1. Thus v = 0 and w = 0. Now, assume that λ = 0 and q ≡ 3 (mod 4). Then λ = 2 cos π 2 = 2 cos q+1 2 π q+1 . By Lemma 9.2, the eigenvector v of M leads to an eigenvector v 0 of C; and there is an additional eigenvector 0 w = 0. Finally, assume that λ = 0 and q ≡ 1 (mod 4). In this case, by Lemma 9.2, it cannot be that v = 0 is an eigenvector of M because this would lead to 0 = 1, v = λw = 0, a contradiction. Thus the only eigenvector is 0 w = 0.

Proof of the Asymptotic Result
Proof of Corollary G We work out the conditions and parameters for using Theorem A.

Joint Spectral Radius.
As all the square matrices A 0 , . . . , A q−1 have a maximum absolute row sum norm equal to 1, the joint spectral radius of these matrices is bounded by 1. Let r ∈ {1, . . . , q − 1}. Then any product with alternating factors A r −1 and A r , i.e., a finite product A r −1 A r A r −1 · · · , has absolute row sum norm at least 1 as the word (r − 1)r (r − 1) . . . is q-esthetic. Therefore the joint spectral radius of A r −1 and A r is at least 1. Consequently, the joint spectral radius of A 0 , . . . , A q−1 equals 1.
We now assume that q is even. In this case, we still have to show that the functions j are actually 1-periodic. We now need to use Theorem B. Let w 1 , w 2 , . . . , w q−1 , w q be the rows of T where the order is chosen in such a way that J = diag 2 cos π q + 1 , . . . , 2 cos qπ q + 1 , 0 .
We write e q+1 = q k=1 c k w k for suitable c k ∈ R. Setting c := c 1 c 2 · · · c q , this means that e q+1 = cT , or equivalently, c = e q+1 T −1 . The columns of T −1 are the right eigenvectors of C described in Proposition 9.4. Then Proposition 9.4 (1) implies that c k = 0 for even k with 1 ≤ k ≤ q. This means that all fluctuations corresponding to eigenvalues 2 cos(kπ/(q + 1)) for even k with 1 ≤ k ≤ q are multiplied by 0 and do not contribute to the result. As |cos( q+1−k q+1 π)| = |cos( k q+1 π)|, but q + 1 − k and k The same argument can be used for the case of odd q, but in this case, q + 1 − k and k have the same parity. So Proposition 6.4 is used for odd k, and fluctuations to both eigenvalues 2 cos(kπ/(q + 1)) and 2 cos((q + 1 − k)π/(q + 1)) vanish for even k.
Fourier Coefficients. We can compute the Fourier coefficients according to Theorem A and Proposition 6.4; see also Part IV.

Pascal's Rhombus
We discuss the asymptotic analysis of odd entries in Pascal's rhombus; see also Sect. 5.5 for an overview.
We consider Pascal's rhombus R which is, for integers i ≥ 0 and j, the array with entries r i, j , where • r 0, j = 0 all j, • r 1,0 = 1 and r 1, j = 0 for all j = 0, • and We are interested in the number of odd entries in the first N rows of this rhombus. In [21] the authors investigate this quantity for N being a power of 2. We again aim for a more precise analysis and asymptotic description.
So, let X be equal to R but with entries taken modulo 2; see also Fig. 3. We partition X into the four sub-arrays • E consisting only of the rows and columns of X with even indices, i.e., the entries r 2i,2 j , • Y consisting only of the rows with odd indices and columns with even indices, i.e., the entries r 2i−1,2 j , • Z consisting only of the rows with even indices and columns with odd indices, i.e., the entries r 2i,2 j−1 , and • N consisting only of the rows and columns with odd indices, i.e., the entries r 2i−1,2 j−1 .

Recurrence Relations and 2-Regular Sequences
Let X (N ), Y (N ) and Z (N ) be the number of ones in the first N rows (starting with row index 1) of X, Y and Z, respectively. Goldwasser, Klostermeyer, Mays and Trapp [21, (12)- (14)] get the recurrence relations Figs. 2

and 3]). Distinguishing between even and odd indices gives
. These x(n), y(n) and z(n) are the number of ones in the nth row of X, Y and Z, respectively, and clearly holds. We obtain Let us write our coefficients as the vector v(n) = x(n), x(n + 1), y(n + 1), z(n), z(n + 1) .

(10.2)
It turns out that the components included into v(n) are sufficient for a self-contained linear representation of v(n). In particular, it is not necessary to include y(n). By using the recurrences (10.1), we find that for all 10  and with v(0) = (0, 1, 1, 0, 2) . Therefore, the sequences x(n), y(n) and z(n) are 2-regular.
Moreover, we can effectively compute the Fourier coefficients of (as explained in Part IV).
We get analogous results for the sequences Y (N ) and Z (N ) (each with its own periodic function , but the same exponent κ). The fluctuation of X (N ) is visualized in Fig. 4 and its first few Fourier coefficients are shown in Table 2.   All stated digits are correct; see also Part IV z(n) as well as x(n + 1) and z(n + 1); both pairs are asymptotically equal in the sense of (10.3). Therefore, we head for an only 3-dimensional functional system of equations for our Dirichlet series of x(n), y(n) and z(n) (instead of a 5-dimensional system).

Proof of (10.3) We use Theorem A.
Joint Spectral Radius. First we compute the joint spectral radius ρ of A 0 and A 1 . Both matrices have a maximum absolute row sum equal to 2, thus ρ ≤ 2, and both matrices have 2 as an eigenvalue. Therefore we obtain ρ = 2. Moreover, the finiteness property of the linear representation is satisfied by considering only products with exactly one matrix factor A 0 or A 1 . Thus, we have R = ρ = 2.
Asymptotics. By using Theorem A, we obtain an asymptotic formula for X (N − 1). Shifting from N − 1 to N does not change this asymptotic formula, as this shift is absorbed by the error term O (N log 2 N ).

Dirichlet Series and Meromorphic Continuation
In the lemma below, we provide the functional equation (10.4) as a system of three equations. This is in contrast to the generic functional equation provided by Theorem D which is a system of five equations.

Lemma 10.1 Set
with the notion of as in Lemma 6.3, provides meromorphic continuations of the Dirichlet series X n 0 (s), Y n 0 (s), and Z n 0 (s) for s > κ 0 = 1 with the only possible poles at κ + χ for ∈ Z, all of which are simple poles.
Proof We split the proof into several steps. 1 2 , Z n 0 ) . It is an entire function.
It is no surprise that the κ of this lemma and the κ in the proof of Corollary H which comes from the 2-linear representation of Sect. 10.1 coincide.
Meromorphic Continuation. Let D n 0 ∈ {X n 0 , Y n 0 , Z n 0 }. The Dirichlet series D n 0 (s) is analytic for s > 2 = log 2 ρ + 1 with ρ = 2 being the joint spectral radius by Theorem D. We use the functional equation (10.4) which provides the continuation, as we write D n 0 (s) in terms of J n 0 (s), K n 0 (s) and L n 0 (s). By Lemma 6.3, these three functions are analytic for s > 1.
The zeros (all are simple zeros) of the denominator (s) are the only possibilities for the poles of D n 0 (s) for s > 1.

Fourier Coefficients
We are now ready to prove the rest of Corollary H.

Proof of Corollary H We verify that we can apply Theorem E.
The steps of this proof in Sect. 10.2 provided us already with an asymptotic expansion (10.3). Lemma 10.1 gives us the meromorphic function for s > κ 0 = 1 which comes from the Dirichlet series X n 0 (s), Y n 0 (s), Z n 0 (s) . It can only have poles (all simple) at s = κ + χ for ∈ Z and satisfies the assumptions in Theorem E by Theorem D and Remark 6.2.
Therefore a computation of the Fourier coefficients via computing residues [see (3.6)] is possible by Theorem E, and this residue may be computed from (10.4) via Cramer's rule.
We refer to Part IV for details on the actual computation of the Fourier coefficients.

Part III: Proofs
Before reading this part on the collected proofs, it is recommended to recall the definitions and notations of Sect. 6.2. Some additional notations which are only used in the proofs are introduced in the following section.

Additional Notations
We use Iverson's convention [expr] = 1 if expr is true and 0 otherwise, which was popularised by Graham, Knuth and Patashnik [26]. We use the notation z := z(z − 1) · · · (z − + 1) for falling factorials. We use n k 1 ,...c,k r for multinomial coefficients. We sometimes write a binomial coefficient n a as n a,b with a + b = n when we want to emphasise the symmetry and analogy to a multinomial coefficient.

Decomposition into Periodic Fluctuations: Proof of Theorem B
We first give an overview over the proof.

Overview of the Proof of Theorem B
The first step will be to express the summatory function F in terms of the matrices C, B r and A r . Essentially, this corresponds to the fact that the summatory function of a q-regular function is again q-regular. This expression of F will consist of two terms: the first is a sum over 0 ≤ j < log q N involving a jth power of C and matrices B r and A r depending on the − j most significant digits of N . The second term is again a sum, but does not depend on the digits of N ; it only encodes the fact that f (0) = A 0 f (0) may not hold. The fact that we are interested in wF(N ) for the generalised left eigenvector w corresponding to the eigenvalue λ allows to express wC j in terms of wλ j (plus some other terms if w is not an eigenvector).
The second term can be disposed of by elementary observations using a geometric series. We reverse the order of summation in the first summand and extend it to an infinite sum. The infinite sum is written in terms of periodic fluctuations; the difference between the infinite sum and the finite sum is absorbed by the error term. In order not to have to deal with ambiguities due to non-unique q-ary expansions of real numbers, we define the fluctuations on an infinite product space instead of the unit interval.

Upper Bound for Eigenvalues of C
We start with an upper bound for the eigenvalues of C in terms of the joint spectral radius.
by (7.1). Taking th roots and the limit → ∞ yields |λ| ≤ q R. This last inequality does not depend on our particular (cf. Sect. 6.2) choice of R > ρ, so the inequality is valid for all R > ρ, and we get the result.

Explicit Expression for the Summatory Function
In this section, we give an explicit formula for F(N ) = 0≤n<N f (n) in terms of the matrices A r , B r and C.

Proof We claim that
holds for non-negative integers N and r with 0 ≤ r < q. We now prove (12.1): Using (6.1) and f (0) = I yields This concludes the proof of (12.1).

Proof of Theorem B For readability, this proof is split into several steps.
Setting. Before starting the actual proof, we introduce the setting using an infinite product space which will be used to define the fluctuations k . We also introduce the maps linking the infinite product space to the unit interval.
We will first introduce functions k defined on the infinite product space We equip it with the metric such that two elements x = x with a common prefix of length j and x j = x j have distance q − j . We consider the map lval : see Fig. 5. By using the assumption that the zeroth component of elements of is assumed to be non-zero, we easily check that lval is Lipschitz-continuous; i.e., for x = x with a common prefix of length j. For y ∈ [0, 1), let reprq(y) be the unique x ∈ with lval(x) = y such that x does not end on infinitely many digits q − 1, i.e., reprq(y) represents a q-ary expansion of q y . This means that lval • reprq is the identity on [0, 1).
From the definition of the metric on , recall that a function : → C d is continuous if and only if for each ε > 0, there is a j such that (x )− (x) < ε holds for all x and x that have a common prefix of length j. Further recall from the universal property of quotients that if such a continuous function satisfies (x) = (x ) Notation. We will deal with the two sums in Lemma 12.2 separately. We will first introduce notations corresponding to this split and to the eigenvector structure. Let N have the q-ary expansion r −1 . . . r 0 and set holds for all j ≥ 0. These vectors are suitable linear combinations of the vectors v 0 , . . . , v m−1 . We note that we have Second Summand. We will now rewrite wF 2 (N ) by evaluating the geometric sum and rewriting it in terms of a fluctuation. We claim that for suitable continuously differentiable functions (2) k on R, 0 ≤ k < m. If R = 0, then O(N log q R ) shall mean that the error vanishes for almost all N .
Consider first the case that λ = 1. Because of wC j = wC j and wT −1 DT = w (see Sect. 7.3) we have If λ = 0, then wC = 0 for almost all . We may set (2) k = 0 for 0 ≤ k < m and (12.5) is shown. Otherwise, as we have − 1 = log q N = log q N − {log q N } and by (12.3), we can rewrite wC as for reals L and ν, i.e., By the binomial theorem, we have This leads to a representation G 2 (L, ν) = 0≤k<m L k (2) k (ν) for continuously differentiable functions (2) for 0 ≤ k < m. As the functions (2) k are continuously differentiable, they are Lipschitz continuous on compact subsets of R. We note that in the case k = m − 1, the only occurring summand is for r = 0, which implies that Rewriting λ log q N as N log q λ and recalling that wϑ m = 0 yields (12.5) for λ = 1.
We now turn to the case λ = 1. We use wC j = 0≤k<m j k v k for j ≥ 0 as above. Thus where the identity [26, (5.10)] ("summation on the upper index") has been used in the last step. Thus wF 2 (N ) is a polynomial in of degree m. By writing = 1 + log q N − {log q N }, we can again rewrite this as a polynomial in log q N whose coefficients The additional factor T −1 (I − D)T in ϑ m has been introduced in order to annihilate generalised eigenvectors to other eigenvalues. By construction of K , we have wK = 0. Thus we have shown (12.5) for λ = 1, too.
Lifting the Second Summand. For later use-at this point, this may seem to be quite artificial-we set k is continuously differentiable, it is Lipschitz continuous on [0, 1]. As lval is also Lipschitz continuous, so is (2) k . First Summand We now turn to wF 1 (N ). To explain our plan, assume that w is in fact an eigenvector. Then wF 1 (N ) = 0≤ j< λ j w B r j A r j+1 . . . A r −1 . For |λ| ≤ R, it will be rather easy to see that the result holds. Otherwise, we will factor out λ and write the sum as wF 1 We will then reverse the order of summation and extend the sum to an infinite sum, which will be represented by periodic fluctuations. The difference between the finite and the infinite sums will be absorbed by the error term. The periodic fluctuations will be defined on the infinite product space .
We now return to the general case of a generalised eigenvector w and the actual proof. If λ = 0, we certainly have |λ| ≤ R and we are in one of the first two cases of this theorem. Furthermore, we have wC j = 0 for j ≥ m, thus by using (7.1). Together with (12.5), the result follows.
From now on, we may assume that λ = 0. By using (12.3), we have We first consider the case that |λ| < R [corresponding to Theorem B, (1)]. We get where (7.1) was used. Together with (12.5), the result follows. Next, we consider the case where |λ| = R [Theorem B, (2)]. In that case, we get Again, the result follows.
From now on, we may assume that |λ| > R. We set Q := |λ|/R and note that 1 < Q ≤ q by assumption and Lemma 12.1. We claim that there are continuous functions (1) k on for 0 ≤ k < m such that and such that when the first j entries of x and x ∈ coincide. Write N = q −1+{log q N } and let x = reprq({log q N }), i.e., x is the q-ary expansion of q {log q N } = N /q −1 ∈ [1, q) ending on infinitely many zeros. This means that x j = r −1− j for 0 ≤ j < and x j = 0 for j ≥ . Reversing the order of summation in (12.7) yields For j ≥ , we have x j = 0 and therefore B x j = 0. Thus we may extend the sum to run over all j ≥ 0, i.e., We insert − 1 = log q N − {log q N } and obtain for L ∈ R and x ∈ . Note that in contrast to G 2 , the second argument of G 1 is an element of instead of R. Collecting G 1 (L, x) by powers of L, we get which are continuously differentiable and therefore Lipschitz continuous on the unit interval. This shows (12.8). For k = m − 1, only summands with r = s = 0 occur, thus |λ| − j j m−1 R j according to (7.1). We now prove (12.9). So let x and x have a common prefix of length i. Consider the summand of (1) k (x) with index j. First consider the case that j < i. For all r , we have due to Lipschitz continuity of ψ kr • lval. As the matrix product A x j−1 . . . A x 0 is the same for x and x , the difference with respect to this summand is bounded by Thus the total contribution of all summands with , which leads to a total contribution of O(i m−1 Q −i ). Adding the two bounds leads to a bound of O(i m−1 Q −i ), as requested.

Descent. As we have defined the periodic fluctuations
(1) k on the infinite product space , we now need to prove that the periodic fluctuation descends to a periodic fluctuation on the unit interval. To do so, we will verify that the values of the fluctuation coincide whenever sequences in the infinite product space correspond to the same real number in the interval. By whenever x and x ∈ have a common prefix of length j. It remains to show that k (x) = k (x ) holds whenever lval(x) = lval(x ) or lval(x) = 0 and lval(x ) = 1.
Choose x and x such that one of the above two conditions on lval holds and such that x j = 0 for j ≥ j 0 and x j = q − 1 for j ≥ j 0 . Be aware that now the prefixes of x and x of length j 0 do not coincide except for the trivial case j 0 = 0.
We estimate n + 1 as n(1 + O(1/n)) and get (12.13) We have w f (n) = O(R j ) = O(R log q n ) = O(n log q R ) by (6.2) and (7.1). By (12.12), which is used below to replace x by x . Inserting these estimates in (12.13) and dividing by n log q λ yields (12.14) Note that k (x ) − k (x) does not depend on j. Now we let j (and therefore n) tend to infinity. We see that (12.14) can only remain true if k (x ) = k (x) for 0 ≤ k < m, which we had set out to show. Therefore, k descends to a continuous function k on [0, 1] with k (0) = k (1); thus k can be extended to a 1-periodic continuous function.
Hölder Continuity. We will now prove Hölder continuity. As the fluctuations have been defined on the infinite product space , we will basically have to prove Hölder continuity there. The difficulty will be that points in the unit interval which are close to each other there may have drastically different q-ary expansions, thus correspond to drastically different points in the infinite product space . To circumvent this problem, the interval between the two points will be split into two parts.
We first claim that for 0 ≤ y < y < 1, we have as y → y. To prove this, let x := reprq(y) and x := reprq(y ). Let be the length of the longest common prefix of x and x and choose j ≥ 0 such that q − j ≤ q y − q y < q − j+1 . We define x and x ∈ such that and set y := lval(x ) and y := lval(x ). As lval(x) = y < y = lval(x ), we have x > x . We conclude that y ≤ y = y ≤ y . Therefore, q y − q y ≤ q y − q y < q − j+1 , so in view of the fact that each entry of x is greater or equal than the corresponding entry of x, the expansions x and x must have a common prefix of length j. Similarly, the expansions x and x must have a common prefix of length j. Thus (12.12) implies that Noting that − j = log q (q y − q y ) + O(1) leads to (12.15).
In order to prove Hölder continuity with exponent α < log q Q, we first note that Lipschitz-continuity of y → q y on the interval [0, 1] shows that (12.15) implies This can then easily be extended to arbitrary reals y < y by periodicity of k because it is sufficient to consider small y − y and the interval may be subdivided at an integer between y and y .

Constant Dominant Fluctuation.
To finally prove the final assertion on constant fluctuations, we will have to inspect the explicit expression for the fluctuations using the additional assumption.
Under the additional assumption that the vector w(C − I ) m−1 = v m−1 is a left eigenvector to all matrices A 0 , . . . , A q−1 associated with the eigenvalue 1, the same holds for v m−1 by (12.4). Then v m−1 is also a left eigenvector of C associated with the eigenvalue q. In particular, λ = q = 1.

Proof of Theorem C
Proof of Theorem C We denote the rows of T as w 1 , . . . , w d and the columns of T −1 by t 1 , . . . , t d . Thus 1≤ j≤d t j w j = I and w j is a generalised left eigenvector of C of some rank m j corresponding to some eigenvalue λ j ∈ σ (C). Theorem B and the fact that there are no eigenvalues of C of absolute value between ρ and R then immediately imply that for some 1-periodic Hölder continuous functions jk with exponent less than log q |λ j |/R. The first summand K as well as the error term already coincide with the result stated in the theorem. From Sect. 7.3 we recall that w j ϑ m j = 0 for λ j = 1. We set Then we still have to account for (1) . (12.16) The factor (C − I ) m(1)−1 in the definition of ϑ m (1) implies that w j ϑ m(1) vanishes unless λ j = 1 and m j = m (1). Therefore, the sum in (12.16) equals ϑ.

Meromorphic Continuation of the Dirichlet Series: Proof of Theorem D
For future use, we state an estimate for the binomial coefficient. Unsurprisingly, it is a consequence of a suitable version of Stirling's formula. Alternatively, it can be seen as the most basic case of Flajolet and Odlyzko's singularity analysis [19,Proposition 1], where uniformity in s is easily checked.
uniformly for s in a compact subset of C and k → ∞.
Proof of Lemma 6. 3 We have for s > log q R + 1. We note that Therefore, and the series converges for s > log q R . As this holds for all R > ρ, we obtain (s, β, D) = O(| s|) as | s| → ∞ uniformly for log q ρ + δ ≤ s ≤ log q ρ + δ + 1. In the language of [27, § 3.3], (s, β, D) has order at most 1 for log q ρ + δ ≤ s ≤ log q ρ + δ + 1. As log q ρ + δ + 1 is larger than the abscissa of absolute convergence of (s, β, D), it is clear that (s, β, D) = O(1) for s = log q ρ + δ + 1, i.e., (s, β, D) has order at most 0 for s = log q ρ + δ + 1. By Lindelöf's theorem (see [27,Theorem 14]), we conclude that (s, β, D) = O | s| μ δ ( s) for log q ρ + δ ≤ s ≤ log q ρ + δ + 1. For s > log q R + 1, we may rewrite (13.2) using the binomial series as Switching the order of summation was legitimate because for s + k > log q R + 1 and Lemma 13.1 imply absolute and uniform convergence for s in a compact set. Noting that the previous arguments hold again for all R > ρ and that the inner sum in (13.3) is D(s + k) completes the proof.
Proof of Theorem D As f (n) = O(R log q n ) = O(n log q R ) by (6.2) and (7.1), the Dirichlet series F n 0 (s) = n≥n 0 n −s f (n) (see Sect. 6.2) converges absolutely and uniformly on compact sets for s > log q R + 1. As this holds for all R > ρ, i.e., does not depend on our particular (cf. Sect. 6.2) choice of R > ρ, this convergence result holds for s > log q ρ + 1. We use (6.1) and Lemma 6.3 (including its notation) to rewrite F n 0 as for s > log q R+1. By Lemma 6.3 we have H n 0 (s) = O | s| μ δ ( s) for log q ρ +δ ≤ s ≤ log q ρ + δ + 1. Rewriting the expression for H n 0 (s) using the binomial series (see Lemma 6.3 again) yields Combining this with (13.4) yields the expression (6.4) for G n 0 . Solving (6.3) for F n 0 yields the meromorphic continuation of F n 0 (s) to s > log q R (and thus to s > log q ρ) with possible poles where q s is an eigenvalue of C. As long as q s keeps a fixed positive distance δ from the eigenvalues, the bound for G n 0 (coming from the bound for H n 0 ) carries over to a bound for F n 0 , i.e., (6.5).
To estimate the order of the poles, let w be generalised left eigenvector of rank m of C corresponding to an eigenvalue λ with |λ| > R. We claim that wF n 0 (s) has a pole of order at most m at s = log q λ + χ k and no other poles for s > log q R. We prove this by induction on m.
Set v := w(C − λI ). By definition, v = 0 or v is a generalised eigenvector of rank m − 1 of C. By induction hypothesis, vF n 0 (s) has a pole of order at most m − 1 at s = log q λ + χ k for k ∈ Z and no other poles for s > log q R.
The right-hand side has a pole of order at most m − 1 at log q λ + χ k for k ∈ Z and 1 − q −s λ has a simple zero at the same places. This proves the claim.

Fourier Coefficients: Proof of Theorem E
In contrast to the rest of this paper, this section does not directly relate to a regular sequence but gives a general method to derive Fourier coefficients of fluctuations.

Pseudo-Tauberian Theorem
In this section, we generalise the pseudo-Tauberian argument by Flajolet, Grabner, Kirschenhofer, Prodinger and Tichy [18,Proposition 6.4]. The basic idea is that for a 1-periodic Hölder-continuous function and γ ∈ C, there is a 1-periodic continuously differentiable function such that 1≤n<N n γ (log q n) = N γ +1 (log q N ) + o(N γ +1 ), and there is a straight-forward relation between the Fourier coefficients of and the Fourier coefficients of . This relation exactly corresponds to the additional factor s +1 when transitioning from the zeroth order Mellin-Perron formula to the first order Mellin-Perron formula.
In contrast to [18,Proposition 6.4], we allow for an additional logarithmic factor, have weaker growth conditions on the Dirichlet series and quantify the error. We also extend the result to all complex γ . The generalisation from q = 2 there to our real q > 1 is trivial.
Denote the Fourier coefficients of j and j by ϕ j := for ∈ Z and Z → 0. If q γ +1 = 1, then −1 vanishes.

Remark 14.2
Note that the constant c is absorbed by the error term if γ + 1 > α, in particular if γ > 0. Therefore, this constant does not occur in the article [18].

Remark 14.3
The factor γ + 1 + 2 πi log q + Z in (14.2) will turn out to correspond exactly to the additional factor s + 1 in the first order Mellin-Perron summation formula with the substitution s = γ + 2 πi log q + Z such that the local expansion around the pole in s = γ + 2 πi log q of the Dirichlet generating function is conveniently written as a Laurent series in Z . See the proof of Theorem E for details.
Before actually proving Proposition 14.1, we give an outline.
Overview of the Proof of Proposition 14. 1 We start with the left-hand side of (14.1) and split the range of summation according to log q n , thereby, in terms of our periodic functions, split after each period. We then use periodicity of the j and collect terms. This results in Riemann sums which converge to the corresponding integrals. Therefore, we can approximate these sums by the integrals.
More rewriting constructs and reveals the functions j [of the right-hand side of (14.1)]: these functions are basically defined via the above mentioned integral. We then show that these functions are indeed periodic and that their Fourier coefficients relate to the Fourier coefficients of the j . The latter is done by a direct computation of the integrals defining these coefficients.
For this proof, we use an approach via exponential generating functions. This reduces the overhead for dealing with the logarithmic factors (log n) k in (14.1) such that we can essentially focus on the case m = 1. The resulting formula (14.1) follows by extracting a suitable coefficient of this power series.
There is another benefit of the generating function approach: this formulation allows to easily translate the relation between the Fourier coefficients here to the additional factors occurring when transitioning to higher order Mellin-Perron summation formulae, in particular the factor s + 1 in the first order Mellin-Perron summation.
Proof of Proposition 14. 1 We split the proof into six parts.
Notations. We start by defining quantities that are used through the whole proof.
Without loss of generality, we assume that q γ +1 = q α : otherwise, we slightly decrease α keeping the inequality β < α intact. We use the abbreviations := log q N , ν := {log q N }, i.e., N = q +ν . We use the generating functions for 0 ≤ u ≤ 1 and 0 < |Z | < 2r where r > 0 is chosen such that r < (α − β)/2 and such that Q(Z ) = 1 and |Q(Z )| = q α for these Z . (The condition Z = 0 is only needed for the case q 1+γ = 1.) We will stick to the above choice of r and restrictions for Z throughout the proof.
It is easily seen that the left-hand side of (14 Approximation of the Sum by an Integral. We will now rewrite L(N , Z ) so that its shape is that of a Riemann sum, therefore enabling us to approximate it by an integral.
Splitting the range of summation with respect to powers of q yields L(N , Z ) = 0≤ p< q p ≤n<q p+1 n γ +Z (log q n, Z ) + q ≤n<q +ν n γ +Z (log q n, Z ).
We write n = q p x (or n = q x for the second sum), use the periodicity of in u and get The inner sums are Riemann sums converging to the corresponding integrals for p → ∞. We set It will be convenient to change variables x = q w in I (u, Z ) to get We define the error ε p (u, Z ) by As the sum and the integral are both analytic in Z , their difference ε p (u, Z ) is analytic in Z , too. We bound ε p (u, Z ) by the difference of upper and lower Darboux sums (step size q − p ) corresponding to the integral I (u, Z ): On each interval of length q − p , the maximum and minimum of a Hölder continuous function can differ by at most O(q −α p ). As the integration interval as well as the range for u and Z are finite, this translates to the bound ε p (u, Z ) = O(q −α p ) as p → ∞ uniformly in 0 ≤ u ≤ 1 and |Z | < 2r . This results in If |Q(Z )|/q α = q γ +1+ Z −α < 1, i.e., γ + Z < α − 1, the second sum involving the integration error converges absolutely and uniformly in Z for → ∞ to some analytic function c (Z ); therefore, we can replace the second sum by c (Z ) + O q ( γ +1+2r −α) = c (Z ) + O N γ +1+2r −α in this case. If γ + Z > α − 1, then the second sum is O q ( γ +2r +1−α) = O N γ +1+2r −α . By our choice of r , the case γ + Z = α − 1 cannot occur. So in any case, we may write the second sum as c (Z ) + O N γ +1−β by our choice of r . The last summand involving ε (ν, Z ) is absorbed by the error term of the second summand. Note that the error term is uniform in Z and, by its construction, analytic in Z .
Thus we end up with It remains to rewrite S(N , Z ) in the form required by (14.1). We emphasise that we will compute S(N , Z ) exactly, i.e., no more asymptotics for N → ∞ will play any rôle.
Construction of . We will now rewrite the expression S(N , Z ) such that the generating function [i.e., the fluctuations of the right-hand side of (14.1)] appears. After this, we will gather properties of including properties of its Fourier coefficients.
We rewrite (14.5) as We replace by log q N − ν and use .
Periodic Extension of . A priori, it is not clear that the function (u, Z ) defined above can be extended to a periodic function (and therefore Fourier coefficients can be computed later on). The aim now is to show that it is possible to do so. It is obvious that (u, Z ) is continuously differentiable in u ∈ [0, 1]. We have because I (0, Z ) = 0 by (14.3). The derivative of (u, Z ) with respect to u is which implies that We can therefore extend (u, Z ) to a 1-periodic continuously differentiable function in u on R.

Fourier Coefficients of
Knowing that is a periodic function, we can now head for its Fourier coefficients and relate them to those of . By using equations (14.7) and (14.3), Q(Z ) = q γ +1+Z , and exp(− 2 πiu) = q −χ u with χ = 2πi log q , we now express the Fourier coefficients of (u, Z ) in terms of those of (u, Z ) by .
The second and third summands cancel, and we get Extracting Coefficients. So far, we have proven everything in terms of generating functions. We now extract the coefficients of these power series which will give us the result claimed in Proposition 14.1. By (14.7), (u, Z ) is analytic in Z for 0 < |Z | < 2r . If q γ +1 = 1, then it is analytic in Z = 0, too. If q γ +1 = 1, then (14.7) implies that (u, Z ) might have a simple pole in Z = 0. Note that all other possible poles have been excluded by our choice of r . For j ≥ −1, we write and use Cauchy's formula to obtain This and the properties of (u, Z ) established above imply that j is a 1-periodic continuously differentiable function.
Inserting (14.6) in (14.4) and extracting the coefficient of Z m−1 using Cauchy's theorem and the analyticity of the error in Z yields (14.1) Rewriting (14.8) in terms of j and j leads to (14.2). Note that we have to add O(Z m ) in (14.2) to compensate the fact that we do not include ψ j for j ≥ m.
We prove a uniqueness result.
Proof If γ < 0 and c = 0, then (14.9) is impossible as the growth of the right-hand side of the equation is larger than that on the left-hand side. So we can exclude this case from further consideration. We proceed indirectly and choose k maximally such that k = k . Dividing (14.9) by (log q N ) k yields for N → ∞. Let 0 < u < 1 and set N j = q j+u . We clearly have lim j→∞ N j = ∞. Then We define ν j := log q N j − j − u and see that ν j = O(q − j ) for j → ∞, i.e., lim j→∞ ν j = 0. This implies that lim j→∞ {log q N j } = u and therefore Setting N = N j in (14.10) and letting j → ∞ shows that If k = 0 or γ > 0, we immediately conclude that k (u) − k (u) = 0. If γ < 0 we have c = 0, which again implies that k (u) − k (u) = 0. Now we assume that γ = 0 and k = 0. We set β := − log q 2πi γ , which implies that N −γ = exp(2πiβ log q N ). We choose sequences (r ) ≥1 and (s ) ≥1 such that lim →∞ s = ∞ and lim →∞ |s β − r | = 0: For rational β = r /s, we simply take r = r and s = s, and for irrational β, we consider the sequence of convergents (r /s ) ≥1 of the continued fraction of β and the required properties follow from the theory of continued fractions; see for example [28,Theorems 155 and 164]. By using log q N j = j + u + ν j , we get These two limits are distinct as β / ∈ Z by assumption. Thus lim j→∞ N −γ j does not exist. Therefore, (14.11) implies that c = 0 and therefore k (u) − k (u) = 0.
We proved that k (u) = k (u) for u / ∈ Z. By continuity, this also follows for all u ∈ R; contradiction.

Proof of Theorem E
We again start with an outline of the proof.

Overview of the Proof of Theorem E
The idea is to compute the repeated summatory function of F twice: On the one hand, we use the pseudo-Tauberian Proposition 14.1 to rewrite the right-hand side of (6.6) in terms of periodic functions aj . On the other hand, we compute it using a higher order Mellin-Perron summation formula, relating it to the singularities of F. More specifically, the expansions at the singularities of F give the Fourier coefficients of aj . The Fourier coefficients of the functions aj are related to those of the functions j via (14.2).
And up next comes the actual proof.

Proof of Theorem E Initial observations and notations.
As j is Hölder continuous, its Fourier series converges by Dini's criterion; see, for example, [40, p. 52].
Asymptotic Summation. We first compute the Ath repeated summatory function S A F of F (i.e., the (A + 1)th repeated summatory function S A+1 f of the function f ) by applying Proposition 14.1 A times. This results in an asymptotic expansion involving new periodic fluctuations while keeping track of the relation between the Fourier coefficients of the original fluctuations and those of the new fluctuations.
A simple induction based on (6.6) and using Proposition 14.1 shows that there exist 1-periodic continuous functions aj for a ≥ 0 and − 1 ≤ j < m and some constants c ab for 0 ≤ b < a such that (14.12) for integers N → ∞. In fact, 0 j = j for 0 ≤ j < m. For a ≥ 1 and − 1 ≤ j < m, aj is continuously differentiable. Note that the case that q γ +a+1 = 1 occurs for at most one 0 ≤ a < A, which implies that the number of non-vanishing fluctuations increases at most once in the application of Proposition 14.1. Also note that the assumption α > γ − γ 0 implies that the error terms arising in the application of Proposition 14.1 are absorbed by the error term stemming from (6.6). We denote the corresponding Fourier coefficients by for 0 ≤ a < A, ∈ Z and Z → 0. Iterating this recurrence yields 0≤ j<m for ∈ Z and Z → 0.
Explicit Summation. We now compute S A+1 f explicitly with the aim of decomposing it into one part which can be computed by the Ath order Mellin-Perron summation formula and another part which is smaller and can be absorbed by an error term. Explicitly, we have Note that we formally write the outer sum over the range 1 ≤ n < N although the inner sum is empty (i.e., equals 0) for n ≥ N − a; this will be useful later on. The inner sum counts the number of selections of a elements out of {n + 1, . . . , N − 1}, thus we have for 0 ≤ a ≤ A and falling factorials z a := z(z − 1) · · · (z − a + 1). The polynomials 1 a! (U − 1) a , 0 ≤ a ≤ A, are clearly a basis of the space of polynomials in U of degree at most A. Thus, there exist rational numbers b 0 , . . . , b A such that Comparing the coefficients of U A shows that b A = 1. Substitution of U by N − n, multiplication by f (n) and summation over 1 ≤ n < N yield by (14.14). When inserting the asymptotic expressions from (14.12), the summands involving fluctuations for 0 ≤ a < A are absorbed by the error term O(N γ 0 +A ) of the summand for a = A because γ − γ 0 < 1. Thus there are some constants for integers N → ∞. If γ + A = b + χ for some 0 ≤ b < A and ∈ Z, then we assume without loss of generality that c b = 0: Otherwise, we replace A(m−1) (u) by A(m−1) (u) + c b exp(− 2 πiu) and c b by 0. Both (14.15) and (14.13) remain intact: the former trivially, the latter because the factor for a = A−b in (14.13) equals γ + A−b−χ + Z = Z which compensates the fact that the Fourier coefficient ψ A(m−1)(− ) is modified.
Mellin-Perron summation. We use the Ath order Mellin-Perron summation formula to write the main contribution of S A+1 f as determined above in terms of new periodic fluctuations j whose Fourier coefficients are expressed in terms of residues of a suitably modified version of the Dirichlet generating function F.
Without loss of generality, we assume that σ abs > 0: the growth condition (6.8) trivially holds with η = 0 on the right of the abscissa of absolute convergence of the Dirichlet series. By the Ath order Mellin-Perron summation formula (see [18,Theorem 2.1]), we have with the arbitrary choice σ abs + 1 > σ abs for the real part of the line of integration.
The growth condition (6.8) allows us to shift the line of integration to the left such that ds.
The summand for a in the second term corresponds to a possible pole at s = −a which is not taken care of in the first sum; note that F(s) is analytic at s = −a in this case by assumption because of γ 0 < −a.
We now compute the residue at s = γ + χ . We use to split up the residue as Res F(s)N s+A s(s + 1) · · · (s + A) , s = γ + χ = N γ +A+χ for j ≥ −1. Note that we allow j = −1 for the case of γ ∈ −a + 2πi log q Z for some 1 ≤ a ≤ A when F(s)/ s · · · (s + A) might have a pole of order m + 1 at s = −a. Using the growth condition (6.8) and the choice of A yields for | s| → ∞ and s which are at least a distance δ away from the poles γ + χ . By writing the residue in (14.16) in terms of an integral over a rectangle around s = γ +χ (distance again at least δ away from γ + χ ), we see that (14.17) implies Thus we proved that where the ξ j are given in (14.16). By (14.18), the Fourier series (14.20) converges uniformly and absolutely. This implies that j is a 1-periodic continuous function.
Fourier Coefficients. We will now compare the two asymptotic expressions for S A+1 f obtained so far to see that the fluctations coincide. We know explicit expressions for the Fourier coefficients of the j in terms of residues, and we know how the Fourier coefficients of the fluctuations of the repeated summatory function are related to the Fourier coefficients of the fluctuations of F. Therefore, we are able to compute the latter. By (14.15), (14.19), elementary asymptotic considerations for the terms N b with b > γ + A, Lemma 14.4 and the fact that c b = 0 if b ∈ γ + A + 2πi log q Z for some 0 ≤ b < A, we see that j = A j for − 1 ≤ j < m. This immediately implies that F(0) = 0 if γ 0 < 0 and γ / ∈ 2πi log q Z. To compute the Fourier coefficients ψ A j = ξ j , we set Z := s − γ − χ to rewrite (14.16) using (6.7) as for − 1 ≤ j < m and ∈ Z. This is equivalent to −1≤ j<m ψ A j Z j = j≥0 ϕ j Z j 1≤a≤A (γ + a + χ + Z ) for ∈ Z and Z → 0. Clearing the denominator and using (14.13)  for ∈ Z and Z → 0. Comparing coefficients shows that ψ 0 j = ϕ j for 0 ≤ j < m and ∈ Z. This proves (6.9). For computing the Fourier coefficients, we denote the rows of T by w 1 , . . . , w d . Thus w a is a generalised left eigenvector of C of some order m a associated to some eigenvalue λ a of C. We can write e 1 = 1≤a≤d c a w a for some suitable constants The reason for incorporating v(0) into the value for n = 1 is that the corresponding Dirichlet series H (a) (s) := n≥1 n −s h a (n) only takes values at n ≥ 1 into account. By definition, we have H (a) (s) = w a v(0) + w a V(s). Taking the linear combination yields 1≤a≤d c a H (a) (s) = x(0) + X (s). We choose γ 0 > log q R such that there are no eigenvalues λ ∈ σ (C) with log q R < log q λ ≤ γ 0 and such that γ 0 / ∈ Z ≤0 . By Theorem B, we have for N → ∞ for suitable 1-periodic Hölder continuous functions ak (which vanish if |λ a | ≤ R). By Theorem D, the Dirichlet series H (a) (s) is meromorphic for s > γ 0 with possible poles at s = log q λ a + χ for ∈ Z.

Proof of Theorem
The sequence h a satisfies the prerequisites of Theorem E, either with γ = log q λ a if log q λ a > γ 0 or with arbitrary real γ > γ 0 and j = 0 for all j if log q λ a ≤ γ 0 . The theorem then implies that H (a) (0) = 0 (15.2) if γ 0 < 0 and λ a = 1.
If |λ a | > R, Theorem E also yields where the ψ ak are given by the singular expansion for s > γ 0 . Note that (15.2) ensures that there is no additional pole at s = 0 when γ 0 < 0 and λ a = 1. Also note that in comparison to Theorem E, m a −1−k there corresponds to ak here. We now have to relate the results obtained for the sequences h a with the results claimed for the original sequence f . For λ ∈ σ (C) with |λ| > R, we have λk (u) = 1≤a≤d λ a =λ c a ak (u).
We denote the Fourier coefficients of λk by ϕ λk for ∈ Z and will show that these Fourier coefficients actually fulfil (3.5). Taking linear combinations of (15.3) shows that Summing over all λ ∈ σ (C) yields (3.5) because summands λ with |λ| ≤ R are analytic for s > γ 0 and do therefore not contribute to the right-hand side.
It might seem to be somewhat artificial that Theorem E is used to prove that H ( j) (0) = 0 in some of the cases above. In fact, this can also be shown directly using the linear representation; we formulate and prove this in the following remark. As R < 1 and λ j = 1, the Dirichlet series H ( j) (s) is analytic in s = 0 by Theorem D. It is therefore legitimate to set s = 0 in the above equation. We use the induction hypothesis that H ( j+1) (0) = 0 as well as the fact that v(n) = A n v(0) (note that v(0) is a right eigenvector of A 0 to the eigenvalue 1; see Sect. 7.1) for 0 ≤ n < q to get
We set thus is a p-periodic function.
For the Fourier series expansion, we get Replacing p + j by leads to the Fourier series claimed in the proposition.

Part IV: Computational Aspects
The basic idea for computing the Fourier coefficients is to use the functional equation in Theorem D. This part describes in detail how this is done. We basically follow an approach found in Grabner and Hwang [25] and Grabner and Heuberger [23], but provide error bounds. An actual implementation is also available; SageMath [38] code can be found at https://gitlab.com/dakrenn/regular-sequence-fluctuations . We use the Arb library [33] (more precisely, its SageMath bindings) for ball arithmetic which keeps track of rounding errors such that we can be sure about the precision and accuracy of our results.
We use the results of this part to compute Fourier coefficients for our examples, in particular for esthetic numbers (Sect. 9) and Pascal's rhombus (Sect. 10).

Strategy for Computing the Fourier Coefficients
The computation of the Fourier coefficients relies on the evaluation of Dirichlet series at certain points s = s 0 . It turns out to be numerically preferable to split up the sum as for some suitable n 0 (see Sect. 18.2), compute the sum of the first n 0 − 1 summands directly and evaluate F n 0 (s 0 ) as it is described in the following.
For actually computing the Fourier coefficients, we use a formulation in terms of a residue; for instance, see (3.6) where this is formulated explicitly in the set-up of Theorem A. As said, we will make use of the functional equation (6.3) for the matrixvalued Dirichlet series F n 0 (s) with its right-hand side, the matrix-valued Dirichlet series G n 0 (s).
Let us make this explicit for a simple eigenvalue λ = 1 of C and a corresponding eigenvector w. Then w(I − q −s C) = w(1 − q −s λ) and (6.3) can be rewritten as Thus, w F 1 (s) has simple poles at s = log q λ + χ for all ∈ Z, where χ = 2 πi log q . By (6.7) and (6.9) of Theorem E (with κ = log q λ and m = 1), the th Fourier coefficient is given by the residue
Note that log q is the derivative of 1 − q −s λ with respect to s evaluated at the pole s = log q λ. By (6.4), G n 0 (log q λ + χ ) is expressed in terms of an infinite sum containing F n 0 (log q λ + χ + k) for k ≥ 1. We truncate this sum and bound the error; this is the aim of Sect. 18.1 and in particular Lemma 18.2. We can iterate the above idea for the shifted Dirichlet series F n 0 (log q λ + χ + k) which leads to a recursive evaluation scheme. Note that once we have computed G n 0 (log q λ + χ + k), we get F n 0 (log q λ + χ + k) by solving a system of linear equations.

Bounding the Error
We need to estimate the approximation error which arises if the infinite sum over k ≥ 1 in (6.4) is replaced by a finite sum. It is clear that for large s and n 0 , the value F n 0 (s) will approximately be of the size of its first summand n −s 0 f (n 0 ). In view of f (n 0 ) = O(ρ log q n 0 ), this will be rather small. We give a precise estimate in a first lemma. For 0 ≤ x < 1, Taylor's theorem (or induction on K ≥ 1 using integration by parts) implies that For 0 ≤ t ≤ x < 1, we can bound |(1 + t) −s−K | from above by 1 since we have assumed that s + K > 0. Thus Thus we obtain the bound Bounding the remaining Dirichlet series by Lemma 18.1 yields the result.

Choices of Parameters
As mentioned at the beginning of this part, we choose the Arb library [33] for reliable numerical ball arithmetic. In our examples (esthetic numbers in Sect. 9 and Pascal's rhombus in Sect. 10), we choose n 0 = 1024 and recursively compute F n 0 (log q λ + χ + k) for k ≥ 1 by (6.4). In each step, we keep adding summands for k ≥ 1 until the bound of the approximation error in Lemma 18.2 is smaller than the smallest increment which can still be represented with the chosen number of bits. For plotting the graphs, we simply took machine precision; for the larger number of significant digits in Table 2, we used 128 bits precision.

Non-vanishing Coefficients
Using reliable numerical arithmetic for the computations (see above) yields small balls in which the true value of the Fourier coefficients is. If such a ball does not contain zero, we know that the Fourier coefficient does not vanish. If the ball contains zero, however, we cannot decide whether the Fourier coefficient vanishes. We can only repeat the computation with higher precision and hope that this will lead to a decision that the coefficient does not vanish, or we can try to find a direct argument why the Fourier coefficient does indeed vanish, for instance using the final statement of Theorem B (3).
Vanishing Fourier coefficients appear in our introductory Example 3.1: In its continuation (Example 3.3) an alternative approach is used to compute these coefficients explicitly symbolically. In this way a decision for them being zero is possible. The same is true for the example of transducers in Sect. 8.
It should also be noted that in the analysis of esthetic numbers (example in Sect. 9) we could have modelled the problem by a complete transducer (by just introducing a sink) and then applied the results of Sect. 8. This would have led to an asymptotic expansion where the fluctuations of the main term (corresponding to the eigenvalue q) would in fact have vanished, but an argument would have been needed. So we chose a different approach in Sect. 9 to avoid this problem. There the eigenvalue q does no longer occur. This implies that the fluctuations for q of the transducer approach vanish. Note also that half of the remaining fluctuations still turn out to vanish: this is shown in the proof of Corollary G.