C++ routine for complex Gamma Function

[Tina]

Dear all,

I'm looking for a routine which calculates the Gamma Function
for a complex valued variable. I'm using

#include <complex>

to work with complex numbers.
I define a complex number as
complex <double> alpha(3.0, 1.0)

So here alpha is a complex number with real part of alpha = 3.0
and imaginary part of alpha = 1.0.

By now, I've found "gammln" from Numerical Recipes. But
"gammln" is only working for double, i.e. real valued variables.

Typing
complex <double> tau(3.0,1.0);
cout<<"Gamma(3,1) "<<gammln(tau)<<endl;

gives
GammaComplexTest.cpp: In function `int main()':
GammaComplexTest.cpp:89: error: cannot convert `std::complex<double>'
to `
double' for argument `1' to `DP NR::gammln(double)'

It would be much appreciated if someone could tell me
where to find a routine which also allows to use complex valued
variables for the Gamma Function.

Many thanks!

Tina

 

[mlimber]

I'm looking for a routine which calculates the Gamma Function
for a complex valued variable. I'm using

#include <complex>

to work with complex numbers.
I define a complex number as
complex <double> alpha(3.0, 1.0)

So here alpha is a complex number with real part of alpha = 3.0
and imaginary part of alpha = 1.0.

By now, I've found "gammln" from Numerical Recipes. But
"gammln" is only working for double, i.e. real valued variables.

Typing
complex <double> tau(3.0,1.0);
cout<<"Gamma(3,1) "<<gammln(tau)<<endl;

gives
GammaComplexTest.cpp: In function `int main()':
GammaComplexTest.cpp:89: error: cannot convert `std::complex<double>'
to `
double' for argument `1' to `DP NR::gammln(double)'

It would be much appreciated if someone could tell me
where to find a routine which also allows to use complex valued
variables for the Gamma Function.

See the websites in this FAQ or ask in a math newsgroup:

http://www.parashift.com/c++-faq-lite/class-libraries.html#faq-37.9

Cheers! --M

 

[loic-dev]

Hello Tina,

Dear all,

I'm looking for a routine which calculates the Gamma Function
for a complex valued variable. I'm using

It would be much appreciated if someone could tell me
where to find a routine which also allows to use complex valued
variables for the Gamma Function.

I am sure that there are libraries available on the net that provide
what you need.

However, in your particular case, the gamma function isn't too
difficult to implement yourself using Lanczos approximation:
http://en.wikipedia.org/wiki/Lanczos_approximation

Translating the Python code to C++ gives the program after the
signature. Without any explicit or implied warranty ;-)

HTH,
Loic.
--
#include <complex>
#include <iostream>

using namespace std;


static const int g=7;
static const double pi =
3.1415926535897932384626433832795028841972 ;
static const double p[g+2] = {0.99999999999980993, 676.5203681218851,
-1259.1392167224028, 771.32342877765313, -176.61502916214059,
12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6,
1.5056327351493116e-7};

complex<double> gamma( complex<double> z)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<double> x=p[0];
for (int i=1; i<g+2; i++) {
x += p/(z+complex<double>(i,0));
}
complex<double> t = z + (g + 0.5);
return sqrt(2*pi) * pow(t,z+0.5) * exp(-t) * x;
}

int
main()
{
cout << gamma(complex<double>(5,0)) << endl; // should be 4!
}

 

[Steve Pope]

I'm looking for a routine which calculates the Gamma Function
for a complex valued variable. I'm using

#include <complex>

to work with complex numbers.
I define a complex number as
complex <double> alpha(3.0, 1.0)

So here alpha is a complex number with real part of alpha = 3.0
and imaginary part of alpha = 1.0.

By now, I've found "gammln" from Numerical Recipes. But
"gammln" is only working for double, i.e. real valued variables.

That's because the gamma (and hence, the gammaln) function is only
implemented for real-valued arguments by Matlab. Google on
"mathworks" and "gammaln". If you need the complex gamma function,
you may have to write it yourself, and there may be significant
numerical analysis work to get it to function properly.

I'd ask over in sci.math if anyone knows of an algorithm.

Steve

 

[Steve Pope]

complex<double> gamma( complex<double> z)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<double> x=p[0];
for (int i=1; i<g+2; i++) {
x += p/(z+complex<double>(i,0));
}
complex<double> t = z + (g + 0.5);
return sqrt(2*pi) * pow(t,z+0.5) * exp(-t) * x;

The above line has a complex<double> argument to pow().
Will that work?

Thanks

Steve

 

[loic-dev]

Hello Steve,

complex<double> gamma( complex<double> z)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<double> x=p[0];
for (int i=1; i<g+2; i++) {
x += p/(z+complex<double>(i,0));
}
complex<double> t = z + (g + 0.5);
return sqrt(2*pi) * pow(t,z+0.5) * exp(-t) * x;

The above line has a complex<double> argument to pow().
Will that work?

I am not a C++ expert, but I believe that <complex> is mandated by the
ISO C++ standard to define pow() for the following cases:
template<class T> complex<T> pow(const complex<T>&, int);
template<class T> complex<T> pow(const complex<T>&, const T&);
template<class T> complex<T> pow(const complex<T>&, const complex<T>&);
template<class T> complex<T> pow(const T&, const complex<T>&);

HTH,
Loic.

 

[loic-dev]

Hello Steve,

complex<double> gamma( complex<double> z)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<double> x=p[0];
for (int i=1; i<g+2; i++) {
x += p/(z+complex<double>(i,0));
}
complex<double> t = z + (g + 0.5);
return sqrt(2*pi) * pow(t,z+0.5) * exp(-t) * x;

The above line has a complex<double> argument to pow().
Will that work?

I am not a C++ expert, but I believe that <complex> is mandated by the
ISO C++ standard to define pow() for the following cases:
template<class T> complex<T> pow(const complex<T>&, int);
template<class T> complex<T> pow(const complex<T>&, const T&);
template<class T> complex<T> pow(const complex<T>&, const complex<T>&);
template<class T> complex<T> pow(const T&, const complex<T>&);

HTH,
Loic.

 

[Steve Pope]

Hello Steve,

I am not a C++ expert, but I believe that <complex> is mandated by the
ISO C++ standard to define pow() for the following cases:
template<class T> complex<T> pow(const complex<T>&, int);
template<class T> complex<T> pow(const complex<T>&, const T&);
template<class T> complex<T> pow(const complex<T>&, const complex<T>&);
template<class T> complex<T> pow(const T&, const complex<T>&);

Thanks, this does seem to work if I include both <complex> and
<cmath>.

Steve