Séminaire Lotharingien de Combinatoire 54 (2006), Article B54c

FUCHSIAN DIFFERENTIAL EQUATION FOR THE PERIMETER GENERATING FUNCTION OF THREE-CHOICE POLYGONS

Anthony J. Guttmann and Iwan Jensen

ARC Centre of Excellence for Mathematics and Statistics of Complex Systems Department of Mathematics and Statistics, The University of Melbourne, Victoria 3010, Australia

Abstract. Using a simple transfer matrix approach we have derived very long series expansions for the perimeter generating function of three-choice polygons. We find that all the terms in the generating function can be reproduced from a linear Fuchsian differential equation of order 8. We perform an analysis of the properties of the differential equation.

1. Introduction

A well-known long standing problem in combinatorics and statistical mechanics is to find the generating function for self-avoiding polygons (or walks) on a two-dimensional lattice, enumerated by perimeter. Recently, we have gained a greater understanding of the difficulty of this problem, as Rechnitzer [13] has proved that the (anisotropic) generating function for square lattice self-avoiding polygons is not differentiably finite [14], as had been previously conjectured, on numerical grounds [6], but not proved. That is to say, it cannot be expressed as the solution of an ordinary differential equation with polynomial coefficients. There are many simplifications of this problem that are solvable [1], but all the simpler models impose an effective directedness or other constraint that reduces the problem, in essence, to a one-dimensional problem.

One model, that of so-called three-choice polygons, has remained unsolved despite the knowledge that its solution must be D-finite. In this paper we report on recent numerical work resulting in an exact differential equation apparently satisfied by the perimeter generating function of three-choice polygons. While our results do not constitute a rigorous mathematical proof the numerical evidence is compelling.

Three-choice self-avoiding walks on the square lattice, ℤ^{2}, were introduced by Manna [12] and can be defined as follows: Starting from the origin one can step in any direction; after a step upward or downward one can head in any direction (except backward); after a step to the left one can only step forward or head downward, and similarly after a step to the right one can continue forward or turn upward. Alternatively put, one cannot make a right-hand turn after a horizontal step. Whittington [16] showed that the growth constant for three-choice walks is exactly 2, so that if w_{n} denotes the number of such walks of n steps on an infinite lattice, equivalent up to a translation, then w_{n} ∼ 2^{n+o(n)}. It is perhaps surprising that the best known result for the sub-dominant term is 2^{o(n)} but attempts to improve on this have not been successful. Even numerically, there is no firmly based Figure 1. Examples of the two types of three-choice polygons. In the right panel we indicate the origin (O) and the direction of the first step (note that rotation by 180 degrees also leads to a valid three-choice polygon).

conjecture for the sub-dominant term, unlike for ordinary self-avoiding walks, for which the sub-dominant term is widely believed to be O(n^{g}).

As usual one can define a polygon version of the walk model by requiring the walk to return to the origin. So a three-choice polygon [8] is simply a three-choice self-avoiding walk which returns to the origin, but has no other self-intersections. There are two distinct classes of three-choice polygons. The three-choice rule either leads to staircase polygons or imperfect staircase polygons [3] as illustrated in figure 1. In the case of staircase polygons any vertex on the perimeter can act as the origin of the three-choice walk (which then proceeds counter-clockwise), while for imperfect staircase polygons there is only one possible origin but the polygon could be rotated by 180 degrees. If we denote by t_{n} the number of three-choice polygons with perimeter 2n then, t_{n} = 2nc_{n} + 2p_{n}, where c_{n} is the number of staircase polygons with perimeter 2n, and p_{n} is the number of imperfect staircase polygons of perimeter 2n. Note that t_{n}, p_{n} and c_{n} all grow like 4^{n} and in particular we recall the well-known result that c_{n+1} = C_{n} = \frac{1}{n+1} \binom{2n}{n} are given by the Catalan numbers C_{n}.

In this paper we report on recent work which has led to an exact Fuchsian [9] linear differential equation of order 8 apparently satisfied by the perimeter generating function, T (x) = ∑_{n≥0} t_{n}x^{n}, for three-choice polygons (that is T (x) is conjectured to be one of the solutions of the ODE, expanded around the origin). The first few terms in the generating function are

T (x) = 4x^{2} + 12x^{3} + 42x^{4} + 152x^{5} + 562x^{6} + · · · .

The generating function for the coefficients p_{n} is no simpler.

If we distinguish between steps in the x and y direction, and let t_{m,n} denote the number of three-choice polygons with 2m horizontal steps and 2n vertical steps, then the anisotropic generating function for T can be written

T (x, y) = ∑_{m,n} t_{m,n}x^{m}y^{n} = ∑_{n} H_{n}(x)y^{n},

Figure 2. Showing the decomposition of an imperfect staircase polygon into a sequence of 2-4-2 non-intersecting walkers, each expressible as a Gessel-Viennot determinant

where H_{n}(x) = \frac{R_{n}(x)}{S_{n}(x)} is the (rational [15]) generating function for three-choice polygons with 2n vertical steps. In earlier, unpublished, numerical work, we found that, for imperfect staircase polygons, the denominators were:

S_{n}(x) = (1 − x)^{2n−1}(1 + x)^{(2n−7)_{+}} n even,

and

S_{n}(x) = (1 − x)^{2n−1}(1 + x)^{(2n−8)_{+}} n odd.

This was subsequently proved by Bousquet-Mélou [2]. Further, Bousquet-Mélou showed that the numerators satisfied:

R_{n}(−1) = \frac{−12(4m)!}{m!(m + 1)!(m + 2)!(m + 3)!} n = 2m + 4,

and

R_{n}(−1) = \frac{−96(4m + 1)!}{m!(m + 1)!(m + 2)!(m + 4)!} n = 2m + 5.

Unfortunately, we still do not have enough information to identify the numerators, though we observe that they are of degree 3n − 7 for n ≥ 4 and n even, and of degree 3n − 8 for n ≥ 5 and n odd.

It is also possible to express the generating function T (x) as a five-fold sum, with one constraint, [2] of 4 × 4 Gessel-Viennot determinants [5]. This is clear from figure 2, where the enumeration of the lattice paths between the dotted lines is just the classical problem of 4 non-intersecting walkers, and these must be joined to two non-intersecting walkers to the left, and to two non-intersecting walkers to the right. Then one must sum over different possible geometries. However the fact that the generating function is so expressible implies that it is differentiably finite [11].

In the following sections we report on our work leading to an ODE for the perimeter generating function of three-choice polygons. We started by generating the counts for three-choice polygons up to half-perimeter 260. Using numerical experimentation we then found what we believe to be the underlying ODE. This calculation required the use of the first 206 coefficients with the resulting ODE then correctly predicting the next 54 coefficients. While the possibility that this ODE is not the correct one is extraordinarily small, our result does of course not constitute a proof. Unfortunately we cannot usefully bound the size of the underlying ODE, otherwise we could use the knowledge of Dfiniteness to provide a proof of our results. That is to say, any bounds that follow from closure theorems [11] are too large to be useful.

2. Computer enumeration

The algorithm we use to count the number of imperfect polygons is a slightly modified version of the algorithm of Conway et al. [3]. Before proceeding to the full problem it is useful to briefly outline the transfer matrix algorithm for enumerating staircase polygons. Recall that a staircase polygon consists of two directed walks starting at the origin, moving only to the right and up, and terminating once the walks join at a vertex. If we look at a diagonal line x + y = k + 1/2 then for any integer k this line will intersect a polygon at 0 (miss the polygon) or 2 edges (intersect the polygon), see figure 1. We start with k = 0 such that the line intersects the first two edges of the staircase polygon. We then move the line upward (increase k by 1) and as we do this we add an edge to each walk. There are only four new configurations corresponding to the four possible steps. We need only keep track of the gap between the two walks, where the gap is the minimal number of iterations required in order to join the two walks. As we move the line the gap is either increased by a unit (the upper walk moves up and the lower walk moves right), decreased by a unit (the upper walk moves right and the lower walk moves up) or remains constant in two possible ways (both walks move up or right). Let C(i, k) be the number of configurations with a gap of i at step k. We then have the following very simple algorithm:

Set C(1, 0) = x (where x is a variable conjugate to the half-perimeter of the polygon). Run through all possible gaps i = 1, …, k + 1 and do the following updates: C(i +

1, k + 1) :=^{+} xC(i, k), C(i − 1, k + 1) :=^{+} xC(i, k) and C(i, k + 1) :=^{+} 2xC(i, k). Here a:=^{+} b is

short-hand for assign to a the value a + b.

Formally we can view the transformation from the set of states C(i, k) to C(j, k + 1) as a matrix multiplication (hence our use of the nomenclature transfer matrix algorithm) with k counting the number of iterations of the transfer matrix algorithm. However, as can be readily seen from the algorithm the transfer matrix is extremely sparse and there is no reason to list it explicitly (it is given implicitly by the updating rules).

The term C(0, k) is the number of staircase polygons of half-perimeter k +1. Note that the use of the variable x is somewhat superfluous in the case of staircase polygons since the generating function at iteration k is just x^{k+1}C(i, k), where C(i, k) is the number of configuration with gap i after k iterations. But it is included here for reasons of generality and in the case of imperfect staircase polygons the generating function will be a (nontrivial) polynomial in x. Naturally we need not actually keep all the entries C(i, k) since only the current and subsequent values are needed for the calculation so we can replace C(i, k) with C(i, k mod 2). We just have to initially set to zero all entries in the next step and keep a running total c(k) of the number of staircase polygons.

