On the Number of Conjugate Classes of Derangements

The number of conjugate classes of derangements of order $n$ is the same as the number $h(n)$ of the restricted partitions with every portion greater than $1$. It is also equal to the number of isotopy classes of $2\times n$ Latin rectangles. Sometimes the exact value is necessary, while sometimes we need the approximation value. In this paper, a recursion formula of $h(n)$ will be obtained, also will some elementary approximation formulae with high accuracy for $h(n)$ be presented. Although we may obtain the value of $h(n)$ 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 an 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.


Introduction
Below n 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 order n 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. ) Let S n 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 S n is also called a permutation (of order n). If σ ∈ S n , σ(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 D n (or !n in some literature). It is mentioned in nearly every combinatorics textbook that, Here x is the floor function, which stands for the maximum integer that is less than or equal to the real x.
For x, y ∈ S n , if ∃z ∈ S n , 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 size k × 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 }.) A 2 × 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 in S n share the same cycle structure. Here we define the cycle structure of a derangement as the sequence in non-decreasing 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 s 1 +s 2 +· · ·+s q = n, (2 s 1 s 2 · · · s q ), (1) where s 1 , s 2 , · · · , s q are unknowns. So the number of conjugate classes of derangements of order n is h(n). Since h(n) is the number of a type of restricted partitions, it is tightly connected with the partition number p(n).
Following the notation of [2], denote by P q (n) the number of integer solutions of equation s 1 +s 2 +· · ·+s q = n, (1 s 1 s 2 · · · s q ) (2) for a fixed q, where 1 q n, and by p(n) the number of all the (unrestricted) partitions of n. It is clear that 1 There is a brief introduction of the important results on the partition number (or partition function) p(n) and P q (n) in reference [2], such as the recursion formula of p(n) and P q (n). More information about the partition number p(n) 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 p(n), which are revised from the Hardy-Ramanujan's asymptotic formula.
In [22], we can find many cases of Restricted Partitions (some of them are introduced in [23], [24] or [25]). One class are concerned on the restriction of the sizes of portions, such as portions restricted to Fibonacci numbers, powers (of 2 or 3), unit, primes, non-primes, composites or non-composites; another class are 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 are about the restrictions for both, for example, the cases that the number of parts is restricted while the parts restricted to powers or primes. But the author has not found too much information on the number h(n), 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] and [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.
Section 2 will deduce the recursion formula for h(n) and will show the relation of h(n) and p(n). Subsection 3.1 will deduce the asymptotic formula of h(n) from Hardy and Ramanujan's Asymptotic formula of p(n) (mentioned in [1]). This new asymptotic formula I g (n) coincides with Ingham's result (refer [29] and [30]). By bringing in two parameters C 1 (n) and C 2 (n) in the new asymptotic formula I g (n), we have reached several estimation formulae for h(n) with high accuracy in subsection 3.3, using the same ideas described in [1]. By fitting h(n) − I g (n), we have another two estimation formulae for h(n) in subsection 3.4. When n < 100, we have a more accurate estimation formula for h(n) in subsection 3.5. The relative errors of these estimation formulae will be presented to shown the accuracy.

Some Formulae for h(n)
In this section, we will derive a recursive formula from the method mentioned in reference [2] (page 53~55).
It is mentioned in [2] (page 52) that in 1942 Auluck gave an estimation of P q (n) by P q (n) ≈ 1 q! n − 1 q − 1 when n is large enough.
By the same method shown in reference [2] (page 53, 57), we can obtain the generation function of h(n): and a formula where h(0) = 1, h(1) = 0, and C is a contour around the original point.
It is difficult to get a simple formula to count the solutions of Equation (1) in general. But for a fixed integer q, the number H q (n) of solutions is 0 (when q > n 2 ) or . It is not difficult to find out that And there is a recursion for P q (n) in reference [2] (page 51) where t = min{q, n − q}, so there is no difficulty to obtain the values of P q (n) and h(n) when n is small.
For the value of p(n) there is a recursion, where and assume that p(0) = 1. (Refer [2], page 55). Here we assume that p(x) = 0 when x < 0.
We can obtain the same recursion for h(n), where k 1 and k 2 are determined by Equation (9) and assume that h(0) = 1, h(k) = 0 when k < 0.
The proof of Equation (10) is easy to understand. Since Compare Equation (11) and Equation (12), we have by assumption h(k) = p(k) = 0 when k < 0. Hence, h(n) = p(n) − p(n − 1), (n = 0, 1, 2, · · · ). (13) By Equation (13), We can easily obtain the solutions of Equation (1) by hand when n < 13. By Equation (10), we can obtain the number h(n) of solutions of Equation (1) by writing a small program in some Computer Algebra System (CAS) softwares such as "maple", "maxima", "axiom" or some other softwares likewise (be aware of that 0 is not a valid index value in some software such like maple).
The value of h(n) when n < 250 are shown on Table  1 (on page 5) and Table 2 (on page 5). Some value of H q (n) are shown on Table 3 (on page 5).
Obviously, h(n) < p(n) holds by definition (when n > 1). As p(n) grows much more slowly than exponential functions, i.e., for any r > 1, p(n) < r n will hold when n is large enough, which means we can not estimate p(n) and h(n) by an exponential function. As p(n) grows faster than any power of n, which means we can not estimate p(n) by a polynomial function. (refer [2], page 53) So, h(n) can not be estimated by a polynomial function, either.

The Estimation of h(n)
The recursion formula Equation (10) for h(n) 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 figure of the data n, ln(h(n)) ( n = 60 + 20k, k = 1, 2, · · · , 397) are shown on 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 on Figure 19 on page 14). But later in Figure 18, the circles distribute uniformly on a curve, it will be difficult to distinguish the circles and the curve passes through the centers of the them.
The author has not found a practical estimation formula with high accuracy of the number h(n) before.
Actually, it is very difficult to find directly a simple function to fit the data on Figure 1 with high accuracy. The main reason is that the fitting functions        obtained by the methods used frequently could not reach satisfying accuracy.
Since we have several accurate estimation formula of p(n) (refer [1]), such as where a 2 = 0.4432884566, b 2 = 0.1325096085, c 2 = 0.274078 and By Equation (13), we can obtain h(n) by and the error of this formula will not exceed twice of the error of R h2 (n) or R h0 (n). Of course, this formula is not as simple as we want, but the accuracy is very good.
In coincidence, half a year after the main results were obtained in this paper, the author found an asymptotic formula in [30]. When a = 1, b = 2, we will have which coincides with the asymptotic formula obtained here.
The formula (17) will be called the Ingham-Meinardus asymptotic formula in this paper, since Daniel mentioned in [30] that this general asymptotic formula (18) 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, √ n will be denoted by I g (n) for short.
It is not satisfying to estimate h(n) by I g (n) when n is small. The relative error of I g (n) to h(n) is greater than 6% as shown on Table 4 (on page 6). The round approximation will not change the accuracy distinctly, as shown on Table 5 (on page 6). So it is necessary to modify the asymptotic formula for better accuracy.

