Closed form solution for $ \int\limits_{0}^{1} x^{x} (1-x)^{1-x} dx $

Question

Does anyone have any hints or ways forward (or even better, a solution) for this integral?

$$\int_0^1 x^x (1-x)^{1-x} dx $$

I've tried multiple contour's but none seem to work. The only path forward I've managed to find after extensive searching is this: \begin{align} I &= \int_0^1 \left( \sum_{n=0}^\infty \frac{(x \ln x)^n}{n!} \right) \left( \sum_{m=0}^\infty \frac{((1-x) \ln(1-x))^m}{m!} \right) ,円 dx \\ &= \sum_{n=0}^\infty \sum_{m=0}^\infty \frac{1}{n! m!} \int_0^1 (x \ln x)^n ((1-x) \ln(1-x))^m ,円 dx \end{align}

But this is just so much more complicated :(. I'm trying to find a way to do this using complex analysis but I'd also love to see other solutions using different methods. The numerical value of the integral according to WA is $$\int_0^1 x^x (1-x)^{1-x} ,円 dx \approx 0.617826...$$

Any help would be greatly appreciated.

• Shaun: Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Mathematics Meta, or in Mathematics Chat. Comments continuing discussion may be removed.
• user1707952: @ТymaGaidash Please do share. I don't really know anything about beta functions but would love to learn :)
• Accelerator: FWIW... I got it to equal $\int_{-\infty}^{\infty}\frac{e^{2u}}{\left(1+e^{u}\right)^{3}}\exp\left(\frac{-u}{1+e^{u}}\right)du$. Using a semicircle in the upper half plane with radius 2ドルn\pi,ドル where $n$ is a large positive integer, Desmos seems to suggest that the sequence of integrals over the semicircles goes to 0ドル$ due to the semicircles' circumferences each having a length that's purposely placed to avoid the essential singularites. The issue is that I can't figure out how to express the integrand as a Laurent series that has a term like $\frac{\text{some coefficient}}{z-(2n+1)\pi i}$.

Answer 1

Too long for a comment.

It could be interesting to find the patterns in the coefficients of the expansion $$x^x,円(1-x)^{1-x}=\sum_{n=0}^\infty \frac{1}{n!} P_n(L) ,円x^n$$ where $L=\log(x)$, the first coefficients being $$\left( \begin{array}{cc} n & P_n(L) \\ 0 & 1 \\ 1 & L-1 \\ 2 & L^2-2 L+2 \\ 3 & L^3-3 L^2+6 L-3 \\ 4 & L^4-4 L^3+12 L^2-12 L+8 \\ 5 & L^5-5 L^4+20 L^3-30 L^2+40 L-10 \\ 6 & L^6-6 L^5+30 L^4-60 L^3+120 L^2-60L+54 \\ 7 & L^7-7 L^6+42 L^5-105 L^4+280 L^3-210 L^2+378 L+42 \\ 8 & L^8-8 L^7+56 L^6-168 L^5+560 L^4-560 L^3+1512 L^2+336 L+944 \\ \end{array} \right)$$ and to use

$$\int_0^1 x^m,円\log^n(x),円dx=(-1)^n \frac {\Gamma (n+1) } {(m+1)^{n+1} }$$

• Travis Willse: Cf. oeis.org/A293472
• user1707952: Hi Claude, I'm actually finding these numerical approximation answers quite nice and am very interested in learning more about this area of math. Do you have any book or online resources for learning more?
• Claude Leibovici: @RimuruTempest. This is a Taylor series. Cheers :-)
• user1707952: @ClaudeLeibovici I am aware of what it is, what I am asking is do you have recommendations for good books to learn numerical approximation techniques.
• Claude Leibovici: @RimuruTempest. I am so old that I do not know the books for studies. Sorry about that.

Answer 2

Too long for a comment. Let $I$ denote the integral. Then by substituting $x = \frac{1-y}{2}$, we get

