On the Unreasonable Difficulty of a Particular Cosine Integral

Consider the following integral, which arises when attempting to normalize the Blinn–Phong BRDF:

\[ I(C,m) := \int_0^{2\pi} \frac{d x}{( C + \cos{x} )^m} \]

This integral certainly doesn't look very intimidating. We can even make some nice assumptions about it—say $C>1$ and maybe $m \in \mathbb{Z}^+$. However, I defy you to try solving it!

Essentially everything you try will fail, no matter how many weeks (I definitely spent more than 1) you spend trying everything. You can try $u$-substitution, Weierstrass substitution, Euler substitution, trigonometric substitution, iterated partial fractions, series expansion, complex contour integration, integration by parts, Feynman integration, reduction, and of course various computer algebra systems—and any sequence of these together you can think of.

This integral, in fact, seems almost certain to be unsolvable in terms of elementary functions. I did find two different solutions in terms of the hypergeometric function, just in case you came here looking for results instead of my rambling philosophy. The better of the two is:

\[ I(C,m) = \frac{2\pi}{C^m} \cdot {_2}F_1\!\left( \frac{m}{2},\frac{m+1}{2} ; 1 ; \frac{1}{C^2} \right) \]

import math
import scipy.special

def I( C, m ):
	return 2*math.pi * C**-m * scipy.special.hyp2f1( 0.5*m,0.5*(m+1), 1, C**-2 )

However, even this requires a very subtle proof: not only does it require a fair amount of Calculus and Analysis sophistication, straying at all from the one path will lead down garden paths of hopeless messes. No CAS can get a result as good, as of this writing^[1].

I'll show the "right" way through the forest below, but that's not really the point of all this. What I want to focus on is the lurking, interconnected, monsters of difficulty, interrelated in ways unimaginable, visible only when we stray from the narrow path to success, that somehow conspire to make us "lose". What are these monsters? Why are they here? Where do they come from? Nobody knows, but by stepping off the path and becoming "lost", as it were, we can catch glimpses of their true forms, eluding while yet defining the shoreline of our knowledge.

Difficulty, Magically Appearing

The mathematics humanity has invented—perhaps that it is even possible to invent—has carved out an island of sanity in an ocean of madness. We think of mathematics as being a powerful tool—and it surely is—but few ever dip their feet into the great waters beyond. Many have forgotten it is there.

When we probe the present frontiers of mathematics, we find ourselves somehow rebuked. Often, we fail entirely. Sometimes, we escape with an answer in terms of a polylogarithm, a Bessel function, a convergent sum, a solution to a special case. Very rarely, we prise from the depths a result in terms of something we can understand (such results are celebrated with the holiest of epithets: "elegant").

But for all this, we do not understand. We have no idea why some problems are solvable and others aren't. We just have to try them and see.

In fact, when we fail to solve a problem, what is peculiar—and in fact astonishing when properly apprehended!—is that we seem to also fail every other way we try too!

In the following, I will take the above integral as an example to show what I mean—sketching out for various methods how absolutely arbitrarily we seem to be blocked from an elementary answer. Naturally, no procedure will give an answer if it doesn't exist—but it is not at all clear that an answer unreachable one way should be unreachable in every way. What possible connection does "not being able to move a limit inside a particular sum" have to do with "having an extra term in the denominator" or "lacking a derivative for a $u$-substitution"? On the flip-side, problems that can be solved one way can usually be solved in many ways. How do these completely different structures seemingly "know" that the problem is "supposed" to be impossible, blocking you at every turn—or else that the problem is supposed to be "easy", granting answers by any method?

Surely this is not a coincidence! I believe that these disparate connections, ineffably linking vastly different areas of mathematical analysis, are really reflections of the same underlying knot—a bramble of ineffable petulance that stymies us whenever we approach. Perhaps the immense abstract formalism we have built so keenly in mathematics itself actually curves obliquely around this sort of complexity, so that we are forever seeing the same demon through a glass darkly, just from different angles. If this seems unlikely, or mayhap overly poetic, how else can we account for the sheer improbability of seeing difficulty reflected in so many completely different ways?