Method A: Modifying the exponent
In this subsection we consider fitting h(n) by I ga = π 12 √ 2n 3 exp 2 3 π n + C 1 (n) , or fitting − n ( n = 60 + 20k, k = 1, 2, · · · , 397) by a function The reason that we fit C 1 (x) by the function in the form displayed in (20) is the same as that described in section 3 of [1] (although the data differs distinctly). Many other types of functions have been tested, but they can not fit these data very well.
But here it is not valid to obtain the constants in f 1 (n) by iteration method described in reference [1].
6 Figure   2: The graph of the data − n ( n = 60 + 20k, k = 1, 2, · · · , 397) by a function in the form That means we have assumed that e 1 = 1/2, temporarily. We will explain the reason in subsection 3.2.3.
The average error of f 1 (x) is where K 1 = 397, n ranges from 80 to 8000, by step 20. Here only a 1 , b 1 , and c 1 are unknown, so we can consider E 1 as a function of the variables (a 1 , b 1 , c 1 ).
We want to find a triple (a 1 , b 1 , c 1 ) such that E 1 reaches its minimum, or to make E 1 as small as possible.
Since a lot of functions have several local minimum points, it is necessary to find out whether E 1 = E 1 (a 1 , b 1 , c 1 ) has more than one local minimum before we start to calculate the minimum point by numeral method. But E 1 = E 1 (a 1 , b 1 , c 1 ) is too complicate, it is very difficult to know all the critical points in the range we are considering.

