# On the Number of Conjugate Classes of Derangements

Authors: Wen-Wei Li; Zhong-Lin Cheng; Jia-Bao Liu
Journal: Journal of Mathematics
 (2021)

---

## Abstract

The number of conjugate classes of derangements of ordern is the same as the number hn of the restricted partitions with every portion greater than 1. It is also equal to the number of isotopy classes of 2×n Latin rectangles. Sometimes the exact value is necessary, while sometimes we need the approximation value. In this paper, a recursion formula of hn will be obtained and also will some elementary approximation formulae with high accuracy for hn be presented. Although we may obtain the value of hn in some computer algebra system, it is still meaningful to find an efficient way to calculate the approximate value, especially in engineering, since most people are familiar with neither programming nor CAS software. This paper is mainly for the readers who need a simple and practical formula to obtain the approximate value (without writing a program) with more accuracy, such as to compute the value in a pocket science calculator without programming function. Some methods used here can also be applied to find the fitting functions for some types of data obtained in experiments.


---

## Body

## 1. Introduction

Belown is a positive integer greater than 1.On some occasions, it is necessary to know the number of conjugate classes of derangements.When generating the representatives of all the isotopy classes of Latin rectangles of ordern by some method, we need to know the number of the isotopy classes of 2×n Latin rectangles for verification. In some cases, we need the approximate value in a simple and efficient method. (When writing a C program to generate the representatives of all the isotopy classes of Latin rectangles of order n, we need to prepare some space in the memory module (RAM) to store the cycle structures of derangements so as to make the program more efficient; otherwise, we have to allocate memory dynamically, which will cost more time in memory addressing when writing and reading data frequently in the particular position in the memory module. So, we need to know the number of the isotopy classes of 2×n Latin rectangles for verification.)LetSn be the symmetry group of the set X = {1, 2, ⋯, n}, i.e., the set (together with the operation of combination) of the bijections from X to itself. An element σ in the symmetry group Sn is also called a permutation (of order n). If σ∈Sn, σi≠i (∀i∈X), σ will be called a derangement of order n. If a permutation σ transforms any element in X to a distinct element, then the sequence σ1,σ2,…,σn will also be called a derangement. The number of derangements of order n is denoted by Dn (or !n in some literatures). It is mentioned in nearly every combinatorics textbook that(1)Dn=n−1Dn−1+Dn−2=n!∑i=0n−1ii!≐n!e+12,n⩾1,where x is the floor function, which stands for the maximum integer that is less than or equal to the real x.Forx,y∈Sn, if ∃z∈Sn, s.t. x=zyz−1, then x and y will be called conjugate and y is called the conjugation of x. Of course the conjugacy relation is an equivalence relation. So, the set of derangements of order n can be divided into some conjugate classes. This paper is mainly concerned on the number of conjugate classes of derangements of order n. The main method is similar to that described in reference [1].A matrix of sizek×n (1≤k≤n−1) with every row being a reordering of a fixed set of n elements and every column being a part of a reordering of the same set of n elements is called a Latin rectangle. Usually, the set of the n elements is assumed to be {1, 2, 3, ⋯, n}. (In some literatures, the members in a Latin rectangle is assumed in the set {0, 1, 2, ⋯, n−1}.)A2×n Latin rectangle with the first row in increasing order could be considered as a derangement. An isotopy class of 2×n Latin rectangles will correspond to a unique conjugate class of derangements. So, the number of isotopy classes of 2×n Latin rectangles is the same as the number of conjugate classes of derangements of order n.All the members in a conjugate class of derangements inSn share the same cycle structure. Here, we define the cycle structure of a derangement as the sequence in nondecreasing order of the lengths (with duplicate entries) of all the cycles in the cycle decomposition of the derangement. A cycle structure of a derangement of order n could be considered as an integer solution of the equation as follows:(2)s1+s2+⋯+sq=n,2≤s1≤s2≤⋯≤sq,where s1, s2, ⋯, sq are unknowns.For a fixedq, designate the number of integer solutions of equation (2) as Hqn where q is less than n/2+1 (otherwise, Hqn is defined by 0), and denote hn the number of all the integer solutions of equation (1) for all possible q, i.e.,(3)hn=∑n/2q=1Hqn.So, the number of conjugate classes of derangements of ordern is hn. Since hn is the number of a type of restricted partitions, it is tightly connected with the partition number pn.Following the notation of [2], denote by Pqn the number of integer solutions of equation(4)s1+s2+⋯+sq=n,1≤s1≤s2≤⋯≤sq,for a fixed q, where 1≤q≤n, and by pn the number of all the (unrestricted) partitions of n. It is clear that (in a lot of articles, pn,q is used instead of Pqn, but in some other literatures, pn,q stands for some other number.)(5)pn=∑q=1nPqn.There is a brief introduction of the important results on the partition number (or partition function)pn and Pqn in reference [2], such as the recursion formula of pn and Pqn. More information about the partition number pn may be found in reference [3]. There is a list of some important papers and book chapters on the partition number in [4] (including the “LINKS” and “REFERENCES”) and [5]. Reference [1] presented some estimation formulae with high accuracy for pn, which are revised from Hardy–Ramanujan’s asymptotic formula.There are also a lot of literatures on the number of some types of restricted partitions ofn (such as [6–13]) or on the congruence properties of (restricted) partition function (such as [14–21]).In [22], we can find many cases of restricted partitions (some of them are introduced in [23, 24] or [25]). One class is concerned on the restriction of the sizes of portions, such as portions restricted to Fibonacci numbers, powers (of 2 or 3), unit, primes, nonprimes, composites, or noncomposites; another class is related to the restriction of the number of portions, such as the cases that the number of parts will not exceed k; the third class is about the restrictions for both, for example, the cases that the number of parts is restricted while the parts restricted to powers or primes. However, the author has not found too much information on the number hn, especially on the approximate calculation, although we can find a lot of information on other restricted partition numbers.The basic method of function fitting (curve fitting) could be found in any textbook of numerical analysis, such as [26, 27]. Some tricks used here may be found in some books of data analysis such as [28]. They may be helpful for understanding the methods described in Section 3.Section2 will deduce the recursion formula for hn and will show the relation of hn and pn. Subsection 3.1 will deduce the asymptotic formula of hn from Hardy–Ramanujan’s asymptotic formula of pn (mentioned in [1]). This new asymptotic formula Ign coincides with Ingham’s result (refer [29, 30]). By bringing in two parameters C1n and C2n in the new asymptotic formula Ign, we have reached several estimation formulae for hn with high accuracy in Subsection 3.3, using the same ideas described in [1]. By fitting hn−Ign, we have another two estimation formulae for hn in Subsection 3.4. When n<100, we have a more accurate estimation formula for hn in Subsection 3.5. The relative errors of these estimation formulae will be presented to shown the accuracy.


## 2. Some Formulae forhn

In this section, we will derive a recursive formula from the method mentioned in reference ([2], p. 53˜55).By definition,hk = 0 when k<2, but here we assume that hk = 0 when k<0, h0=1, and h1=0, for convenience.It is mentioned in reference ([2], p. 52) that in 1942, Auluck gave an estimation of Pqn by Pqn≈1/q!n−1q−1 when n is large enough.By the same method shown in reference ([2], p. 53, 57), we can obtain the generation function of hn:(6)Gx=∑n=0∞hnxn=∏i=2∞1−xi−1=11−x211−x311−x4⋯11−xi⋯⋯,and a formula(7)hn=12πi∮CGxxn+1dx,where h0=1, h1=0, and C is a contour around the original point.It is difficult to get a simple formula to count the solutions of equation (2) in general. However, for a fixed integer q, the number Hqn of solutions is 0 (when q>n/2) or(8)∑s1=2n/q∑s2=s1n−s1/q−1⋯∑sq−1=sq−2n−s1−s2⋯−sq−2/21=∑s1=2n/q∑s2=s1n−s1/q−1⋯∑sq−2=sq−3n−s1−s2⋯−sq−3/3n−s1−s2⋯−sq−22−sq−2+1=Pqn−qwhenq≤n2.It is not difficult to find out that(9)hn=∑q=1n/2Hqn=∑q=1n/2Pqn−q.and there is a recursion for Pqn in reference ([2], p. 51) as follows:(10)Pqn=∑j=1tPjn−q,where t=minq,n−q, so there is no difficulty to obtain the values of Pqn and hn when n is small.For the value ofpn, there is a recursion as follows:(11)pn=pn−1+pn−2−pn−5−pn−7+⋯+−1k−1pn−3k2±k2+⋯⋯,=∑k=1k1−1k−1pn−3k2+k2+∑k=1k2−1k−1pn−3k2−k2,where(12)k1=24n+1−16,k2=24n+1+16,and assume that p0=1 (refer [2], p. 55). Here, we assume that px = 0 when x<0.We can obtain the same recursion forhn as follows:(13)hn=hn−1+hn−2−hn−5−hn−7+⋯+−1k−1hn−3k2±k2+⋯⋯=∑k=1k1−1k−1hn−3k2+k2+∑k=1k2−1k−1hn−3k2−k2,where k1 and k2 are determined by equation (12), and assume that h0=1 and hk = 0 when k<0.The proof of equation (13) is easy to understand.By equation (6), we have(14)∑n=0∞hnxn∏i=2∞1−xi=1.SinceFx=∑n=0∞pnxn=∏i=1∞1−xi−1 (where p0=1), ∑n=0∞pnxn∏i=1∞1−xi=1 or(15)∑n=0∞pnxn1−x∏i=2∞1−xi=1.Comparing equations (14) and (15), we have(16)∑n=0∞hnxn=∑n=0∞pnxn1−x=∑n=0∞pn−pn−1xn,by assumption hk=pk=0 when k<0. Hence,(17)hn=pn−pn−1,n=0,1,2,⋯.By equation (17),(18)hn=hn−1+hn−2−hn−5−hn−7+⋯+−1k−1hn−3k2±k2+⋯⋯=∑k=1k1−1k−1hn−3k2+k2+∑k=1k2−1k−1hn−3k2−k2.We can easily obtain the solutions of equation (2) by hand when n < 13. By equation (13), we can obtain the number hn of solutions of equation (2) by writing a small program in some Computer Algebra System (CAS) software such as “maple,” “maxima,” and “axiom” or some other software likewise (be aware of that 0 is not a valid index value in some software such like maple).The value ofhn when n < 250 is shown in Table 1 (on page 5) and Table 2 (on page 5). Some values of Hqn are shown on Table 3 (on page 5).Table 1
The value ofhn when 1≤n≤100.
nh(n)nh(n)nh(n)nh(n)nh(n)1021165417245611550388122078512122210428591621786518225019283123253431008763205343832832214422432044119146423613184320519152253834513959652709288536236976426478461642466310962864095605742757447191966735616987462471187287084822519684080468852204369829847492625269466610895887816101230103950307017053362390663824811143112385135717716092379174781861221321507524164672695578928421448132433179453483427379290693947637014343421675456224749038119410659543154135257355651217510287649511981699165536309456755477611708279613462885176637366057873317713307729715116626188838437858101066781512301981696720619105395170591166007917164869919031739201374061536013464780194782610021339417Table 2
The value ofhn when 101≤n≤250.
nh(n)nh(n)nh(n)nh(n)10123911834116124763797131593224104146260819459010226784253117138801828132656291385147287161937910329983571118154364067133725798623148316074751910433552415119171594522134802411183149347793570310537524344120190680895135886795381150382588011310641950627121211798491136979745604160977543091110746873053122235172861137108206333617024329692015108523534551232610173291381194696815180591106378161095844339612428960225913913186080641901404538044681106521750612532118685214014549282402003269265972631117273945712635609534014116048110732107465212729801128109895312739464160314217696041122201674422848222113903744721284372143051431950689437230369330486166511410067403712948419327014421496716882408019313019148115112093786130536043530145236820356425017156634544056Table 3
The number of solutions of equation (1) for different q.
nh(n)H1 (n)H2 (n)H3 (n)H4 (n)H5 (n)H6 (n)H7 (n)421152116412174121871321981331101214421111414531122115752113241586311434161095211541161211731Obviously,hn < pn holds by definition (when n>1). As pn grows much more slowly than exponential functions, i.e., for any r>1, pn<rn will hold when n is large enough, which means we cannot estimate pn and hn by an exponential function. As pn grows faster than any power of n, it means we cannot estimate pn by a polynomial function (refer [2], p. 53). So, hn cannot be estimated by a polynomial function, either.