To me, this is the most interesting reason to be a mathematician. What good is it to discover new structure when we have no idea how these deep connections, the foundation, arise? When we slam headfirst into an insurmountable obstacle, we wonder why this, right here and now is the wall. Why this? Why now? This is why being a mathematician is awesome! We have the tools to seek out the higher-dimensional manifold of allowable and unallowable problems, to touch the boundary between civilization and madness. We can probe the curiously sharp line between friendly polynomials, trigonometry, and the like, against the holomorphic jungle, where prime numbers hunt zeros in the Riemann zeta function, exponents twist space in circles, and some answers can be glimpsed in reflections, but never reached: a mirage bespeaking a truth that never was.

Maybe the deep structure is an illusion. Maybe every problem that can be solved has been rendered easy in our mathematics, with mature tactics in every field, while every problem that seems impossible is actually truly impossible. If we can't find an elementary answer, it's because it doesn't exist. The frontier isn't a chaotic ocean after all, but the fences built around the truly impossible holes in reality . . . but I don't believe in such a worldview, and I don't think many other mathematicians do either. There are so many open questions—and every. now. and then, the frontier is pushed back, fundamentally, our shoreline expanding.

I don't have the answers to these questions, and I don't have much to offer besides this pseudo-philosophical sense-poetry. I don't think anyone does. I don't think very many people are even working on such problems which, in my view, is a shame. What I present now is a path skirting along the frontier, playing in the surf before the endless ocean of our ignorance. It is for you to decide what it means.

Stuff That Doesn't Work

$u$-Substitution and Integration by Parts

The obvious starting point is to take a $u$-substitution. With any such substitution, the challenge is to pick a useful one—and in this case, the first time we see our great wall manifest, this turns out to be impossible.

$u:=\cos x$ or $u:=C\!+\!\cos x$ is immediately problematic, since $d u/d x = -\sin x$ and there's no $\sin x$ anywhere. One could try to use the pythagorean identity to get it, in which case, we get (showing for the latter substitution, since it seems more useful to get $C$ outside of the exponent):

\[ u:=C\!+\!\cos x ,\hspace{1cm} \frac{d u}{d x} = -\sin x, \hspace{1cm} \sin x = \sqrt{1-\cos^2 x} = \sqrt{1-(u-C)^2} \] \[ 2\int_0^\pi \frac{d x}{( C + \cos{x} )^m} = 2 \int_{C+1}^{C-1} \frac{d u}{u^m (-\sqrt{1-(u-C)^2})} = \int_{C-1}^{C+1} \frac{2~d u}{u^m \sqrt{1-(u-C)^2}} \]

You might think you can integrate this by parts, but you can't. The obvious part $u^{-m}$ will immediately make you integrate arcsine in the other part, while choosing the square-root just makes the next integral worse^[2].

If either term had been absent, the integral would be easy, as would the square-root being just a little different. The complexity demon has shown itself.

Weierstrass $t$-Substitution

Maybe we should be more aggressive in removing the trigonometry. The Weierstrass substitution does this:

\[ t := \tan(x/2) ,\hspace{1cm} \cos x = \frac{1-t^2}{1+t^2} ,\hspace{1cm} \frac{d x}{d t} = \frac{2}{1+t^2} \] \begin{align*} 2\int_0^\pi \frac{d x}{( C + \cos{x} )^m} &= 2 \int_0^1 \frac{1}{\left( C + \frac{1-t^2}{1+t^2} \right)^m} \frac{2}{1+t^2} d t \\ &= \int_0^1 \frac{4 ~(1+t^2)^{m-1} ~d t}{ \left( C (1+t^2) + 1-t^2 \right)^m } = \int_0^1 \frac{4 ~(1+t^2)^{m-1} ~d t}{ \left(~ (C\!-\!1)t^2 + (C\!+\!1) ~\right)^m } \end{align*}

It's not clear where we should go from here. Probably back to the drawing board. Perhaps most remarkably, the complexity demon has made that last constant term in the denominator $(C\!+\!1)$. If only it had been $(C\!-\!1)$ instead, then the integral would be trivial!

Feynman Integration

Feynman Integration: How-To

As a basic explanation^[3], let's say we want to know the following integral:

\[ \int_a^b f(x) ~ d x \]

First, we define some $g(x,t)$ such that for some $t=t_f$ we have $g(x,t_f)=f(x)$. We also define a "generalized" integral in terms of $g$ instead of $f$:

\[ G(t) := \int_a^b g(x,t) ~ d x \]

Notice that if we can get $G(t_f)$ we solve our problem!

We start with computing $G(t_0)=G_0$ at any convenient $t=t_0$ that makes the integral easy.

Then, we differentiate $G(t)$. The idea is that this differentiation makes the integrand simpler, allowing us to solve it:

\[ \frac{d}{d t} G(t) = \int_a^b \left( \frac{\partial}{\partial t} g(x,t) \right) d x = h(t) \]

Finally, we get our answer by integrating this out to our desired $t=t_f$. This is just the fundamental theorem of calculus:

\[ G(t_f) = G_0 + \int_{t_0}^{t_f} h(t) ~ dt \]

The best demonstration I've seen is for $f(x):=(x^2-1)/\ln(x)$. Many methods will prove difficult, but the Feynman integration generalization $g(x,t):=(x^t-1)/\ln(x)$ makes the solution simple.

Since it's a definite integral^[4], we can try Feynman integration. You basically "generalize" the integrand to be a function of some new parameter $t$ (that is, $f(x)$ becomes $g(x,t)$ for some particular value $t_f$ of $t$). Then you differentiate with respect to $t$, solve the now easier integral, and integrate with respect to $t$ to get back where you started, but with the original integral solved. It's simpler than it sounds or is usually explained (see sidebar).

It's a neat trick, to be sure. The core of the process is trying to find an appropriate generalization. You're looking for an integrand with a parameter that reproduces the original integral but can also be set to make integration easier—while also ensuring that the whole thing also differentiates to an easy integrand. Finding such an generalization relies on creativity, although various heuristics can of course help.

I tried several generalizations by hand, but when these failed for one reason or another, I resorted to computer search. In fact, you can think up many generalizations, none of which will work:

\begin{alignat*}{2} & g(x,t) := \frac{1-t\sin(x)}{(C+\cos(x))^m} && g(x,t) := \frac{1}{(t+\cos(x))^m}\\ & g(x,t) := \frac{1}{(C+t\cos(x))^m} && g(x,t) := \frac{1}{(C+\cos(x)^t)^m}\\ & g(x,t) := \frac{1}{(C+\cos(x)+t\sin(x))^m}\hspace{1cm} && g(x,t) := \frac{1}{(C+\cos(x)+\sin(x)^t)^m}\\ & g(x,t) := \frac{1}{(C+\cos(t x))^m} && g(x,t) := \frac{1}{(C+\cos(x^t))^m}\\ & g(x,t) := \frac{1}{(C+\cos(t^x))^m} && g(x,t) := (C+\cos(x))^t\\ & g(x,t) := \frac{e^{-t\sin(x)}}{(C+\cos(x))^m} && g(x,t) := \frac{\sin^t(x)}{(C+\cos(x))^m}\\ & g(x,t) := \frac{\cos^t(x)}{(C+\cos(x))^m} && \cdots \end{alignat*}

You will always find that you can't get $f(x)$, or you can't find a parameter value that will make integration easy, or the derivative won't be integrable. It's always something that won't work.

And that is the key here—nothing will work. The ineffable difficulty of this problem has manifested itself in our inability to find such a function! Nevermind that we're going to be differentiating it along a different dimension first. Somehow, we are blocked from finding such a function at all!

Getting an Actual Answer