Preparation work
For a given pair (a 1 , c 1 ), by the property of the arithmetic mean, 2 it is clear that E 1 reaches its minimum when where Let be the average error of the the fitting function Here only a 1 and b 1 are unknowns. G 1 could be considered as a function of (a 1 , c 1 ). In order to find the minimum point of G 1 , we can draw the figure of the function G 1 = G 1 (a 1 , c 1 ), (In a cube coordinate system with axis a 1 , c 1 and G 1 ) as shown on Figure  3, Figure 4 and Figure 7. Figure 5 and Figure 6 are the projection of the graph of (a 1 , c 1 , G 1 ) (when −100 a 1 100, −50 c 1 100, which is a part of a surface) on the a 1 − G 1 plane (spanned by the axis a 1 and G 1 ) and c 1 − G plane, respectively.
From these figures, we can find out that the influence of c 1 to G 1 is much less than that of a 1 . In Figure  9, we find that when G 1 reaches its minimum, a 1 is between 0.50 and 0.53, but there is not a definite range for c 1 .
It is possible that for different range of c 1 , the range of a 1 when G 1 reaches its minimum will be different. But considering that a1 √ n+c1 + b 1 is a real, c 1 should be greater than −1 in theory. (For the fitting data used here, c 1 should be greater than −80.) From Figure 5, we can see that G 1 touches its bottom when −15 a 1 15. Although we can not see clearly the exact value of of a 1 in the minimum points, we can draw another figure of (a 1 , c 1 , G 1 ) when −15 a 1 15 and −1 c 1 100 to observe more details (the figure is not presented here), then we will find that 2 For some given data x i (i = 1, 2, · · · , k; x i ∈ R), the the more detailed range of a 1 for the minimum points is [−3, 3] in the new figure (not presented), next, draw the figures of (a 1 , c 1 , G 1 ) when −3 a 1 3, 0 a 1 1, 0.2 a 1 0.8 or 0.45 a 1 0.6, respectively, while −1 c 1 100, we will find the range of a 1 of the minimum points of G 1 is about [0.50, 0.53]. The projections of the figure 7 of (a 1 , c 1 , G 1 ) when 0.45 a 1 0.6 and −1 c 1 100 is shown on Figure 8 and Figure 9.
In Figure 6, for the curves on the bottom, G 1 decrease with c 1 at first then increases with c 1 , but it is difficult to find the critical point in the figure, since different curves have different critical points.
Although we can not find an satisfying value of a 1 or c 1 from the figures above to construct a fitting function f 1 , these pictures show that the figure of G 1 = (a 1 , c 1 , G 1 ) has only one bottom in the domain we are considering, unlike the figure of another function shown on Figure 10, so the existence of the minimum point is almost confirmed, therefore we are confident to find the value of a 1 or c 1 in the minimum point by numerical method. This guarantees the validity of the numerical calculation by loop in next step.

Find c 1
On the other hand, by the least square method, to fit the data (x k , y k ) (k = 1, 2, · · · , K 1 ) by a linear function y = a × x + b, the result is that where is the square of the average value of x n . So, by the least square method, a and b are uniquely determined by the given data (x k , y k ) (k = 1, 2, · · · , K 1 ).
For every given value of c 1 (greater than −80), we can fit (n, C 1 (n)) (n = 60 + 20k, k = 1, 2, · · · , 397) by a function f 1 (x) . = a1 √ x+c1 + b 1 by the least square method, if we consider 1 √ 60+20k+c1 and C 1 (60 + 20k) as x k and y k , respectively. Then So a 1 and b 1 could both be considered as functions of c 1 , denoted by a 1 = a 1 (c 1 ), b 1 = b 1 (c 1 ), since they are uniquely determined by c 1 with the given data. Then It will cost some time to plot the figure of the function G 2 = G 2 (c 1 ) in a CAS software.
If we plot the figure of the function G 2 = G 2 (c 1 ) on the coordinates (as shown on Figure 11, Figure  12 and Figure 13), we will find that G 2 reaches its minimum when c 1 ≈ −3.2594807. On Figure 13, we find that the curve of G 2 = G 2 (c 1 ) 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 13 will be more smooth, at the cost of much more time. By writing a small program (since the default function to find the minimum provide by the software Maple 18 are unable to deal with such a complicated function G 2 = G 2 (c 1 ) involving so much data), we can obtain a more accurate value of the critical point When the value of c 1 is obtained, we can find the value of a 1 and b 1 by the least square method without difficulty, i.e., But here c 1 is less than −1, so the estimation formula for h(n) constructed from these coefficients is invalid when n < 4.