Imperfect staircase polygons start out as ordinary staircase polygons (see figure 2) Then at some vertex two additional directed walks (sharing the same starting point) are inserted between the two original walks (at the first dashed line marked with a ‘4’ in figure 2). The diagonal line will thus intersect these polygon configurations at four edges. Imperfect staircase polygons are created by connecting the first two walks and the last two walks (as illustrated at the last dashed line marked with a ‘4’ in figure 2). With four walks we need to retain three pieces of information, namely, the three gaps l, m, and n between consecutive walks. Each existing configuration can produce 16 new configurations as each walk is extended by a step either up or to the right. The resulting updating is easily worked out [3]. Let G(l, m, n) be the generating function (a polynomial in the variable x) for partially completed polygons at a given diagonal k. As we proceed to the next diagonal k + 1 we add x^{2}G(l, m, n) (the factor x^{2} arise because we extend all walks by a step) to G(l, m, n) (twice), G(l + 1, m, n), G(l, m + 1, n), G(l, m, n + 1), G(l − 1, m, n), G(l, m − 1, n), G(l, m, n − 1), G(l + 1, m − 1, n), G(l + 1, m, n − 1), G(l − 1, m + 1, n), G(l, m + 1, n − 1), G(l − 1, m, n + 1), G(l, m − 1, n + 1), G(l − 1, m + 1, n − 1) and G(l + 1, m − 1, n + 1). Any update resulting in G(l, 0, n) has to be rejected because it corresponds to a configuration in which we have joined the two middle walks in and this can never lead to an imperfect staircase polygon. Obviously once any two walks have been connected the remaining walks follow the usual staircase polygon updating rules.

The configurations with two walks already connected can also be encoded by the G functions. We simply let G(l, 0, 0) be the generating function for partial polygons with two walks already connected (note that if the boundary line intersects four edges m > 0). So in the updating of imperfect staircase polygons we can set G(0, m, n) (we connect the two lower-most walks) to G(n, 0, 0) (this case is illustrated at the last dashed line marked with a ‘4’ in figure 2). Likewise we can set G(l, m, 0) (we connect the two uppermost walks) to G(l, 0, 0). The condition for the formation of a valid polygon is l = n = 0 (note that we can’t demand m = 0 as well, since we could connect both the lower- and uppermost walks simultaneously).

The ‘creation’ of a configuration with three gaps, alternatively, one in which a diagonal line intersects four edges of an imperfect staircase polygon, which we refer to as a Gtype configuration, is also very simple (see the first dashed line marked with a ‘4’ in figure 2). We start with a staircase type configuration C(i) and from this we can create four G-type configurations by assigning the value x^{2}C(i) to G(j, 1, i−j), G(j −1, 1, i−j), G(j, 1, i−j −1) and G(j −1, 1, i−j −1), where 1 ≤ j ≤ i−1 (the factor x^{2} arises because we extend the outer walks by a step and insert two new walks each containing a single step).

The algorithm outlined above is already very efficient, but it can be further enhanced by the following simple observation. If we wish to calculate the number of polygons up to a given maximal half-perimeter length N, we need not consider all possible gaps since some configurations can only lead to polygons of a size exceeding N. First of all since gaps only increase or decrease by one at each iteration we need never consider configurations with gaps exceeding N/2. Furthermore, any G-type configuration with m > 0 must have half-perimeter at least k + m. Here we get the contribution k from the outermost walks (k is the number of forward steps or iterations taken) and the contribution m from the innermost walks (a gap m requires at least m steps). In order to produce an imperfect staircase polygon we have to add at least l + n additional steps (we have to join both the two upper-most and two lower-most walks), so if M = k + l + m + n > N we can discard this configuration. Not only can we thus discard some configurations when M > N but we can also further decrease the memory use since rather than storing N terms per configuration we only need to store N − M terms.

We calculated the number of imperfect staircase polygons up to perimeter 520. The integer coefficients become very large so the calculation was performed using modular arithmetic [10]. This involves performing the calculation modulo various prime numbers p_{i} and then reconstructing the full integer coefficients at the end. We used primes of the form p_{i} = 2^{30} − r_{i}, where r_{i} are small positive integers, less than 1000, chosen so that p_{i} is prime, and p_{i} ≠ p_{j} unless i = j. 18 primes were needed to represent the coefficients correctly. The calculation for each prime used about 250Mb of memory and about 18 minutes of CPU time on a 2.8 GHz Xeon processor. Naturally we could have carried the calculation much further but as we shall demonstrate in the next section this more than sufficed to identify an exact differential equation satisfied by T (x).

3. The Fuchsian differential equation

In recent papers Zenine et al. [17, 18, 19] obtained the linear differential equations whose solutions give the 3- and 4-particle contributions χ^{(3)} and χ^{(4)} to the Ising model susceptibility. In this paper we use their method to find a linear differential equation which has as a solution the generating function T (x) for three-choice polygons. We briefly outline the method here. Let us assume we have a function F (x) with a singularity at x = x_{c} = 1/μ. Starting from a (long) series expansion for the function F (x) we look for a linear differential equation of order m of the form

(1)  ∑_{k=0}^{m} P_{k}(x)\frac{d^{k}}{dx^{k}} F (x) = 0,

