Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Trying to compute this limit, things seem to get frozen, even for hours. What would you
recommend to fastly compute it?

Limit[(Integrate[Product[x^k (1 - x^k), {k, 1, n}], {x, 0, 1}])^(1/n),n -> Infinity]
share|improve this question
up vote 20 down vote accepted

I very much suspect that the limit is not $1/4$, but rather $$\exp \left( \max_{y>0} \int_{t=0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt \right) \approx 0.185155.$$

If I were writing this up on math.SE, I'd start right in on a proof of this, but on this site I think that it would be more welcome to show some of the Mathematica techniques I used to come up with this.

The first step in estimating an integral is to see where it is large.

integrand[x_, n_] := Product[x^k (1 - x^k), {k, 1, n}]
xMax[n_] := xMax[n] = x /. Last[NMaximize[{integrand[x, n], 0 < x < 1}, x]]
yMax[n_] := yMax[n] = First[NMaximize[{integrand[x, n], 0 < x < 1}, x]]
Do[{xMax[n], yMax[n]}, {n,5,30}]

Here xMax is where the function is largest, and yMax is its largest value. The NMaximize's take time, so I am computing the first bunch once and for all.

My first goal is to get an idea of the shape of the function. Right now, I don't care so much about its absolute size, I want to know what it looks like. Here is integrand, renormalized to take values from $0$ to $1$.

Plot[Table[integrand[x, 5*m]/yMax[5*m], {m, 1, 5}] , {x, 0, 1},PlotRange -> {0, 1}]

enter image description here

So we're talking about a spike near $1$, which gets narrower as $n \to \infty$. The width of that spike looks like $1/n^c$ for some small constant $c$ -- in any case not exponential -- so $\int_{x=0}^1 (\mbox{your function}) dx$ should be roughly $n^{-c} \max_{0 < x < 1} (\mbox{your function})$ and we should have $$\lim_{n \to \infty} \left( \int_{x=0}^1 (\mbox{your function}) dx \right)^{1/n} \approx \lim_{n \to \infty} \left( \max_{0 < x <1 } (\mbox{your function}) dx \right)^{1/n}.$$ This will be the first of many steps I won't make rigorous.

So, let's figure out how tall that spike is and, while we're at it, where it is. Here is a semi-log plot of the maximum value:

ListPlot[Table[Log[yMax[n]], {n, 5, 30}]]

enter image description here

Straight as a pin! So yMax is dropping exponentially, and the constant in that exponential is the constant you want. For the fun of it, let's do a best fit line:

Fit[Table[{n, Log[yMax[n]]}, {n, 5, 30}], {n, 1}, n]
(* Output 0.795824 - 1.6565 n *)

That $-1.66$ is pretty good: I'll later calculate that it is actually $-1.69$. For now, I'll continue the mathematical analysis.

The next step is to figure out where that maximum is. We saw above that xMax is approaching $1$, so the most useful number is 1-xMax[n]. Here is a log-log plot. (Disclaimer: I actually figured out what the answer should be with pencil and paper first, and did the plot to check my work.)

ListPlot[Table[{Log[n], Log[1 - xMax[n]]}, {n, 10, 30}]]

enter image description here

Looks linear again, suggesting that xMax is roughly $1-a n^{-b}$. Let's find out what $b$ is:

Fit[Table[{n, Log[1 - xMax[n]]}, {n, 10, 30}], {Log[n], 1}, n]
(* Output -0.0213984 - 0.908425 Log[n] *)

That $0.9$ is close enough to $1$ to suggest that, probably, xMax is acting like $1-a n^{-1}$. (I'm cheating: I already got an exponent of $1$ in my hand computation, so this is just confirming evidence.)

It is usually a good idea to make a change of variables which puts the maximum of the integrand in a position independent of $n$. Usually, I'd use $x=1-y/n$, but in this case $x$ is being raised to a lot of powers, so I went with $x=e^{-y/n}$. So our integrand is $$\prod_k e^{-k y/n} (1-e^{-ky/n}) = \exp \left( \sum_k - y (k/n) + \log(1-e^{-y k/n}) \right).$$ That sum looks like a Riemann sum, so (nonrigorously) $$\prod_k e^{-k y/n} (1-e^{-ky/n}) \approx \exp \left( n \int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right)dt \right).$$

How good is this approximation?

    {integrand[E^(-y/(5*m)), 5*m], 
     E^(5*m*NIntegrate[-y*t + Log[1 - E^(-y*t)], {t, 0, 1}])}, 
    {y, 0, 3}, PlotRange -> {0, yMax[5*m]}]  , {m, 1, 5}]

