CoCalc -- Examples.ipynb

CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual users to large groups and classes!

| Download

Project: FPSAC - Algebraico-logarithmic

Path: Examples.ipynb

Views: ²³

Visibility: Unlisted (only visible to those who know the link)

Image: ubuntu2204

Kernel: SageMath 10.2

Theorem 1:

Let $H(x, y)$ be an analytic function near the origin whose power series expansion at $(0, 0)$ has non-negative coefficients. Define $\mathcal{V} = \{(x, y) : H(x, y) = 0\}$ . Assume that there is a single smooth strictly minimal critical point of $\mathcal{V}$ at $(p, q)$ within the domain of analyticity of $H$ where $p$ and $q$ are real and positive. Let $\lambda = \frac{r + O(1)}{s}$ as $r,s \to \infty$ with $r$ and $s$ integers. Define the following quantities: $\begin{align*} \chi_1 &= \frac{H_y(p,q)}{H_x(p,q)} = \frac{p}{\lambda q}, \\ \chi_2 &= \left.\frac{1}{2H_x}\left(\chi_1^2H_{xx} - 2\chi_1H_{xy} + H_{yy}\right)\right|_{(x,y)=(p,q)}, \\ M &= -\frac{2\chi_2}{p} - \frac{\chi_1^2}{p^2} - \frac{1}{\lambda q^2}. \end{align*}$ Assume that $H_x(p,q)$ and $M$ are nonzero. Fix $\alpha \in \mathbb{R}$ where $\alpha \not \in \mathbb{Z}_{\leq 0}$ and $\beta \in \mathbb{Z}_{\geq 0}$ . Then, the following expression holds as $r,s \to \infty$ : $\begin{align*} [x^ry^s]H(x,y)^{-\alpha}\log^\beta(H(x,y)) \sim (-1)^\beta\frac{(-pH_x(p,q))^{-\alpha}r^{\alpha-1}}{\Gamma(\alpha)\sqrt{-2\pi q^2Mr}}p^{-r}q^{-s}\log^\beta r\left[1+\sum_{j\geq1}\frac{\mathcal{E}_j}{\log^jr}\right], \end{align*}$ where $\mathcal{E}_j = \sum_{k=0}^j\binom{\beta}{j}\binom{j}{k}\log^k\left(\frac{-1}{pH_x(p,q)}\right)\Gamma(\alpha)\left.\frac{d^{j-k}}{dt^{j-k}}\frac{1}{\Gamma(t)}\right|_{t=\alpha}$ .

In the case that $\alpha = 0$ and $\beta \in \mathbb{Z}_{>0}$ , the asymptotic formula of theorem 1 becomes $[x^ry^s]\log^\beta H(x,y) \sim -\frac{\beta r^{-3/2}p^{-r}q^{-s}}{\sqrt{-2\pi q^2M}}\log^{\beta-1}r.$ In particular, for $\beta = 1$ we have $[x^ry^s]\log H(x,y) \sim -\frac{r^{-3/2}p^{-r}q^{-s}}{\sqrt{-2\pi q^2M}}.$

In [1]:

var('x,y,z,u,v,t,ell,m,n,p,r'); # Declare variables

Example 1: (Cyclical Interlaced Permutations)

C(x,y) = \log\left(\frac{1}{1-x-y}\right) = \sum_{m+n \geq 1} \frac{1}{n+m}\binom{n+m}{n}x^ny^m

In [2]:

C = log(1/(1-x-y));
HC = 1-x-y;

HCx = -1;
HCy = -1
HCxx = 0;
HCxy = 0;
HCyy = 0;

Identify smooth minimal critical point in the direction $(1,\ell)$ :

In [3]:

# Check for non smooth critical points
solve([HC==0, HCx==0, HCy==0],x,y)

[]

In [4]:

# Check for minimal smooth critical points
solve([HC==0, ell*x*HCx-y*HCy==0], x, y)

[[x == (1/(ell + 1)), y == ell/(ell + 1)]]

In [5]:

# The minimal smooth critical point is (P_C,Q_C), where
P_C = 1/(1+ell);
Q_C = ell/(1+ell);