such that F (x) is a solution to this homogeneous linear differential equation, where the P_{k}(x) are polynomials. In order to make it as simple as possible we start by searching for a Fuchsian [9] equation. Such equations have only regular singular points. There are several reasons for searching for a Fuchsian equation, rather than a more general D-finite equation. Computationally the Fuchsian assumption simplifies the search for a solution. One may also argue, less precisely, that for “sensible” combinatorial models one would expect Fuchsian equations, as irregular singular points are characterized by explosive, super-exponential behaviour. Such behaviour is not normally characteristic of combinatorial problems. (The point at infinity may be an exception to this somewhat imprecise observation). One may also ask the question whether most of the problems in combinatorics with D-finite solutions have Fuchsian solutions? While we have not made an exhaustive study, we know of no counter-example to this suggestion.

From the general theory of Fuchsian [9] equations it follows that the degree of P_{k}(x) is at most n − m + k where n is the degree of P_{m}(x). To simplify matters (reduce the order of the unknown polynomials) it is often advantageous to explicitly assume that the origin and x = x_{c} are regular singular points and to set P_{k}(x) = Q_{k}(x)S(x)^{k}, where S(x) = xR(x) and R(x) is a polynomial of minimal degree having x_{c} as a root (in our case we have R(x) = 1 − 4x). S(x) could be generalised to include more regular singular points if some were known from other methods of analysis, but we have not found this to be advantageous. Thus when searching for a solution of Fuchsian type there are only two parameters: namely the order m of the differential equation and the degree q_{m} of the polynomial Q_{m}(x). Let ρ be the degree of S(x) (2 in our case), then for given m and q_{m} there are L = (m + 1)(q_{m} + 1) + ρm(m + 1)/2 − 1 unknown coefficients, where we have assumed without loss of generality that the leading order coefficient in P_{m}(x) = Q_{m}(x)S(x)^{m} is 1. We can then search systematically for solutions by varying m and q_{m}. In this way we first found a solution with m = 10 and q_{m} = 12, which required the determination of L = 206 unknown coefficients. We have 260 terms in the half-perimeter series and thus have more than 50 additional terms with which to check the correctness of our solution. Having found this conjectured solution we then turned the ODE into a recurrence relation and used this to generate more series terms in order to search for a lower order Fuchsian equation. The lowest order equation we found was eighth order and with q_{m} = 30, which requires the determination of L = 321 unknown coefficients. Thus from our original 260 term series we could not have found this 8^{th} order solution since we did not have enough terms to determine all the unknown coefficients in the ODE. This raises the question as to whether perhaps there is an ODE of lower order than 8 that generates the coefficients? The short answer to this is no. Further study [7] of our differential operator revealed that it can be factorised. In fact we found a factorization into three first-order linear operators, a second order and a third order. The generating function is a solution of the 8^{th} order operator, not of any of the smaller factors.

So the (half)-perimeter generating function T (x) for three-choice polygons is conjectured to be a solution of the linear differential equation of order 8

(2)  ∑_{k=0}^{8} P_{k}(x)\frac{d^{k}}{dx^{k}} F (x) = 0

with

P_{8}(x) = x^{3}(1 − 4x)^{4}(1 + 4x)(1 + 4x^{2})(1 + x + 7x^{2})Q_{8}(x),
P_{7}(x) = x^{2}(1 − 4x)^{3}Q_{7}(x), P_{6}(x) = 2x(1 − 4x)^{2}Q_{6}(x),

P_{5}(x) = 6(1 − 4x)Q_{5}(x), P_{4}(x) = 24Q_{4}(x), (3)

P_{3}(x) = 24Q_{3}(x), P_{2}(x) = 144x(1 − 2x)Q_{2}(x),
P_{1}(x) = 144(1 − 4x)Q_{1}(x), P_{0}(x) = 576Q_{0}(x),

where Q_{8}(x), Q_{7}(x), …, Q_{0}(x), are polynomials of degree 25, 31, 32, 33, 33, 32, 29, 29, and 29, respectively. The polynomials are listed in Appendix A (note that the polynomials do not factorise).

The singular points of the differential equation are given by the roots of P_{8}(x). One can easily check that all the singularities (including x = ∞) are regular singular points so equation (2) is indeed of the Fuchsian type. It is thus possible, using the method of Frobenius, to obtain from the indicial equation the critical exponents at the singular points. These are listed in Table 1.

For equations of the Fuchsian type the critical exponents satisfy a simple Fuchsian summation relation, which we now take the opportunity to confirm in our case. Let x_{1}, x_{2}, … x_{n}, x_{n+1} = ∞ be the regular singular points of a Fuchsian type equation of order m and α_{j,1}, … α_{j,m} (j = 1, … n + 1) the m exponents determined from the roots of the indicial equation for each regular singular point, x_{j}, then the following Fuchsian Table 1. Critical exponents for the regular singular points of the Fuchsian differential equation satisfied by T (x).