Deriving the solution given above is relatively tricky, but should be understandable for anyone with a solid grasp of Calculus and basic Analysis (if you've made it this far . . .).

The first thing we have to do is factor out $C$ so that we can expand the integrand in terms of a binomial series:

\begin{align*} I(C,m) := \int_0^{2\pi} \frac{d x}{( C + \cos{x} )^m} &= \int_0^{2\pi} C^{-m} \left( 1 + \frac{1}{C}\cos(x) \right)^{-m} d x\\ &= C^{-m} \int_0^{2\pi} \left[ \sum_{k=0}^{\infty} \binom{-m}{k} \left( \frac{1}{C}\cos(x) \right)^k \right] d x \end{align*}

Since $C>1$, then $|\cos(x)/C|<1$ and the series is absolutely convergent. Therefore^[5], we can bring the integral inside the summation:

\[ I(C,m) = C^{-m} \sum_{k=0}^{\infty} \left[ \binom{-m}{k} ~ C^{-k} \int_0^{2\pi} \cos^k(x) ~ d x \right] \]

Now, the integral of a power of cosine is a classic problem for demonstrating integration by parts and expression as a reduction. For $k>0$, there is a reasonably well-known result that (less-well-known) comes from integration by parts^[6]:

\[ \int \cos^k(a x) ~ d x = \frac{\sin(a x) \cos^{k-1}(a x)}{a k} + \frac{k-1}{k} \int \cos^{k-2}(a x) ~ d x \text{, }\hspace{1cm}k>0 \]

For our purposes, we don't need the factor on $x$ and we're interested in the definite integral. This makes the recurrence simplify significantly:

\[ \int_0^{2\pi} \cos^k(x) ~ d x = \frac{k-1}{k} \int_0^{2\pi} \cos^{k-2}(x) ~ d x \text{, }\hspace{1cm}k>0 \]

We can't use a recurrence directly; we have to flatten it. The base case $k=0$ is just the integral of $1$ around the circle, which is $2\pi$. For $k$ odd, the result is always zero since half the integral cancels half the rest. For $k$ even, the important bit is the factors of $(k-1)/k$ that pile up in front of the base case. A table helps to make the patterns clearer:

$k$	$0$	$2$	$4$	$6$	$8$	$1$	$3$	$5$	$7$	$9$
$\int_0^{2\pi} \cos^k(x) ~ d x$	$2\pi$	$\frac{1}{2} \cdot 2\pi$	$\frac{3}{4} \cdot \frac{1}{2} \cdot 2\pi$	$\frac{5}{6} \cdot \frac{3}{4} \cdot \frac{1}{2} \cdot 2\pi$	$\frac{7}{8} \cdot \frac{5}{6} \cdot \frac{3}{4} \cdot \frac{1}{2} \cdot 2\pi$	$0$	$0$	$0$	$0$	$0$

For $k≥0$ even, we can summarize the result as:

\[ \int_0^{2\pi} \cos^k(x) ~ d x = 2\pi \prod_{j=1}^{k/2} \frac{2 j - 1}{2 j} \text{, }\hspace{1cm}k≥0\text{ even} \]

We substitute this result back into the series. Again, if $k$ is even, then the substitution should be the product above, while if $k$ is odd, the substitution should be zero. We can accomplish this by skipping the "odd" terms, which we can do by just doubling $k$ everywhere:

\[ I(C,m) = \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \binom{-m}{2 k} ~ C^{-2 k} \prod_{j=1}^k \frac{2 j - 1}{2 j} \right] \]

Now, we rearrange the series slightly. Note that in the following, $\Gamma(z)$ is the gamma function, an (the) extension of the (shifted) factorial function to the complex plane. The notation $(z)_n=\Gamma(z\!+\!n)/\Gamma(z)$ is the Pochhammer symbol, standard in advanced series manipulations.

We make the binomial coefficient positive by applying an identity. Using the gamma function, we also expand out the product^[7] and the binomial:

\begin{align*} \newcommand{\G}[1]{\Gamma\!\left(#1\right)} I(C,m) &= \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \binom{ m\!+\!2k\!-\!1 }{ 2k } \cdot C^{-2k} \cdot \prod_{j=1}^k \frac{ j \!-\! \tfrac{1}{2} }{j} \right] \\ &= \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \frac{ (m\!+\!2k\!-\!1)! }{ (m\!+\!2k\!-\!1\!-\!2k)! ~(2k)! } \cdot C^{-2k} \cdot \frac{ \G{k\!+\!\tfrac{1}{2}} ~/~ \G{\tfrac{1}{2}} }{ \G{k\!+\!1} } \right] \\ &= \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \frac{ \G{m\!+\!2k} }{ \G{m} ~\G{2k\!+\!1} } \cdot C^{-2 k} \cdot \frac{ \G{k\!+\!\tfrac{1}{2} } ~/ \sqrt{\pi} }{ \G{k\!+\!1} } \right] \end{align*}

By now we have a proper disaster! We can't sum or simplify, but we can rewrite the series in the form of the hypergeometric function, which will allow us to efficiently compute it.

To do this, we rearrange^[8] and then apply the Legendre duplication formula three times to split the arguments of certain gamma functions. When the dust settles, we rewrite them in terms of Pochhammer symbols and apply the definition of the hypergeometric function:

\begin{align*} I(C,m) &= \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \frac{ \G{m\!+\!2k} }{ \G{m} } \cdot \frac{ \G{k\!+\!\tfrac{1}{2}} ~\G{1} }{ \G{2k\!+\!1} } \cdot \frac{ C^{-2k} }{ \sqrt{\pi} ~\G{k\!+\!1} } \right] \\ & \begin{split} = \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \frac { \G{ \tfrac{m}{2} \!+\! k } ~\G{ \tfrac{m}{2} \!+\! k \!+\! \tfrac{1}{2} } ~/~ ( 2^{\cancel{1-m} \cancel{-2k} } \cancel{\sqrt{\pi}} ) } { \G{ \tfrac{m}{2} } ~\G{ \tfrac{m}{2} \!+\! \tfrac{1}{2} } ~/~ ( \cancel{2^{1-m}} \cancel{\sqrt{\pi}} ) } \cdot \right.\hspace{1cm} \\ \left. \frac { \cancel{\G{ k \!+\! \tfrac{1}{2} }} ~\G{1} } { \cancel{\G{ k \!+\! \tfrac{1}{2} }} ~\G{ k \!+\! 1 } ~/~ ( \cancel{2^{-2k}} \cancel{\sqrt{\pi}} ) } \cdot \frac{ C^{-2k} }{ \cancel{\sqrt{\pi}} ~\G{k\!+\!1} } \right] \end{split} \\ &= \frac{2\pi}{C^m} \sum_{k=0}^{\infty} \left[ \frac { \left(\tfrac{m}{2}\right)_k ~\left(\tfrac{m}{2}+\tfrac{1}{2}\right)_k } { (1)_k } \cdot \frac{ (1/C^2)^k }{ k! } \right] \\ &= \frac{2\pi}{C^m} \cdot {_2}F_1\!\left( \frac{m}{2}, \frac{m+1}{2} ; 1 ; \frac{1}{C^2} \right) \end{align*}

This is the result given above. I derived all this, and to my knowledge, this constitutes the first-ever solution of this problem^[9].

A professional mathematician with a background in analysis would likely have no trouble following this^[10]. However, each step—and particularly the first—seemingly must be followed.

Notice that even though we got a solution, we were still blocked from getting an answer in closed-form. Just because we "succeeded", we still haven't escaped non-elementary-ness. The great chaotic ocean still blocks us from understanding this problem in ways ineffable, even if we escaped back to dry land with a formula we can use.

Another Solution: Contour Integration

Another approach is complex contour integration, by making the substitution $z:=\exp(i x)$, with the following implications:

\[ z:=e^{i x}, \hspace{1cm} \frac{d z}{d x}=i~z, \hspace{1cm} \cos(x)=\frac{1}{2}\left( z + \frac{1}{z} \right) \]

Our integral then becomes a contour integral over the unit circle in the complex plane, and we can get to hacking, pulling out constants and factoring the roots towards a solution by residues, since a $u$-substitution won't work:

\begin{align*} I(C,m) := \int_0^{2\pi} \frac{d x}{( C + \cos{x} )^m} &= \oint_{|z|=1} \frac{d z}{ i ~z ~(~ C + (z+1/z)/2 ~)^m }\\ &= -i~2^m \oint_{|z|=1} \frac{ z^m ~d z }{ z ~(~ 2 C z + z^2 + 1 ~)^m }\\ &= -i~2^m \oint_{|z|=1} \frac{ z^{m-1} ~d z }{ (z-a)^m ~(z-b)^m } ,\hspace{1cm}\begin{matrix} a := -C - \sqrt{C^2-1} \\ b := -C + \sqrt{C^2-1} \\ \end{matrix} \end{align*}

The negative root $a$ is outside the contour (recall that we can assume $C>1$, so $-C$ minus anything will be less than $-1$), so we can ignore its residue. We are left, then, with computing the residue at the positive root $b$. That is:

\begin{align*} I(C,m) &= -i ~2^m \oint_{|z|=1} f(z)~d z = (-i ~2^m) (2 \pi i) \cdot \Res(f,b) ,\hspace{1cm} f(z) :=\frac{z^{m-1}}{(z-a)^m ~ (z-b)^m} \\ &= 2^{m+1}~\pi \cdot \Res(f,b) \end{align*}

If $m=0$, then $f(z)=z^{-1}$ and it is well-known that $\Res(f,b)=1$. For $m>0$, the pole is order $m$, and we need a slightly arcane formula:

\begin{align*} \Res(f,b) &= \frac{1}{(m-1)!} \lim_{z\rightarrow b}\frac{d^{m-1}}{d z^{m-1}} (z-b)^m f(z) ,\hspace{1cm} m>0 \\ &= \frac{1}{(m-1)!} \lim_{z\rightarrow b}\frac{d^{m-1}}{d z^{m-1}} \frac{ z^{m-1} }{ (z-a)^m }\\ &= \frac{1}{(m-1)!} \lim_{z\rightarrow b}\frac{d^{m-1}}{d z^{m-1}} \sum_{k=0}^\infty (-a)^k \binom{-m}{k} \frac{1}{z^{1+k}}\\ &= \frac{1}{(m-1)!} \lim_{z\rightarrow b} \sum_{k=0}^\infty (-a)^k \binom{-m}{k} (-k\!-\!1)(-k\!-\!2)(-k\!-\!3)\cdots(-k\!-\!m\!-\!1) z^{-k-m}\\ &= \frac{1}{(m-1)!} \lim_{z\rightarrow b} \sum_{k=0}^\infty (-a)^k \binom{-m}{k} (-1)^{m-1} \frac{(m+k-1)!}{k!} z^{-k-m}\\ &= \frac{ (-1)^{m-1} }{ b^m } \sum_{k=0}^\infty \binom{-m}{k} \frac{ (m+k-1)! }{ (m-1)! ~k! } \left(\frac{a}{b}\right)^k\\ &= \frac{ (-1)^{m-1} }{ b^m } {_2}F_1\!\left( m,m; 1; \frac{a}{b} \right) \end{align*}

Thus:

\begin{align*} I(C,m) &= \begin{cases} (2^{m+1}~\pi) ~1 & \text{ if } m=0 \\ (2^{m+1}~\pi) \frac{ (-1)^{m-1} }{ b^m } {_2}F_1\!\left( m,m; 1; \frac{a}{b} \right) & \text{ if } m>0 \end{cases} \\ &= \begin{cases} 2~\pi & \text{ if } m=0 \\ \frac{ (-2)^{m+1}~\pi }{ (-C + \sqrt{C^2-1})^m } {_2}F_1\!\left( m,m; 1; \frac{-C-\sqrt{C^2-1}}{-C+\sqrt{C^2-1}} \right) & \text{ if } m>0 \end{cases} \end{align*}

This result also works, but appears to be strictly worse. The biggest problem is that it doesn't generalize to real-valued $m$. It's also more complicated, and the factorial terms in the series are larger (and so, slower and/or less-accurate to compute).

Some other obvious substitutions: $z:=\exp(2 i x)$ results in a bunch of square roots, while $z:=\exp(i x/2)$ is hard to make go around the pole (mainly because of the differential $d z/dx=i z/2$, without a $z^2$).

Ancillary Thoughts: Elementary-ness

Maybe you think the lack of an elementary solution is obvious—after all, don't computer algebra systems use the Risch Algorithm, which is exhaustive? Moreover, who cares about integrals with non-elementary solutions? We should only care about integrals that we can compute as so are actually useful!

First, as a point of order, the Risch Algorithm is only exhaustive over a limited class of functions. It was made significantly more powerful by later developments, and that's what CASes work at implementing, but there are still limits. In fact, there are significant theoretical limitations on any kind of symbolic integration. For example, allowing antiderivatives with absolute-value makes the problem undecidable in general^[11], which implies that no algorithm can solve all such cases.

Anyway, one issue is that the Risch algorithm is also immensely complicated—on the order of hundreds of pages of mathematics or many thousands of lines of code. The Risch algorithm has never actually been fully implemented, and even so the chance of errors in the partial implementations is high. For example, in Sympy, I discovered that even the $m=1$ special-case exposed a bug. One might expect that paid implementations are higher-quality, but in the past I have found bugs in e.g. Mathematica too—simply being (ridiculously^[12]) expensive is no guarantee of perfection.