enter image description here enter image description here

Not great. Here I show the $n=5$ and $n=25$ plots, omitting the intermediate ones. It looks like I missed a constant factor somewhere, or else the convergence is slow. But notice that it is only off by a factor of $2$ or $3$, even though the value of the function is changing by many orders of magnitude, and I nailed the basic shape of the function. Any constant factor will be swamped by the $1/n$ in the exponent when I take the limit.

In short, I believe that the integrand is approximately $\exp \left( n \int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right)dt \right)$. I also pick up a $(1/n) e^{-y/n} dy$ in place of $dx$, but that won't effect the asymtotics on this crude level. The largest value of the integrand should then be roughly $\exp \left( n \max_{y>0} \int_{t=0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt \right)$ and the $n$-th root of the integral should approach $\exp \left( \max_{y>0} \int_{t=0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt \right)$.

Before we compute that maximum, let's take a look at the function we are maximizing:

Plot[NIntegrate[-y*t + Log[1 - E^(-y*t)], {t, 0, 1}], {y, 0, 3}]

enter image description here

Good, nothing funny about it. A nice smooth peak. Incidently, this helps us analyze the width of the spike from back at the start: For $\exp \left( n \int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right)dt \right)$ to be above $1/e$ of it's maximum value, $\int_{t=0}^1 \left(-yt + \log(1-e^{yt}) \right) dt$ must be no less than $1/n$ below its maximum. Since that curve looks like it is smooth with a nonzero second derivative at it's max, this means that $y$ must be within $\pm 1/n^{1/2}$ of its optimum, so $x$ must be within $\pm 1/n^{3/2}$ of its optimum. In short, the width of the spikes from my first figure is probably $1/n^{3/2}$; a careful proof would have to establish this in order to justify approximating the integral by the maximum of its integrand.

There is only one thing left to do.

maxIntegral = NMaximize[{NIntegrate[-y*t + Log[1 - E^(-y*t)], {t, 0, 1}], 0 < y}, y]
(*  {-1.68656, {y -> 1.40505}} *)

(* 0.185155 *)
share|improve this answer
So (BesselI[0, 1] + Log@2)^Catalan? – Dr. belisarius Dec 17 '13 at 4:05
Very nice. And an upvote (I mention this in part to surpass the 15 character minimum). – Daniel Lichtblau Dec 17 '13 at 15:26
@belisarius How did you come to that? – ziyuang Dec 17 '13 at 20:52
@ziyuang‌​ – Dr. belisarius Dec 24 '13 at 0:17
In the formula after "That sum looks like a Riemann sum..." is a typo, inside the integral must be e^(-yt), not e^(yt) – Vaclav Kotesovec Jun 9 '15 at 16:42

Even computing this numerically for large n is a bit of work..

integrand[x_?NumericQ, n_?IntegerQ] := Product[x^k (1 - x^k), {k, 1, n}]
       NIntegrate[integrand[x, n], {x, 0, 1}, WorkingPrecision -> 100, 
       MaxRecursion -> 20]^(1/n) }, {n, {1, 2, 3, 8, 10, 20, 50, 100, 1000}}], Joined -> True]

Looks to be approaching ~0.1841

Edit..perhpas useful.. We can series expand the product:

Series[ Product[x^k (1 - x^k), {k, 1, n}] , {n, Infinity, 1}]

This allows a somewhat faster numerical integration

n = 10000; 
NIntegrate[x^(n/2 + n^2/2) QPochhammer[x, x, n] , {x, 0, 1}, 
     WorkingPrecision -> 100, MaxRecursion -> 20]^(1/n)

-> 0.1850069

Also, since QPochhammer[x, x, large_n] is bounded between 0:1 we can do this:

 Limit[Simplify[ Integrate[x^(n/2 + n^2/2) (1)  , {x, 0, 1} ]^(1/n) , 
         Assumptions -> {n > 0}] , n -> Infinity] -> 1

Which I'd take as evidence that the original expression is convergent to something less than 1.