Singularity  Exponents
x = 0  −1, 0, 0, 0, 1, 2, 3, 4
x = 1/4  −1/2, −1/2, 0, 1/2, 1, 3/2, 2, 3
x = −1/4  0, 1, 2, 3, 4, 5, 6, 13/2
x = ± i/2  0, 1, 2, 3, 4, 5, 6, 13/2
1 + x + 7x^{2} = 0 0, 1, 2, 2, 3, 4, 5, 6
x = ∞  −2, −3/2, −1, −1, −1/2, 1/2, 3/2, 5/2
Q_{8}(x) = 0  0, 1, 2, 3, 4, 5, 6, 8

relation holds:

(4)  ∑_{j=1}^{n+1} ∑_{k=1}^{m} α_{j,k} = \frac{(n − 1)m(m − 1)}{2}.

In this case the number of regular singular points is m + 1 = 33, namely the 25 roots of Q_{8}(x), the two roots of 1 + x + 7x^{2}, x = ± i/2, x = ±1/4, x = 0 and x = ∞. It is easy to verify that the Fuchsian relation is satisfied with m = 8, n = 32, and all the exponents α_{j,k} summing to 868, which is a useful check on our results.

We shall now consider the local solutions of the differential equation around each singularity. Recall that in general it is known [4, 9] that if the indicial equation yields k critical exponents which differ by an integer, then the local solutions may contain logarithmic terms up to log^{k−1}. However, for the Fuchsian equation (2) only multiple roots of the indicial equation give rise to logarithmic terms in the local solution around a given singularity, so that a root of multiplicity k gives rise to logarithmic term up to log^{k−1}.

In particular this means that near any of the 25 roots of Q_{8}(x) the local solutions have no logarithmic terms and the solutions are thus analytic since all the exponents are positive integers. The roots of Q_{8}(x) are thus apparent singularities [4, 9] of the Fuchsian equation (2). There are methods for distinguishing real and apparent singularities (see, e.g, [4] §45) and in principle one should check that the roots of Q_{8}(x) satisfy the conditions for being apparent singularities. However, this theoretical method is quite cumbersome. An easier numerical way to see that the roots of Q_{8}(x) must be apparent singularities is as follows. We already found a 10th order Fuchsian equation for which the polynomial P_{10}(x) was of a form similar to P_{8}(x) as listed in equation (3), but with the degree of Q_{10}(x) being only 7. That is all the singularities as tabulated in Table 1 also appear in this higher order equation with the exception of the 25 roots of Q_{8}(x) (at most 7 of these could appear in the order 10 Fuchsian equation). In fact we can find a solution of order 14 of the same form as above but with Q_{14}(x) being just a constant. So at this order none of the roots of Q_{8}(x) appear. Clearly any real singularity of the system cannot be made to vanish and we conclude that the 25 roots of Q_{8}(x) must indeed be apparent singularities.

Assuming that only repeated roots give rise to logarithmic terms, and thus that a sequence of positive integers give rise to analytic terms, then near the physical critical point x = x_{c} = 1/4 we expect the singular behaviour

(5)  T (x) ∼ A(x)(1 − 4x)^{−1/2} + B(x)(1 − 4x)^{−1/2} log(1 − 4x),

where A(x) and B(x) are analytic in the neighbourhood of x_{c}. Note that the terms associated with the exponents 1/2 and 3/2 become part of the analytic correction to the (1 − 4x)^{−1/2} term. Near the singularity on the negative x-axis, x = x_{−} = −1/4 we expect the singular behaviour

(6)  T (x) ∼ C(x)(1 + 4x)^{13/2},

where again C(x) is analytic near x_{−}. We expect similar behaviour near the pair of singularities x = ±i/2, and finally at the roots of 1 + x + 7x^{2} we expect the behaviour T (x) ∼ D(x)(1 + x + 7x^{2})^{2} log(1 + x + 7x^{2}).

Next we turn our attention to the asymptotic behaviour of the coefficients of T (x). To standardise our analysis, we assume that the critical point is at 1. The growth constant of staircase and imperfect staircase polygons is 4, so we normalise the series by considering a new series with coefficients r_{n}, defined by r_{n} = t_{n+2}/4^{n}. Thus the generating function we study is R(y) = ∑_{n≥0} r_{n}y^{n} = 4 + 3y + 2.625y^{2} + · · ·. From equations (5) and (6) it follows that the asymptotic form of the coefficients is

(7)  [y^{n}]R(y) = r_{n} = \frac{1}{\sqrt{n}} ∑_{i≥0} (\frac{a_{i} log n + b_{i}}{n^{i}} + (−1)^{n} (\frac{c_{i}}{n^{7+i}} )) + O(λ^{−n}).

The last term includes the effect of other singularities, further from the origin than the dominant singularities. These will decay exponentially since λ > 1 in the scaled variable y = x/4.

