To have more insight, let us compute explicitly the involved polynomials.
This is not needed, but the ambitious reader may observe the structure of the polynomial $P,Q,f$ and start the own path to a solution. Here and below, we omit the index $n$, considered to be fixed in the context.
The following sage code defines a function fPQ of a variable $n=1,2,3,\dots$ (in computational reach) which builds by definition and returns the polynomials f,P,Q:
def fPQ(n):
R0 = [+1, -1]
R = [R0 for counter in range(n)]
F.<z> = CyclotomicField(2*n)
S.<X> = PolynomialRing(F)
P = prod([X - sum([mu[k-1]*z^k for k in [1..n]])
for mu in cartesian_product(R) if prod(mu) == +1])
Q = prod([X - sum([mu[k-1]*z^k for k in [1..n]])
for mu in cartesian_product(R) if prod(mu) == -1])
f = (P+Q)/2 - X^(2^(n-1))
return f, P, Q
for n in [1..7]:
f, P, Q = fPQ(n)
print(f"n = {n}\nP = {P}\nQ = {Q}\nf = {f}\n")
for n in [1..7]:
f, P, Q = fPQ(n)
print(f"n = {n} :: f = {f}")
And we obtain from the prints in the last loop:
n = 1
P = X + 1
Q = X - 1
f = 0
n = 2
P = X^2 + 2z
Q = X^2 - 2z
f = 0
n = 3
P = X^4 - 8X
Q = X^4 + 8X
f = 0
n = 4
P = X^8 - 48z^2X^4 - 64
Q = X^8 + 48z^2X^4 - 64
f = -64
n = 5
P = X^16 + 384X^11 + 10240X^6 - 32768X
Q = X^16 - 384X^11 + 10240X^6 + 32768X
f = 10240*X^6
n = 6
P = X^32 + 3840z^3X^26 - 1712128X^20 - 17825792z^3X^14 + 6996099072X^8 + 8589934592z^3X^2
Q = X^32 - 3840z^3X^26 - 1712128X^20 + 17825792z^3X^14 + 6996099072X^8 - 8589934592z^3X^2
f = -1712128X^20 + 6996099072X^8
n = 7
P = X^64 - 46080X^57 + 345899008X^50 - 394096803840X^43 + 300517251088384X^36 + 341689456127901696X^29 - 38826236075002822656X^22 - 1189544776776623849472X^15 + 3201446845511100268544X^8 - 1180591620717411303424X
Q = X^64 + 46080X^57 + 345899008X^50 + 394096803840X^43 + 300517251088384X^36 - 341689456127901696X^29 - 38826236075002822656X^22 + 1189544776776623849472X^15 + 3201446845511100268544X^8 + 1180591620717411303424X
f = 345899008X^50 + 300517251088384X^36 - 38826236075002822656X^22 + 3201446845511100268544X^8
Based on the above data we can formulate the following guess:
Claim:
The polynomials $P, Q$ live in monomial degrees that differ from $2^{n-1}=\deg P=\deg Q$ by a multiple of $n$. (All other coefficients of $P,Q$ are zero.) In case of an even multiple of $n$ the corresponding coefficients in $P,Q$ coincide, for an odd multiple of $n$ they have opposite sign.
Proof of the claim:
Let $F=\Bbb Q(z)$ be the cyclotomic field of degree $2n$, where $z=\zeta_{2n}$ is (chosen to be) $\exp\frac{2\pi i}{2n}$, "the" standard choice of the primitive root $\zeta_{2n}$ of order $2n$ of the unit. The full list of the roots of order $2n$ is
$$
\pm 1\ ,\ \pm z\ ,\ \pm z^2\ ,\ \dots\ , \pm z^{n-1}\ .
$$
They form a group $(G,\cdot)$ with $2n$ elements, isomorphic to $(\Bbb Z/(2n), +)$.
The sums involved in the roots of $P,Q$ use the system of representatives
$S=\{z_1,z_2,\dots,z_n\}$, where
$$
z_1=z\ ,\ z_2=z^2\ ,\ \dots\ ,\ z_n=z^n=-1
$$ for $G$ modulo $\pm 1$. The
Galois structure of $F$ is abelian, and the
action of the Frobenius automorphism is
$$
\Phi:F\to F\ ,\qquad\Bbb Q(z)\overset \Phi\longrightarrow \Bbb Q(z)\ ,\qquad
z\to z^g\ ,
$$
for a suitable generator $g$ of $((\Bbb Z/2n)^\times,\cdot)$. In particular, $g$ is odd.
We apply the Frobenius substitution $\Phi$ on the representatives in $S$,
elementwise, and obtain an other system of representatives, $\Phi(S)$. (Because $\Phi(\pm 1)=\pm 1$.) The representatives differ each by a sign factor $\pm$, and we ask for the product of these factor differences. Well, it is
$$
\frac{\prod \Phi(S)}{\prod S}
=
\frac
{\Phi(z_1z_2\dots z_n)}
{z_1z_2\dots z_n}
=
\frac
{\Phi\left(z^{n(n+1)/2}\right)}
{z^{n(n+1)/2}}
=
\frac
{z^{gn(n+1)/2}}
{z^{n(n+1)/2}}
= z^{(g-1)n(n+1)/2}
\ .
$$
For an odd $n$ we have a power $(g-1)n(n+1)/2$ which is an integer, moreover a multiple of $2n$ (since $(g-1)$ and $(n+1)$ are both even), so the last expression is one.
For an even $n$ both $(g-1)/2$ and $(n+1)$ are odd, so we have minus one instead. But applying $\Phi$ twice, $\Phi^2(z)=(z^g)^g=z^{g^2}$, we obtain always a multiple of $2n$, because $(g^2-1)(n+1)/2$ is an even integer.
This implies that $P,Q,f$ have coefficients
- in the fixed field of $\Phi$ which is $\Bbb Q$ in case of an odd $n$,
- respectively in the fixed field of $\Phi^2$ which is $\Bbb Q(i)$ where $
i=\sqrt{-1}=z^{n/2}$ in case of an even $n$.
Because the roots of $P$, and $Q$ is permuted by the named first, respectively second power of $\Phi$.
(This explains in the case of $n=6$ the presence of $z^3=\zeta_{12}^3$ in some of the coefficients of $P,Q$.)
We have done some progress, so that we can formulate:
Refined claim: In case of an odd $n$ we have $P(X),Q(X)\in\Bbb Z[X]$,
$\deg P=\deg Q=2^{n-1}=:d$, and the coefficients of $P,Q$ live only in degrees of the shape $d-jn$ with $j=0,1,2,3,\dots$. For an even $j$ the coefficients of $P,Q$ in degree $d-jn$ coincide, for an odd $j$ they differ by a $(-1)$ factor (sign difference).
In case of an even $n$ we have $P(X),Q(X)\in\Bbb Z[i][X]$,
$\deg P=\deg Q=2^{n-1}=:d$, and the coefficients of $P,Q$ live only in degrees of the shape $d-jn$ with $j=0,1,2,3,\dots$. For an even $j$ the coefficients of $P,Q$ in degree $d-jn$ coincide and are in $\Bbb Z$, for an odd $j$ they differ by a $(-1)$ factor (sign difference) and are in $i\Bbb Z$.
Continuation of the proof: The information that the coefficients of $P,Q$ are algebraic integers follows from the fact that their roots are algebraic integers in the cyclotomic field $F$.
We attack now the property in the OP focus, the vanishing of coefficients in degrees not of the shape $d-jn$.
Some notations let us type simpler formulas. Let $I(P)$, respectively $I(Q)$ be the index set of all $\mu\in\{+1,-1\}^n$ with product of the components $\Pi \mu$ equal to $+1$, respectively $-1$.
I have a problem with the notation $z$ for the $\zeta_{2n}$, and for the indexed version of other units $z_1,z_2,\dots,z_n$, which are $z^1,z^2,\dots,z^n$. It is the same $z$ used as a letter in the notation.
Instead of the sum $\sum\mu_j z_j$ with $1\le j\le n$ i will use
either the compact notation $z(\mu)$, or an explicit scalar product notation, $\langle \mu_1,\mu_2,\dots,\mu_n\mid z_1,z_2,\dots,z_n\rangle$. (The notation $\mu\cdot z$ would be misleading. I may have used $w$ for the representatives, but this is not so good at other places.) So we compute with this long scalar product notation:
$$
\begin{aligned}
P(X)
&=
\prod_{\mu \in I(P)}
\Big(X-z(\mu)\Big)
=
\prod_{\mu \in I(P)}
\Big(X-\langle \mu_1,\mu_2,\dots,\mu_n\mid z_1,z_2,\dots,z_n\rangle\Big)
\ ,\\
Q(X)
&=
\prod_{\mu \in I(Q)}
\Big(X-z(\mu)\Big)
=
\prod_{\mu \in I(Q)}
\Big(X-\langle \mu_1,\mu_2,\dots,\mu_n\mid z_1,z_2,\dots,z_n\rangle\Big)
\ ,\\[2mm]
P(zX)
&=
\prod_{\mu \in I(P)}
\Big(zX-\langle \mu_1,\mu_2,\dots,\mu_n\mid z_1,z_2,\dots,z_n\rangle\Big)
\\
&=
z^{|I(P)|}
\prod_{\mu \in I(P)}
\Big(X-\left\langle \mu_1,\mu_2,\dots,\mu_n\mid \frac{z_1}z,\frac{z_2}z,\dots,\frac{z_n}z\right\rangle\Big)
\\
&=z^d
\prod_{\mu \in I(P)}
\Big(X-\langle \mu_1,\mu_2,\dots,\mu_n\mid 1,z_1,\dots,z_{n-1}
\rangle\Big)
\\
&=z^d
\prod_{\mu \in I(Q)}
\Big(X-\langle \mu_1,\mu_2,\dots,\mu_n\mid -1,z_1,\dots,z_{n-1}
\rangle\Big)
\\
&=z^d
\prod_{\mu \in I(Q)}
\Big(X-\langle \mu_1,\mu_2,\dots,\mu_n\mid z_1,\dots,z_{n-1},z_n
\rangle\Big)
\\
&=z^d\;Q(X)\ ,\\[2mm]
&\qquad\text{ and in a similar manner }
\\[2mm]
Q(zX) &= z^d\;P(X)\ .
\end{aligned}
$$
Let us use the relation $P(z^2X) = z^{2d}\;P(X)$ and see which are the consequences for the coefficient $c$ in degree $d-j$ of $P$.
In $P(z^2X)$ the coefficient in degree $d-j$ is $c\cdot(z^2)^{d-j}$, and on the other side, in $z^{2d}\;P(X)$ it is $c\cdot z^{2d}$. They differ by the factor $z^{2j}$, which is not one if $j$ is not a multiple of $n$, so in this case we have $c=0$.
So we expect non-zero coefficients in $P,Q$ exactly in degrees $d,d-n,d-2n,d-3n,\dots$ and so on.
What happens in degree $d'=d-n$? Denote the corresponding coefficient of $P$ by $c'$. From $P(zX)=z^d Q(X)$ we obtain in $P(zX)$ the the $d'$-term $c'(zX)^{d-n} =z^d(c'z^{-n})X^{d-n}=z^d(-c')X^{d-n}$, so the corresponding coefficient in $Q$ is $-c'$. This shows that in $P+Q$ there is no coefficient in degree $d'$. The degree of $f$ thus drops from this expectation and is at most $d''=d-2n$.
This concludes the part of the claim insuring $\deg f\le d''$.
As mentioned, the first (odd $n$) or the second power (even $n$) of the Frobenius automorphism $\Phi$ stabilizes $P,Q$, so we have integral coefficients, in $\Bbb Z$ in the first case, respectively in $\Bbb Z[i]$ in the second case for both $P,Q$.
Let $c''$ be the coefficient in degree $d''=d-2n$ in both $P$ and $Q$. They coincide by the same argument, looking at the terms in degree $d''$ in both $P,Q$. Note that $c''$ is an element in $\Bbb Z$, because it is a symmetric elementary polynomial (of even degree) in $\{z(\mu)\ :\ \mu\in I(P)\}$
and/or in the same time in $\{z(\mu)\ :\ \mu\in I(Q)\}$, and applying $\Phi$
we either map the first set in itself, or in the other set. This gives $\Phi(c'')=c''$, i.e. $c''\in \Bbb Z$. The same argument shows that we have in degree $d'=d-n$ - in a first case with an odd $n$ the relation $\Phi(c')=c'$, so $c'\in\Bbb Z$, in a second case with an even $n$ the relation $\Phi(c')=-c'$ and $\Phi^2(c')$, so $c'\in \Bbb Z[i]$ which is the fixed field of $\Phi^2$, and even more special $c'\in i\Bbb Z$ so that the complex conjugation brings $c'\to -c'$.
It remains to show $c''\ne 0$, and to look at the divisibility property: $c''$
is zero modulo $64$. These are claims of a different kind, compared to what we did before. The divisibility is an information on the $2$-valuation
$$\nu=\nu_2$$
of the coefficient $c''$. The best tool to deal with the situation is to use the
Newton polygon
of $P$. To see how to use it, here is an example that has nothing to do with our problem, but illustrates the way we trace back divisibility information on the coefficients by knowing it for the roots. Let us consider the polynomial $g$ in $\Bbb Z[X]$ defined as follows:
$$
\small
\begin{aligned}
g
&=(X+2)^2(X-4)^3(X-12)(X-22)(X+32)
\\
&=
X^8 - 10X^7 - 804X^6 + 15112X^5 - 71104X^4 - 32256X^3
\\
&\qquad\qquad + 729088X^2 - 329728X - 2162688\ ,
\\
&=
g_8X^8 +
g_7X^7 +
g_6X^6 +
g_5X^5 +
g_4X^4 +
g_3X^3 +
g_2X^2 +
g_1X +
g_0\ ,
\\[3mm]
\nu(g_0) &= \nu(2162688) = \nu\left(2^{16} \cdot 3 \cdot 11\right) = 16 \ , &P_0&=(0,16)\ , \\
\nu(g_1) &= \nu(329728) = \nu\left(2^{11} \cdot 7 \cdot 23\right) = 11 \ , &P_1&=(1,11)\ , \\
\nu(g_2) &= \nu(729088) = \nu\left(2^{13} \cdot 89\right) = 13 \ , &P_2&=(2,13)\ , \\
\nu(g_3) &= \nu(32256) = \nu\left(2^{9} \cdot 3^{2} \cdot 7\right) = 9 \ , &P_3&=(3,9)\ , \\
\nu(g_4) &= \nu(71104) = \nu\left(2^{6} \cdot 11 \cdot 101\right) = 6 \ , &P_4&=(4,6)\ , \\
\nu(g_5) &= \nu(15112) = \nu\left(2^{3} \cdot 1889\right) = 3 \ , &P_5&=(5,3)\ , \\
\nu(g_6) &= \nu(804) = \nu\left(2^{2} \cdot 3 \cdot 67\right) = 2 \ , &P_6&=(6,2)\ , \\
\nu(g_7) &= \nu(10) = \nu\left(2 \cdot 5\right) = 1 \ , &P_7&=(7,1)\ , \\
\nu(g_8) &= \nu(1) = 0 \ , &P_8&=(8,0) \ .
\end{aligned}
$$
So we draw the following polygon:

The slopes and the widths match the roots $2$-adic information:
- the slope $-5$ has width one, there is exactly one root, $-32=-2^5$, with valuation five,
- the slope $-2$ has width four, there are exactly four roots (counting multiplicities), $4=2^2$ three times and $12=2^2\cdot 3$, with valuation two,
- the slope $-1$ has width four, there are exactly three roots (counting multiplicities), $-2=-2^1$ two times and $22=2^1\cdot 11$, with valuation one.
In our case, we will have for $P$ also a Newton polygon. It is monic, so the right most point is $(d,0)$, on the bottom of the first quadrant.
So the Newton polygon of $P$ (which is the same for $Q$, since coefficients differ up to units $\pm 1$) is "above" the line with the last slope.
$P$ is "lacunary" (many coefficients vanish), and instead of plotting them with height $\infty$ we omit them. Such zero terms are also "at the beginning".
For instance, in the case $n=7$ we have for $z=\zeta_{14}$ the minimal polynomial
(cyclotomic polynomial) $Z^6 - Z^5 + Z^4 - Z^3 + Z^2 - Z + 1$.
So
$z-z^2+z^3-z^4+z^5-z^6+z^7$ is zero, which makes $\mu_0(z)=0$ for $\mu_0
=(1,-1,1,-1,1,-1,1)$. One of $\pm \mu_0$ is in $I(P)$, one in $I(Q)$, so $P,Q$ have the factor $X$.
We neglect such factors equal to $X$ to some power when building the Newton polynomial.
So the information on the "last slope" can be obtained (estimated) from the valuations (or their estimations) of the roots of $P$.
So let us consider these roots.
We show that $2$ divides each $\mu(z)$. It is enough to restrict to $\mu\in I(P)$.
Let $\mu_0=(1,1,\dots,1)\in I(P)$.
For any $\mu\in I(P)$ the difference $\mu(z)-\mu_0(z)$ is a multiple of two.
So we will show that $2$ divides
$$\mu_0(z)=z+z^2+\dots+z^n
=z\cdot\frac {1-z^n}{1-z}
\sim
\frac {1-z^n}{1-z}
=\frac {1-(-1)}{1-z}
=\frac 2{1-z}\ .
$$
Claim A: Assume that $n$ is not a power of two. Then $(1-z)$ is a unit.
Proof of the claim:
Below, till we address Claim B, we assume that $n$ is not a power of two.
Let us write $n=t\cdot (2k+1)$ with a two-power $t$, multiplied
by the odd factor $2k+1$, $k>0$. Let $w$ be for short $w=z^t$. Then we have:
$$
(1-w)\ \cdot\ \color{blue}{w\frac{1-w^{2k}}{1-w^2}}
=\frac{w(1-w^{2k})}{1+w}
=\frac{w-w^{2k+1}}{1+w}
=\frac{w-z^n}{1+w}
=\frac{w-(-1)}{1+w}
=1\ .
$$
The blue element is $w(1+w^2+w^4+\dots)$, a cyclotomic integer. So $1-w=1-z^t=(1-z)(\dots)$ is a unit. So $1-z$ is a unit.
$\square$
So $\mu_0(z)=2\cdot (1-z)^{-1}$ is divisible by $2$ in the ring of cyclotomic integers $\mathcal O_F$ of the cyclotomic field $F=\Bbb Q(z)$ of order $2n$.
This shows that the slopes of the Newton polygon of $P$ lives above the
line through $(d,0)$ with slope $-1$.
In particular,
- in degree $d'=d-n$ we have a coefficient $c'$ divisible by $2^n$,
- in degree $d''=d-2n$ we have a coefficient $c''$ divisible by $2^{2n}$,
and so on.
This is enough to conclude the wanted divisibility by $64$, and gives bonus information.
Claim B: Assume that $n=2^r$ is a power of two.
Then $(1-z)$ has the valuation $\displaystyle\nu(1-z)=\frac 1{2^r}$.
Proof of Claim B: The primitive root $z$ has order $2n$, so
it is a root of the cyclotomic polynomial of order $2n$, which is
$X^n+1$, a polynomial of degree $\varphi(2n)=2^r=n$. So we can compute the norm of $z-1$ (associated to $1-z$), knowing that
its minimal polynomial is $(X+1)^n+1$. Its free part is $1^n+1=2$. This is the sum of the valuations of all $n$ conjugates of $z$, which are all equal in
valuation.
From $n\nu(1-z)=\nu(2)=1$ we obtain $\nu(1-z)=1/n$ as claimed.
$\square$
In this case the minimal slope
"slightly goes down" (in modulus), but the information is enough to conclude the wanted divisibility.
For instance, for $n=4$ we have the slope $-3/4$, for $n=8$ we have $-7/8$, for $n=16$ we have $-15/16$ and so on.
The case $n=8$ is plotted below.
Let us illustrate this in a table in some cases with small $n$.
Here $d=2^{n-1}$, $d'=d-n$ and $d''=d-2n$ are the "first" degrees where we expect non-zero coefficients, and these are $1,c,c'$, listed below, together with their $2$-valuations $\nu(c'),\nu(c'')$.
For the OP question we need only $\nu(c'')\ge 6$,
but we even have
- $\nu(c'')\ge 2n$ in case $n$ is not a power of two, respectively
- $\nu(c'')\ge \left(1-\frac 1n\right)\cdot 2n=2(n-1)$ when $n$ is a $2$-power.
(Thanks to hbghlyj for the comment, the above it is edited and corrected now.)
$$
\small
\begin{array}{|r|r|r|r|r|r|r|r|}
\hline
n & d & d' & c' & \nu(c') & d'' & c'' & \nu(c'') \\\hline
3 & 4 & 1 & -8 & 3 & -2 & 0 & \infty \\\hline
4 & 8 & 4 & -48i & 4 & 0 & -64 & 6 \\\hline
5 & 16 & 11 & 384 & 7 & 6 & 10240 & 11 \\\hline
6 & 32 & 26 & 3840i & 8 & 20 & -1712128 & 13 \\\hline
7 & 64 & 57 & -46080 & 10 & 50 & 345899008 & 17 \\\hline
8 & 128 & 120 & -645120i & 11 & 112 & -86249439232 & 17 \\\hline
9 & 256 & 247 & 10321920 & 15 & 238 & 26408252342272 & 23 \\\hline
10 & 512 & 502 & 185794560i & 16 & 492 & -9813054673911808 & 25 \\\hline
11 & 1024 & 1013 & -3715891200 & 18 & 1002 & 4368756701919182848 & 29 \\\hline
12 & 2048 & 2036 & -81749606400i & 19 & 2024 & -2301705013494443671552 & 30 \\\hline
13 & 4096 & 4083 & 1961990553600 & 22 & 4070 & 1419008274778271286034432 & 35 \\\hline
14 & 8192 & 8178 & 51011754393600i & 23 & 8164 & -1013348982986729992159756288 & 37 \\\hline
15 & 16384 & 16369 & -1428329123020800 & 25 & 16354 & 830669684968147584741265113088 & 41 \\\hline
\end{array}
$$
I have no way to exhibit an argument, an idea to attack $c''\ne 0$, but i think this is not so important in the given context. (If it is, then let me know - also the motivation for this non-vanishing, so that a partial progress can be done in the good direction.)
I will finish with the picture of two Newton polygons for the polynomial $P$, and compare it with the computed valuations of its roots.
n = 8 Newton polygon plot:
For the value $n=8$ first, a two-power. The degree of $P$ is $d=2^{n-1}=2^7=128$, so the plot will have a chance of visibility.
This polynomial has the following coefficients, that are not ignored in the plot of the vertices of the Newton polygon:
$$
\small
\begin{array}{|r|r|r|r|r|r|r|r|}
\hline
\text{degree} & \text{coefficient }c & \nu(c) \\\hline
0 & 36220203446805248187373993663980270398734336 & 112\\\hline
8 & 35827446505073955166559408040216869783539286016i & 109\\\hline
16 & -38131836494277129848049667299302474170834616320 & 101\\\hline
24 & -6246341177070973482312350675389078531781689344i & 95\\\hline
32 & -137038837715960541148251713274135582640963584 & 86\\\hline
40 & -161205944795104966960364917291191883530240i & 81\\\hline
48 & -18486565864433401437058515697614222000128 & 73\\\hline
56 & -38028552829022542432817495343762505728i & 67\\\hline
64 & -1261645642172270042112757092843520 & 57\\\hline
72 & 31499645556987584648598515613696i & 53\\\hline
80 & -5152933192025219317681356800 & 45\\\hline
88 & 149888122836514731196416i & 39\\\hline
96 & 37667807777003667456 & 30\\\hline
104 & 2900646999097344i & 25\\\hline
112 & -86249439232 & 17\\\hline
120 & -645120i & 11\\\hline
128 & 1 & 0\\\hline
\end{array}
$$

There is only one slope, $-7/8$, and we expect that all roots of $P$ have valuation $7/8$. This is the case, as it can be easily checked by code:
sage: n = 8
sage: f, P, Q = fPQ(n)
sage: F = P.parent().base_ring()
sage: set(valuation( ZZ(root.norm()) , 2) / F.degree() for root in P.roots(multiplicities=False))
{7/8}
n = 7 Newton polygon plot:
We have the following non-zero coefficients with corresponding valuations in $P$:
$$
\small
\begin{array}{|r|r|r|}
\hline
\text{degree} & \text{coefficient }c & \nu(c) \\\hline
1 & -1180591620717411303424 & 70\\\hline
8 & 3201446845511100268544 & 56\\\hline
15 & -1189544776776623849472 & 54\\\hline
22 & -38826236075002822656 & 47\\\hline
29 & 341689456127901696 & 38\\\hline
36 & 300517251088384 & 29\\\hline
43 & -394096803840 & 25\\\hline
50 & 345899008 & 17\\\hline
57 & -46080 & 10\\\hline
64 & 1 & 0\\\hline
\end{array}
$$

sage: n = 7
sage: f, P, Q = fPQ(n)
sage: ZT.<T> = ZZ[]
sage: PT = ZT( sum(ZZ(P[k]) * T^k for k in [0..P.degree()]) / T )
sage: set(PT.newton_slopes(2))
{1, 2}