Compute asymptotics in direction $(1,\ell)$ :

In [6]:

# From Theorem 2, we have
chi1 = 1;
chi2 = 0;
M = -(1+ell)^3/ell;

In [7]:

# The asymptotics estimates of [x^n y^ln]C(x,y) are then given by 
EstC = 1/sqrt(2*pi*ell*(ell+1))*n^(-3/2)*ell^(-n*ell)*(1+ell)^((1+ell)*n);

Test results:

In the direction $(1,4)$ , the coefficients of $C(x,y)$ are given by $[x^ny^{4n}]C(x,y) = \frac{1}{5n}\binom{5n}{n}$

In [8]:

# Actual value
ActC = 1/(5*n)*binomial(5*n,n);
N(ActC.subs(n=500))

1.61024605591596e538

In [9]:

# Asymptotic estimate of Theorem 2
N(EstC.subs(ell=4,n=500))

1.61052787359801e538

In [10]:

# Comparing values
N(ActC.subs(n=500)/EstC.subs(ell=4,n=500))

0.999825015333995

Example 2: (Necklaces)

N(x,y) = \sum_{k\geq 1}\frac{\varphi(k)}{k}\log\left(\frac{1-x^k}{1-x^k-y^kx^{2k}}\right),

where $\varphi$ is Eulers totient function.

Checking for non-smooth critical points:

For each $k\geq1$ , let $H_k=1-x^k-y^kx^{2k}$ . Then $(H_k)_x = -kx^{k-1}-2ky^kx^{2k-1}$ , and $(H_k)_y = -ky^{k-1}x^{2k}$ . We can write $H_k = 1+\frac{x}{k}(H_k)_x-\frac{y}{k}(H_k)_y,$ which implies the system $[H_k=0, (H_k)_x=0, (H_k)_y=0]$ has no solution. So there are no non-smooth critical points.

Asymptotics are determined by the first term:

We consider the direction $(\ell,1)$ , with $\ell>2$ . For each $k\geq1$ , the $k^{\text{th}}$ term of $N$ contributes smooth critical points given by solutions to the system $[H_k=0,x(H_k)_x-\ell y(H_k)_y=0].$ It is easy to show that solutions to this system satisfy $x^k=\frac{\ell-2}{\ell-1}, y^k=\frac{\ell-1}{(\ell-2)^2}.$ For any such solution, the exponential growth is $x^{-\ell}y^{-1} = \left[(x^k)^{-\ell}(y^k)^{-1}\right]^{1/k} = \left[\frac{(\ell-1)^{\ell-1}}{(\ell-2)^{\ell-2}}\right]^{1/k},$ which is maximized by $k=1$ . So the minimal smooth critical point is $(p,q) = \left(\frac{\ell-2}{\ell-1},\frac{\ell-1}{(\ell-2)^2)}\right)$ . That is, the asymptotics are determined by the first term of $N$ .

To further justify that we only need to consider the first term of $N(x,y)$ , note that $[x^{rk}y^{sk}]N_k(x,y) = \frac{\varphi(k)}{k}[x^ry^s]N_1(x,y),$ where $N_k(x,y) = \frac{\varphi(k)}{k}\log\left(\frac{1-x^k}{1-x^k-y^kx^{2k}}\right).$

In [11]:

# Computing the asymptotics in the direction (ell,1) using the first term of $N$
N1 = log((1-x)/(1-x-y*x^2));

H = 1-x-y*x^2;
Hx = diff(H,x);
Hy = diff(H,y);
Hxx = diff(Hx,x);
Hxy = diff(Hx,y);
Hyy = diff(Hy,y);

# Smooth critical point
P = (ell-2)/(ell-1);
Q = (ell-1)/(ell-2)^2;

# For the values of Theorem 2, we have
Chi_1 = (ell-2)^3/(ell*(ell-1)^2);
Chi_2 = -(2*ell-1)*(ell-2)^5/(ell^3*(ell-1)^3);
M = -(ell-2)^5/(ell^3*(ell-1));

# Using these values, the estimate of Theorem 2 gives
EstN1 = (ell*n)^(-3/2)*P^(-ell*n)*Q^(-n)/sqrt(-2*pi*Q^2*M)