Using the recurrence relations for t_{n} (derived from the ODE) it is easy and fast to generate many more terms r_{n}. We generated the first 100000 terms and saved them as floats with 500 digit accuracy (this calculation took less than 15 minutes). With such a long series it is possible to obtain accurate numerical estimates of the first 20 amplitudes a_{i}, b_{i}, c_{i} for i ≤ 19 with precision of more than 100 digits for the dominant amplitudes, shrinking to 10–20 digits for the the case when i = 18, or 19. In making these estimates we have ignored the exponentially decaying term, which is the last term in eq.(6). In this way we confirmed an earlier conjecture [3] that a_{0} = \frac{3\sqrt{3}}{π^{3/2}} , where we have taken into account the different normalisation, as discussed in the introduction. We also find that b_{0} = 3.173275384589898481765 … and c_{0} = \frac{−24}{π^{3/2}} , though we have not been able to identify b_{0}. However, we have successfully identified further sub-dominant amplitudes, and find a_{1} = \frac{−89}{8\sqrt{3}π^{3/2}} , a_{2} = \frac{1019}{384\sqrt{3}π^{3/2}} , and a_{3} = \frac{−10484935}{248832\sqrt{3}π^{3/2}} , and c_{1} = \frac{225}{π^{3/2}} , c_{2} = \frac{−16575}{16π^{3/2}}, and c_{3} = \frac{389295}{128π^{3/2}} . It seems possible that the amplitudes π^{3/2}\sqrt{3}a_{i} and π^{3/2}c_{i} are rational.

Estimates for the amplitudes were obtained by fitting r_{n} to the form given above using an increasing number of amplitudes. ‘Experimentally’ we find we need about the same total number of terms at x_{c} and −x_{c} = x_{−}.

So in the fits we used the terms with amplitudes a_{i}, and b_{i}, i = 0, …, K and c_{i}, i = 0, …, 2K. Going only to i = K with the c_{i} amplitudes results in much poorer convergence and going beyond 2K leads to no improvement. For a given K we thus have to estimate 4K + 3 unknown amplitudes. So we use the last 4K + 3 terms r_{n} with n ranging from 100000 to 100000 − 4K − 2 and solve the resulting 4K + 3 system of linear equations. We find that the amplitudes are fairly stable up to around 2K/3. We observed this by doing the calculation with K = 30 and K = 40 and then looking at the difference in the amplitude estimates. For a_{0} and b_{0} the difference is less than 10^{−131}, while for c_{0} the difference is less than 10^{−123}. Each time we increase the amplitude index by 1 we loose around 10^{6} in accuracy. With i = 20 the differences are respectively around 10^{−16} and 10^{−8}.

The excellent convergence is solid evidence (though naturally not a proof) that the assumptions leading to equation (7) are correct. Further evidence was obtained as follows. We can add extra terms to the asymptotic form and check what happens to the amplitudes of the new terms. If the amplitudes are very small it is highly likely that the terms are not truly present (if the calculation could be done exactly these amplitudes would be zero). One possibility is that our assumption about integer exponents leading only to analytic terms is incorrect. To test this we fitted to the form

\frac{1}{\sqrt{n}} ∑_{i≥0} (\frac{ã_{i} log n + b̃_{i}}{n^{i/2}} + (−1)^{n} (\frac{c̃_{i}}{n^{7+i}} )) + O(λ^{−n}),

(as above, in making these estimates we have ignored the exponentially decaying term, which is the last term in the above equation.) With K = 30 we found that the amplitudes ã_{1} and b̃_{1} of the terms log n/n and 1/n, respectively, were less than 10^{−60}, while the amplitudes ã_{3} and b̃_{3} were less than 10^{−50}. We think we can safely say that all the additional terms we just added are not present. We found similar results if we added terms like log^{2} n or additional log n terms at y = −1. That is, we found that those terms were not present. So this fitting procedure provides convincing evidence that the asymptotic form (7), and thus the assumption leading to this formula, is correct.

4. Conclusion

We have developed an improved algorithm for enumerating three choice polygons. The extended series, coupled with a search program that assumes the solution is a Fuchsian ODE, enabled us to discover the underlying ODE, which was of 10^{th} order. We did this without using more than 50 of the coefficients that we had generated. That is to say, 50 of the known coefficients were unused, and so their value provided a check on the solution found. This leads us to believe that we have found the correct ODE, as it reproduces the known, unused coefficients, though we have not provided a proof. Further refinement allowed us to find an 8^{th} order ODE.

A numerical technique we have developed specifically for such problems then allowed us to find accurate numerical estimates for the amplitudes of the first several terms in the asymptotic form for the coefficients.

We have also initiated an investigation of the area generating function. We expect this to involve q-series, and thus far our investigations only lead us to believe that the area generating function A(q) is of the form

A(q) = (G(q) + H(q)/\sqrt{1 − q/η})/[J_{0}(1, 1, q)^{2}],

where J_{0}(x, y, q) is a q-generalisation of the Bessel function, and occurs, for example, in the solution of the problem of staircase polygons enumerated by area [1]. Here q = η is the first zero of J_{0}(1, 1, q), and G and H are regular in the neighbourhood of q = η. The coefficients thus behave asymptotically as

a_{n} = [q^{n}]A(q) ∼ const.η^{−n}n^{3/2}.