## 3. The Estimation ofhn

The recursion formula equation (13) for hn is not convenient in practical for a lot of people who do not want to write programs. Sometimes we need the approximation value, such as the cases mentioned in [1], so an estimation formula is necessary.The graph of datan,lnhn (n = 60+20k, k = 1, 2, ⋯, 397) is shown in Figure 1 on page 5. Here, the data points are displayed by small hollow circles, and the circles are very crowded that we may believe that the circles themselves be a very thick curve if we notice only the right-hand part. In this figure, the data points in the lower left part are sparse (compared with the points in the upper right part), and we may find some hollow circles easily. If there is a curve passes through these hollow circles, we will notice it (as shown in Figure 2 on page 14). However, later in Figure 3, the circles distribute uniformly on a curve, it will be difficult to distinguish the circles, and the curve passes through the centers of them.Figure 1
The graph of the datan,lnhn.Figure 2
The graph of the datan,lnhn and the fitting curve lnIg1n.Figure 3
The graph of the datan,πexpπ2/3n/122hn and the fitting curve.The author has not found a practical estimation formula with high accuracy of the numberhn before.Actually, it is very difficult to find directly a simple function to fit the data on Figure1 with high accuracy. The main reason is that the fitting functions obtained by the methods used frequently could not reach satisfying accuracy.We have several accurate estimation formula ofpn (refer [1]), such as(19)Rh2′n=exp2/3πn43n+a2n+c2+b2+12,n≥80,and(20)Rh0′n=exp2/3πn43n+C2′n+12,1≤n≤100,where a2=0.4432884566, b2=0.1325096085, c2=0.274078, and(21)C2′n=0.4527092482×n+4.35278−0.05498719946,n=3,5,7,…,99;0.4412187317×n−2.01699+0.2102618735,n=4,6,8⋯,100.By (17), we can obtain hn by(22)h1n=Rh0′n−Rh0′n−1,2≤n≤80,Rh2′n−Rh2′n−1,n>80.and the error of this formula will not exceed twice of the error of Rh2′n or Rh0′n. Of course, this formula is not as simple as we want, but the accuracy is very good.


### 3.1. Asymptotic Formula

Ashn=pn−pn−1, by Hardy–Ramanujan’s asymptotic formula,(23)pn∼14n3exp23πn,(refer [1, 3, 31–35]), we assume that, when n≫1, hn∼exp2/3πn/43n−exp2/3πn−1/43n−1. So,(24)hn∼exp2/3πn−143expπ2/3n−n−1n−1n−1=exp2/3πn−143expπ2/3/n+n−1n−1n−1∼exp2/3πn43expπ2/3/2nn−1n−1∼exp2/3πn431+π/6nn−1n−1∼exp2/3πn43π6n3=πexp2/3πn122n3.So,(25)hn∼π122n3exp2/3πn.In coincidence, half a year after the main results were obtained in this paper, the author found an asymptotic formula as follows:(26)Pa,bn∼Γb/aπb/a−12−3/2−b/2a3−b/2aa−1/2+b/2an−a+b/2aeπ2n/3a,in [30]. When a = 1 and b = 2, we will have(27)P1,2n∼π122n3expπ2/3n,which coincides with the asymptotic formula obtained here.Formula (25) will be called the Ingham–Meinardus asymptotic formula in this paper since Daniel mentioned in [30] that this general asymptotic formula (26) was first given by A. E. Ingham in [29] and the proof was refined by G. Meinardus later in another two papers written in German.Later in this paper,π/122n3exp2/3πn will be denoted by Ign for short.It is not satisfying to estimatehn by Ign when n is small. The relative error of Ign to hn is greater than 6% as shown in Table 4 (on page 6). The round approximation(28)Ig′n=Ign+12,will not change the accuracy distinctly, as shown in Table 5 (on page 6). So, it is necessary to modify the asymptotic formula for better accuracy.Table 4
The relative error ofIgn to hn when n≤1000.
nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)11650.304032.1022013.105208.392146.241756.825028.6024012.505408.233202.891846.696025.9026012.005608.08495.591952.757023.9028011.505807.935156.432044.948022.3030011.106007.79668.622148.489020.9032010.806407.547121.382243.4710019.8034010.406807.31865.432346.0011018.8036010.107207.10988.382441.0912018.003809.867606.911062.582543.6813017.204009.608006.731179.472639.9314016.604209.368406.561253.292741.2715016.004409.148806.411370.982838.5016015.404608.939206.271453.122939.7018014.504808.749606.131560.353037.0020013.705008.5610006.01Table 5
The relative error ofIgn+1/2 to hn when n≤30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—7125%1370.83%1952.38%2543.60%2100%871.43%1452.94%2045.26%2639.96%3200%987.50%1560.98%2148.48%2741.29%4100%1066.67%1650.91%2243.33%2838.56%5150%1178.57%1757.58%2345.85%2939.67%675%1252.38%1846.59%2440.94%3037.05%


### 3.2. Method A: Modifying the Exponent

In this subsection, we consider fittinghn by Iga=π/122n3exp2/3πn+C1n or fitting.n,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) by a function(29)f1x≐a1x+c1e1+b1.LetC1n = 3/2π2ln122n3hn/π2−n. The reason that we fit C1x by the function in the form displayed in (29) is the same as that described in Section 3 of [1] (although the data differ distinctly). Many other types of functions have been tested, but they cannot fit these data very well.However, here it is not valid to obtain the constants inf1n by the iteration method described in reference [1].The graph of the datan,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) is shown in Figure 4 on page 7.Figure 4
The graph of the datan,3/2π2ln122n3hn/π2−n.First, we try to fitn,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) by a function in the form as follows:(30)f1x≐a1x+c1+b1.That means we have assumed thate1=1/2, temporarily. We will explain the reason in Subsection 3.2.3.The average error off1x is(31)E1=1K1∑nC1n−f1n2=1K1∑nC1n−a1n+c1−b12E1=1K1∑k=1K1C160+20k−a160+20k+c1−b12,where K1 = 397 and n ranges from 80 to 8000, by step 20. Here, only a1, b1, and c1 are unknown, so we can consider E1 as a function of the variables a1,b1,c1.We want to find a triplea1,b1,c1 such that E1 reaches its minimum or to make E1 as small as possible.Since a lot of functions have several local minimum points, it is necessary to find out whetherE1 = E1a1,b1,c1 has more than one local minimum before we start to calculate the minimum point by the numeral method. However, E1 = E1a1,b1,c1 is too complicate, and it is very difficult to know all the critical points in the range we are considering.


#### 3.2.1. Preparation Work

For a given paira1,c1, by the property of the arithmetic mean (for some given data xi (i = 1, 2, ⋯, k; xi∈ℝ), the function st=1/k∑i=1kxi−t2 reaches its minimum at t = 1/k∑i=1kxi), it is clear that E1 reaches its minimum when(32)b1=1K1∑nC1n−a1n+c1,=C1¯−1K1∑n′a1n′+c1,where C1¯ = 1/K1∑n′C1n′.LetE1=1/K1∑nC1n−a1/n+c1−C1¯+1/K1∑n′a1/n′+c12 be the average error of the fitting function f1x≐a1/x+c1+C1¯−1/K1∑n′a1/n′+c1. (Here, a1 and b1 are undetermined coefficients.)LetG1=E12=1/K1∑nC1n−a1/n+c1−C1¯+1/K1∑n′a1/n′+c12. Here, only a1 and b1 are unknowns. G1 could be considered as a function of a1,c1. In order to find the minimum point of G1, we can draw the figure of the function G1=G1a1,c1 (in a cube coordinate system with axis a1,c1, and G1) as shown in Figures 5–7. Figures 8 and 9 are the projection of the graph of a1,c1,G1 (when −100≤a1≤ 100 and −50≤c1≤ 100, which is a part of a surface) on the a1−G1 plane (spanned by the axis a1 and G1) and c1−G plane, respectively.Figure 5
The graph of the dataa1,c1,G1 when −100≤a1≤ 100 and −50≤c1≤ 100.Figure 6
The graph of the dataa1,c1,G1 when −100≤a1≤ 100 and −50≤c1≤ 100 version 2.Figure 7
The graph of the dataa1,c1,G1 when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.Figure 8
The projection of the graph of the dataa1,c1,G1 on the a1−G1 plane when −100≤a1≤ 100 and −50≤c1≤ 100.Figure 9
The projection of the graph of the dataa1,c1,G1 on the c1−G1 plane when −100≤a1≤ 100 and −50≤c1≤ 100.From these figures, we can find out that the influence ofc1 to G1 is much less than that of a1. In Figure 10, we find that when G1 reaches its minimum, a1 is between 0.50 and 0.53, but there is not a definite range for c1.Figure 10
The projection of the graph of the dataa1,c1,G1 on the c1−G1 plane when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.It is possible that for different range ofc1, the range of a1 when G1 reaches its minimum will be different. However, considering that a1/n+c1+b1 is a real, c1 should be greater than −1 in theory. (For the fitting data used here, c1 should be greater than −80.) From Figure 8, we can see that G1 touches its bottom when −15≤a1≤ 15. Although we cannot see clearly the exact value of of a1 in the minimum points, we can draw another figure of a1,c1,G1 when −15≤a1≤ 15 and −1≤c1≤ 100 to observe more details (the figure is not presented here), and then we will find that the more detailed range of a1 for the minimum points is −3,3 in the new figure (not presented), and next we will draw the figures of a1,c1,G1 when −3≤a1≤ 3, 0≤a1≤ 1, 0.2≤a1≤ 0.8, or 0.45≤a1≤ 0.6, respectively, while −1≤c1≤ 100, we will find the range of a1 of the minimum points of G1 is about [0.50, 0.53]. The projections of Figure 7 of a1,c1,G1 when 0.45≤a1≤ 0.6 and −1≤c1≤ 100 are shown in Figures 11 and 10.Figure 11
The projection of the graph of the dataa1,c1,G1 on the a1−G1 plane when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.In Figure9, for the curves on the bottom, G1 decreases with c1 at first then increases with c1, but it is difficult to find the critical point in the figure since different curves have different critical points.Although we cannot find a satisfying value ofa1 or c1 from the figures to construct a fitting function f1, these pictures show that the figure of G1=a1,c1,G1 has only one bottom in the domain we are considering, unlike the figure of another function shown in Figure 12, so the existence of the minimum point is almost confirmed; therefore, we are confident to find the value of a1 or c1 in the minimum point by the numerical method. This guarantees the validity of the numerical calculation by loop in the next step.Figure 12
Example of a surface with more than one bottom.


#### 3.2.2. Findc1

