469,939 Members | 2,345 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 469,939 developers. It's quick & easy.

C++ routine for complex Gamma Function

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

Dec 13 '06 #1
7 8052
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
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-lit....html#faq-37.9

Cheers! --M

Dec 13 '06 #2
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.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<doublegamma( complex<doublez)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<doublex=p[0];
for (int i=1; i<g+2; i++) {
x += p[i]/(z+complex<double>(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<double>(5,0)) << endl; // should be 4!
}

Dec 13 '06 #3
Tina <ti*********@yahoo.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
Dec 14 '06 #4
<lo******@gmx.netwrote:
>complex<doublegamma( complex<doublez)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<doublex=p[0];
for (int i=1; i<g+2; i++) {
x += p[i]/(z+complex<double>(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<doubleargument to pow().
Will that work?

Thanks

Steve
Dec 14 '06 #5
Hello Steve,
complex<doublegamma( complex<doublez)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<doublex=p[0];
for (int i=1; i<g+2; i++) {
x += p[i]/(z+complex<double>(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<doubleargument 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(const complex<T>&, int);
template<class Tcomplex<Tpow(const complex<T>&, const T&);
template<class Tcomplex<Tpow(const complex<T>&, const complex<T>&);
template<class Tcomplex<Tpow(const T&, const complex<T>&);

HTH,
Loic.

Dec 15 '06 #6
Hello Steve,
complex<doublegamma( complex<doublez)
{

if ( real(z)<0.5 ) {
return pi / (sin(pi*z)*gamma(1.0-z));
}
z -= 1.0;
complex<doublex=p[0];
for (int i=1; i<g+2; i++) {
x += p[i]/(z+complex<double>(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<doubleargument 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(const complex<T>&, int);
template<class Tcomplex<Tpow(const complex<T>&, const T&);
template<class Tcomplex<Tpow(const complex<T>&, const complex<T>&);
template<class Tcomplex<Tpow(const T&, const complex<T>&);

HTH,
Loic.

Dec 15 '06 #7
<lo******@gmx.netwrote:
>Hello Steve,
return sqrt(2*pi) * pow(t,z+0.5) * exp(-t) * x;

The above line has a complex<doubleargument 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(const complex<T>&, int);
template<class Tcomplex<Tpow(const complex<T>&, const T&);
template<class Tcomplex<Tpow(const complex<T>&, const complex<T>&);
template<class Tcomplex<Tpow(const T&, const complex<T>&);
Thanks, this does seem to work if I include both <complexand
<cmath>.

Steve
Dec 15 '06 #8

This discussion thread is closed

Replies have been disabled for this discussion.

Similar topics

2 posts views Thread by Scott | last post: by
8 posts views Thread by Roberto Dias | last post: by
45 posts views Thread by OcelotIguana | last post: by
reply views Thread by james.duckworthy | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.