In our case, the integrand actually is integrable in closed-form . . . for a particular value of $m$. In this case, the integrand will be a rational function (e.g. see the Weierstrass substitution section above) and the vanilla Risch algorithm can solve it easily. However, when we insist on generalizing it, the complexity demon strikes and makes the problem almost intractable.

The good news is that the non-elementary-ness of a solution does not necessarily pose an obstacle to computation. For example, in a recent paper, I provided a non-elementary solution for normalizing the Phong BRDFs. This was the first complete solution, as well as the fastest by far—it is in fact constant-time! Our integral here actually derives from the study of a related problem, normalizing the Blinn–Phong BRDF.

Ancillary Thoughts: Lessons Learned

It seems that if we can't solve a problem one way, we can't solve it any way. This is non-obvious and subtle, and the researcher would confront it. The pragmatist should count on it. That is, if you can't reach an elementary solution from a few directions, you should assume that an elementary solution does not exist.
Similarly, beyond a certain level of mathematical competence, if you can't get an answer with some technique, it's a good hint that the problem isn't you. It's that an elementary answer does not exist.
If an elementary solution does not exist, embrace it early, breaking the problem up immediately into small problems that we can solve. In our example above, we finally got an answer with a series expansion straightaway, rather than going down successive rabbitholes. Not only does this save time cutting off avenues that only fail many steps deep, it embraces the non-elementary-ness of the problem, allowing us to tackle it better.

Notes

[1]	Some CASes are (just barely) able to get a result in terms of the Appell hypergeometric series—a strictly worse result.
[2]	The arcsine can be integrated by parts, which is of course pointless because this leads you back to the original integral. It's possible there's a recurrence here, which could be equivalent to the non-closed-form solution I give below, with the second way, though (I haven't checked).
[3]	(Glossing over a few required technical restrictions. Fortunately, these rarely are of any concern and are relatively obvious.)
[4]	This seems to be necessary. Refer to the underlying formalism.
[5]	This is a little involved to prove. It can be shown from the Fubini–Tonelli Theorem or applying the Absolute Convergence Theorem to the partial sum. See this higher-level discussion.
[6]	The exact substitution is slightly tricky to discover. Use $u:=\cos^{k-1}(a x)$ and $v:=\sin(a x)/a$. Then proceed as in the link: the right-hand side will turn into an integral for the power $k-2$ and a factor times the original integral. Factor the original integral onto the left-hand side and divide through to get the recurrence.
[7]	This is slightly tricky to see if you haven't before. Try writing out the individual terms, and then as a factorial first. Then convert with $\Gamma(z)=(z-1)!$. Also note that the numerator would commonly be seen as $(\cdots)!!$ the double-factorial/semifactorial function, which much not be confused with the factorial function applied twice $((\cdots)!)!$.
[8]	Note that in the first line we also added a factor of $\G{1}=(1-1)!=0!=1$.
[9]	Since this writing, CASes since have been able to return a worse version in terms of the Appell series, as mentioned.
[10]	All they'd need to reproduce this would be something like: "Take the binomial series expansion, integrate each term, then rewrite as gammas and sum as a hypergeometric function.". It would take them a while to rederive it, but it would be straightforward.
[11]	See Richardson's Theorem.
[12]	Seriously Mathematica, $1 800 per year for professionals outside a university? Do you think mathematicians are made out of money?

Return to the list of blog posts.

\(k\)	\(0\)	\(2\)	\(4\)	\(6\)	\(8\)	\(1\)	\(3\)	\(5\)	\(7\)	\(9\)
\(\int_0^{2\pi} \cos^k(x) ~ d x\)	\(2\pi\)	\(\frac{1}{2} \cdot 2\pi\)	\(\frac{3}{4} \cdot \frac{1}{2} \cdot 2\pi\)	\(\frac{5}{6} \cdot \frac{3}{4} \cdot \frac{1}{2} \cdot 2\pi\)	\(\frac{7}{8} \cdot \frac{5}{6} \cdot \frac{3}{4} \cdot \frac{1}{2} \cdot 2\pi\)	\(0\)	\(0\)	\(0\)	\(0\)	\(0\)