How can I compute the factorial of a complex number?

Question

To work with complex numbers, I created a struct Cmplx that stores the real and imaginary parts in two separate variables.

struct Cmplx
{
    long double real, imag;
};

I'm having difficulty computing the factorial of a Cmplx, i.e., (a + bi)!. I could use the gamma function, but I don't know how to extend it to a struct. I tried to calculate it using powers, but it required log((a+bi)!) which leads me back to the problem of computing the factorial.

Answer 1

You can use Lanczos approximation for estimating gamma:

#include <iostream>
#include <cmath>
#include <complex>

struct Cmplx
{
    long double real, imag;

    Cmplx(long double r = 0, long double i = 0) : real(r), imag(i) {}

    static const std::complex<long double> gamma(const std::complex<long double> &z)
    {
        static const long double g = 8.0;
        static const long double p[] = {
            0.9999999999999999298L,
            1975.3739023578852322L,
            -4397.3823927922428918L,
            3462.6328459862717019L,
            -1156.9851431631167820L,
            154.53815050252775060L,
            -6.2536716123689161798L,
            0.034642762454736807441L,
            -7.4776171974442977377e-7L,
            6.3041253821852264261e-8L,
            -2.7405717035683877489e-8L,
            4.0486948817567609101e-9L};

        std::complex<long double> x = p[0];
        for (int i = 1; i < sizeof(p) / sizeof(p[0]); ++i)
            x += p[i] / (z + static_cast<long double>(i));

        std::complex<long double> t = z + g - 0.5L;
        return std::pow(2 * M_PI, 0.5L) * std::pow(t, z - 0.5L) * std::exp(-t) * x;
    }

    Cmplx factorial() const
    {
        std::complex<long double> z(real, imag);
        std::complex<long double> res = gamma(z + 1.0L);
        return Cmplx(res.real(), res.imag());
    }
};

int main()
{
    Cmplx z(1.5L, 0.5L);
    Cmplx res = z.factorial();
    std::cout << z.real << " + " << z.imag << "i\n";
    std::cout << res.real << " + " << res.imag << "i\n";
}

Prints

1.5 + 0.5i

0.579116 + 0.294109i