Confirm e 1
In [1] we fit C 1 (n) = 3 2 · (ln(4n √ 3p(n)))  the iteration method does not work well, so we fit √ n+c1 + b 1 directly, which means that we have assumed that e 1 = 1 2 . Here we may doubt that whether e 1 = 0.5 is the best option for us?
Here we use the same idea described in subsection 3.2.1.
For every pair (e 1 , c 1 ), we can obtain corresponding a 1 and b 1 by the least square method, just like (24) and (25), except that here So the square of the average error     After that, by written another program, we can obtain the approximate value of (e 1 , c 1 ) where G 3 touches the bottom, i.e. e 1 ≈ 0.494, c 1 ≈ −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 can not 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 e 1 = 0.50 while it differs very little with 0.494. Another reason is that we prefer simple exponent, as the time spend on computing a square root is much less than that to compute a power with exponent 0.494 in general. Here the value of c 1 ≈ −4.85 is obviously different from the value obtained at the end of subsection (3.2.2), because of the little difference on e 1 . Therefore, it will be fine to use the result in subsection (3.2.2).

The Result
By Then we could fit h(n) by The figure of the function ln (I ga (x)) together with the figure of the data (n, ln (h(n))) ( n = 60 + 20k, k = 1, 2, · · · , 397) are shown on Figure 16. It seems that I ga (n) fits h(n) very well.
The relative error of I ga is shown on Table 6 (when n 1000 ) and Figure 17 (when 1000 < n 10000).
When n < 20, the relative error of I ga is still greater than 2%. Although it is much better than the error of I g , it is not as good as expect when n < 40. If we take the round approximation by the relative error will be obviously smaller with a few exceptions, as shown on Table 7.
Later we will find out that it is obviously greater than the relative error of I g1 and I g2 obtained in the next subsection by modifying the denominator part; when 4000 < n < 10000, the relative error of I ga is about 1000 times of that of I g2 .     Table 7: The relative error of I ga (n) + 1 2 to h(n) when n 30. When 30 < n 1000, the relative error of I ga (n) + 1 2 is closed to that of I ga (n).

Method B: Modifying the Denominator
n , (i.e., fit by a function C 3 (n)), where C 3 (x) is a cubic function or a function like ax 3 + bx 2.5 + cx 2 + dx 1.5 + ex + f x 0.5 + g.
But the results are worse, as the relative errors are obviously much greater than the relative error of I g (n) when n < 350. by a function The result is very good. The figure of the data n, and the fitting curve C 4 (n) are shown on Figure 18 on page 12. 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.     Table 9: The relative error of I g1 (n) + 1 2 to h(n) when n 30. The value of a 4 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 We may call it the Ingham-Meinardus revised estimation formula 1. The graph of ln (I g1 (n)) is shown on Figure 19 on page 14, together with the data points of (n, ln h(n)). 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 on Figure 20 on page 14), and less than 3 when n 30 (as shown on Table 8 on page 13). The relative error of the round approximation I g1 (n) = I g1 (n) + 1 2 is shown on Table 9 on page 13.
But Equation (29) 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 that h(n) ∼ by a function C 4 (n) shown in Equation (28), the coefficient a 4 should be exactly 1, hence we should fit by a function C 4 (n) = n 3/2 + b 5 n + c 5 n 1/2 + d 5 , or fit − n 3/2 by a function The figure of the data n, shown on Figure 21 on page 14 (together with the figure of the fitting function C 5 (n) generated by the least square method).
The values of the coefficients in Equation (30)      We may call it the Ingham-Meinardus revised estimation formula 2. The graph of ln (I g2 (n)) is nearly the same as that of ln (I g1 (n)) shown on Figure 19 on page 14. 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 on Figure 22 on page 14), about 1 500 of the relative error of I g1 (n). When n < 10, the relative error is also distinctly less than that of I g1 (n) (as shown on Table  10 on page 15). The relative error of the round approximation I g2 (n) = I g2 (n) + 1 2 is shown on Table  11 (on page 15).
It should be mentioned that in Figure 21 on page 14, the graph of the data points lie in a line, so we might be willing to fit this line by a first order equation. The result is C 5 (n) = 1.873818457 × n + 27.08318017.
If we use this fitting function instead of C 5 (n) generated above, the relative error to fit h(n) will be about 10000 times more, that is about 20 times more than that of I g1 (n). So we did not use linear function to fit the data n, − n 3/2 before.