In [12]:

# As a test, we look in the direction (3,1) at the coefficient of x^225*y^75
N = sum(euler_phi(k)/k*log((1-x^k)/(1-x^k-y^k*x^(2*k))) for k in range(1,75)) 
    # Only the first 75 terms of N(x,y) can contribute to [x^225y^75]N(x,y), so we omit the rest.
T = N.taylor((x,0),(y,0),300)

In [13]:

# For this coefficient, we obtain:
ActVal = T.list()[225].list()[75]
EstVal = EstN1.subs(ell=3,n=75)

print("The true value: ", float(ActVal))
print("The approximation: ", float(EstVal))
print("The ratio: ", float(EstVal/ActVal))

The true value:  6.18840464911392e+41
The approximation:  6.1987271801224774e+41
The ratio:  1.0016680439618693

Example 3: (Logarithm of Narayana Numbers)

The Narayana number $a_{n, s}$ counts Dyck paths of length $n$ with $s$ peaks, and $a_{n, s} = [z^nt^s]N(z, t)$ with N defined by

N(z, t) = \frac{1 + z - tz - \sqrt{(1 + z - tz)^2 - 4z}}{2z}.

Here, we look for coefficients of $\log^r N(z, t)$ .

In [14]:

N = (1 + z - t*z - sqrt((1 + z - t*z)^2 - 4*z))/(2*z)
# We check that the singularity at z = 0 is always removable:
N.limit(z=0)

1

In [15]:

# Therefore, the singularities are defined by the zero set of the function underneath the square root.
H = (1 + z - t*z)^2 - 4*z
# We initialize derivatives:
Hz = H.diff(z)
Ht = H.diff(t)
Hzz = Hz.diff(z)
Htt = Ht.diff(t)
Hzt = Hz.diff(t)

In [16]:

# First, we check when N(z, t) = 0, which could create singularities within the logarithm:
solve((1 + z - t*z)^2 == ((1 + z - t*z)^2 - 4*z), [z, t])

[[z == 0, t == r1]]

In [17]:

# This solution can be thrown out because we already saw the limit as z tends to zero is 1.
# Now, we solve for critical points in the [z, t] = [ell, 1] direction:
solve([H, 1*z*Hz - ell*t*Ht], [z, t])

[[z == (ell^2 - 2*ell + 1)/ell^2, t == (1/(ell^2 - 2*ell + 1))], [z == 1, t == 0]]

In [18]:

# By definition, critical points have non-zero coordinates, which leaves us with the single critical point
p = ((ell^2 - 2*ell + 1)/ell^2).factor()
q = (1/(ell^2 - 2*ell + 1)).factor()
print('p = ', p)
print('q = ', q)

p =  (ell - 1)^2/ell^2
q =  (ell - 1)^(-2)

In [21]:

# Following the series expansion in Example 4 of our abstract, we have the following approximation:
Approx = r/(2*pi)*(ln(ell/(ell-1)))^(r-1)*1/n^2*(ell-1)^(-2*n*(ell-1)-1)*ell^(2*ell*n-1)

# We'll check for r = 3, ell = 2, and the coefficient z^80 t^40.
TestF = ln(N)^3
# print(TestF)
TestSeries = TestF.taylor([z, t], 0, 121)
trueCoeff = TestSeries.coefficient(z^(80)).coefficient(t^(40))

approxCoeff = Approx.subs({r:3, ell:2, n:40})
print("The true value: ", float(trueCoeff))
print("The approximation: ", float(approxCoeff))
print("The ratio: ", float(approxCoeff/trueCoeff))

The true value:  9.594835943046605e+43
The approximation:  1.0477113202334314e+44
The ratio:  1.0919533449581384

In [22]:

# We also check for r = 1, ell = 3, and the coefficient z^90 t^30.
TestF = ln(N)
# print(TestF)
TestSeries = TestF.taylor([z, t], 0, 121)
trueCoeff = TestSeries.coefficient(z^(90)).coefficient(t^(30))

approxCoeff = Approx.subs({r:1, ell:3, n:30})
print("The true value: ", float(trueCoeff))
print("The approximation: ", float(approxCoeff))
print("The ratio: ", float(approxCoeff/trueCoeff))

