For application purposes, I would like to compute the Fourier transform of a power law distribution, $$f(t) = b a^b (a + t)^{-b-1}, \qquad a, b, t > 0,$$ which is the integral given by \begin{align} \tilde f(\xi) &= b a^b \int_0^{+\infty} (a+t)^{-b-1} e^{-i \xi t} \mathrm{d}t\\ &= b e^{i \xi a} \int_1^{+\infty} u^{-b-1}e^{-i \xi a u} \mathrm{d}u. \end{align} Note that the Fourier transform is computed on $[0, +\infty]$ specifically, as $f(t) = 0$ for all $t < 0$.
This problem can be reduced to finding a good approximation for the incomplete Gamma function with pure imaginary argument (see below): $$\Gamma(\nu, ix) = \int_1^{+\infty} t^{\nu - 1} e^{-i x t} \mathrm{d}t,\qquad 0 < \nu < 1.$$
So I'm looking to either find an explicit formulation or an approximation for the Fourier transform of the power-law distribution, or for the incomplete Gamma function with pure imaginary argument. Ideally, I would like this Fourier transform to be easy to calculate in C++, as it is the programming language I am currently using for a statistical package I am developing (so any answer involving existing C++ libraries is welcome!).
What I have achieved so far: when $b$ is an integer, using successive integrations by part, I come up with: $$\tilde f(\xi) = 1 + \sum_{k=1}^{b-1} \frac{(-i \xi a)^k}{(b-1)\cdots(b-k)} + \frac{(-i \xi a)^b}{(b-1)!} e^{i \xi a} E_1(i \xi a),\tag{1}\label{eq1}$$ where $E_1(z)$ denotes the exponential integral $$E_1(z) = \int_1^{+\infty} t^{-1}e^{- tz} \mathrm{d}t.$$ It turns out the exponential integral with pure imaginary argument can be related to the trigonometric integrals $\mathrm{Si}$ and $\mathrm{Ci}$ by the relation (https://en.wikipedia.org/wiki/Exponential_integral): $$E_1(ix) = i \left[ -\frac{1}{2} \pi + \mathrm{Si}(x) \right] - \mathrm{Ci}(x), \qquad x > 0,$$ and that $\mathrm{Si}$ and $\mathrm{Ci}$ have readily available Padé approximants (i.e. can be approximated by rational functions) which are accurate to $10^{-16}$ (https://en.wikipedia.org/wiki/Trigonometric_integral). These steps allow fast computation of the Fourier transform of the power law distribution, provided that $b$ is an integer.
What I struggle with: when $b$ is not an integer, with similar successive integrations by part than before, I end up with $\Gamma(b - \lfloor b \rfloor, i \xi a)$ instead of $E_1(i \xi a)$ in \eqref{eq1}, where $\Gamma(\nu, z)$ is the incomplete Gamma function: $$\Gamma(\nu, z) = \int_1^{+\infty} t^{\nu - 1} e^{-z t} \mathrm{d}t,\qquad 0 < \nu < 1.$$ However, I have not found any explicit method to calculate a good approximation of the incomplete Gamma function for $z$ pure imaginary: the C++ libraries I found only work for real arguments, and I could not find good approximations like the ones I did for the exponential integral.
Another stackexchange post (Computing Fourier transform of power law) has already discussed this problem, but with a symmetric power-law; I suspect the symmetry facilitates the resolution in this case.
Using the residue theorem, I tried relating $\Gamma(\nu, ix)$ to $\Gamma(\nu, x)$, but this leaves me with an intractable integral along a quarter circle between the real and the imaginary axis; something along the line of $$\int_0^{\pi/2} e^{i\theta \nu} e^{-\cos \theta - i \sin \theta} \mathrm{d}\theta.$$
Ending thoughts: Since the power-law distribution is quite a common function to study, I am surprised that I was not able to find an explicit formulation, or at least an approximation method for its Fourier transform. I am sure that I missed a seminal work on this matter, and if one indeed exists, I would be grateful if you could direct me to it.
It's a strange thing to take a Fourier transform of as it's not really periodic. But perhaps it is a characteristic function or similar.
The Mathematica Command
FourierTransform[b a^b (a+t)^(-b-1),t,x]gives the output $$ \frac{\sqrt{\frac{\pi }{2}} b \left(-\frac{1}{a}\right)^{-b} \csc (\pi b) | x| ^{b-1} }{\Gamma (b+1)} \left(\left(-\frac{1}{a}\right)^b a^b \left(-| x| \cos \left(a | x| +\frac{\pi b}{2}\right)+i x \sin \left(a | x| +\frac{\pi b}{2}\right)\right)+i x \sin \left(\frac{\pi b}{2}-a | x| \right)+| x| \cos \left(\frac{\pi b}{2}-a | x| \right)\right) $$ which seems to be composed of standard library functions.This simplifies down to $$ \frac{i \sqrt{\frac{\pi }{2}} a^b (x-| x| ) | x| ^{b-1} e^{i a | x| -\frac{i \pi b}{2}}}{\Gamma (b)} $$ with $a,b>0$. Let me know if this is not useful, or not what you were after.