Method C: Fittong I g (n) − h(n)
We wonder whether we can fit I g (n) − h(n) by a function r(n), then estimate h(n) by I g (n) − r(n) which may be believed more accurate than I g2 (n) at the price of being more complicated.
By the same tricks used at the beginning of this subsection, we will have √ n where C 6 (n) is a quadratic function or a function like an 2 + bn 1.5 + cn + dn 0.5 + e.
That means, we can fit 3(Ig(n)−h(n)) by a function C 6 (n). But the result is useless. Although C 6 (n) will fit the data  h(n) is much greater than that of I g1 (n) or I g2 (n), and the relative error differs very little with that of I g (n) when n is small. Besides, the formula I g (n) − π 2 24 √ 3C6(n) exp 2 3 π √ n are much more complicated than I g1 (n) and I g2 (n).
If we use the trick mentioned in 3.3, to fit I g (n) − h(n) by n where E 6 (n) is a function like bn 1.5 + cn + dn 0.5 + e, as we already know that the coefficient of n 2 should be 1 in theory. The result will be a little better, but useless too. The accuracy is not as good as that of I g0 (n).
Then we consider fitting 3n 2 (Ig(n)−h(n)) by a function C 7 (n). If C 7 (n) is in the form a n +b or a n + b n 2 +c, the result is useless, too. If C 7 (n) is in the form a n 0.5 +b, it will be barely satisfactory. If C 7 (n) is in the form a n 0.5 + b n + c n 1.5 + d n 2 + e or a n 0.5 + b n + c n 1.5 + e, the result will be much better than the previous forms, but the accuracy (when estimating h(n)) is not as good as that of I g1 (n) and I g2 (n).
The result of C 7 (n) is The relative error of and to h(n) when 1000 n 10000 are shown on Figure 23 and Figure 24 (page 16), respectively. In this interval (1000, 10000), F 7a (n) is obviously more accurate than F 7b (n). When n 1000 the relative error of F 7a (n) + 1 2 and F 7b (n) + 1 2 are shown on Table   12 (page 15) and Table 13 ( page 15). In this case, F 7b (n) is better than F 7a (n). But neither of them is as good as I g1 (n) or I g2 (n), although they are more complicated than I g1 (n) and I g2 (n).

Estimate h(n) When n 100
All the estimation function for h(n) 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 I g1 (n) and I g2 (n) are better than others, the relative error are still greater than 1 when n < 40.      − n 3/2 ( n = 3, 4, · · · , 100) is not so complicated (as shown on 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 below: .291226268, n = 3, 5, 7, · · · , 99; 1.803056782 × x + 2.356539877 × √ x −6.043824511, n = 4, 6, 8 · · · , 100.
Consider that h(n) is an integer, we can take the round approximation of Equation (35), Here n begins from 3, not 1 or 2, because is meaningless since h(1) = 0, and I g0 (2) differs from h(2) a lot. Besides, the value of h(1) and h (2) are clear by definition, so there is no need to use a complicated formula to estimate them.
The relative error of I g0 (n) (or I g0 (n)) to h(n) are shown on Table 14 (or Table 15) on page 17. Compared them with Table 11 on page 15, we will find that when n 80, I g2 (n) is more accurate than I g0 (n); when n < 80, I g0 (n) is better.

Conclusion
We have presented a recursion formula and several practical estimation formulae with high accuracy to calculated the number h(n) 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 of h(n), we can use the recursion formula (10) 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 of h(n) with high accuracy, we can use the formulae (31), (36), (29), etc.
With the asymptotic formula (18) 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].