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 <doublealpha(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 <doubletau(3.0, 1.0);
cout<<"Gamma(3, 1) "<<gammln(tau)< <endl;

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

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]

Tina wrote:

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 <doublealpha(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 <doubletau(3.0, 1.0);
cout<<"Gamma(3, 1) "<<gammln(tau)< <endl;

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

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-lit....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

<snip>

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.1415926535897 932384626433832 795028841972 ;
static const double p[g+2] = {0.999999999999 80993, 676.52036812188 51,
-1259.1392167224 028, 771.32342877765 313, -176.61502916214 059,
12.507343278686 905, -0.1385710952657 2012, 9.9843695780195 716e-6,
1.5056327351493 116e-7};

complex<doubleg amma( complex<doublez )
{

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

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

[Steve Pope]

Tina <ti*********@ya hoo.cawrote:

>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 <doublealpha(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]

<lo******@gmx.n etwrote:

>complex<double gamma( complex<doublez )
{

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

The above line has a complex<doublea rgument to pow().
Will that work?

Thanks

Steve

[loic-dev]

Hello Steve,

complex<doubleg amma( complex<doublez )
{

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

The above line has a complex<doublea rgument to pow().
Will that work?

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

HTH,
Loic.

[loic-dev]

Hello Steve,

complex<doubleg amma( complex<doublez )
{

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

The above line has a complex<doublea rgument to pow().
Will that work?

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

HTH,
Loic.

[Steve Pope]

<lo******@gmx.n etwrote:

>Hello Steve,

return sqrt(2*pi) * pow(t,z+0.5) * exp(-t) * x;

The above line has a complex<doublea rgument to pow().
Will that work?

>I am not a C++ expert, but I believe that <complexis mandated by the
ISO C++ standard to define pow() for the following cases:
template<cla ss Tcomplex<Tpow(c onst complex<T>&, int);
template<cla ss Tcomplex<Tpow(c onst complex<T>&, const T&);
template<cla ss Tcomplex<Tpow(c onst complex<T>&, const complex<T>&);
template<cla ss Tcomplex<Tpow(c onst T&, const complex<T>&);

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

Steve