Alright, as I've mentioned in the comments, I've already had a half typed answer so it feels wasteful to just delete it. I'll add a generalization to the integral in question as well as a few insights to make the answer more meaningful.
(And i have to thank @Mark Viola for pointing out the phase shift formula for the trig function lol I was wondering why my result wasn't in the correct format )
Let $n=p-1$.
$$f(z)=\frac{z^n}{z^2+1}$$
Simple poles from denominator at $x=\pm i$, and since $n\neq\mathbb{Z}$ we have a branch point at $0$ and complex infinity as well as a branch cut connecting them. We'll define $$0\le\arg(z)<2\pi$$ and have our branch cut be on the positive real axis.
Define a contour as follows

$$\oint_\mathcal{C}f(z)\text{ d}z=\int_T+\int_B+\int_\psi+\int_\Gamma f(z)\text{ d}z=2\pi i\sum\text{Residue at }\pm i$$
Decay rate of function has $\int_\Gamma\to 0$. You can prove this rigorously as an exercise by first paramterizing it(choose some dummy angle say $\delta$ as part of your bounds since the outer contour does not circle the full $2\pi$) and then using typical inequalities(CS, ML, triangle etc). For example, here are some $$\left|\int f(z)\right|\le\int|f(z)|,\quad\int\frac{1}{|a-b|}\le\int\frac{1}{||a|-|b||}\text{, and so on}$$ once you get the limit at the end just note that our restrictions will guarentee it going to $0$.
Similar situation for the branch point integral, $\int_\psi\to 0$. Prove this via parametrization the same way as before with a dummy angle and use inequalities/take limit to show it goes to $0$. I recommend trying these are exercises but if you want just tell me and I'll type out the process.
Paths about $T$ and $B$ can be parameterized as follows.
$$T: z=x+i\epsilon,\quad dz=dx,\quad \epsilon\le z\le R$$
$$B: z=x-i\epsilon,\quad dz=dx,\quad R\le z\le\epsilon$$
\begin{align*}
\lim_{\epsilon\to0,R\to\infty}\int_T&=\int^\infty_0\lim_{\epsilon\to0} \frac{(x+i\epsilon)^n}{(x+i\epsilon)^2+1}\\
\end{align*}
Justify the limit swap using DCT. Let a function $\int^R_0f(z)$ and induce a limit inside first. Bound the integrand by a constant since $0\le x<R$ and use DCT. Then take the $R$ limit. Same thing for the integral about $B$.
Expand out the numerator and denominator with $e^{\ln}$ and rewrite the log into the complex logarithm definition. Split the terms and factor out the argument. In general, we can see that $$\lim_{\epsilon\to0}\int_{\text{path}}f(z)\text{ d}z=\exp\left( i\arg(f(z)\text{ in specified direction}) \right)\int_{\text{path}}\lim_{\epsilon\to0}|f(z)|\text{ d}z$$
Our function is always positive, so we ignore the absolute value in this case.
$T$ path has an argument of $0$ so $\int_T$ is just our integral we seek.
$B$ path has argument $2\pi$, so $\int_B$ is $-e^{2\pi i n}$ times our integral($n$ comes from numerator power)
The residues can be found by just using the derivative definition, which gives(be careful of your argument $i=e^{\pi i/2}$, $-i=e^{3\pi i/2}$)
$$2\pi i\sum\text{Residue at }\pm i=2\pi i\left(\frac{i^n}{2i}+\frac{(-i)^n}{-2i}\right)=2\pi i\left(\frac{e^{i\pi n/2}}{2i}+\frac{e^{3i\pi n/2}}{-2i}\right)=-2\pi ie^{i\pi n}\left(\frac{e^{i\pi n/2}-e^{-i\pi n/2}}{2i}\right)=-2\pi ie^{i\pi n}\sin\left(\frac{n\pi}2\right)$$
Hence we have $$I(1-e^{2\pi i n})=-2\pi ie^{i\pi n}\sin\left(\frac{n\pi}2\right)$$
$$\implies I=\color{blue}{-2}\pi \color{blue}{i}\color{red}{e^{i\pi n}}\sin\left(\frac{n\pi}2\right)\cdot\frac{\color{blue}{1}}{\color{red}{e^{\pi i n}}\color{blue}{\left(e^{-\pi i n}-e^{\pi i n}\right)}}=\pi\cdot\sin\left(\frac{p\pi}2-\frac{\pi}2\right)\cdot\color{blue}{\csc(p\pi-\pi)}$$
$$\implies \pi\cdot\color{red}{-\cos\left(\frac{p\pi}2\right)}\cdot\frac{\color{red}{-}1}{2\sin\left(\frac{p\pi}2\right)\color{red}{\cos\left(\frac{p\pi}2\right)}}=\frac{2}{\pi}\csc\left(\frac{\pi p}{2}\right)$$
Now to make this more useful, consider the more general integral $n,\,m\in\mathbb{R}$, $n\ge0$ and $n+1<m$
$$I=\int_0^\infty {\frac{{{x^n}}}{{1 + {x^m}}}\text{ d}x }$$
Let $f(z)$ be the integrand. If $n,\,m\notin\mathbb{Z}$, the function has a branch cut from the origin out to complex infinity. Set this on the negative real axis. The poles of our function lie at all the places where $1+z^m=0$, which is along the unit circle. For an integer power $m$, we will get the $m$ roots of unity reflected across the imaginary axis as our poles.
For each noninteger $m$, the number of poles will be the nearest even integer. Note the fact that the only time we will get a pole on the negative real axis is when $m$ is an odd integer, but this also means there is no branch cut.
We'll use a pizza slice contour here as follows

Here, $\theta=\frac{2\pi}m$, and the contour path $U$ ends along the ray $e^{\frac{2\pi i}{m}}$
We have $$\oint_{\mathcal{C}}f(z)\text{ d}z=\int_B+\int_{\Gamma}+\int_U+\int_{\gamma}f(z)\text{ d}z=2\pi i\mathop{\mathrm{Res}}_{z = \exp\left(\frac{\pi i}m\right)}f(z)$$
Residue first this time. $$\mathop{\mathrm{Res}}_{z = \exp\left(\frac{\pi i}m\right)}\frac{z^n}{z^m+1}=\frac{\exp\left(\frac{\pi i}m\right)^n}{\dfrac{\text{d}}{\text{d}z}\left[z^m+1\right]\Big|_{z=\exp\left(\frac{\pi i}m\right)}}=\frac{\exp\left(\frac{\pi in}m\right)}{-m\exp\left(-\frac{\pi i}m\right)}=\frac{\exp\left(\frac{\pi i}m(n+1)\right)}{-m}$$
Before we move on, I'd like to note a small problem when we have $1<m<2$: the path $U$ may actually lie on a pole or enclose another pole and $\Gamma,\,\gamma$ crosses the branch cut. I'll deal with this at the end.
If we parameterize our paths we get \begin{alignat*}{5}
B&:\text{ }z=x,\qquad &\text{d}z&=\text{d}x,\qquad &x&\in[\,\epsilon, R\,]\\
\Gamma &:\text{ }z=Re^{i\theta},\qquad &\text{d}z&=iRe^{i\theta}\text{ d}\theta,\qquad &\theta&\in\left[0, \frac{2\pi}{m}\right]\\
U&:\text{ }z=re^{\frac{2\pi i}m},\qquad &\text{d}z&=e^{\frac{2\pi i}m}\text{d}r,\qquad &r&\in[\,R, \epsilon\,]\\
\gamma &:\text{ }z=\epsilon e^{i\theta},\qquad &\text{d}z&=i\epsilon e^{i\theta}\text{ d}\theta,\qquad &\theta&\in\left[\frac{2\pi}{m},0\right]
\end{alignat*}
With the integrals along the two rays being really nice.
$$\lim_{R\to+\infty,\,\epsilon\to+0}\int_Bf(z)\text{ d}z=\lim_{R\to+\infty,\,\epsilon\to+0}\int^{R}_{\epsilon}\frac{x^n\text{ d}x}{1+x^m}=I$$
and
$$\lim_{R\to+\infty,\,\epsilon\to+0}\int_Uf(z)\text{ d}z=\lim_{R\to+\infty,\,\epsilon\to+0}\int^{\epsilon}_R\frac{r^ne^{\frac{2\pi i n}{m}}}{1+r^m\color{red}{e^{\frac{2\pi i m}{m}}}}\cdot e^{\frac{2\pi i}m}\text{ d}r=-e^{\frac{2\pi i}m(n+1)}I$$
Both arcs go to $0$, but I will show this i guess, so we have
\begin{align}
\left|\int_{\Gamma}f(z)\text{ d}z\right|&\le\int^{\frac{2\pi}m}_0\frac{\left|R^n\right|\color{red}{\left|e^{in\theta}\right|}}{\left|\displaystyle R^m e^{im\theta}+1 \right|}\cdot\color{red}{|i|}|R|\color{red}{\left|e^{i\theta}\right|}\text{ d}\theta\\
&\le\int^{\frac{2\pi}m}_0\frac{R^{n+1}}{\left|\left|\displaystyle R^m\right|\color{red}{\left|e^{im\theta}\right|}-|-1|\right|}\text{ d}\theta\\
&\le\int^{\frac{2\pi}m}_0\frac{R^{n+1}\text{ d}\theta}{\displaystyle R^m-1}=\frac{2\pi}{m}\cdot\frac{R^{n+1}}{\displaystyle R^m-1}
\end{align}
and the limit
$$\frac{2\pi}{m}\lim_{R\to+\infty}\frac{R^{n+1}}{\displaystyle R^m-1}=0$$
only happens if $m>0$, $m>n+1$, and $n>-1$, which is satisfied since we have set the following restrictions where we have $m>2$, $m>n+1$, and $n\ge0$.
next, \begin{align}
\left|\int_{\gamma}f(z)\text{ d}z\right|&\le\int_{\frac{2\pi}m}^0\frac{\left|\epsilon^n\right|\color{red}{\left|e^{in\theta}\right|}}{\left|\displaystyle \epsilon^m e^{im\theta}+1 \right|}\cdot\color{red}{|i|}|\epsilon|\color{red}{\left|e^{i\theta}\right|}\text{ d}\theta\\
&\le\int_{\frac{2\pi}m}^0\frac{\epsilon^{n+1}}{\left|\left|\displaystyle \epsilon^m\right|\color{red}{\left|e^{im\theta}\right|}-|-1|\right|}\text{ d}\theta\\
&\le\int_{\frac{2\pi}m}^0\frac{\epsilon^{n+1}\text{ d}\theta}{\displaystyle \epsilon^m-1}=-\frac{2\pi}{m}\cdot\frac{\epsilon^{n+1}}{\displaystyle \epsilon^m-1}
\end{align}
the limit will become $$-\frac{2\pi}{m}\lim_{\epsilon\to+0}\frac{\epsilon^{n+1}}{\displaystyle \epsilon^m-1}=0$$
with laxer restrictions, with us needing only $m>0$ and $n>-1$, which are both satisfied.
Hence, we have $$\oint_{\mathcal{C}}f(z)\text{ d}z=\int_B+\int_U f(z)\text{ d}z=I-e^{\frac{2\pi i}{m}(n+1)}I=2\pi i\cdot\frac{-1}{m}e^{\frac{\pi i}{m}(n+1)}$$
So we have
$$I=\frac{-\frac{2\pi i}{m}\exp\left(\frac{\pi i}{m}(n+1)\right)}{1-\exp\left(\frac{2\pi i}{m}(n+1)\right)}=\frac{\pi}{m}\cdot\frac{2i}{e^{\frac{\pi i}m(n+1)}-e^{-\frac{\pi i}m(n+1)}}=\frac{\pi}{m}\csc\left(\frac{\pi(n+1)}{m}\right),\qquad m>2,\,m>n+1,\,n\ge0$$
Okay, now let's revisit the discrepancy at the beginning and attempt to get rid of the $m>2$ restriction. The region of integration has $z = re^{i \theta}$ where $\epsilon \leq r \leq R$, $0 \leq \theta \leq \frac{2\pi}{m}$. So to get a function mostly continuous in the region, we can take the branch cut along the ray bisecting the remaining angle between $\frac{2\pi}{m}$ and $2\pi$. With $r,\, \theta \in \mathbb{R}$ and $r \geq 0$, this branch is given by
$$ (r e^{i\theta})^m = r^m \exp\left(i \theta m - 2im\pi \left\lfloor \frac{\theta}{2\pi} + \frac{m-1}{2m} \right\rfloor\right) $$
This new definition utilizes the number of complete counterclockwise rotations relative to the chosen branch cut in the direction of $-e^{\frac{i \pi}m}$. Using this definition of $z^m$, the integrand $f(z) = \frac{z^n}{1+z^m}$ has exactly one pole inside the region of integration, at $z=e^{\frac{i \pi}m}$.
For example, suppose $m = \frac{4}{3}$. The region has $0 \leq \theta \leq \frac{3\pi}{2}$, and $f(z)$'s only pole is at $z=e^{\frac{3i \pi}4}$. Although some branches of $z^{\frac43}$, including the principal branch, have a pole at $z = e^{\frac{-3i \pi}4}$, this new function $f$ does not, as it's defined in terms of a "$z^{\frac43}$" branch which gives $e^{\frac{-3i\pi}4} \mapsto \exp(\frac{4}{3}\cdot\frac{5i \pi}{4}) = e^{\frac{-i \pi}3}$ so there is no division by zero.
We also see that along the straight line segment edges of the path of integration where $\theta = 0$ or $\theta = \frac{2\pi}{m}$, or alternatively $B$ or $U$ respectively, the chosen branch gives $z^m = 1$, so $f(z) = \frac{z^n}{2}$, and the symmetry argument is allowed since $f(z)$ has different values on $U$ if the principal branch of $z^m$ is used with $m<2$. This justifies using our pizza slice contour for $1<m<2$. Thus, we must have it such that
$$I=\frac{\pi}{m}\csc\left(\frac{\pi(n+1)}{m}\right),\qquad m>n+1,\,n\ge0$$