The true value:  1.678177782844727e+45
The approximation:  1.6890898091302997e+45
The ratio:  1.0065023064880978

Example 4: (An Artificial Example)

F(x,y) = \frac{1}{\sqrt{1-x-y}}\log\left(\frac{1}{1-x-y}\right)

In [6]:

F = (1-x-y)^(-1/2)*log(1/(1-x-y));
H = 1-x-y;

Hx = -1;
Hy = -1;
Hxx = 0;
Hxy = 0;
Hyy = 0;

Identify smooth critical point in direction $(1,\ell)$

From example 1, the smooth minimal point in the direction $(1,\ell)$ is $(P_A,Q_A) = (\frac{1}{1+\ell},\frac{\ell}{1+\ell})$ . So we again have

In [0]:

chi1 = 1;
chi2 = 0;
M = -(1+ell)^3/ell;

Compute asymptotics in the direction $(1,\ell)$

In [7]:

# The enumeration formula of Theorem 2, taken up to j=1, gives
Est = (1/(n*pi*sqrt(2*ell)))*((1+ell)^((1+ell)*n)/(ell^(ell*n)))*(log(n)+euler_gamma+2*log(2)+log(1+ell))

Test results:

In [9]:

T_F = F.taylor((x,0),(y,0),100)

In [10]:

# With ell=1
N(Est.subs(n=50,ell=1))/N((T_F.list()[50]).list()[50])

1.00375635630651

In [11]:

# With ell=2
N(Est.subs(n=33,ell=2))/N((T_F.list()[33]).list()[66])

1.00421688271818

In [12]:

# With ell=4
N(Est.subs(n=20,ell=4))/N((T_F.list()[20]).list()[80])

1.00563985529602

Example 5: (Cyclic Partitions of an $n$ -set) [Gao and Richmond 1992]

F(z,x) = \log\frac{1}{1-x(e^z-1)}

From Gao and Richmond's "Central and local limit theorems applied to asymptotic enumeration IV: multivariate generating functions"

https://www.sciencedirect.com/science/article/pii/037704279290247U

$[z^nx^k]F(z, x) = a_{n, k}/n!$ , with $a_{n, k}$ equal to the number of cyclic partitions of a set of size n with k blocks. In other words, the blocks of the partition are arranged in a cyclical order.

In [3]:

# When searching for critical points in the direction [z, x] = [1-p, p], we first verify there are no non-smooth points:
H = 1 - x*(exp(z) - 1)
Hx = H.diff(x)
Hz = H.diff(z)
print(solve([H, Hx, Hz], [x, z]))
# (This check is trivial by hand.)

[

]

In [4]:

# Then, have the following smooth critical point equations:
print(H, "= 0")
print((1-p)*x*Hx, " = ", p*z*Hz)

-x*(e^z - 1) + 1 = 0
(p - 1)*x*(e^z - 1)  =  -p*x*z*e^z

We can rewrite the second equation as follows:

-\frac{1-p}{p} e^{-\frac{1-p}{p}} = \left(z - \frac{1-p}{p}\right)e^{\left(z - \frac{1-p}{p}\right)}

Then, by definition of the Lambert W function $W_k(x)$ ,

z = \frac{1-p}{p} + W_k\left(\frac{p-1}{p} e^{\frac{p-1}{p}}\right)

The Lambert W function is real exactly when k = -1 or 0. Additionally, we trivially have

\frac{1-p}{p} + W_{-1}\left(\frac{p-1}{p}e^{\frac{p-1}{p}}\right) = 0

and can verify by sight that z = 0 is a solution for all values of p. The solution of interest is thus when k = 0, which yields the positive z critical point coordinate

v = \frac{1-p}{p} + W_{0}\left(\frac{p-1}{p}e^{\frac{p-1}{p}}\right)

from which we also have the positive x critical point coordinate

w = \frac{1}{e^v - 1}

The form of the GF makes it clear this is a positive minimal critical point, and hence it contributes to the asymptotics.

In [5]:

# We check the accuracy of this estimate for a specific direction.  For the combinatorially relevant range, p < 1/2.  We choose p = 1/8.  Then,
v = N(((1-p)/p + lambert_w((p-1)/p*exp((p-1)/p))).subs(p=1/8))
w = N(1/(exp(v) - 1))
print("The critical point is approximately (", v, ", ", w, ")")

The critical point is approximately ( 6.99357568672816 ,  0.000918602094210073 )

In [12]:

# Now, substituting into the theorem yields:
Hzz = Hz.diff(z)
Hzx = Hz.diff(x)
Hxx = Hx.diff(x)
lam = ((1-p)/p).subs(p=1/8)
chi1 = (Hx/Hz).subs(z=v, x=w)
chi2 = (1/(2*Hz)*(chi1^2*Hzz - 2*chi1*Hzx + Hxx)).subs(p=1/8, z=v, x=w)
M = (-2*chi2/v - chi1^2/v^2 - 1/(lam*w^2))
var('s')
Approx = factorial(lam*s)*((lam*s)^(-3/2)*v^(-(lam*s))*w^(-s))/sqrt(-2*pi*w^2*M)

In [22]:

# This approximation for s = 10, corresponding to a_{70, 10} (since the ratio lambda = (1-p)/p = 7):
approx70 = N(Approx.subs(s=10))
# The exact value:
F = log(1/(1 - x*(exp(z)-1)))
Fseries = F.taylor([x, z], 0, 81)
true70 = N(factorial(70)*Fseries.coefficient(z^70).coefficient(x^10))
print("The true value: ", true70)
print("The approximation: ", approx70)
print("The ratio: ", approx70/true70)

The true value:  9.93741615566349e68
The approximation:  9.95271893468655e68
The ratio:  1.00153991528415

In [24]:

# This approximation for s = 50, corresponding to a_{105, 15}:
approx70 = N(Approx.subs(s=15))
# The exact value:
F = log(1/(1 - x*(exp(z)-1)))
Fseries = F.taylor([x, z], 0, 121)
true70 = N(factorial(105)*Fseries.coefficient(z^105).coefficient(x^15))
print("The true value: ", true70)
print("The approximation: ", approx70)
print("The ratio: ", approx70/true70)

The true value:  2.03622855313551e122
The approximation:  2.03832272703266e122
The ratio:  1.00102845719059

In [0]:

CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual users to large groups and classes!

Theorem 1:

Example 1: (Cyclical Interlaced Permutations)

Identify smooth minimal critical point in the direction $(1,\ell)$ :

Compute asymptotics in direction $(1,\ell)$ :

Test results:

Example 2: (Necklaces)

Checking for non-smooth critical points:

Asymptotics are determined by the first term:

Example 3: (Logarithm of Narayana Numbers)

Example 4: (An Artificial Example)

Identify smooth critical point in direction $(1,\ell)$

Compute asymptotics in the direction $(1,\ell)$

Test results:

Example 5: (Cyclic Partitions of an $n$ -set) [Gao and Richmond 1992]

Product

Resources

Company

CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual users to large groups and classes!

Theorem 1:

Example 1: (Cyclical Interlaced Permutations)

Identify smooth minimal critical point in the direction (1,ℓ)(1,\ell)(1,ℓ):

Compute asymptotics in direction (1,ℓ)(1,\ell)(1,ℓ):

Test results:

Example 2: (Necklaces)

Checking for non-smooth critical points:

Asymptotics are determined by the first term:

Example 3: (Logarithm of Narayana Numbers)

Example 4: (An Artificial Example)

Identify smooth critical point in direction (1,ℓ)(1,\ell)(1,ℓ)

Compute asymptotics in the direction (1,ℓ)(1,\ell)(1,ℓ)

Test results:

Example 5: (Cyclic Partitions of an nnn-set) [Gao and Richmond 1992]

Identify smooth minimal critical point in the direction $(1,\ell)$ :

Compute asymptotics in direction $(1,\ell)$ :

Identify smooth critical point in direction $(1,\ell)$

Compute asymptotics in the direction $(1,\ell)$

Example 5: (Cyclic Partitions of an $n$ -set) [Gao and Richmond 1992]