$$\begin{align*} I &= \frac{1}{2} \int_{0}^{1} \sqrt{1-y^2} \exp(y\operatorname{artanh}y) ,円 \mathrm{d}y \\ &= \frac{1}{2} \int_{0}^{1} \exp\biggl( \sum_{k=1}^{\infty} \frac{x^{2k}}{(2k-1)(2k)} \biggr) ,円 \mathrm{d}y. \end{align*} $$

The last integrand can be expanded using integer partitions. Indeed,

$$ \exp\biggl( \sum_{k=1}^{\infty} \frac{x^{2k}}{(2k-1)(2k)} \biggr) = \sum_{n=0}^{\infty} \sum_{1^{m_1}2^{m_2}\cdots \vdash n} \frac{y^{2n}}{\prod_k m_k! [(2k-1)(2k)]^{m_k}}, $$

where the inner sum is taken for all integer partitions 1ドル^{m_1}2^{m_2}\cdots$ of $n$. Therefore $I$ can be expanded as

$$ I = \frac{1}{2} \sum_{n=0}^{\infty} \frac{1}{2n+1} \sum_{1^{m_1}2^{m_2}\cdots \vdash n} \frac{1}{\prod_k m_k! [(2k-1)(2k)]^{m_k}}. $$

Answer 3

You can start with a series expansion, where near $x=0$, the function approaches 0ドル$ instead of 1ドル$, but more terms make the series closer to the function and error go to 0ドル$ graphically and extract derivatives using $\ln^n(x)=\left.\frac{d^a}{da^n}x^a\right|_{a=0}$:

$$\int_0^1 x^x (1-x)^{1-x} dx=\sum_{n=0}^\infty\frac1{n!}\int_0^1x\ln^n\left(\frac x{1-x}\right)dx=\sum_{n=0}^\infty\frac{(-1)^n}{n!}\frac{d^n}{da^n}\int_0^1x^{a+1}(1-x)^{n-a}dx\biggl|_{a=0}=\sum_{n=0}^\infty\frac{(-1)^n}{n!}\frac{d^n}{da^n}\operatorname B(a+2,n-a+1)\biggl|_{a=0}=\frac12+0+\frac{\pi^2-6}{72}+\frac{\pi^2-6}{144}+\frac{14\pi^4-75\pi^2-270}{21600}+\dots$$

The convergence is slow, but is numerically shown here. The reflection formula can be used to give Pochhammer symbols with a finite series like this in powers of $a-u$ times $\csc(\pi a)$. Afterwards, the derivatives can be taken.

• PM 2Ring: I explored your solution in Sage, but those beta function derivatives quickly involve a lot of terms.
• Тyma Gaidash: Solution without derivatives with the input sum[sum[1/(n+2)!π^(n-k+1)(StirlingS1[1+n,k-1]+StirlingS1[1+n,k])((-1)^((n-k+1)/2-1)(2^(n-k+1)-2)BernoulliB[n-k+1])/(n-k+1)!,{k,1,20}],{n,0, 20}]

Answer 4

This answer is derived from the $\log(x)$ series described in Claude's answer:

$$\sum_{n=0}^\infty \frac{1}{n!} P_n(L) ,円x^n$$

As Travis Willse mentions in a comment, the coefficients for that log series are related to A293472, which gives "coefficients of polynomials in t = log(x) of the n-th derivative of x^x, evaluated at x = 1". Sadly, there's no magic formula for the coefficients, we have to get them the hard way, by computing derivatives. (Note that for Claude's series, we actually need the sign-alternating versions of those coefficients).

Let $$z = x^x (1-x)^{1-x}$$

That log series nicely approximates $z$ for small $x$, but it takes a lot of terms for it to get close to $z$ for $x>1/2$.

Here are the coefficients for 8 terms, where each term corresponds to a power of $x$.