On the other hand, by the least square method, to fit the dataxk,yk (k = 1, 2, ⋯, K1) by a linear function y=a×x+b, the result is that(33)a=xn⋅yn¯−xn¯⋅yn¯xn2¯−xn¯2,(34)b=xn¯⋅xn⋅yn¯−xn2¯⋅yn¯xn¯2−xn2¯=yn¯−xn¯⋅a,where xn¯ = 1/K1∑i=1K1xi is the average value of xi, yn¯ = 1/K1∑i=1K1yi, xi2¯ = 1/K1∑i=1K1xi2, and xn⋅yn¯ = 1/K1∑i=1K1xi⋅yi. By definition, xn¯2 = 1/K1∑i=1K1xi2 is the square of the average value of xn. So, by the least square method, a and b are uniquely determined by the given data xk,yk (k = 1, 2, ⋯, K1).For every given value ofc1 (greater than −80), we can fit n,C1n (n = 60+20k, k = 1, 2, ⋯, 397) by a function f1x≐a1/x+c1+b1 by the least square method if we consider 1/60+20k+c1 and C160+20k as xk and yk, respectively. Then,(35)xn¯=1K1∑k=1K1xk=1K1∑k=1K1160+20k+c1,yn¯=1K1∑k=1K1yk=1K1∑k=1K1C160+20k,xn2¯=1K1∑k=1K1xk2=1K1∑k=1K1160+20k+c1,xn⋅yn¯=1K1∑k=1K1xk⋅yk=1K1∑k=1K1C160+20k60+20k+c1,a1=xn⋅yn¯−xn¯⋅yn¯xn2¯−xn¯2,b1=yn¯−xn¯⋅a1.So,a1 and b1 could both be considered as functions of c1, denoted by a1=a1c1 and b1=b1c1, since they are uniquely determined by c1 with the given data.Then,G2=E12=1/K1∑k=1K1C160+20k−a1c1/60+20k+c1−b1c12 is a function of c1.It will cost some time to plot the figure of the functionG2 = G2c1 in a CAS software.If we plot the figure of the functionG2 = G2c1 on the coordinates (as shown in Figures 13–15), we will find that G2 reaches its minimum when c1≈−3.2594807. In Figure 15, we find that the curve of G2 = G2c1 is not so smooth. The reason is that we hold up 18 significant digits in the process. If we compute more significant digits in the process, the curve on Figure 15 will be more smooth, at the cost of much more time. By writing a small program (since the default function to find the minimum provided by the software Maple 18 is unable to deal with such a complicated function G2 = G2c1 involving so much data), we can obtain a more accurate value of the critical point as follows:(36)c1=−3.259480684.Figure 13
The graph of the functionG2 = G2c1 when −50≤c1≤ 100.Figure 14
The graph of the functionG2 = G2c1 when −4≤c1≤ −2.Figure 15
The graph of the functionG2 = G2c1 when −3.259483≤c1≤−3.259478.When the value ofc1 is obtained, we can find the value of a1 and b1 by the least square method without difficulty, i.e.,(37)a1=0.5097429624,b1=−1.453552800.However, herec1 is less than −1, so the estimation formula for hn constructed from these coefficients is invalid when n<4.


#### 3.2.3. Confirme1

In [1], we fit C1n = 3/2⋅ln4n3pn2/π2−n by a function f1x≐a1/n+c1e1+b1 when estimating pn and found that e1≈0.50 by iteration. Here, the iteration method does not work well, so we fit C1n = 3/2π2ln122n3hn/π2−n by a function f1x≐a1/n+c1+b1 directly, which means that we have assumed that e1=1/2. Here, we may doubt that whether e1=0.5 is the best option for us?Here, we use the same idea described in Subsection3.2.1.For every paire1,c1, we can obtain corresponding a1 and b1 by the least square method, just like (33) and (34), except that here xk=1/60+20k+c1e1.So, the square of the average error(38)G3=E12=1K1∑nC1n−a1e1,c1n+c1e1−b1e1,c12,could be considered as a function of e1 and c1 as both a1=a1e1,c1 and b1=b1e1,c1 could be expressed by certain elementary functions of e1 and c1.If we draw the figure of the functionG3 = G3e1,c1, we will find that the surface has only one bottom when 0.1 ≤e1≤0.9 and −50≤c1≤100, as shown in Figure 16. However, the process to draw the figure is time-consuming. It costs more than 5 hours on a notebook (ThinkPad E40 Edge, with 6 GB RAM and AMD P360 Dual-Core Processor 2.30 GHz) by Maple 18 in Ubuntu 14.04.1 system.Figure 16
The graph of the functionG3 = G3e1,c1 when 0.1≤e1≤0.9 and −50≤c1≤100.After that, by writing another program, we can obtain the approximate value ofe1,c1 where G3 touches the bottom, i.e., e1≈0.494 and c1≈−4.85, when 18 significant digits are involved in the process, which still costs tens of minutes. Considering that we have used only a small part of data, we cannot afford the time for computing more significant digits in process, and the computing is so complicated; hence, error accumulation effect is considerable, so we choose e1=0.50 while it differs very little with 0.494. Another reason is that we prefer simple exponent, as the time spent on computing a square root is much less than that to compute a power with exponent 0.494 in general. Here, the value of c1≈−4.85 is obviously different from the value obtained at the end of Subsection 3.2.2 because of the little difference on e1. Therefore, it will be fine to use the result in Subsection 3.2.2).In this figure, the points are shown as small circles which are very close to each other. These crowded circles seem like a thick curve. A fitting curve passes though the center of these circles. The fitting curve might not be found in reduce printing. That means the curve fits the points (displayed as circles) very well.


#### 3.2.4. The Result

By the valuea1 = 0.5097429624, b1 = −1.453552800, and c1 = −3.259480684 obtained in Subsection 3.2.2, we will have a fitting function f1x≐a1/x+c1+b1.The graph off1x (when −80≤x≤ 8000) and the graph of the data n,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) are shown in Figure 17. It shows that f1n fits 3/2π2ln122n3hn/π2−n very well.Figure 17
The graph of the datan,3/2π2ln122n3hn/π2−n and the fitting curve.Then, we could fithn by(39)Igan=πexp2/3πn+a1/x+c1+b1122n3,n≥4.The graph of the functionlnIgax and the graph of the data n,lnhn (n = 60+20k, k = 1, 2, ⋯, 397) is shown in Figure 18. It seems that Igan fits hn very well.Figure 18
The graph of the datan,lnhn and the fitting curve.The relative error ofIga is shown in Table 6 (when n≤1000) and Figure 19 (when 1000 <n≤10000).Table 6
The relative error ofIgan to hn when n≤1000.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—16−2.16%40−0.18%2202.62E − 06520−8.93E − 072—173.29%50−7.12E − 042401.60E − 06540−9.04E − 073—18−2.34%60−3.12E − 042608.76E − 07560−8.75E − 0748.89%192.71%70−1.01E − 042803.54E − 07580−8.50E − 07534.03%20−1.64%80−3.21E − 053003.00E − 09600−8.46E − 076−10.34%211.64%90−2.28E − 06320−2.87E − 07640−8.04E − 07721.05%22−1.00%1001.16E − 05340−4.83E − 07680−7.54E − 078−6.99%231.51%1101.51E − 05360−6.16E − 07720−7.43E − 0798.67%24−1.20%1201.51E − 05380−7.26E − 07760−6.98E − 0710−4.02%251.29%1301.40E − 05400−8.05E − 07800−6.69E − 07118.17%26−0.72%1401.22E − 05420−8.47E − 07840−6.06E − 0712−5.88%270.83%1501.04E − 05440−8.98E − 07880−5.64E − 07136.76%28−0.58%1608.81E − 06460−8.97E − 07920−5.65E − 0714−2.93%290.80%1806.04E − 06480−8.94E − 07960−5.17E − 07153.07%30−0.59%2004.02E − 06500−9.00E − 071000−4.75E − 07Figure 19
The relative error ofIgan when 1000 ≤n≤ 10000, step 300.Whenn<20, the relative error of Iga is still greater than 2%. Although it is much better than the error of Ig, it is not as good as expected when n<40. If we take the round approximation by(40)Iga′n=πexp2/3πn+a1/x+c1+b1122n3+12,n≥4,the relative error will be obviously smaller with a few exceptions, as shown in Table 7.Table 7
The relative error ofIgan+1/2 to hn when n≤30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—725%138.33%192.86%251.31%2—8014−2.94%20−1.46%26−0.63%3—912.5%152.44%211.82%270.87%4010016−1.82%22−0.95%28−0.56%550%117.14%173.03%231.58%290.83%6012−4.76%18−2.27%24−1.25%30−0.58%Later, we will find out that it is obviously greater than the relative error ofIg1 and Ig2 obtained in the next subsection by modifying the denominator part; when 4000<n<10000, the relative error of Iga is about 1000 times of that of Ig2.When 30<n≤ 1000, the relative error of Igan+1/2 is close to that of Igan.


### 3.3. Method B: Modifying the Denominator

Sincehn∼π/122n3exp2/3πn, we consider estimating hn byπ/122C3nexp2/3πn (i.e., fit π2exp2π2/3n/288h2n by a function C3n), where C3x is a cubic function or a function like(41)ax3+bx2.5+cx2+dx1.5+ex+fx0.5+g.However, the results are worse as the relative errors are obviously much greater than the relative error ofIgn when n<350.Then, we consider estimatinghn by π/122C4nexp2/3πn or fit πexpπ2/3n/122hn by a function as follows:(42)C4n=a4n1.5+b4n+c4n0.5+d4.The result is very good. The graph of the datan,πexpπ2/3n/122hn and the fitting curve C4n is shown in Figure 3. Here, the fitting curve is displayed by a thick continuous curve, which lies in the middle of the area the circles occupied. Since the circles are too crowded, the circles themselves look like a very thick curve.The values of the coefficients in the expression ofC4n are as follows:(43)a4=1.000010809,b4=1.862505234,c4=1.169930087,d4=−0.7005460222.The value ofa4 is very close to 1, which means that this fitting function coincides with the Ingham–Meinardus asymptotic formula very well.So, we have an estimation formula as follows:(44)hn∼Ig1n=π122C4nexp23πn.We may call it the Ingham–Meinardus revised estimation formula (1). The graph oflnIg1n is shown in Figure 2, together with the data points of n,lnhn. This revised estimation formula is much more accurate than the asymptotic formula. The relative error is less than 1×10−6 when n > 2000 (as shown in Figure 20 on page 14) and less than 3 when n≥30 (as shown in Table 8 on page 14). The relative error of the round approximation Ig1′n=Ig1n+1/2 is shown in Table 9 on page 13.Figure 20
The relative error ofIg1n when 1000 ≤n≤ 10000, step 300.Table 8
The relative error ofIg1n to hn when n≤1000.
nRel-Err (%)nRel-Err (%)nRel-ErrnRel-ErrnRel-Err1—16−1.6340−2.21E − 052207.23E − 055201.04E − 062−7.23173.82505.35E − 042405.74E − 055403.01E − 07329.9718−1.87606.16E − 042604.59E − 05560−3.53E − 074−8.44193.18706.15E − 042803.68E − 05580−8.95E − 07527.9420−1.21805.35E − 043002.97E − 05600−1.37E − 066−11.61212.06904.56E − 043202.40E − 05640−2.10E − 06720.7622−0.611003.89E − 043401.93E − 05680−2.64E − 068−6.74231.891103.30E − 043601.55E − 05720−2.97E − 0699.2124−0.851202.80E − 043801.24E − 05760−3.20E − 0610−3.44251.631302.40E − 044009.79E − 06800−3.36E − 06118.8526−0.401402.06E − 044207.63E − 06840−3.43E − 0612−5.28271.141501.78E − 044405.82E − 06880−3.51E − 06137.4228−0.291601.55E − 044604.32E − 06920−3.49E − 0614−2.35291.081801.18E − 044803.04E − 06960−3.43E − 06153.6630−0.322009.20E − 055001.97E − 061000−3.37E − 06Table 9
The relative error ofIg1n+1/2 to hn when n⩽30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—725%138.33%192.86%251.57%208014−2.94%20−1.46%26−0.42%30912.5%154.88%211.82%271.22%4010016−1.82%22−0.48%28−0.28%550%117.14%174.55%231.98%291.06%6012−4.76%18−2.27%24−0.94%30−0.29%However, equation (44) is not so satisfying when n<30, especially when n<15 as the relative error is not negligible for some value of n.As we already know thathn∼π/122n3exp2/3πn, or n3/2∼π/122hnexp2/3πn, which means that when fitting πexpπ2/3n/122hn by a function C4n shown in equation (42), the coefficient a4 should be exactly 1; hence, we should fit πexpπ2/3n/122hn by a function C4′n=n3/2+b5n+c5n1/2+d5 or fit πexpπ2/3n/122hn−n3/2 by a function(45)C5n=b5n+c5n1/2+d5.The graph of the datan,πexpπ2/3n/122hn−n3/2 is shown in Figure 21 on page 14 (together with the figure of the fitting function C5n generated by the least square method).Figure 21
The graph of the datan,πexpπ2/3n/122hn−n3/2 and the fitting curve C5n.The values of the coefficients in equation (45) are as follows:(46)b5=1.864260743,c5=1.084436400,d5=0.4754177757.So, we have another estimation formula forhn as follows:(47)hn∼Ig2n=πexp2/3πn122n3/2+C5n.We may call it the Ingham–Meinardus revised estimation formula (2). The graph oflnIg2n is nearly the same as that of lnIg1n shown in Figure 2. The second revised estimation formula is much more accurate than the first one. The relative error is less than 2×10−9 when n > 3000 (as shown in Figure 22 on page 14), about 1/500 of the relative error of Ig1n. When n < 10, the relative error is also distinctly less than that of Ig1n (as shown in Table 10 on page 15). The relative error of the round approximation Ig2′n=Ig2n+1/2 is shown in Table 11 (on page 15).Figure 22
The relative error ofIg2n when 1000 ≤n≤ 10000, step 300.Table 10
The relative error ofIg2n to hn when n≤1000.
NRel-Err (%)nRel-Err (%)nRel-ErrnRel-ErrnRel-Err116−2.4940−0.21%2202.15E − 065201.78E − 072−18.69172.9950−9.06E − 042401.76E − 065401.83E − 07319.7518−2.5960−4.33E − 042601.46E − 065601.59E − 074−13.56192.4870−1.80E − 042801.20E − 065801.51E − 07522.5120−1.8380−8.70E − 053009.94E − 076001.32E − 076−14.58211.4690−4.13E − 053208.51E − 076401.03E − 07717.4322−1.15100−1.68E − 053407.18E − 076805.80E − 088−8.89231.37110−5.98E − 063606.14E − 077206.80E − 0897.0624−1.32120−7.10E − 073805.23E − 077606.70E − 0810−5.09251.181302.07E − 064004.61E − 078005.10E − 08117.2326−0.821403.17E − 064203.90E − 078405.40E − 0812−6.53270.741503.54E − 064403.34E − 07880−4.30E − 09136.1628−0.661603.59E − 064602.96E − 079207.00E − 0914−3.38290.721803.16E − 064802.50E − 079602.80E − 08152.6730−0.662002.64E − 065002.22E − 0710003.30E − 08Table 11
The relative error ofIg2n+1/2 to hn when n≤30.
nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)11650.304032.1022013.105208.392146.241756.825028.6024012.505408.233202.891846.696025.9026012.005608.08495.591952.757023.9028011.505807.935156.432044.948022.3030011.106007.79668.622148.489020.9032010.806407.547121.382243.4710019.8034010.406807.31865.432346.0011018.8036010.107207.10988.382441.0912018.003809.867606.911062.582543.6813017.204009.608006.731179.472639.9314016.604209.368406.561253.292741.2715016.004409.148806.411370.982838.5016015.404608.939206.271453.122939.7018014.504808.749606.131560.353037.0020013.705008.5610006.01It should be mentioned that in Figure21 on page 14, the graph of the data points lies in a line, so we might be willing to fit this line by a first-order equation. The result is(48)C5′n=1.873818457×n+27.08318017.If we use this fitting function instead ofC5n generated before, the relative error to fit hn will be about 10000 times more, that is, about 20 times more than that of Ig1n. So, we did not use linear function to fit the data n,πexpπ2/3n/122hn−n3/2 before.


