SM-7-Partitions

4137 days ago by MathFest

First, let us set up our precision, and some variables we will use throughout the worksheet.

R=RealField(1000); var('t,c,x,y') 
       
(t, c, x, y)
(t, c, x, y)

In this worksheet, we will investigate some questions regarding variants of the partition function, and in particular, how we used Sage to discover and prove that the partition function is almost log concave.  More precisely, defining $p(n)$ to be the number of partitions of $n$, we show that for $n>26$

\[

p(n)^2 > p(n-1)p(n+1)

\]

Prior to this, the inequality is true only for even $n$.

 

First, let's just look at the first 10 partition numbers: Sage uses number_of_partitions(n) to denote this.

for n in range(0,10): number_of_partitions(n) 
       
1
1
2
3
5
7
11
15
22
30
1
1
2
3
5
7
11
15
22
30

We can also, say, list all partitions of $n$, or all partitions of $n$ with a given number of parts:  one nice feature of partitions is that we have a bijection between partitions with exactly $k$ parts, and partitions with largest part equal to $k$.

Partitions(7) 
       
Partitions of the integer 7
Partitions of the integer 7

This suggests that Sage really doesn't want us to look at this list without telling it explicitly:

use list(partitions(7)) to view the elements:

list(Partitions(7)) 
       