0 [1]
1 [-1, 1]
2 [2, -2, 1]
3 [-3, 6, -3, 1]
4 [8, -12, 12, -4, 1]
5 [-10, 40, -30, 20, -5, 1]
6 [54, -60, 120, -60, 30, -6, 1]
7 [42, 378, -210, 280, -105, 42, -7, 1]

Here's a plot showing the resulting function $s$ (in green) and $z$ (in blue).

Plot of z & s, 8 terms

And here's a plot showing the error $s-z$ over $[0, 1/2]$ plot of s-z, 8 terms

So if we exploit the symmetry around $x=1/2$, we can get much faster convergence by evaluating the integral over $[0, 1/2]$ and doubling it.

Of course, we end up dealing with messier integrals for each term in the $s$ series, but at least they have a closed form, using the upper incomplete gamma function: $$\Gamma(s, x) = \int_x^\infty t^{s-1}e^{-t},円dt$$

In place of Claude's compact equation, we get

$$\int_0^{1/2}x^a \log^b(x),円dx =$$

$$(-1)^b \left(\frac{2^{a + 1} b \Gamma(b, (a + 1)\log(2)) + (a + 1)^b \log^b(2)} {2^{a+1}(a + 1)^{b+1}}\right)$$

I have implemented this procedure in Sage.

""" x^x * (1 - x)^(1 - x) integral on [0, 1].
 https://math.stackexchange.com/a/5106309
 Also see https://oeis.org/A293472
 Written by PM 2Ring 2025年11月08日
"""
RF = RealField(prec=64)
# integral(x^a * log(x)^b, x, 0, 1/2)
def ff(a, b):
 aa = a + 1
 return (-1)^b * (b * gamma(b, aa * log(2)) + aa^b * log(2)^b / 2^aa) / aa^(b + 1) #
var('x,t')
y = x^x
d = {log(x): -t}
s = ifac = 1
for i in srange(1, 12):
 y = y.diff(x)
 q = y.subs(d).expand()(x=1)
 ifac *= -i
 s += 2 * sum(RF(k * ff(i, b))
 for k,b in q.coefficients()) / ifac
 print(f"{i:>2}:", s)

Here's its output for 12 terms

 1: 0.451713204860013673
 2: 0.669118804205858581
 3: 0.603691681804339028
 4: 0.620885573710924924
 5: 0.617151698289766411
 6: 0.617929873049376159
 7: 0.617802795217129002
 8: 0.617828024247258519
 9: 0.617825822660470872
10: 0.617826732883800971
11: 0.617826843140138036

Here's an interactive live version running on the SageMathCell server. You can select the number of terms and the precision (in bits). It takes 44 terms to get 64 bit precision: 0.617826964729011473

Here's the interactive script which prints the coefficients, and does the plots.

---

Sage can easily print $s$ in symbolic form, but the expression is rather large for large numbers of terms. Here's what it looks like in Sage notation for 6 terms (which evaluates to 0.617151698289766):

-1/23040*log(2)^5 + 7/23040*log(2)^4 - 23/5760*log(2)^3 + 119/3840*log(2)^2 - 1/559872*gamma(5, 6*log(2)) - 1/23328*gamma(4, 6*log(2)) + 1/9375*gamma(4, 5*log(2)) - 1/1296*gamma(3, 6*log(2)) + 1/625*gamma(3, 5*log(2)) - 1/256*gamma(3, 4*log(2)) - 1/216*gamma(2, 6*log(2)) + 2/125*gamma(2, 5*log(2)) - 1/32*gamma(2, 4*log(2)) + 2/27*gamma(2, 3*log(2)) - 557/2880*log(2) + 123961/172800

• PM 2Ring: The best value I can get before the server kicks me off is 0.61782696472901147331619643118758832
• user1707952: I managed to get: 0.61782696472901147331602683571104152694
• PM 2Ring: @RimuruTempest It depends on how busy the server is. If you're lucky, you can get 95 terms or so (at 128 bit precision).
• user1707952: Yep, it takes 106 terms to get to 128 bit precision. This method is pretty neat
• user1707952: How did you learn Sage?