### 3.4. Method C: FittingIgn−hn

We wonder whether we can fitIgn−hn by a function rn and then estimate hn by Ign−rn which may be believed more accurate than Ig2n at the price of being more complicated.By the same tricks used at the beginning of this subsection, we will have(49)Ign−Ign−t∼tπ2243n2exp23πn.t≪n,So, we may fitIgn−hn by π2/243C6nexp2/3πn where C6n is a quadratic function or a function like(50)an2+bn1.5+cn+dn0.5+e.That means, we can fitπ2exp2/3πn/243Ign−hn by a function C6n. However, the result is useless. Although C6n will fit the data π2exp2/3πn/243Ign−hn very well, the relative error of Ign−π2/243C6nexp2/3πn to hn is much greater than that of Ig1n or Ig2n, and the relative error differs very little with that of Ign when n is small. Besides, the formula Ign−π2/243C6nexp2/3πn is much more complicated than Ig1n and Ig2n.If we fitIgn−hn by π2/243n2+E6nexp2/3πn with the trick described in subsection 3.3, where E6n is a function in the form as follows:(51)bn1.5+cn+dn0.5+e,as we already know the coefficient of n2 should be 1 in theory. The result will be a little better, but useless too. The accuracy is not as good as that of Ig0n.Then, we consider fittingπ2exp2/3πn/243n2Ign−hn by a function C7n. If C7n is in the form a/n+b or a/n+b/n2+c, the result is useless, too. If C7n is in the form a/n0.5+b, it will be barely satisfactory. If C7n is in the form a/n0.5+b/n+c/n1.5+d/n2+e or a/n0.5+b/n+c/n1.5+e, the result will be much better than the previous forms, but the accuracy (when estimating hn) is not as good as that of Ig1n and Ig2n.The result ofC7n is(52)C7an=0.8782296151n0.5+0.2567016063n−3.580442785n1.5+21.28305831n2+0.6879945549,or(53)C7bn=0.8861039149n0.5−0.05719053203n+0.9843423289n1.5+0.6879343652.The relative errors of(54)F7an=Ign−π2exp2/3πn243n2C7an,and(55)F7bn=Ign−π2exp2/3πn243n2C7bn.to hn when 1000 ≤n≤ 10000 are shown in Figures 23 and 24 (page 16), respectively. In this interval (1000, 10000), F7an is obviously more accurate than F7bn. When n≤ 1000, the relative error of F7an+1/2 and F7bn+1/2 is shown in Table 12 (page 15) and Table 13 (page 15). In this case, F7bn is better than F7an. However, neither of them is as good as Ig1n or Ig2n although they are more complicated than Ig1n and Ig2n.Figure 23
The relative error ofF7an when 1000 ≤n≤ 10000, step 300.Figure 24
The relative error ofF7bn when 1000 ≤n≤ 10000, step 300.Table 12
The relative error ofF7an+1/2 to hn when n≤1000.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—16040−8.13E − 04220−1.74E − 065203.10E − 072100%174.55%50−3.26E − 04240−1.51E − 065403.29E − 073100%18−1.14%60−1.49E − 04260−1.25E − 065603.42E − 07450%193.81%70−2.81E − 05280−9.91E − 075803.51E − 07550%20−0.73%80−4.62E − 06300−7.56E − 076003.56E − 0760212.42%903.46E − 06320−5.49E − 076403.57E − 07750%22−0.48%1006.70E − 06340−3.72E − 076803.49E − 0780231.98%1105.27E − 06360−2.23E − 077203.36E − 07912.5%24−0.63%1203.37E − 06380−9.87E − 087603.19E − 07100251.83%1301.93E − 064003.54E − 098003.00E − 071114.29%26−0.21%1405.77E − 074208.70E − 088402.80E − 07120271.22%150−4.01E − 074401.55E − 078802.59E − 07138.33%28−0.28%160−1.04E − 064602.09E − 079202.39E − 07140291.06%180−1.72E − 064802.51E − 079602.19E − 07154.88%30−0.29%200−1.86E − 065002.85E − 0710001.99E − 07Table 13
The relative error ofF7bn+1/2 to hn when n⩽1000.
nRel-ErrnRel-Err (%)nRel-ErrnRel-ErrnRel-Err1—16−1.8240−0.16%220−1.21E − 06520−1.11E − 0620173.0350−6.19E − 04240−1.77E − 06540−1.01E − 063018−2.2760−2.60E − 04260−2.08E − 06560−9.24E − 0740192.8670−7.50E − 05280−2.21E − 06580−8.41E − 07550%20−1.4680−1.85E − 05300−2.24E − 06600−7.64E − 0760211.82903.62E − 06320−2.20E − 06640−6.25E − 07725%22−0.951001.30E − 05340−2.12E − 06680−5.05E − 0780231.581101.37E − 05360−2.02E − 06720−4.00E − 07912.5%24−1.251201.21E − 05380−1.90E − 06760−3.10E − 07100251.311301.01E − 05400−1.78E − 06800−2.31E − 07117.14%26−0.631407.74E − 06420−1.66E − 06840−1.63E − 0712−4.76%270.871505.68E − 06440−1.54E − 06880−1.04E − 07138.33%28−0.561603.97E − 06460−1.42E − 06920−5.24E − 0814−2.94%290.831801.40E − 06480−1.31E − 06960−7.83E − 09152.44%30−0.58200−2.22E − 07500−1.21E − 0610003.08E − 08


### 3.5. Estimatehn When n≤100

All the estimation functions forhn found now are with very good accuracy when n is greater than 1000, but they are not so accurate when n<50, especially when n<25. Although Ig1′n and Ig2′n are better than others, the relative error is still greater than 1 when n<40.Whenn<40, it is too difficult to fit πexpπ2/3n/122hn−n3/2 by a simple smooth function with high accuracy, as shown in Figure 25 (on page 16). The points n,πexpπ2/3n/122hn−n3/2 (n = 3, 4, ⋯, 100) are not so complicated (as shown in Figure 25). It seems that we can fit them by a simple piecewise function with 2 pieces, as the even points (where n is even) lie roughly on a smooth curve, so do the odd points. If we try to fit them, respectively, we will have the fitting function as follows:(56)C8x=1.942141112×x−0.4796781366×x+8.291226268,n=3,5,7,…,99;1.803056782×x+2.356539877×x−6.043824511,n=4,6,8,…,100.Figure 25
The graph of the datan,C8n.Hence, we can calculatehn (3≤n≤100) by(57)hn∼Ig0n=πexp2/3πn122n3/2+C8n,3≤n≤100.Considering thathn is an integer, we can take the round approximation of equation (57):(58)Ig0′n=πexp2/3πn122n3/2+C8n+12,3≤n≤100.Here,n begins from 3, not 1 or 2, because Ig0′1−h1/h1 is meaningless since h1=0, and Ig0′2 differs from h2 a lot. Besides, the value of h1 and h2 is clear by definition, so there is no need to use a complicated formula to estimate them.The relative error ofIg0n (or Ig0′n) to hn is shown in Table 14 (or Table 15) on page 17. Compared them with Table 11 on page 15, we will find that when n≥80, Ig2′n is more accurate than Ig0′n; when n<80, Ig0′n is better.Table 14
The relative error ofIg0n to hn when n≤100.
nRel-Err (%)nRel-ErrnRel-ErrnRel-ErrnRel-Err1—21−0.16%41−5.72E − 0461−1.04E − 04817.19E − 052—220.12%423.75E − 04621.43E − 0482−5.18E − 053−14.85230.05%43−4.68E − 0463−1.04E − 04838.12E − 05412.7224−0.28%444.91E − 04641.18E − 0484−6.50E − 0551.99258.37E − 0445−6.29E − 0465−5.63E − 05858.69E − 056−1.83264.75E − 04465.45E − 04667.10E − 0586−6.90E − 0574.7627−0.17%47−4.45E − 0467−2.40E − 05878.65E − 058−0.64286.69E − 04483.17E − 04685.74E − 0588−7.40E − 059−0.9229−3.78E − 0449−3.42E − 0469−1.67E − 05899.01E − 05100.6930−4.43E − 04503.79E − 04704.05E − 0590−8.07E − 05111.4431−1.98E − 0451−3.98E − 04711.30E − 05919.09E − 0512−2.46324.21E − 04523.55E − 04726.81E − 0692−8.24E − 05131.8633−0.12%53−2.81E − 04733.41E − 05938.80E − 0514−0.24349.27E − 04542.47E − 0474−1.74E − 0694−8.37E − 0515−0.5435−7.86E − 0455−2.21E − 04753.84E − 05958.70E − 0516−0.04361.82E − 04562.44E − 0476−1.56E − 0596−8.65E − 05170.4637−4.80E − 0457−2.25E − 04775.70E − 05978.44E − 0518−0.67386.53E − 04582.29E − 0478−3.58E − 0598−8.56E − 05190.4739−9.11E − 0459−1.55E − 04796.90E − 05997.92E − 0520−0.28406.34E − 04601.44E − 0480−4.41E − 05100−8.48E − 05Table 15
The relative error ofIg0′n to hn when n≤100.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—9017025033−0.11%2—10018−1.14%260349.23E-043011019027−0.17%35−7.77E-044012−4.76%200280363.23E-045013021029037−5.46E-0460140220300386.85E-047015023031039−9.67E-048016024−0.31%320.07%406.50E-04