In a subsequent publication [7] we propose to investigate the area generating function more fully, and to say more about the properties of the ODE we have found for the perimeter generating function. In particular, we discuss some simple solutions of the ODE, and ask what these can tell us about the full solution.

Acknowledgments

We would like to thank N Zenine and J-M Maillard for access to their Mathematica routines for identifying differential equations and useful advice about their use. We would also like to thank M. Bousquet-Mélou for communicating her unpublished results on this problem, which we mentioned in the introduction. We gratefully acknowledge financial support from the Australian Research Council.

Appendix A: Polynomials Q_{n}(x)

Q_{8}(x) = −180 + 4005x − 45340x^{2} + 352567x^{3} − 2100653x^{4} + 8247059x^{5} + 1869782x^{6}

−198745492x^{7} + 222232422x^{8} + 7981490552x^{9} − 58454247760x^{10} +223070561538x^{11} − 653903984242x^{12} + 1691567153918x^{13} −3628069390936x^{14} + 9508812403200x^{15} − 42130737708796x^{16} +151950842991736x^{17} − 347187650580720x^{18} + 558723092175488x^{19} −722483977609792x^{20} + 551434913787008x^{21} + 82817126361856x^{22} −426478334005248x^{23} + 279157576126464x^{24} + 2780644737024x^{25}

Q_{7}(x) = −3420 + 82530x − 926615x^{2} + 6866662x^{3} − 37878392x^{4} + 131975108x^{5}

+198512462x^{6} − 5322566116x^{7} + 16816064102x^{8} + 88956629348x^{9} −872972184658x^{10} + 3395585125316x^{11} − 8662194926872x^{12} +2179593948608x^{13} + 130585482759744x^{14} − 698610495175368x^{15} +2229946022661696x^{16} − 6216128747042864x^{17} +15724091332879132x^{18} − 38607908490402392x^{19} +128963713249678592x^{20} − 464640056155209952x^{21} +1296873363475699328x^{22} − 2966555758830491904x^{23} +5741739615453110784x^{24} − 7824348079140616704x^{25} +8096625038421797888x^{26} − 6327622359115208704x^{27} −663175049190105088x^{28} + 4390942020748738560x^{29} −3449431865352388608x^{30} − 33011814317948928x^{31}

Q_{6}(x) = −9180 + 310275x − 4493475x^{2} + 40204094x^{3} − 262917778x^{4}

+1302960911x^{5} − 3743237840x^{6} − 8573351756x^{7} + 140454430666x^{8} −409322626730x^{9} − 1570504457342x^{10} + 18303308342032x^{11} −89658228463172x^{12} + 259420736216632x^{13} + 26862202296376x^{14} −4190021721023184x^{15} + 21897720821926584x^{16} −75837533674259508x^{17} + 212508813586272428x^{18} −476010656497826944x^{19} + 1034090705056496672x^{20} −3196181326637410304x^{21} + 10833216991064882848x^{22} −30172750280212408832x^{23} + 70340668591569812736x^{24} −132300506186507025408x^{25} + 177280513560453634560x^{26} −184990945124657242112x^{27} + 135828858351882342400x^{28} +12754320650381836288x^{29} − 85576383794502107136x^{30} +61165017902554546176x^{31} + 565560894352785408x^{32}

Q_{5}(x) = −4500 + 244800x − 4876845x^{2} + 55164150x^{3} − 438701640x^{4}

+2758453094x^{5} − 13804842198x^{6} + 45370091528x^{7} − 3608230380x^{8} −892524490064x^{9} + 4421327158154x^{10} − 2297315126532x^{11} −103201897035096x^{12} + 748998082407080x^{13} − 2329708885595260x^{14} −457382726191024x^{15} + 35817660448173240x^{16} −188156345496838984x^{17} + 677783573996257364x^{18} −1904649390940935752x^{19} + 4199594693024922016x^{20} −8814226144821806432x^{21} + 23568486792872894272x^{22} −70089404940793421632x^{23} + 188311273940137111552x^{24} −435002993494719438848x^{25} + 791152555777632593920x^{26} −1045593345640931730432x^{27} + 1096015208846337957888x^{28} −774016903940080771072x^{29} − 37178029375778357248x^{30} +412071049964952354816x^{31} − 275345921075326746624x^{32} −2464051649845395456x^{33}

Q_{4}(x) = 31500 − 1114080x + 17560755x^{2} − 178469565x^{3} + 1412918104x^{4}

−9431590849x^{5} + 52336335969x^{6} − 220707961458x^{7} + 525965711332x^{8} +915935968370x^{9} − 13996439933349x^{10} + 35303141246088x^{11} +231992664240696x^{12} − 2180352456480752x^{13} + 6859298731027888x^{14} −1272161420555012x^{15} − 75338205421491734x^{16} +406836568590013948x^{17} − 1513874477368697252x^{18} +4439738234446975124x^{19} − 10514406278248398472x^{20} +22797584086521040520x^{21} − 52624647215757093584x^{22} +130673617185226821792x^{23} − 324680301683722155712x^{24} +724969297042825531136x^{25} − 1271869215082051692800x^{26} +1661614177465373698560x^{27} − 1744486537247742479360x^{28} +1209724637255295010816x^{29} − 32142663152460406784x^{30} −498040622799430975488x^{31} + 321317702703841148928x^{32} +2787318284392857600x^{33}

