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.

Latashialatch

std::complex already exists. – Posey

Instead of pow three times, look up function hypot. – Brewis

There's the std::tgamma function for buit-in numeric types. But no counterpart for std::complex. I guess the mathematical evaluation of the integral in complex domain is too complex. – Teal

Factorials of complex numbers are rather esoteric. If your question is "How is factorial defined for complex numbers (in terms of real and imaginary components)?" that's a question for math.stackexchange.com. If you're having difficulty converting that pre-existing formula into code, that's a slightly different question than you've posted here. There's also the big question of why you're attempting this in the first place. I've a feeling this may be an XY problem, where the thing you're really trying to address doesn't actually need the computation. – Cherey

The best I got was gsl_sf_lngamma_complex_e: gnu.org/software/gsl/doc/html/… – Teal

Re "I'm having difficulty computing the factorial": Can you be more specific? E.g., what were the symptoms? Please respond by editing (changing) your question, not here in comments (but *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** without *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** "Edit:", "Update:", or similar - the question should appear as if it was written right now). – Koloski

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

Spue