## 3.1. Asymptotic Formula

Ashn=pn−pn−1, by Hardy–Ramanujan’s asymptotic formula,(23)pn∼14n3exp23πn,(refer [1, 3, 31–35]), we assume that, when n≫1, hn∼exp2/3πn/43n−exp2/3πn−1/43n−1. So,(24)hn∼exp2/3πn−143expπ2/3n−n−1n−1n−1=exp2/3πn−143expπ2/3/n+n−1n−1n−1∼exp2/3πn43expπ2/3/2nn−1n−1∼exp2/3πn431+π/6nn−1n−1∼exp2/3πn43π6n3=πexp2/3πn122n3.So,(25)hn∼π122n3exp2/3πn.In coincidence, half a year after the main results were obtained in this paper, the author found an asymptotic formula as follows:(26)Pa,bn∼Γb/aπb/a−12−3/2−b/2a3−b/2aa−1/2+b/2an−a+b/2aeπ2n/3a,in [30]. When a = 1 and b = 2, we will have(27)P1,2n∼π122n3expπ2/3n,which coincides with the asymptotic formula obtained here.Formula (25) will be called the Ingham–Meinardus asymptotic formula in this paper since Daniel mentioned in [30] that this general asymptotic formula (26) was first given by A. E. Ingham in [29] and the proof was refined by G. Meinardus later in another two papers written in German.Later in this paper,π/122n3exp2/3πn will be denoted by Ign for short.It is not satisfying to estimatehn by Ign when n is small. The relative error of Ign to hn is greater than 6% as shown in Table 4 (on page 6). The round approximation(28)Ig′n=Ign+12,will not change the accuracy distinctly, as shown in Table 5 (on page 6). So, it is necessary to modify the asymptotic formula for better accuracy.Table 4
The relative error ofIgn to hn when n≤1000.
nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)11650.304032.1022013.105208.392146.241756.825028.6024012.505408.233202.891846.696025.9026012.005608.08495.591952.757023.9028011.505807.935156.432044.948022.3030011.106007.79668.622148.489020.9032010.806407.547121.382243.4710019.8034010.406807.31865.432346.0011018.8036010.107207.10988.382441.0912018.003809.867606.911062.582543.6813017.204009.608006.731179.472639.9314016.604209.368406.561253.292741.2715016.004409.148806.411370.982838.5016015.404608.939206.271453.122939.7018014.504808.749606.131560.353037.0020013.705008.5610006.01Table 5
The relative error ofIgn+1/2 to hn when n≤30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—7125%1370.83%1952.38%2543.60%2100%871.43%1452.94%2045.26%2639.96%3200%987.50%1560.98%2148.48%2741.29%4100%1066.67%1650.91%2243.33%2838.56%5150%1178.57%1757.58%2345.85%2939.67%675%1252.38%1846.59%2440.94%3037.05%


## 3.2. Method A: Modifying the Exponent

In this subsection, we consider fittinghn by Iga=π/122n3exp2/3πn+C1n or fitting.n,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) by a function(29)f1x≐a1x+c1e1+b1.LetC1n = 3/2π2ln122n3hn/π2−n. The reason that we fit C1x by the function in the form displayed in (29) is the same as that described in Section 3 of [1] (although the data differ distinctly). Many other types of functions have been tested, but they cannot fit these data very well.However, here it is not valid to obtain the constants inf1n by the iteration method described in reference [1].The graph of the datan,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) is shown in Figure 4 on page 7.Figure 4
The graph of the datan,3/2π2ln122n3hn/π2−n.First, we try to fitn,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) by a function in the form as follows:(30)f1x≐a1x+c1+b1.That means we have assumed thate1=1/2, temporarily. We will explain the reason in Subsection 3.2.3.The average error off1x is(31)E1=1K1∑nC1n−f1n2=1K1∑nC1n−a1n+c1−b12E1=1K1∑k=1K1C160+20k−a160+20k+c1−b12,where K1 = 397 and n ranges from 80 to 8000, by step 20. Here, only a1, b1, and c1 are unknown, so we can consider E1 as a function of the variables a1,b1,c1.We want to find a triplea1,b1,c1 such that E1 reaches its minimum or to make E1 as small as possible.Since a lot of functions have several local minimum points, it is necessary to find out whetherE1 = E1a1,b1,c1 has more than one local minimum before we start to calculate the minimum point by the numeral method. However, E1 = E1a1,b1,c1 is too complicate, and it is very difficult to know all the critical points in the range we are considering.


### 3.2.1. Preparation Work

For a given paira1,c1, by the property of the arithmetic mean (for some given data xi (i = 1, 2, ⋯, k; xi∈ℝ), the function st=1/k∑i=1kxi−t2 reaches its minimum at t = 1/k∑i=1kxi), it is clear that E1 reaches its minimum when(32)b1=1K1∑nC1n−a1n+c1,=C1¯−1K1∑n′a1n′+c1,where C1¯ = 1/K1∑n′C1n′.LetE1=1/K1∑nC1n−a1/n+c1−C1¯+1/K1∑n′a1/n′+c12 be the average error of the fitting function f1x≐a1/x+c1+C1¯−1/K1∑n′a1/n′+c1. (Here, a1 and b1 are undetermined coefficients.)LetG1=E12=1/K1∑nC1n−a1/n+c1−C1¯+1/K1∑n′a1/n′+c12. Here, only a1 and b1 are unknowns. G1 could be considered as a function of a1,c1. In order to find the minimum point of G1, we can draw the figure of the function G1=G1a1,c1 (in a cube coordinate system with axis a1,c1, and G1) as shown in Figures 5–7. Figures 8 and 9 are the projection of the graph of a1,c1,G1 (when −100≤a1≤ 100 and −50≤c1≤ 100, which is a part of a surface) on the a1−G1 plane (spanned by the axis a1 and G1) and c1−G plane, respectively.Figure 5
The graph of the dataa1,c1,G1 when −100≤a1≤ 100 and −50≤c1≤ 100.Figure 6
The graph of the dataa1,c1,G1 when −100≤a1≤ 100 and −50≤c1≤ 100 version 2.Figure 7
The graph of the dataa1,c1,G1 when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.Figure 8
The projection of the graph of the dataa1,c1,G1 on the a1−G1 plane when −100≤a1≤ 100 and −50≤c1≤ 100.Figure 9
The projection of the graph of the dataa1,c1,G1 on the c1−G1 plane when −100≤a1≤ 100 and −50≤c1≤ 100.From these figures, we can find out that the influence ofc1 to G1 is much less than that of a1. In Figure 10, we find that when G1 reaches its minimum, a1 is between 0.50 and 0.53, but there is not a definite range for c1.Figure 10
The projection of the graph of the dataa1,c1,G1 on the c1−G1 plane when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.It is possible that for different range ofc1, the range of a1 when G1 reaches its minimum will be different. However, considering that a1/n+c1+b1 is a real, c1 should be greater than −1 in theory. (For the fitting data used here, c1 should be greater than −80.) From Figure 8, we can see that G1 touches its bottom when −15≤a1≤ 15. Although we cannot see clearly the exact value of of a1 in the minimum points, we can draw another figure of a1,c1,G1 when −15≤a1≤ 15 and −1≤c1≤ 100 to observe more details (the figure is not presented here), and then we will find that the more detailed range of a1 for the minimum points is −3,3 in the new figure (not presented), and next we will draw the figures of a1,c1,G1 when −3≤a1≤ 3, 0≤a1≤ 1, 0.2≤a1≤ 0.8, or 0.45≤a1≤ 0.6, respectively, while −1≤c1≤ 100, we will find the range of a1 of the minimum points of G1 is about [0.50, 0.53]. The projections of Figure 7 of a1,c1,G1 when 0.45≤a1≤ 0.6 and −1≤c1≤ 100 are shown in Figures 11 and 10.Figure 11
The projection of the graph of the dataa1,c1,G1 on the a1−G1 plane when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.In Figure9, for the curves on the bottom, G1 decreases with c1 at first then increases with c1, but it is difficult to find the critical point in the figure since different curves have different critical points.Although we cannot find a satisfying value ofa1 or c1 from the figures to construct a fitting function f1, these pictures show that the figure of G1=a1,c1,G1 has only one bottom in the domain we are considering, unlike the figure of another function shown in Figure 12, so the existence of the minimum point is almost confirmed; therefore, we are confident to find the value of a1 or c1 in the minimum point by the numerical method. This guarantees the validity of the numerical calculation by loop in the next step.Figure 12
Example of a surface with more than one bottom.


### 3.2.2. Findc1

On the other hand, by the least square method, to fit the dataxk,yk (k = 1, 2, ⋯, K1) by a linear function y=a×x+b, the result is that(33)a=xn⋅yn¯−xn¯⋅yn¯xn2¯−xn¯2,(34)b=xn¯⋅xn⋅yn¯−xn2¯⋅yn¯xn¯2−xn2¯=yn¯−xn¯⋅a,where xn¯ = 1/K1∑i=1K1xi is the average value of xi, yn¯ = 1/K1∑i=1K1yi, xi2¯ = 1/K1∑i=1K1xi2, and xn⋅yn¯ = 1/K1∑i=1K1xi⋅yi. By definition, xn¯2 = 1/K1∑i=1K1xi2 is the square of the average value of xn. So, by the least square method, a and b are uniquely determined by the given data xk,yk (k = 1, 2, ⋯, K1).For every given value ofc1 (greater than −80), we can fit n,C1n (n = 60+20k, k = 1, 2, ⋯, 397) by a function f1x≐a1/x+c1+b1 by the least square method if we consider 1/60+20k+c1 and C160+20k as xk and yk, respectively. Then,(35)xn¯=1K1∑k=1K1xk=1K1∑k=1K1160+20k+c1,yn¯=1K1∑k=1K1yk=1K1∑k=1K1C160+20k,xn2¯=1K1∑k=1K1xk2=1K1∑k=1K1160+20k+c1,xn⋅yn¯=1K1∑k=1K1xk⋅yk=1K1∑k=1K1C160+20k60+20k+c1,a1=xn⋅yn¯−xn¯⋅yn¯xn2¯−xn¯2,b1=yn¯−xn¯⋅a1.So,a1 and b1 could both be considered as functions of c1, denoted by a1=a1c1 and b1=b1c1, since they are uniquely determined by c1 with the given data.Then,G2=E12=1/K1∑k=1K1C160+20k−a1c1/60+20k+c1−b1c12 is a function of c1.It will cost some time to plot the figure of the functionG2 = G2c1 in a CAS software.If we plot the figure of the functionG2 = G2c1 on the coordinates (as shown in Figures 13–15), we will find that G2 reaches its minimum when c1≈−3.2594807. In Figure 15, we find that the curve of G2 = G2c1 is not so smooth. The reason is that we hold up 18 significant digits in the process. If we compute more significant digits in the process, the curve on Figure 15 will be more smooth, at the cost of much more time. By writing a small program (since the default function to find the minimum provided by the software Maple 18 is unable to deal with such a complicated function G2 = G2c1 involving so much data), we can obtain a more accurate value of the critical point as follows:(36)c1=−3.259480684.Figure 13
The graph of the functionG2 = G2c1 when −50≤c1≤ 100.Figure 14
The graph of the functionG2 = G2c1 when −4≤c1≤ −2.Figure 15
The graph of the functionG2 = G2c1 when −3.259483≤c1≤−3.259478.When the value ofc1 is obtained, we can find the value of a1 and b1 by the least square method without difficulty, i.e.,(37)a1=0.5097429624,b1=−1.453552800.However, herec1 is less than −1, so the estimation formula for hn constructed from these coefficients is invalid when n<4.


### 3.2.3. Confirme1