Q_{3}(x) = −156000 + 3778920x − 46727325x^{2} + 457371630x^{3} − 3919246431x^{4}

+27446185200x^{5} − 152613919692x^{6} + 659637747242x^{7} − 1723470963068x^{8} −1667066145852x^{9} + 27889854017778x^{10} + 15933308039400x^{11} −972460279627326x^{12} + 4552136023731292x^{13} − 7976188460233924x^{14} −4422880527966948x^{15} + 63325989574562728x^{16} −287206984863975352x^{17} + 1115308575007981980x^{18} −3508943115779966584x^{19} + 8987842561562515768x^{20} −19184807012355087408x^{21} + 37821550927408731776x^{22} −83609060792238083072x^{23} + 194683017390969665280x^{24} −366258975512082319872x^{25} + 498254429378056694784x^{26} −558421919820222289920x^{27} + 441211762632912959488x^{28} −80404063142199537664x^{29} − 110342796490113417216x^{30} +85904442856027127808x^{31} + 720965567415582720x^{32}

Q_{2}(x) = 102000 − 1245240x + 445275x^{2} + 77507430x^{3} − 505005638x^{4} + 674357270x^{5}
+7410398802x^{6} − 50751541108x^{7} + 109730141494x^{8} + 263567061768x^{9}
−2398666258514x^{10} + 4447124418524x^{11} + 33544348232760x^{12}
−341405641395740x^{13} + 1843130781900080x^{14} − 7441271357292384x^{15}
+23827305830694324x^{16} − 59142500096057112x^{17}
+113845825936073424x^{18} − 169659492965796928x^{19}
+190085091157739584x^{20} − 160391840217609984x^{21}
+95477320250924800x^{22} − 21461546279272960x^{23}
−73590898428536832x^{24} + 43442402559821824x^{25}
+129164030193680384x^{26} − 136460131311329280x^{27}
+54532752690511872x^{28} + 389290263183360x^{29}

Q_{0}(x) = Q_{1}(x) = −Q_{2}(x)

References

[1] Bousquet-Mélou M 1996 A method for the enumeration of various classes of column-convex polygons Disc. Math. 154 1–25

[2] Bousquet-Mélou M Private communication.

[3] Conway A R, Guttmann A J and Delest M 1997 The number of three-choice polygons Mathl. Comput. Modelling 26 51–58

[4] Forsyth A R 1902 Part III. Ordinary linear equations vol. IV of Theory of differential equations. (Cambridge: Cambridge University Press)

[5] Gessel I and Viennot X G 1989 Determinants, paths and plane partitions preprint

[6] Guttmann A J and Conway A R 2001 Square lattice self-avoiding walks and polygons Ann. Comb. 5 319–345

[7] Guttmann A J and Jensen I Properties of Fuchsian differential equations for polygon enumeration problems In preparation

[8] Guttmann A J, Prellberg T and Owczarek A L 1993 On the symmetry classes of planar self-avoiding walks J. Phys. A: Math. Gen. 26 6615–6623

[9] Ince E L 1927 Ordinary differential equations (London: Longmans, Green and Co. Ltd.)

[10] Knuth D E 1969 Seminumerical Algorithms. The Art of Computer Programming, Vol 2. (Reading, Mass: Addison Wesley)

[11] Lipshitz L 1989 D-finite power series, J. Algebra 122 353–373

[12] Manna S S 1984 Critical behaviour of anisotropic spiral self-avoiding walks J. Phys. A: Math. Gen. 17 L899–L903

[13] Rechnitzer A 2003 Haruspicy and anisotropic generating functions Advances in Applied Mathematics 30 228–257

[14] Stanley R P 1980 Differentiably finite power series European J. Combin. 1 175–188

[15] Stanley R P 1999 Enumerative Combinatorics vol. 2 (Cambridge: Cambridge University Press)

[16] Whittington S G 1985 Anisotropic spiral self-avoiding walks J. Phys. A: Math. Gen. 18 L67–L69 [17] Zenine N, Boukraa S, Hassani S and Maillard J M 2004 The Fuchsian differential equation of the

square lattice Ising model χ^{(3)} susceptibility J. Phys. A: Math. Gen. 37 9651–9668

[18] Zenine N, Boukraa S, Hassani S and Maillard J M 2005 Square lattice Ising model susceptibility: series expansion method and differential equation for χ^{(3)} J. Phys. A: Math. Gen. 38 1875–1899

[19] Zenine N, Boukraa S, Hassani S and Maillard J M 2005 Ising model susceptibility: the Fuchsian differential equation for χ^{(4)} and its factorization properties J. Phys. A: Math. Gen. 38 4149–4173