share|improve this answer
Yeah, it seems the limit exceeds $1/(2e)$. +1 for showing that. Finally, the limit should hit $1/4$ as expected. – A_math_ninja Dec 16 '13 at 20:21

Just a note in the integration that takes a bit of CPU time. As the integrand is a polynomial its integral 0->1 can be obtained replacing $x^k$ by $1/(k+1)$ throughout. This will compute the integral by first computing the coefficient of the polynomial :

F[n_] := Module[{coefList},
   coefList = CoefficientList[Product[x^k (1 - x^k), {k, 1, n}], x];
   (Total[coefList[[#]]/# & /@ Range@Length@c[n] ])^(1/n) // N
F[20] // Timing 
(* {0.007211, 0.16846} *)
F[500] // Timing 
(* {78.636256, 0.183303} *)

whereas Integrate is way slower :

integrand[x_, n_] := Product[x^k (1 - x^k), {k, 1, n}]
With[{n = 20}, (Integrate[integrand[x, n], {x, 0, 1}])^(1/n) // N] // Timing
(* {9.616129, 0.16846} *)

F also appears to converge :

DiscretePlot[F[k], {k, 10, 200, 10}, Joined -> True]

BTW, the coefficients are rather neat, look at

   With[{n = 100},
     ListPlot[CoefficientList[Product[x^k (1 - x^k), {k, 1, n}], x], 
     PlotRange -> Full, ImageSize -> Large]]

(I updated this plot with <code>PlotRange->Full</code>

share|improve this answer

Products are polynomials in x:

Table[Expand[Product[x^k*(1-x^k), {k, 1, n}]], {n, 1, 3}]
(* {x - x^2, x^3 - x^4 - x^5 + x^6, x^6 - x^7 - x^8 + x^10 + x^11 - x^12} *)

Integration of polynomials is easy:

Table[Integrate[Expand[Product[x^k*(1-x^k), {k, 1, n}]], x], {n, 1, 3}]
(* {x^2/2 - x^3/3, x^4/4 - x^5/5 - x^6/6 + x^7/7, x^7/7 - x^8/8 - x^9/9 + x^11/11 + x^12/12 - x^13/13} *)

For integral x=(0..1) set x=1:

n=1 1/2 - 1/3 = 1/6
n=2 1/4 - 1/5 - 1/6 + 1/7 = 11/420
n=3 1/7 - 1/8 - 1/9 + 1/11 + 1/12 - 1/13 = 293/72072

This efficient program generate a partial values of the integral (example is for the first six terms):

nmax=6; p=1; Table[p=Expand[p*x^n*(1-x^n)]; Total[CoefficientList[p,x]/Range[1,Exponent[p,x]+1]],{n,1,nmax}]
(* {1/6, 11/420, 293/72072, 487/760760, 129952159/1266697832400, 13084761625/783333734619744} *)

1000 terms took 23 minutes

nmax=1000; p=1; pint=Table[PrintTemporary[n]; p=Expand[p*x^n*(1-x^n)]; (Total[CoefficientList[p,x]/Range[1,Exponent[p,x]+1]]),{n,1,nmax}];

The convergence is slow

N[Table[pintn[[n]]^(1/n), {n, 995, 1000}], 40]

Now I used the Richardson extrapolation and get a numerical values of the limit


Summary: I confirmed that $\lim\limits_{n\to\infty} \left(\int_0^1 \prod\limits_{k=1}^n x^k (1-x^k)\, \mathrm dx\right)^{1/n} = 0.18515\dots$

I added a new sequences to the OEIS, a values of the integral (numerators) (denominators)

The new results from 2000 terms (4 hours of CPU, 2 GB of RAM)


Richardson extrapolation:


Wynn epsilon:

(* 0.1851544656 *)

Conjecture: The limit is equal to

= 1 / 5.400871904118154152466091119104270052029...

see OEIS

A similar limit: $\lim\limits_{n\to\infty} \left(\int_0^1 \prod\limits_{k=1}^n (1-x^k)\, \mathrm dx\right) = 0.368412535931433652321316597327851\dots$

after the Euler pentagonal theorem is the limit equal to

exact solution

I created a new sequences (numerators), (denominators).

share|improve this answer
An alternative to Richardson extrapolation would be to use the Wynn epsilon algorithm. In Mathematica, it's as easy as SequenceLimit[pintn^(1/Range[Length[pintn]])]. – J. M. May 23 '15 at 16:08
Yes, I used also extrapolation epsilon with result 0.185148345990579. But the results from Richardson extrapolation are in this case better. – Vaclav Kotesovec May 23 '15 at 16:12
I posted Mathematica routines for Richardson and Wynn here, in case you might be interested. – J. M. May 23 '15 at 16:32

I don't know how to show the limit is 1/4. Ithink it can be bounded above by 1/4 as follows. Note that there is a gap in the argument below.

(1) Replace factors 1-x^k by 1-x^(n/2).

nn = Integrate[
  Product[x^k, {k, 1, n}]*Product[(1 - x^(n/2)), {k, 1, n}],
{x, 0, 1}, Assumptions -> n > 1000]

(* Out[147]= (2 Gamma[n] Gamma[1 + 2/n + n])/Gamma[2 (1 + 1/n + n)] *)

Limit[nn^(1/n), n -> Infinity]

(* Out[148]= 1/4 *)

(2) Investigate the error factor. (Also, check to make sure I didn't mess this up in some woeful manner.)

factor = 
 Limit[Product[(1 - x^(n/2))/(1 - x^k), {k, 1, n}], n -> Infinity, 
  Assumptions -> 0 < x < 1]

(* Out[152]= ((-1 + x)^2 (1 + x))/(QPochhammer[x, x] QPochhammer[x, x, 2]) *)

f2 = Together[FunctionExpand[factor, Assumptions -> 0 < x < 1]]

(* Out[176]= 1/QPochhammer[x, x] *)

(3) Show that this factor is bounded below by 1. This can be done by grouping the denominator in pairs of the form (1-x^k)*(1-x^(n-k)) and showing these are smaller than (1-x^(n/2))^2. Here I gotta run so this is a detail to check.

share|improve this answer
We can use AM-GM inequality to show the upper bound $$x^k (1-x^k) \le \left(\frac{x^k+1-x^k}{2}\right)^2=\frac{1}{2^2}$$ – A_math_ninja Dec 17 '13 at 8:39

A better way how compute the maximum of integral

$$\int_{0}^1 \left( -yt + \log(1-e^{-yt}) \right) dt $$

is following. Derivative must be equal to zero at the maximum.

FullSimplify[D[Integrate[-y*t + Log[1 - Exp[-y*t]], t], y]]

$$\frac{\text{Li}_2\left(e^{t y}\right)+t y \left(\log \left(1-e^{t y}\right)-t y\right)}{y^2}$$

For t=1 we have

(t y (-t y + Log[1-E^(t y)]) + PolyLog[2,E^(t y)])/y^2/.t->1

$$\frac{\text{Li}_2\left(e^y\right)+y \left(\log \left(1-e^y\right)-y\right)}{y^2}$$

For t=0 we must use a limit

Limit[(t y (-t y + Log[1-E^(t y)]) + PolyLog[2,E^(t y)])/y^2, t->0]

$$\frac{\pi ^2}{6 y^2}$$

Maximum is at the point y, where y is the root of the equation $$-\text{Li}_2\left(e^y\right)+y^2-y \log \left(1-e^y\right)+\frac{\pi ^2}{6}=0$$

We have

FindRoot[Pi^2/6 + y^2 - y*Log[1-E^y] - PolyLog[2,E^y] == 0, {y,1}, WorkingPrecision->100]
(* {y->1.405054216578477533211873399082087670321440412190343393755288523553823298202008968431771181607706332} *)

The final value of the limit is

Re[Exp[Integrate[-y*t+Log[1-Exp[-y*t]] /.FindRoot[Pi^2/6 + y^2 - y*Log[1-E^y] - PolyLog[2,E^y]==0, {y,1}, WorkingPrecision->100], {t,0,1}]]]
(* 0.1851552893223595946473132111979542852738223644590750839856055303609697606866994281773438299262 *)

Or simplified

(r-1)/r^2 /.FindRoot[-PolyLog[2, 1-r] == Log[r]^2, {r, E}, WorkingPrecision->100]
(* 0.185155289322359594647313211197954285273822364459075083985605530360969760686699428177343829926166036 *)

$\lim\limits_{n\to\infty} \left(\int_0^1 \prod\limits_{k=1}^n x^k (1-x^k)\, \mathrm dx\right)^{1/n} = 0.18515528932235959464731321119795428527\dots$

This result agree with my previous conjecture, see

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.