In [1], we fit C1n = 3/2⋅ln4n3pn2/π2−n by a function f1x≐a1/n+c1e1+b1 when estimating pn and found that e1≈0.50 by iteration. Here, the iteration method does not work well, so we fit C1n = 3/2π2ln122n3hn/π2−n by a function f1x≐a1/n+c1+b1 directly, which means that we have assumed that e1=1/2. Here, we may doubt that whether e1=0.5 is the best option for us?Here, we use the same idea described in Subsection3.2.1.For every paire1,c1, we can obtain corresponding a1 and b1 by the least square method, just like (33) and (34), except that here xk=1/60+20k+c1e1.So, the square of the average error(38)G3=E12=1K1∑nC1n−a1e1,c1n+c1e1−b1e1,c12,could be considered as a function of e1 and c1 as both a1=a1e1,c1 and b1=b1e1,c1 could be expressed by certain elementary functions of e1 and c1.If we draw the figure of the functionG3 = G3e1,c1, we will find that the surface has only one bottom when 0.1 ≤e1≤0.9 and −50≤c1≤100, as shown in Figure 16. However, the process to draw the figure is time-consuming. It costs more than 5 hours on a notebook (ThinkPad E40 Edge, with 6 GB RAM and AMD P360 Dual-Core Processor 2.30 GHz) by Maple 18 in Ubuntu 14.04.1 system.Figure 16
The graph of the functionG3 = G3e1,c1 when 0.1≤e1≤0.9 and −50≤c1≤100.After that, by writing another program, we can obtain the approximate value ofe1,c1 where G3 touches the bottom, i.e., e1≈0.494 and c1≈−4.85, when 18 significant digits are involved in the process, which still costs tens of minutes. Considering that we have used only a small part of data, we cannot afford the time for computing more significant digits in process, and the computing is so complicated; hence, error accumulation effect is considerable, so we choose e1=0.50 while it differs very little with 0.494. Another reason is that we prefer simple exponent, as the time spent on computing a square root is much less than that to compute a power with exponent 0.494 in general. Here, the value of c1≈−4.85 is obviously different from the value obtained at the end of Subsection 3.2.2 because of the little difference on e1. Therefore, it will be fine to use the result in Subsection 3.2.2).In this figure, the points are shown as small circles which are very close to each other. These crowded circles seem like a thick curve. A fitting curve passes though the center of these circles. The fitting curve might not be found in reduce printing. That means the curve fits the points (displayed as circles) very well.


### 3.2.4. The Result

By the valuea1 = 0.5097429624, b1 = −1.453552800, and c1 = −3.259480684 obtained in Subsection 3.2.2, we will have a fitting function f1x≐a1/x+c1+b1.The graph off1x (when −80≤x≤ 8000) and the graph of the data n,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) are shown in Figure 17. It shows that f1n fits 3/2π2ln122n3hn/π2−n very well.Figure 17
The graph of the datan,3/2π2ln122n3hn/π2−n and the fitting curve.Then, we could fithn by(39)Igan=πexp2/3πn+a1/x+c1+b1122n3,n≥4.The graph of the functionlnIgax and the graph of the data n,lnhn (n = 60+20k, k = 1, 2, ⋯, 397) is shown in Figure 18. It seems that Igan fits hn very well.Figure 18
The graph of the datan,lnhn and the fitting curve.The relative error ofIga is shown in Table 6 (when n≤1000) and Figure 19 (when 1000 <n≤10000).Table 6
The relative error ofIgan to hn when n≤1000.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—16−2.16%40−0.18%2202.62E − 06520−8.93E − 072—173.29%50−7.12E − 042401.60E − 06540−9.04E − 073—18−2.34%60−3.12E − 042608.76E − 07560−8.75E − 0748.89%192.71%70−1.01E − 042803.54E − 07580−8.50E − 07534.03%20−1.64%80−3.21E − 053003.00E − 09600−8.46E − 076−10.34%211.64%90−2.28E − 06320−2.87E − 07640−8.04E − 07721.05%22−1.00%1001.16E − 05340−4.83E − 07680−7.54E − 078−6.99%231.51%1101.51E − 05360−6.16E − 07720−7.43E − 0798.67%24−1.20%1201.51E − 05380−7.26E − 07760−6.98E − 0710−4.02%251.29%1301.40E − 05400−8.05E − 07800−6.69E − 07118.17%26−0.72%1401.22E − 05420−8.47E − 07840−6.06E − 0712−5.88%270.83%1501.04E − 05440−8.98E − 07880−5.64E − 07136.76%28−0.58%1608.81E − 06460−8.97E − 07920−5.65E − 0714−2.93%290.80%1806.04E − 06480−8.94E − 07960−5.17E − 07153.07%30−0.59%2004.02E − 06500−9.00E − 071000−4.75E − 07Figure 19
The relative error ofIgan when 1000 ≤n≤ 10000, step 300.Whenn<20, the relative error of Iga is still greater than 2%. Although it is much better than the error of Ig, it is not as good as expected when n<40. If we take the round approximation by(40)Iga′n=πexp2/3πn+a1/x+c1+b1122n3+12,n≥4,the relative error will be obviously smaller with a few exceptions, as shown in Table 7.Table 7
The relative error ofIgan+1/2 to hn when n≤30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—725%138.33%192.86%251.31%2—8014−2.94%20−1.46%26−0.63%3—912.5%152.44%211.82%270.87%4010016−1.82%22−0.95%28−0.56%550%117.14%173.03%231.58%290.83%6012−4.76%18−2.27%24−1.25%30−0.58%Later, we will find out that it is obviously greater than the relative error ofIg1 and Ig2 obtained in the next subsection by modifying the denominator part; when 4000<n<10000, the relative error of Iga is about 1000 times of that of Ig2.When 30<n≤ 1000, the relative error of Igan+1/2 is close to that of Igan.


## 3.2.1. Preparation Work

For a given paira1,c1, by the property of the arithmetic mean (for some given data xi (i = 1, 2, ⋯, k; xi∈ℝ), the function st=1/k∑i=1kxi−t2 reaches its minimum at t = 1/k∑i=1kxi), it is clear that E1 reaches its minimum when(32)b1=1K1∑nC1n−a1n+c1,=C1¯−1K1∑n′a1n′+c1,where C1¯ = 1/K1∑n′C1n′.LetE1=1/K1∑nC1n−a1/n+c1−C1¯+1/K1∑n′a1/n′+c12 be the average error of the fitting function f1x≐a1/x+c1+C1¯−1/K1∑n′a1/n′+c1. (Here, a1 and b1 are undetermined coefficients.)LetG1=E12=1/K1∑nC1n−a1/n+c1−C1¯+1/K1∑n′a1/n′+c12. Here, only a1 and b1 are unknowns. G1 could be considered as a function of a1,c1. In order to find the minimum point of G1, we can draw the figure of the function G1=G1a1,c1 (in a cube coordinate system with axis a1,c1, and G1) as shown in Figures 5–7. Figures 8 and 9 are the projection of the graph of a1,c1,G1 (when −100≤a1≤ 100 and −50≤c1≤ 100, which is a part of a surface) on the a1−G1 plane (spanned by the axis a1 and G1) and c1−G plane, respectively.Figure 5
The graph of the dataa1,c1,G1 when −100≤a1≤ 100 and −50≤c1≤ 100.Figure 6
The graph of the dataa1,c1,G1 when −100≤a1≤ 100 and −50≤c1≤ 100 version 2.Figure 7
The graph of the dataa1,c1,G1 when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.Figure 8
The projection of the graph of the dataa1,c1,G1 on the a1−G1 plane when −100≤a1≤ 100 and −50≤c1≤ 100.Figure 9
The projection of the graph of the dataa1,c1,G1 on the c1−G1 plane when −100≤a1≤ 100 and −50≤c1≤ 100.From these figures, we can find out that the influence ofc1 to G1 is much less than that of a1. In Figure 10, we find that when G1 reaches its minimum, a1 is between 0.50 and 0.53, but there is not a definite range for c1.Figure 10
The projection of the graph of the dataa1,c1,G1 on the c1−G1 plane when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.It is possible that for different range ofc1, the range of a1 when G1 reaches its minimum will be different. However, considering that a1/n+c1+b1 is a real, c1 should be greater than −1 in theory. (For the fitting data used here, c1 should be greater than −80.) From Figure 8, we can see that G1 touches its bottom when −15≤a1≤ 15. Although we cannot see clearly the exact value of of a1 in the minimum points, we can draw another figure of a1,c1,G1 when −15≤a1≤ 15 and −1≤c1≤ 100 to observe more details (the figure is not presented here), and then we will find that the more detailed range of a1 for the minimum points is −3,3 in the new figure (not presented), and next we will draw the figures of a1,c1,G1 when −3≤a1≤ 3, 0≤a1≤ 1, 0.2≤a1≤ 0.8, or 0.45≤a1≤ 0.6, respectively, while −1≤c1≤ 100, we will find the range of a1 of the minimum points of G1 is about [0.50, 0.53]. The projections of Figure 7 of a1,c1,G1 when 0.45≤a1≤ 0.6 and −1≤c1≤ 100 are shown in Figures 11 and 10.Figure 11
The projection of the graph of the dataa1,c1,G1 on the a1−G1 plane when 0.45≤a1≤ 0.60 and −1≤c1≤ 100.In Figure9, for the curves on the bottom, G1 decreases with c1 at first then increases with c1, but it is difficult to find the critical point in the figure since different curves have different critical points.Although we cannot find a satisfying value ofa1 or c1 from the figures to construct a fitting function f1, these pictures show that the figure of G1=a1,c1,G1 has only one bottom in the domain we are considering, unlike the figure of another function shown in Figure 12, so the existence of the minimum point is almost confirmed; therefore, we are confident to find the value of a1 or c1 in the minimum point by the numerical method. This guarantees the validity of the numerical calculation by loop in the next step.Figure 12
Example of a surface with more than one bottom.


## 3.2.2. Findc1

On the other hand, by the least square method, to fit the dataxk,yk (k = 1, 2, ⋯, K1) by a linear function y=a×x+b, the result is that(33)a=xn⋅yn¯−xn¯⋅yn¯xn2¯−xn¯2,(34)b=xn¯⋅xn⋅yn¯−xn2¯⋅yn¯xn¯2−xn2¯=yn¯−xn¯⋅a,where xn¯ = 1/K1∑i=1K1xi is the average value of xi, yn¯ = 1/K1∑i=1K1yi, xi2¯ = 1/K1∑i=1K1xi2, and xn⋅yn¯ = 1/K1∑i=1K1xi⋅yi. By definition, xn¯2 = 1/K1∑i=1K1xi2 is the square of the average value of xn. So, by the least square method, a and b are uniquely determined by the given data xk,yk (k = 1, 2, ⋯, K1).For every given value ofc1 (greater than −80), we can fit n,C1n (n = 60+20k, k = 1, 2, ⋯, 397) by a function f1x≐a1/x+c1+b1 by the least square method if we consider 1/60+20k+c1 and C160+20k as xk and yk, respectively. Then,(35)xn¯=1K1∑k=1K1xk=1K1∑k=1K1160+20k+c1,yn¯=1K1∑k=1K1yk=1K1∑k=1K1C160+20k,xn2¯=1K1∑k=1K1xk2=1K1∑k=1K1160+20k+c1,xn⋅yn¯=1K1∑k=1K1xk⋅yk=1K1∑k=1K1C160+20k60+20k+c1,a1=xn⋅yn¯−xn¯⋅yn¯xn2¯−xn¯2,b1=yn¯−xn¯⋅a1.So,a1 and b1 could both be considered as functions of c1, denoted by a1=a1c1 and b1=b1c1, since they are uniquely determined by c1 with the given data.Then,G2=E12=1/K1∑k=1K1C160+20k−a1c1/60+20k+c1−b1c12 is a function of c1.It will cost some time to plot the figure of the functionG2 = G2c1 in a CAS software.If we plot the figure of the functionG2 = G2c1 on the coordinates (as shown in Figures 13–15), we will find that G2 reaches its minimum when c1≈−3.2594807. In Figure 15, we find that the curve of G2 = G2c1 is not so smooth. The reason is that we hold up 18 significant digits in the process. If we compute more significant digits in the process, the curve on Figure 15 will be more smooth, at the cost of much more time. By writing a small program (since the default function to find the minimum provided by the software Maple 18 is unable to deal with such a complicated function G2 = G2c1 involving so much data), we can obtain a more accurate value of the critical point as follows:(36)c1=−3.259480684.Figure 13
The graph of the functionG2 = G2c1 when −50≤c1≤ 100.Figure 14
The graph of the functionG2 = G2c1 when −4≤c1≤ −2.Figure 15
The graph of the functionG2 = G2c1 when −3.259483≤c1≤−3.259478.When the value ofc1 is obtained, we can find the value of a1 and b1 by the least square method without difficulty, i.e.,(37)a1=0.5097429624,b1=−1.453552800.However, herec1 is less than −1, so the estimation formula for hn constructed from these coefficients is invalid when n<4.