[[7], [6, 1], [5, 2], [5, 1, 1], [4, 3], [4, 2, 1], [4, 1, 1, 1], [3, 3,
1], [3, 2, 2], [3, 2, 1, 1], [3, 1, 1, 1, 1], [2, 2, 2, 1], [2, 2, 1, 1,
1], [2, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]
[[7], [6, 1], [5, 2], [5, 1, 1], [4, 3], [4, 2, 1], [4, 1, 1, 1], [3, 3, 1], [3, 2, 2], [3, 2, 1, 1], [3, 1, 1, 1, 1], [2, 2, 2, 1], [2, 2, 1, 1, 1], [2, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]

The more approved way of doing this is with Partitions(7).list

Partitions(7).list() 
       
[[7], [6, 1], [5, 2], [5, 1, 1], [4, 3], [4, 2, 1], [4, 1, 1, 1], [3, 3,
1], [3, 2, 2], [3, 2, 1, 1], [3, 1, 1, 1, 1], [2, 2, 2, 1], [2, 2, 1, 1,
1], [2, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]
[[7], [6, 1], [5, 2], [5, 1, 1], [4, 3], [4, 2, 1], [4, 1, 1, 1], [3, 3, 1], [3, 2, 2], [3, 2, 1, 1], [3, 1, 1, 1, 1], [2, 2, 2, 1], [2, 2, 1, 1, 1], [2, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]

In order to investigate what methods can be applied to a list of partitions, give this a name, part7,

and type

   part7.

followed by tab

part17=Partitions(17) 
       
part17 
       
Partitions of the integer 17
Partitions of the integer 17

We can play with each partition in this list:

p1=part17[135]; p1 
       
[6, 6, 4, 1]
[6, 6, 4, 1]

Again, to see which methods are available, type

  p1.

followed by tab

p1. 
       
Traceback (click to the left of this block for traceback)
...
SyntaxError: invalid syntax
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_11.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cDEu"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpu9Wz9C/___code___.py", line 2
    p1.
      ^
SyntaxError: invalid syntax
fd1=p1.ferrers_diagram();fd1 
       
'******\n******\n****\n*'
'******\n******\n****\n*'

One way to display this nicely is to paste it into the text window, replacing the line feeds by new lines.

I don't know how to get Sage to do it nicely.

******
******
****
*

There are other things we can do to p1.

p2=p1.conjugate();p2 
       
[4, 3, 3, 3, 2, 2]
[4, 3, 3, 3, 2, 2]

This is the partition created by taking the Ferrers Diagram and reflecting it in the line $y=-x$.

fd2=p2.ferrers_diagram();fd2 
       
'****\n***\n***\n***\n**\n**'
'****\n***\n***\n***\n**\n**'

With the line feeds shown, this becomes

****
***
***
***
**
**

We can also list partitions by the size of the largest part:

PartitionsGreatestEQ(10,4).list() 
       
[[4, 4, 2], [4, 4, 1, 1], [4, 3, 3], [4, 3, 2, 1], [4, 3, 1, 1, 1], [4,
2, 2, 2], [4, 2, 2, 1, 1], [4, 2, 1, 1, 1, 1], [4, 1, 1, 1, 1, 1, 1]]
[[4, 4, 2], [4, 4, 1, 1], [4, 3, 3], [4, 3, 2, 1], [4, 3, 1, 1, 1], [4, 2, 2, 2], [4, 2, 2, 1, 1], [4, 2, 1, 1, 1, 1], [4, 1, 1, 1, 1, 1, 1]]

and we can pick out one to work with

PartitionsGreatestEQ(10,4)[5] 
       
[4, 2, 2, 2]
[4, 2, 2, 2]
PartitionsGreatestEQ(10,4)[5].conjugate() 
       
[4, 4, 1, 1]
[4, 4, 1, 1]

We can even use the map function to apply conjugate to each element in the list

map(conjugate,PartitionsGreatestEQ(10,4)) 
       
[[3, 3, 2, 2], [4, 2, 2, 2], [3, 3, 3, 1], [4, 3, 2, 1], [5, 2, 2, 1],
[4, 4, 1, 1], [5, 3, 1, 1], [6, 2, 1, 1], [7, 1, 1, 1]]
[[3, 3, 2, 2], [4, 2, 2, 2], [3, 3, 3, 1], [4, 3, 2, 1], [5, 2, 2, 1], [4, 4, 1, 1], [5, 3, 1, 1], [6, 2, 1, 1], [7, 1, 1, 1]]

We see that each of these is the conjugate of a partition with largest part equal to 4, and hence is a partition with 4 parts.

Let's list the number of partitions of $n$ with exactly $k$ parts, as $k$ ranges from 1 to $n$:

for k in range(1,18): PartitionsGreatestEQ(18,k).cardinality() 
       
1
9
27
47
57
58
49
40
30
22
15
11
7
5
3
2
1
1
9
27
47
57
58
49
40
30
22
15
11
7
5
3
2
1

It appears that this list is unimodal: that is, it has a unique (possibly repeated) maximum.

This was shown in 1952, for $n$ sufficiently large, by Szekeres, and a few years ago, for $n \leq 2000$, by Canfield.

We've extended this computation out to $n\leq 25 000$, and are planning to continue further: note that since

$p(n)$ grows like $e^{c\sqrt{n}}$, to do this to $n$ requires about $n^{5/2}$ bits of storage.

We are simultaneously trying to make Szekeres' paper explicit, that is, give an explicit $n$ which is "large enough".

 
       

One alternative method to show that a sequence (finite or infinite) is unimodal is to show something rather stronger:

show that it is log-concave.  A sequence $a_n$ is log-concave if for every $n$,

\[

a_n^2 \geq a_{n+1} a_{n-1}.

\]

Let's see if the partitions of $n$ into $k$ parts ($n$ fixed, $k$ varying) is log concave:

def p(n,k): return PartitionsGreatestEQ(n,k).cardinality() 
       
for k in range(2,8): k,(p(8,k)^2>= p(8,k+1)*p(8,k-1)) 
       
(2, True)
(3, True)
(4, True)
(5, False)
(6, True)
(7, False)
(2, True)
(3, True)
(4, True)
(5, False)
(6, True)
(7, False)

At this point it is clear that the sequence is not log-concave, and we give up.

However, in actuality, in our research, a happy accident occurred.  Several of us

were checking this, and one of us (me) made a mistake, and instead of computing

$p(n,k)$, instead computed $P(n,k)$.

def P(n,k): return PartitionsGreatestLE(n,k).cardinality() 
       
for k in range(2,8): k,(P(8,k)^2 >= P(8,k+1)*P(8,k-1)) 
       
(2, True)
(3, True)
(4, True)
(5, True)
(6, True)
(7, True)
(2, True)
(3, True)
(4, True)
(5, True)
(6, True)
(7, True)

Now, this led us to consider more carefully what happened with $p(n,k)$.

for k in range(2,28): k,(p(28,k)^2>= p(28,k+1)*p(28,k-1)) 
       
(2, True)
(3, True)
(4, True)
(5, True)
(6, True)
(7, True)
(8, True)
(9, True)
(10, True)
(11, True)
(12, True)
(13, True)
(14, True)
(15, False)
(16, True)
(17, False)
(18, True)
(19, False)
(20, True)
(21, False)
(22, True)
(23, False)
(24, True)
(25, False)
(26, True)
(27, False)
(2, True)
(3, True)
(4, True)
(5, True)
(6, True)
(7, True)
(8, True)
(9, True)
(10, True)
(11, True)
(12, True)
(13, True)
(14, True)
(15, False)
(16, True)
(17, False)
(18, True)
(19, False)
(20, True)
(21, False)
(22, True)
(23, False)
(24, True)
(25, False)
(26, True)
(27, False)

Two observations come to mind: first, it appears to start out satisfying

the log-concave inequality and then alternately satisfying and not, 

and secondly, it's getting slow.

def p(n,k): if n>=0 and k>=1 and k<=n: if n<2*k: return(number_of_partitions(n-k)) return(number_of_partitions(n-1,k-1)+number_of_partitions(n-k,k)) 
       
for k in range(2,28): p(28,k) 
       
14
65
169
291
391
436
434
393
340
278
224
174
135
101
77
56
42
30
22
15
11
7
5
3
2
1
14
65
169
291
391
436
434
393
340
278
224
174
135
101
77
56
42
30
22
15
11
7
5
3
2
1
for k in range(2,28): k,(p(28,k)^2>= p(28,k+1)*p(28,k-1)) 
       
(2, True)
(3, True)
(4, True)
(5, True)
(6, True)
(7, True)
(8, True)
(9, True)
(10, True)
(11, True)
(12, True)
(13, True)
(14, True)
(15, False)
(16, True)
(17, False)
(18, True)
(19, False)
(20, True)
(21, False)
(22, True)
(23, False)
(24, True)
(25, False)
(26, True)
(27, False)
(2, True)
(3, True)
(4, True)
(5, True)
(6, True)
(7, True)
(8, True)
(9, True)
(10, True)
(11, True)
(12, True)
(13, True)
(14, True)
(15, False)
(16, True)
(17, False)
(18, True)
(19, False)
(20, True)
(21, False)
(22, True)
(23, False)
(24, True)
(25, False)
(26, True)
(27, False)
for k in range(2,28): if (p(28,k)^2 < p(28,k+1)*p(28,k-1)): print(k) 
       
15
17
19
21
23
25
27
15
17
19
21
23
25
27
def checklc(n): for k in range(2,n): if p(n,k)^2 < p(n,k+1)*p(n,k-1): print(k) return 
       
checklc(28) 
       
15
17
19
21
23
25
27
15
17
19
21
23
25
27
checklc(100) 
       
75
77
79
81
83
85
87
89
91
93
95
97
99
75
77
79
81
83
85
87
89
91
93
95
97
99
def checklc(n): for k in range(2,n): if p(n,k)^2 < p(n,k+1)*p(n,k-1): print(n-k) return 
       
checklc(100) 
       
25
23
21
19
17
15
13
11
9
7
5
3
1
25
23
21
19
17
15
13
11
9
7
5
3
1
checklc(200) 
       
25
23
21
19
17
15
13
11
9
7
5
3
1
25
23
21
19
17
15
13
11
9
7
5
3
1
checklc(300) 
       
25
23
21
19
17
15
13
11
9
7
5
3
1
25
23
21
19
17
15
13
11
9
7
5
3
1
checklc(301) 
       
25
23
21
19
17
15
13
11
9
7
5
3
1
25
23
21
19
17
15
13
11
9
7
5
3
1
 
       

We now have a clearer idea of what is going on: for $n$ sufficiently large, it appears that $p(n,k)$ satisfies the log-concavity inequality

except when $n-k$ is odd and at most 25.

Now, we know that the second half of the sequence $p(n,k)$ is given by $p(n-k)$, that is, it is the reverse of the first $n/2$ terms of 

the partition function, $p(n)$.  So, our first conjecture ought to be the following:

Conjecture: For $n\geq 26$, $p(n)^2 \geq p(n+1)p(n-1)$.  

In what follows, we'll prove this conjecture.  We are close to completing a proof of the remainder of the conjecture for $p(n,k)$ for 

sufficiently large $n$ (and, we hope, all $n$).

We will use the Rademacher series for $p(n)$ to obtain a asymptotic expression, and with this we will fully explain the above behaviour.

Switch to beamer presentation 

 

 

 

 

 

Now that we've deduced our asymptotic expression, with $t=n-1/24$:

\[ p(n)  = \frac{1}{4\pi \sqrt{2}t^{3/2}}  
\left(  c\sqrt{t} -1 \right) \exp\left(c\sqrt{t}\right) +
\frac{(-1)^n}{4 \pi t^{3/2}}
\left(  \frac{c\sqrt{t}}{2} -1 \right) \exp\left(\frac{c\sqrt{t}}{2}\right)  + O \left( t \exp\left(\frac{c\sqrt{t}}{3}\right) \right)
\]

\[

= \frac{\exp(c\sqrt{t})}{4\pi \sqrt{2}t^{3/2}}   \left(  
     \left(  c\sqrt{t} -1 \right) +
(-1)^n\sqrt{2}
\left(  \frac{c\sqrt{t}}{2} -1 \right) \exp\left(\frac{-c\sqrt{t}}{2}\right)  + O \left( t \exp\left(\frac{-2c\sqrt{t}}{3}\right) \right)  \right)

\]

\[

= \exp(c\sqrt{t})\frac{\left(  c\sqrt{t} -1 \right)}{4\pi \sqrt{2}t^{3/2}}    
 \left(    1 + \frac{(-1)^n}{\sqrt{2}}\frac{c\sqrt{t}-2}{c\sqrt{t}-1}\exp\left(\frac{-c\sqrt{t}}{2}\right)   

+ O \left(  \exp\left(\frac{-2c\sqrt{t}}{3}\right)    \right) \right)

\]

\[

= \exp(c\sqrt{t})\frac{\left(  c\sqrt{t} -1 \right)}{4\pi \sqrt{2}t^{3/2}}    
 \left(    1 + \frac{(-1)^n}{\sqrt{2}}\frac{c\sqrt{t}-2}{c\sqrt{t}-1}\exp\left(\frac{-c\sqrt{t}}{2}\right)   \right)  \left(

1 + O \left(  \exp\left(\frac{-2c\sqrt{t}}{3}\right)    \right) \right)

\]

We really ought to check that this looks accurate:

def pasymp1(t): return exp(c* sqrt(t))*(c*sqrt(t)-1)/(4*pi*sqrt(2)*t^(3/2)) 
       
pasymp1(1000-1/24) 
       
18/575952001*(sqrt(23999/6)*c -
2)*sqrt(2)*sqrt(23999/6)*e^(1/2*sqrt(23999/6)*c)/pi
18/575952001*(sqrt(23999/6)*c - 2)*sqrt(2)*sqrt(23999/6)*e^(1/2*sqrt(23999/6)*c)/pi
_.substitute(c=pi*sqrt(2/3)) 
       
18/575952001*(pi*sqrt(2/3)*sqrt(23999/6) -
2)*sqrt(2)*sqrt(23999/6)*e^(1/2*pi*sqrt(2/3)*sqrt(23999/6))/pi
18/575952001*(pi*sqrt(2/3)*sqrt(23999/6) - 2)*sqrt(2)*sqrt(23999/6)*e^(1/2*pi*sqrt(2/3)*sqrt(23999/6))/pi
_.n(digits=40) 
       
2.406146786403262243279475037438734476345e31
2.406146786403262243279475037438734476345e31
_.trunc() 
       
24061467864032622432794750374387
24061467864032622432794750374387
number_of_partitions(1000) 
       
24061467864032622473692149727991
24061467864032622473692149727991

Let's put in the term involving $\exp(c\sqrt{t}/2)$ as well:

def pasymp2(t): return exp(c*sqrt(t))*(c*sqrt(t)-1)/(4*pi*sqrt(2)*t^(3/2))*(1 + (-1)^(t+1/24)/sqrt(2)*(c*sqrt(t)-2)/(c*sqrt(t)-1)*exp(-c*sqrt(t)/2) ) 
       
pasymp2(1000-1/24) 
       
9/575952001*(sqrt(23999/6)*c - 2)*((sqrt(23999/6)*c -
4)*sqrt(2)*e^(-1/4*sqrt(23999/6)*c)/(sqrt(23999/6)*c - 2) +
2)*sqrt(2)*sqrt(23999/6)*e^(1/2*sqrt(23999/6)*c)/pi
9/575952001*(sqrt(23999/6)*c - 2)*((sqrt(23999/6)*c - 4)*sqrt(2)*e^(-1/4*sqrt(23999/6)*c)/(sqrt(23999/6)*c - 2) + 2)*sqrt(2)*sqrt(23999/6)*e^(1/2*sqrt(23999/6)*c)/pi
_.substitute(c=pi*sqrt(2/3)) 
       
9/575952001*((pi*sqrt(2/3)*sqrt(23999/6) -
4)*sqrt(2)*e^(-1/4*pi*sqrt(2/3)*sqrt(23999/6))/(pi*sqrt(2/3)*sqrt(23999/\
6) - 2) + 2)*(pi*sqrt(2/3)*sqrt(23999/6) -
2)*sqrt(2)*sqrt(23999/6)*e^(1/2*pi*sqrt(2/3)*sqrt(23999/6))/pi
9/575952001*((pi*sqrt(2/3)*sqrt(23999/6) - 4)*sqrt(2)*e^(-1/4*pi*sqrt(2/3)*sqrt(23999/6))/(pi*sqrt(2/3)*sqrt(23999/6) - 2) + 2)*(pi*sqrt(2/3)*sqrt(23999/6) - 2)*sqrt(2)*sqrt(23999/6)*e^(1/2*pi*sqrt(2/3)*sqrt(23999/6))/pi
_.n(digits=50) 
       
2.4061467864032622473692179982529763078831257878542e31
2.4061467864032622473692179982529763078831257878542e31
_.trunc() 
       
24061467864032622473692179982529
24061467864032622473692179982529
number_of_partitions(1000) 
       
24061467864032622473692149727991
24061467864032622473692149727991

And we see, as expected, that about two thirds of the digits are correct now.  We thus have some confidence that we've correctly interpreted Rademacher's formula, at least for the first two terms!

We now wish to consider the ratio

\[

\frac{p(n)p(n)}{p(n+1)p(n+1)}

\]

We will first consider the ratio of the exponential parts, then the ratio of the rational

function of $\sqrt{t}$ and finally the contribution for small $n$ of the terms involving $\exp(-c\sqrt{t}/2)$.

Let's start by naively trying to simplify the ratio of the exponential parts of our asymptotic expansion for $p(n)$:

\[

\frac{\exp(c\sqrt{t})\exp(c\sqrt{t})}{\exp(c\sqrt{t+1})\exp(c\sqrt{t-1})}

\]

exppart=exp(2*c*sqrt(t)-c*sqrt(t-1)+c*sqrt(t+1)) 
       
simplify(exppart) 
       
e^(-sqrt(t - 1)*c + sqrt(t + 1)*c + 2*c*sqrt(t))
e^(-sqrt(t - 1)*c + sqrt(t + 1)*c + 2*c*sqrt(t))
taylor(exppart, t,0,1) 
       
1/2*(4*c^2*e^c + (I + 1)*c*e^c)*t*e^(-I*c) + 2*c*sqrt(t)*e^((-I + 1)*c)
+ e^((-I + 1)*c)
1/2*(4*c^2*e^c + (I + 1)*c*e^c)*t*e^(-I*c) + 2*c*sqrt(t)*e^((-I + 1)*c) + e^((-I + 1)*c)

Clearly this is not going to give us what we want.  Let's look instead at the exponent

\[

\sqrt{t} -\sqrt{t-1} = \frac{1}{\sqrt{t} + \sqrt{t-1}}

\]

and

\[

\sqrt{t} - \sqrt{t+1} = \frac{-1}{\sqrt{t}+\sqrt{t+1}}

\]

Hence

\[ 2\sqrt{t} -\sqrt{t+1}-\sqrt{-1} = \frac{\sqrt{t+1}- \sqrt{t-1}}{(\sqrt{t}+\sqrt{t+1})(\sqrt{t}+\sqrt{t-1})}

= \frac{2}{(\sqrt{t}+\sqrt{t+1})(\sqrt{t}+\sqrt{t-1})(\sqrt{t+1}+\sqrt{t-1})} 

\]

\[

= \frac{1}{4 t^{3/2} }(1+O(1/t))

\]

simpexp=2/((sqrt(t)+sqrt(t+1))*(sqrt(t)+sqrt(t-1))*(sqrt(t+1)+sqrt(t-1))) 
       
simpexp2=simpexp.substitute(t=1/x); simpexp2 
       
2/((sqrt(1/x + 1) + sqrt(1/x))*(sqrt(1/x - 1) + sqrt(1/x))*(sqrt(1/x -
1) + sqrt(1/x + 1)))
2/((sqrt(1/x + 1) + sqrt(1/x))*(sqrt(1/x - 1) + sqrt(1/x))*(sqrt(1/x - 1) + sqrt(1/x + 1)))
taylor(simpexp2,x,0,4) 
       
5/64*x^(7/2) + 1/4*x^(3/2)
5/64*x^(7/2) + 1/4*x^(3/2)

So we now see that our exponential contribution is approximately

\[

\frac{\exp(c\sqrt{t})\exp(c\sqrt{t})}{\exp(c\sqrt{t+1})\exp(c\sqrt{t-1})} =

\exp\left( \frac{c}{4t^{3/2}} \left( 1 + \frac{5}{16t^2} + O(1/t^4) \right) \right)

\]

 
       

Let us now consider the rational function of $\sqrt{t}$.

def rat(t): return (c*sqrt(t)-1)/t^(3/2) 
       

We'll compute the ratio

\[

\frac{rat(t)rat(t)}{rat(t+1)rat(t-1)}

\]

ratterm= rat(t)*rat(t)/(rat(t+1)*rat(t-1)); ratterm 
       
(t - 1)^(3/2)*(t + 1)^(3/2)*(c*sqrt(t) - 1)^2/((sqrt(t + 1)*c -
1)*(sqrt(t - 1)*c - 1)*t^3)
(t - 1)^(3/2)*(t + 1)^(3/2)*(c*sqrt(t) - 1)^2/((sqrt(t + 1)*c - 1)*(sqrt(t - 1)*c - 1)*t^3)

In order to take the Taylor series of this, we replace $t$ by $1/x$.  Observe the use of the

ratterm.substitute functionality: the function is a method which can be used on ratterm,

rather than a function which is applied to an argument.

rattermx = ratterm.substitute(t=1/x); rattermx 
       
(1/x - 1)^(3/2)*(1/x + 1)^(3/2)*(c*sqrt(1/x) - 1)^2*x^3/((sqrt(1/x +
1)*c - 1)*(sqrt(1/x - 1)*c - 1))
(1/x - 1)^(3/2)*(1/x + 1)^(3/2)*(c*sqrt(1/x) - 1)^2*x^3/((sqrt(1/x + 1)*c - 1)*(sqrt(1/x - 1)*c - 1))

Now, there is a built-in function "taylor", which we can use to expand this

taylor(rattermx,x,0,4) 
       
-x^2 + 3/4*x^(5/2)/c + x^3/c^2 + 5/4*x^(7/2)/c^3 + 3/2*x^4/c^4 + 1
-x^2 + 3/4*x^(5/2)/c + x^3/c^2 + 5/4*x^(7/2)/c^3 + 3/2*x^4/c^4 + 1

Or we can use the method rattermx.taylor instead:

rattaylorx = rattermx.taylor(x,0,4); rattaylorx 
       
-x^2 + 3/4*x^(5/2)/c + x^3/c^2 + 5/4*x^(7/2)/c^3 + 3/2*x^4/c^4 + 1
-x^2 + 3/4*x^(5/2)/c + x^3/c^2 + 5/4*x^(7/2)/c^3 + 3/2*x^4/c^4 + 1

Since this has fractional powers of $x$ it is not a polynomial, and hence we can't

apply ``collect'' to it directly: we'll substitute in $y^2$ for $x$ to get around this.

assume(y>0); rattaylory= rattaylorx.substitute(x=y^2);rattaylory 
       
y^6/c^2 + 3/2*y^8/c^4 - y^4 + 3/4*(y^2)^(5/2)/c + 5/4*(y^2)^(7/2)/c^3 +
1
y^6/c^2 + 3/2*y^8/c^4 - y^4 + 3/4*(y^2)^(5/2)/c + 5/4*(y^2)^(7/2)/c^3 + 1

Notice that even though we've assumed that $y$ is positive, it doesn't simplify $(y^2)^{5/2}$ to $y^5$ yet.

We will apply the simplify method, telling it to simplify radicals.  (If we hadn't assumed $y>0$ here it

would include a term abs($y$) as well.)

rattaylory=rattaylory.simplify_radical();rattaylory 
       
-1/4*(4*c^4*y^4 - 3*c^3*y^5 - 4*c^2*y^6 - 5*c*y^7 - 4*c^4 - 6*y^8)/c^4
-1/4*(4*c^4*y^4 - 3*c^3*y^5 - 4*c^2*y^6 - 5*c*y^7 - 4*c^4 - 6*y^8)/c^4

Now we can collect all terms in y (this has the effect of sorting the terms into an easier to read order).

rattaylorysorted = rattaylory.collect(y); rattaylorysorted 
       
3/4*y^5/c + y^6/c^2 + 5/4*y^7/c^3 + 3/2*y^8/c^4 - y^4 + 1
3/4*y^5/c + y^6/c^2 + 5/4*y^7/c^3 + 3/2*y^8/c^4 - y^4 + 1

Alternatively, we could find a series expansion for rattaylory:  this has the nice feature that it sorts our terms into increasing rather that decreasing powers of $y$.

rattaylorseries=rattaylory.series(y,9);rattaylorseries; 
       
1 + (-1)*y^4 + (3/4/c)*y^5 + (c^(-2))*y^6 + (5/4/c^3)*y^7 +
(3/2/c^4)*y^8
1 + (-1)*y^4 + (3/4/c)*y^5 + (c^(-2))*y^6 + (5/4/c^3)*y^7 + (3/2/c^4)*y^8

Finally we substitute back $y=\sqrt{t}$ into each of the above, to see

that they both have the annoying feature that they reorder the terms again.

ratty=rattaylorysorted.substitute(y=1/sqrt(t));ratty 
       
-1/t^2 + 3/4/(c*t^(5/2)) + 1/(c^2*t^3) + 5/4/(c^3*t^(7/2)) +
3/2/(c^4*t^4) + 1
-1/t^2 + 3/4/(c*t^(5/2)) + 1/(c^2*t^3) + 5/4/(c^3*t^(7/2)) + 3/2/(c^4*t^4) + 1
rattaylorseries.substitute(y=1/sqrt(t)) 
       
-1/t^2 + 3/4/(c*t^(5/2)) + 1/(c^2*t^3) + 5/4/(c^3*t^(7/2)) +
3/2/(c^4*t^4) + 1
-1/t^2 + 3/4/(c*t^(5/2)) + 1/(c^2*t^3) + 5/4/(c^3*t^(7/2)) + 3/2/(c^4*t^4) + 1

However, by looking at the series expansion in $y$ and working by hand, we can see that

\[

\frac{rat(t)rat(t)}{rat(t+1)rat(t+1)} = 1 - \frac{1}{t^2} + \frac{3}{4ct^{5/2}} + O\left( \frac{1}{t^3} \right)

\]

 
       

Hence we can put together the exponential and rational parts and obtain, as $n \rightarrow \infty$

\[

\frac{p(n)p(n)}{p(n+1)p(n-1)} \sim \exp\left( \frac{c}{4t^{3/2}} \left( 1 + \frac{5}{16t^2} + O(1/t^4) \right) \right)  \left(  1 - \frac{1}{t^2} + \frac{3}{4ct^{5/2}} + O\left( \frac{1}{t^3} \right) \right)

\]

Then, since $\exp(a)>1+a$, provided

\[

\left(1+\frac{c}{4t^{3/2}} \right) \left( 1 - \frac{1}{t^2} \right)

\]

and the other terms are negligible, we will have the ratio is greater than 1.

def ratplot(t): return 1+ pi*sqrt(2/3)/4/t^(3/2)*(1-1/t^2) 
       
plot(ratplot,(1,100)) 
       
plot(ratplot, (10,100)) 
       

This is enough to show that the partition function is log-concave for $n$ sufficiently large.

It now remains to explain the alternating behaviour for small $n$.  To do this, we'll consider the contribution of the term $k=2$.

 
       

It now remains to explain the alternating behaviour for small $n$.  To do this, we'll consider the contribution of the term $k=2$.

The contribution of the term $k =2$, as we saw above, is a multiplicative factor
\[

m(n)=(  1+ \sqrt{2}(−1)^n \frac{(c \sqrt{t}−1)}{c\sqrt{t}−2)}\exp(−\frac{c\sqrt{t}}{2}) )
\]


For even $n$ $m(n)>1$: for odd $n$ $m(n)<1$.  Clearly, then for even $n$, this only makes $p(n)2/(p(n+1)(p(n−1))$ even greater.

However when $n$ is odd, then $m(n)^2/(m(n+1)m(n−1))<1$, and so we need to be more careful.

In this case, as usual, with $t=n−1/24$,
\[

\frac{m(n)^2}{m(n+1)m(n−1)} =

\frac{(1−\frac{1}{\sqrt{2}}\frac{c\sqrt{t}−2)}{c\sqrt{t}−1}

\exp(−\frac{c\sqrt{t}}{2}))^2 }{ (1+\frac{1}{\sqrt{2}}\frac{(c\sqrt{t+1}−2)}{(c\sqrt{t+1}−1)}\exp(−\frac{c\sqrt{t+1}}{2}))  

(1+ \frac{1}{\sqrt{2}}\frac{c\sqrt{t−1}−2}{c\sqrt{t−1}−1}\exp(−\frac{c\sqrt{t−1}}{2}))}

\]

 
       

Now, since each of the terms in the products is of the form $(1 + o(1))$, we have 

\[

\frac{m(n)^2}{m(n+1)m(n−1)} \sim  1 - 4  \frac{1}{\sqrt{2}}\exp\left( -\frac{c\sqrt{t}}{2}   \right)

\]

 
       

Thus, if

\[

2\sqrt{2} \exp\left( -\frac{c\sqrt{t}}{2} \right)  \simeq \frac{c}{4t^{3/2}}

\]

def fn(t): return 2*sqrt(2)*exp(-pi*sqrt(2/3)*sqrt(t)/2)-pi*sqrt(2/3)/(4*t^(3/2)) 
       
fn(t) 
       
2*sqrt(2)*e^(-1/2*pi*sqrt(2/3)*sqrt(t)) - 1/4*pi*sqrt(2/3)/t^(3/2)
2*sqrt(2)*e^(-1/2*pi*sqrt(2/3)*sqrt(t)) - 1/4*pi*sqrt(2/3)/t^(3/2)
plot(fn,(10,30)) 
       

Bingo!

Even with this first approximation we see that the inequality should start to hold for odd $n$ somewhere around $n=25$.  Since we haven't taken into account the $1/t^2$ term, it seems likely that this will account for the remaining error.

def fn2(t): return (1-0.91*2*sqrt(2)*exp(-pi*sqrt(2/3)*sqrt(t)/2))*exp(pi*sqrt(2/3)/(4*t^(3/2)))*(1 - 1/t^2+3/(4*pi*sqrt(2/3)*t^(5/2))) 
       
fn2(t) 
       
-1/8*(-1.82000000000000*sqrt(2)*e^(-1/2*pi*sqrt(2/3)*sqrt(t)) +
1)*(8/t^2 - 9*sqrt(2/3)/(pi*t^(5/2)) - 8)*e^(1/4*pi*sqrt(2/3)/t^(3/2))
-1/8*(-1.82000000000000*sqrt(2)*e^(-1/2*pi*sqrt(2/3)*sqrt(t)) + 1)*(8/t^2 - 9*sqrt(2/3)/(pi*t^(5/2)) - 8)*e^(1/4*pi*sqrt(2/3)/t^(3/2))
plot(log(fn2(t)),25,30) 
       
def f3(t): return (pi*sqrt(2/3)*sqrt(t)-2)/(pi*sqrt(2/3)*sqrt(t)-1) 
       
plot(f3,(20,30))