Fourier coefficient of $\frac{1}{5 + 4\cos(x)}$ wrong

Question

Bug introduced in 12.3.1 or earlier and persisting through 13.2.0 or later.

---

I've been trying to use Mathematica to compute the Fourier coefficient of

$$ f(x) = \frac{1}{5 + 4\cos(x)}, $$ which is a $2\pi$ periodic function. Inputting

FullSimplify[
 FourierCoefficient[1/(5 + 4 Cos[x]), x, n, 
  FourierParameters -> {1, 1}, 
  Assumptions -> {x \[Element] Reals, n \[Element] Integers}]]

I get as an output

-(1/3) (-2)^Abs[n]

This is clearly wrong as the Fourier coefficient of $f$, it being continuous and periodic, should decay with $n$, which the output from Mathematica does not give. Not only that, the Fourier coefficients should be bounded by 1 via $$ \left|\frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)e^{inx} dx\right| \leq \sup_{x\in [-\pi,\pi)}|f(x)|=1 \quad \forall n\in\mathbb{Z}. $$ Any insight into why this is happening would be appreciated.

Edit

Thank you all for your answers. To summarise, it seems that Mathematica evaluates the above correctly if instead you use either of the following:

FullSimplify[FourierCosCoefficient[1/(5 + 4 Cos[x]), x, n], 
 Assumptions -> {n \[Element] PositiveIntegers}]

OR

Simplify[Integrate[Cos[n*x]/(5 + 4 Cos[x]), {x, -Pi, Pi}], 
 Assumptions -> {n \[Element] PositiveIntegers}]

OR similar variations of the above.

The bug is strange, and hopefully it shall be fixed soon.

Answer 1

There is a workaround by making use of the definition:

FullSimplify[Integrate[1/(5 + 4 Cos[x])*Cos[n*x], {x, -Pi, Pi}, 
Assumptions -> n \[Element] PositiveIntegers]/(2*Pi),  
Assumptions -> n \[Element] PositiveIntegers]

1/3 (-(1/2))^n

It's clear that the sine-coefficients equal zero.

Answer 2

Here's an algebraic way to get both series, the divergent one and the convergent one, which are each geometric series. The equation in SolveAlways assumes an even function of x that can be expanded in a geometric series in $z=e^{ix}$.

fz = 1/(5 + 4 Cos[x]) /. x -> Log[z]/I
sol = SolveAlways[1/(a/z + b) + c + 1/(a z + b) - fz == 0 && a != 0, z]
(*
  1/(5 + (2 (1 + z^2))/z)
  {{a -> -6, b -> -3, c -> 1/3}, {a -> 3/2, b -> 3, c -> -(1/3)}}
*)

We develop the positive and negative powers of $z=e^{ix}$ separately:

cplus = Hold[SeriesCoefficient[1/(a z + b), {z, 0, n}]] /. sol // 
  ReleaseHold

cminus = 
 Hold[SeriesCoefficient[1/(a/ z + b), {z, Infinity, n}]] /. sol // 
   ReleaseHold // ReplaceAll[n -> -n]

cn = cplus + (Piecewise[{{c, n == 0}}] /. sol) + cminus // 
   Map[PiecewiseExpand[#, Assumptions -> n \[Element] Integers] &] //
    Simplify

Simplify did not combine ±n into Abs[n] — oh, well. Check series:

Sum[# Exp[I n x], {n, -Infinity, Infinity}] & /@
  cn /.
 Sum -> Inactive[Sum]
% // Together // ExpandAll // Together // ExpToTrig // Simplify