## 3.2.3. Confirme1

In [1], we fit C1n = 3/2⋅ln4n3pn2/π2−n by a function f1x≐a1/n+c1e1+b1 when estimating pn and found that e1≈0.50 by iteration. Here, the iteration method does not work well, so we fit C1n = 3/2π2ln122n3hn/π2−n by a function f1x≐a1/n+c1+b1 directly, which means that we have assumed that e1=1/2. Here, we may doubt that whether e1=0.5 is the best option for us?Here, we use the same idea described in Subsection3.2.1.For every paire1,c1, we can obtain corresponding a1 and b1 by the least square method, just like (33) and (34), except that here xk=1/60+20k+c1e1.So, the square of the average error(38)G3=E12=1K1∑nC1n−a1e1,c1n+c1e1−b1e1,c12,could be considered as a function of e1 and c1 as both a1=a1e1,c1 and b1=b1e1,c1 could be expressed by certain elementary functions of e1 and c1.If we draw the figure of the functionG3 = G3e1,c1, we will find that the surface has only one bottom when 0.1 ≤e1≤0.9 and −50≤c1≤100, as shown in Figure 16. However, the process to draw the figure is time-consuming. It costs more than 5 hours on a notebook (ThinkPad E40 Edge, with 6 GB RAM and AMD P360 Dual-Core Processor 2.30 GHz) by Maple 18 in Ubuntu 14.04.1 system.Figure 16
The graph of the functionG3 = G3e1,c1 when 0.1≤e1≤0.9 and −50≤c1≤100.After that, by writing another program, we can obtain the approximate value ofe1,c1 where G3 touches the bottom, i.e., e1≈0.494 and c1≈−4.85, when 18 significant digits are involved in the process, which still costs tens of minutes. Considering that we have used only a small part of data, we cannot afford the time for computing more significant digits in process, and the computing is so complicated; hence, error accumulation effect is considerable, so we choose e1=0.50 while it differs very little with 0.494. Another reason is that we prefer simple exponent, as the time spent on computing a square root is much less than that to compute a power with exponent 0.494 in general. Here, the value of c1≈−4.85 is obviously different from the value obtained at the end of Subsection 3.2.2 because of the little difference on e1. Therefore, it will be fine to use the result in Subsection 3.2.2).In this figure, the points are shown as small circles which are very close to each other. These crowded circles seem like a thick curve. A fitting curve passes though the center of these circles. The fitting curve might not be found in reduce printing. That means the curve fits the points (displayed as circles) very well.


## 3.2.4. The Result

By the valuea1 = 0.5097429624, b1 = −1.453552800, and c1 = −3.259480684 obtained in Subsection 3.2.2, we will have a fitting function f1x≐a1/x+c1+b1.The graph off1x (when −80≤x≤ 8000) and the graph of the data n,3/2π2ln122n3hn/π2−n (n = 60+20k, k = 1, 2, ⋯, 397) are shown in Figure 17. It shows that f1n fits 3/2π2ln122n3hn/π2−n very well.Figure 17
The graph of the datan,3/2π2ln122n3hn/π2−n and the fitting curve.Then, we could fithn by(39)Igan=πexp2/3πn+a1/x+c1+b1122n3,n≥4.The graph of the functionlnIgax and the graph of the data n,lnhn (n = 60+20k, k = 1, 2, ⋯, 397) is shown in Figure 18. It seems that Igan fits hn very well.Figure 18
The graph of the datan,lnhn and the fitting curve.The relative error ofIga is shown in Table 6 (when n≤1000) and Figure 19 (when 1000 <n≤10000).Table 6
The relative error ofIgan to hn when n≤1000.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—16−2.16%40−0.18%2202.62E − 06520−8.93E − 072—173.29%50−7.12E − 042401.60E − 06540−9.04E − 073—18−2.34%60−3.12E − 042608.76E − 07560−8.75E − 0748.89%192.71%70−1.01E − 042803.54E − 07580−8.50E − 07534.03%20−1.64%80−3.21E − 053003.00E − 09600−8.46E − 076−10.34%211.64%90−2.28E − 06320−2.87E − 07640−8.04E − 07721.05%22−1.00%1001.16E − 05340−4.83E − 07680−7.54E − 078−6.99%231.51%1101.51E − 05360−6.16E − 07720−7.43E − 0798.67%24−1.20%1201.51E − 05380−7.26E − 07760−6.98E − 0710−4.02%251.29%1301.40E − 05400−8.05E − 07800−6.69E − 07118.17%26−0.72%1401.22E − 05420−8.47E − 07840−6.06E − 0712−5.88%270.83%1501.04E − 05440−8.98E − 07880−5.64E − 07136.76%28−0.58%1608.81E − 06460−8.97E − 07920−5.65E − 0714−2.93%290.80%1806.04E − 06480−8.94E − 07960−5.17E − 07153.07%30−0.59%2004.02E − 06500−9.00E − 071000−4.75E − 07Figure 19
The relative error ofIgan when 1000 ≤n≤ 10000, step 300.Whenn<20, the relative error of Iga is still greater than 2%. Although it is much better than the error of Ig, it is not as good as expected when n<40. If we take the round approximation by(40)Iga′n=πexp2/3πn+a1/x+c1+b1122n3+12,n≥4,the relative error will be obviously smaller with a few exceptions, as shown in Table 7.Table 7
The relative error ofIgan+1/2 to hn when n≤30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—725%138.33%192.86%251.31%2—8014−2.94%20−1.46%26−0.63%3—912.5%152.44%211.82%270.87%4010016−1.82%22−0.95%28−0.56%550%117.14%173.03%231.58%290.83%6012−4.76%18−2.27%24−1.25%30−0.58%Later, we will find out that it is obviously greater than the relative error ofIg1 and Ig2 obtained in the next subsection by modifying the denominator part; when 4000<n<10000, the relative error of Iga is about 1000 times of that of Ig2.When 30<n≤ 1000, the relative error of Igan+1/2 is close to that of Igan.


## 3.3. Method B: Modifying the Denominator

Sincehn∼π/122n3exp2/3πn, we consider estimating hn byπ/122C3nexp2/3πn (i.e., fit π2exp2π2/3n/288h2n by a function C3n), where C3x is a cubic function or a function like(41)ax3+bx2.5+cx2+dx1.5+ex+fx0.5+g.However, the results are worse as the relative errors are obviously much greater than the relative error ofIgn when n<350.Then, we consider estimatinghn by π/122C4nexp2/3πn or fit πexpπ2/3n/122hn by a function as follows:(42)C4n=a4n1.5+b4n+c4n0.5+d4.The result is very good. The graph of the datan,πexpπ2/3n/122hn and the fitting curve C4n is shown in Figure 3. Here, the fitting curve is displayed by a thick continuous curve, which lies in the middle of the area the circles occupied. Since the circles are too crowded, the circles themselves look like a very thick curve.The values of the coefficients in the expression ofC4n are as follows:(43)a4=1.000010809,b4=1.862505234,c4=1.169930087,d4=−0.7005460222.The value ofa4 is very close to 1, which means that this fitting function coincides with the Ingham–Meinardus asymptotic formula very well.So, we have an estimation formula as follows:(44)hn∼Ig1n=π122C4nexp23πn.We may call it the Ingham–Meinardus revised estimation formula (1). The graph oflnIg1n is shown in Figure 2, together with the data points of n,lnhn. This revised estimation formula is much more accurate than the asymptotic formula. The relative error is less than 1×10−6 when n > 2000 (as shown in Figure 20 on page 14) and less than 3 when n≥30 (as shown in Table 8 on page 14). The relative error of the round approximation Ig1′n=Ig1n+1/2 is shown in Table 9 on page 13.Figure 20
The relative error ofIg1n when 1000 ≤n≤ 10000, step 300.Table 8
The relative error ofIg1n to hn when n≤1000.
nRel-Err (%)nRel-Err (%)nRel-ErrnRel-ErrnRel-Err1—16−1.6340−2.21E − 052207.23E − 055201.04E − 062−7.23173.82505.35E − 042405.74E − 055403.01E − 07329.9718−1.87606.16E − 042604.59E − 05560−3.53E − 074−8.44193.18706.15E − 042803.68E − 05580−8.95E − 07527.9420−1.21805.35E − 043002.97E − 05600−1.37E − 066−11.61212.06904.56E − 043202.40E − 05640−2.10E − 06720.7622−0.611003.89E − 043401.93E − 05680−2.64E − 068−6.74231.891103.30E − 043601.55E − 05720−2.97E − 0699.2124−0.851202.80E − 043801.24E − 05760−3.20E − 0610−3.44251.631302.40E − 044009.79E − 06800−3.36E − 06118.8526−0.401402.06E − 044207.63E − 06840−3.43E − 0612−5.28271.141501.78E − 044405.82E − 06880−3.51E − 06137.4228−0.291601.55E − 044604.32E − 06920−3.49E − 0614−2.35291.081801.18E − 044803.04E − 06960−3.43E − 06153.6630−0.322009.20E − 055001.97E − 061000−3.37E − 06Table 9
The relative error ofIg1n+1/2 to hn when n⩽30.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—725%138.33%192.86%251.57%208014−2.94%20−1.46%26−0.42%30912.5%154.88%211.82%271.22%4010016−1.82%22−0.48%28−0.28%550%117.14%174.55%231.98%291.06%6012−4.76%18−2.27%24−0.94%30−0.29%However, equation (44) is not so satisfying when n<30, especially when n<15 as the relative error is not negligible for some value of n.As we already know thathn∼π/122n3exp2/3πn, or n3/2∼π/122hnexp2/3πn, which means that when fitting πexpπ2/3n/122hn by a function C4n shown in equation (42), the coefficient a4 should be exactly 1; hence, we should fit πexpπ2/3n/122hn by a function C4′n=n3/2+b5n+c5n1/2+d5 or fit πexpπ2/3n/122hn−n3/2 by a function(45)C5n=b5n+c5n1/2+d5.The graph of the datan,πexpπ2/3n/122hn−n3/2 is shown in Figure 21 on page 14 (together with the figure of the fitting function C5n generated by the least square method).Figure 21
The graph of the datan,πexpπ2/3n/122hn−n3/2 and the fitting curve C5n.The values of the coefficients in equation (45) are as follows:(46)b5=1.864260743,c5=1.084436400,d5=0.4754177757.So, we have another estimation formula forhn as follows:(47)hn∼Ig2n=πexp2/3πn122n3/2+C5n.We may call it the Ingham–Meinardus revised estimation formula (2). The graph oflnIg2n is nearly the same as that of lnIg1n shown in Figure 2. The second revised estimation formula is much more accurate than the first one. The relative error is less than 2×10−9 when n > 3000 (as shown in Figure 22 on page 14), about 1/500 of the relative error of Ig1n. When n < 10, the relative error is also distinctly less than that of Ig1n (as shown in Table 10 on page 15). The relative error of the round approximation Ig2′n=Ig2n+1/2 is shown in Table 11 (on page 15).Figure 22
The relative error ofIg2n when 1000 ≤n≤ 10000, step 300.Table 10
The relative error ofIg2n to hn when n≤1000.
NRel-Err (%)nRel-Err (%)nRel-ErrnRel-ErrnRel-Err116−2.4940−0.21%2202.15E − 065201.78E − 072−18.69172.9950−9.06E − 042401.76E − 065401.83E − 07319.7518−2.5960−4.33E − 042601.46E − 065601.59E − 074−13.56192.4870−1.80E − 042801.20E − 065801.51E − 07522.5120−1.8380−8.70E − 053009.94E − 076001.32E − 076−14.58211.4690−4.13E − 053208.51E − 076401.03E − 07717.4322−1.15100−1.68E − 053407.18E − 076805.80E − 088−8.89231.37110−5.98E − 063606.14E − 077206.80E − 0897.0624−1.32120−7.10E − 073805.23E − 077606.70E − 0810−5.09251.181302.07E − 064004.61E − 078005.10E − 08117.2326−0.821403.17E − 064203.90E − 078405.40E − 0812−6.53270.741503.54E − 064403.34E − 07880−4.30E − 09136.1628−0.661603.59E − 064602.96E − 079207.00E − 0914−3.38290.721803.16E − 064802.50E − 079602.80E − 08152.6730−0.662002.64E − 065002.22E − 0710003.30E − 08Table 11
The relative error ofIg2n+1/2 to hn when n≤30.
nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)nRel-Err (%)11650.304032.1022013.105208.392146.241756.825028.6024012.505408.233202.891846.696025.9026012.005608.08495.591952.757023.9028011.505807.935156.432044.948022.3030011.106007.79668.622148.489020.9032010.806407.547121.382243.4710019.8034010.406807.31865.432346.0011018.8036010.107207.10988.382441.0912018.003809.867606.911062.582543.6813017.204009.608006.731179.472639.9314016.604209.368406.561253.292741.2715016.004409.148806.411370.982838.5016015.404608.939206.271453.122939.7018014.504808.749606.131560.353037.0020013.705008.5610006.01It should be mentioned that in Figure21 on page 14, the graph of the data points lies in a line, so we might be willing to fit this line by a first-order equation. The result is(48)C5′n=1.873818457×n+27.08318017.If we use this fitting function instead ofC5n generated before, the relative error to fit hn will be about 10000 times more, that is, about 20 times more than that of Ig1n. So, we did not use linear function to fit the data n,πexpπ2/3n/122hn−n3/2 before.


## 3.4. Method C: FittingIgn−hn

We wonder whether we can fitIgn−hn by a function rn and then estimate hn by Ign−rn which may be believed more accurate than Ig2n at the price of being more complicated.By the same tricks used at the beginning of this subsection, we will have(49)Ign−Ign−t∼tπ2243n2exp23πn.t≪n,So, we may fitIgn−hn by π2/243C6nexp2/3πn where C6n is a quadratic function or a function like(50)an2+bn1.5+cn+dn0.5+e.That means, we can fitπ2exp2/3πn/243Ign−hn by a function C6n. However, the result is useless. Although C6n will fit the data π2exp2/3πn/243Ign−hn very well, the relative error of Ign−π2/243C6nexp2/3πn to hn is much greater than that of Ig1n or Ig2n, and the relative error differs very little with that of Ign when n is small. Besides, the formula Ign−π2/243C6nexp2/3πn is much more complicated than Ig1n and Ig2n.If we fitIgn−hn by π2/243n2+E6nexp2/3πn with the trick described in subsection 3.3, where E6n is a function in the form as follows:(51)bn1.5+cn+dn0.5+e,as we already know the coefficient of n2 should be 1 in theory. The result will be a little better, but useless too. The accuracy is not as good as that of Ig0n.Then, we consider fittingπ2exp2/3πn/243n2Ign−hn by a function C7n. If C7n is in the form a/n+b or a/n+b/n2+c, the result is useless, too. If C7n is in the form a/n0.5+b, it will be barely satisfactory. If C7n is in the form a/n0.5+b/n+c/n1.5+d/n2+e or a/n0.5+b/n+c/n1.5+e, the result will be much better than the previous forms, but the accuracy (when estimating hn) is not as good as that of Ig1n and Ig2n.The result ofC7n is(52)C7an=0.8782296151n0.5+0.2567016063n−3.580442785n1.5+21.28305831n2+0.6879945549,or(53)C7bn=0.8861039149n0.5−0.05719053203n+0.9843423289n1.5+0.6879343652.The relative errors of(54)F7an=Ign−π2exp2/3πn243n2C7an,and(55)F7bn=Ign−π2exp2/3πn243n2C7bn.to hn when 1000 ≤n≤ 10000 are shown in Figures 23 and 24 (page 16), respectively. In this interval (1000, 10000), F7an is obviously more accurate than F7bn. When n≤ 1000, the relative error of F7an+1/2 and F7bn+1/2 is shown in Table 12 (page 15) and Table 13 (page 15). In this case, F7bn is better than F7an. However, neither of them is as good as Ig1n or Ig2n although they are more complicated than Ig1n and Ig2n.Figure 23
The relative error ofF7an when 1000 ≤n≤ 10000, step 300.Figure 24
The relative error ofF7bn when 1000 ≤n≤ 10000, step 300.Table 12
The relative error ofF7an+1/2 to hn when n≤1000.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—16040−8.13E − 04220−1.74E − 065203.10E − 072100%174.55%50−3.26E − 04240−1.51E − 065403.29E − 073100%18−1.14%60−1.49E − 04260−1.25E − 065603.42E − 07450%193.81%70−2.81E − 05280−9.91E − 075803.51E − 07550%20−0.73%80−4.62E − 06300−7.56E − 076003.56E − 0760212.42%903.46E − 06320−5.49E − 076403.57E − 07750%22−0.48%1006.70E − 06340−3.72E − 076803.49E − 0780231.98%1105.27E − 06360−2.23E − 077203.36E − 07912.5%24−0.63%1203.37E − 06380−9.87E − 087603.19E − 07100251.83%1301.93E − 064003.54E − 098003.00E − 071114.29%26−0.21%1405.77E − 074208.70E − 088402.80E − 07120271.22%150−4.01E − 074401.55E − 078802.59E − 07138.33%28−0.28%160−1.04E − 064602.09E − 079202.39E − 07140291.06%180−1.72E − 064802.51E − 079602.19E − 07154.88%30−0.29%200−1.86E − 065002.85E − 0710001.99E − 07Table 13
The relative error ofF7bn+1/2 to hn when n⩽1000.
nRel-ErrnRel-Err (%)nRel-ErrnRel-ErrnRel-Err1—16−1.8240−0.16%220−1.21E − 06520−1.11E − 0620173.0350−6.19E − 04240−1.77E − 06540−1.01E − 063018−2.2760−2.60E − 04260−2.08E − 06560−9.24E − 0740192.8670−7.50E − 05280−2.21E − 06580−8.41E − 07550%20−1.4680−1.85E − 05300−2.24E − 06600−7.64E − 0760211.82903.62E − 06320−2.20E − 06640−6.25E − 07725%22−0.951001.30E − 05340−2.12E − 06680−5.05E − 0780231.581101.37E − 05360−2.02E − 06720−4.00E − 07912.5%24−1.251201.21E − 05380−1.90E − 06760−3.10E − 07100251.311301.01E − 05400−1.78E − 06800−2.31E − 07117.14%26−0.631407.74E − 06420−1.66E − 06840−1.63E − 0712−4.76%270.871505.68E − 06440−1.54E − 06880−1.04E − 07138.33%28−0.561603.97E − 06460−1.42E − 06920−5.24E − 0814−2.94%290.831801.40E − 06480−1.31E − 06960−7.83E − 09152.44%30−0.58200−2.22E − 07500−1.21E − 0610003.08E − 08


## 3.5. Estimatehn When n≤100

All the estimation functions forhn found now are with very good accuracy when n is greater than 1000, but they are not so accurate when n<50, especially when n<25. Although Ig1′n and Ig2′n are better than others, the relative error is still greater than 1 when n<40.Whenn<40, it is too difficult to fit πexpπ2/3n/122hn−n3/2 by a simple smooth function with high accuracy, as shown in Figure 25 (on page 16). The points n,πexpπ2/3n/122hn−n3/2 (n = 3, 4, ⋯, 100) are not so complicated (as shown in Figure 25). It seems that we can fit them by a simple piecewise function with 2 pieces, as the even points (where n is even) lie roughly on a smooth curve, so do the odd points. If we try to fit them, respectively, we will have the fitting function as follows:(56)C8x=1.942141112×x−0.4796781366×x+8.291226268,n=3,5,7,…,99;1.803056782×x+2.356539877×x−6.043824511,n=4,6,8,…,100.Figure 25
The graph of the datan,C8n.Hence, we can calculatehn (3≤n≤100) by(57)hn∼Ig0n=πexp2/3πn122n3/2+C8n,3≤n≤100.Considering thathn is an integer, we can take the round approximation of equation (57):(58)Ig0′n=πexp2/3πn122n3/2+C8n+12,3≤n≤100.Here,n begins from 3, not 1 or 2, because Ig0′1−h1/h1 is meaningless since h1=0, and Ig0′2 differs from h2 a lot. Besides, the value of h1 and h2 is clear by definition, so there is no need to use a complicated formula to estimate them.The relative error ofIg0n (or Ig0′n) to hn is shown in Table 14 (or Table 15) on page 17. Compared them with Table 11 on page 15, we will find that when n≥80, Ig2′n is more accurate than Ig0′n; when n<80, Ig0′n is better.Table 14
The relative error ofIg0n to hn when n≤100.
nRel-Err (%)nRel-ErrnRel-ErrnRel-ErrnRel-Err1—21−0.16%41−5.72E − 0461−1.04E − 04817.19E − 052—220.12%423.75E − 04621.43E − 0482−5.18E − 053−14.85230.05%43−4.68E − 0463−1.04E − 04838.12E − 05412.7224−0.28%444.91E − 04641.18E − 0484−6.50E − 0551.99258.37E − 0445−6.29E − 0465−5.63E − 05858.69E − 056−1.83264.75E − 04465.45E − 04667.10E − 0586−6.90E − 0574.7627−0.17%47−4.45E − 0467−2.40E − 05878.65E − 058−0.64286.69E − 04483.17E − 04685.74E − 0588−7.40E − 059−0.9229−3.78E − 0449−3.42E − 0469−1.67E − 05899.01E − 05100.6930−4.43E − 04503.79E − 04704.05E − 0590−8.07E − 05111.4431−1.98E − 0451−3.98E − 04711.30E − 05919.09E − 0512−2.46324.21E − 04523.55E − 04726.81E − 0692−8.24E − 05131.8633−0.12%53−2.81E − 04733.41E − 05938.80E − 0514−0.24349.27E − 04542.47E − 0474−1.74E − 0694−8.37E − 0515−0.5435−7.86E − 0455−2.21E − 04753.84E − 05958.70E − 0516−0.04361.82E − 04562.44E − 0476−1.56E − 0596−8.65E − 05170.4637−4.80E − 0457−2.25E − 04775.70E − 05978.44E − 0518−0.67386.53E − 04582.29E − 0478−3.58E − 0598−8.56E − 05190.4739−9.11E − 0459−1.55E − 04796.90E − 05997.92E − 0520−0.28406.34E − 04601.44E − 0480−4.41E − 05100−8.48E − 05Table 15
The relative error ofIg0′n to hn when n≤100.
nRel-ErrnRel-ErrnRel-ErrnRel-ErrnRel-Err1—9017025033−0.11%2—10018−1.14%260349.23E-043011019027−0.17%35−7.77E-044012−4.76%200280363.23E-045013021029037−5.46E-0460140220300386.85E-047015023031039−9.67E-048016024−0.31%320.07%406.50E-04


## 4. Conclusion

We have presented a recursion formula and several practical estimation formulae with high accuracy to calculate the numberhn of conjugate classes of derangements of order n or the number of isotopy classes of 2×n Latin rectangles.If we want to obtain the accurate value ofhn, we can use the recursion formula (13) and write a program based on it, while sometimes we need to know the estimation value in a program for technique reason especially when we use a general programming language.If we want to obtain the approximation value ofhn with high accuracy, we can use formulae (29), (31), (36), and so on.When2≤n≤80, we can use Ig0′n (equation (58)), with a relative error less than 0.11% (while 32≤n≤80) or mainly 0 with very few exceptions (while 2≤n≤31); when n>80, we can use Ig2′n (equation (47)).Whenn≥100, formulae Iga′n (equation (40)), Ig1′n (equation (44), F7an (equation (54)), and F7bn (equation (55)) are also very accurate although they are not as good as equation (47).With the asymptotic formula (26) described in [30], we can obtain some estimation formulae with high accuracy for some other types of restricted partition numbers by the methods mentioned in this paper or in [1]. The updated version of this paper is shown on https://arxiv.org/abs/1612.08186.



---

*Source: 6023081-2021-10-28.xml*