complex numbers

[Blair]

could someone PLEASE tell me why this doesn't work...

-----------------------------------------
#include <complex>
using namespace std;
typedef complex<long double> cld;

void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}
--------------------------------------------
the square root of -1 is j (or i) but why I debug this code temp is given as
(1,0).

or

if I do this...
void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}

I get the same results. temp=(1, 0) when the real is 1 and the imaginary is
0.

So, what is wrong with the pow function? Why can't I use fractions as the
exponent? (BTW, I need to find cube-roots. That's why I'm not using
sqrt(...) )


Thanks, Blair

 

[Dan Cernat]

could someone PLEASE tell me why this doesn't work...

-----------------------------------------
#include <complex>
using namespace std;
typedef complex<long double> cld;

void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}
--------------------------------------------
the square root of -1 is j (or i) but why I debug this code temp is given as
(1,0).

or

if I do this...
void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}

I get the same results. temp=(1, 0) when the real is 1 and the imaginary is
0.

So, what is wrong with the pow function? Why can't I use fractions as the
exponent? (BTW, I need to find cube-roots. That's why I'm not using
sqrt(...) )


Thanks, Blair

my first guess is that pow function (as is declared in <cmath>) gets the
first argument converted to a double or float.

dan

 

[Cy Edmunds]

could someone PLEASE tell me why this doesn't work...

-----------------------------------------
#include <complex>
using namespace std;
typedef complex<long double> cld;

void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}
--------------------------------------------
the square root of -1 is j (or i) but why I debug this code temp is given as
(1,0).

or

if I do this...
void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}

I get the same results. temp=(1, 0) when the real is 1 and the imaginary is
0.

So, what is wrong with the pow function? Why can't I use fractions as the
exponent? (BTW, I need to find cube-roots. That's why I'm not using
sqrt(...) )


Thanks, Blair

Try

temp = pow(cmplx, 0.5L);

 

[Victor Bazarov]

given

my first guess is that pow function (as is declared in <cmath>) gets the
first argument converted to a double or float.

<cmath> is not included and <complex> declares its own 'std:ow' (and
actually the whole four of them).

Got a second guess?

My first guess is that 'void main' causes it.

Just kidding. Actually, as I found by trying, changing it to

temp = pow(cmplx, cld(0.5, 0));

gives the right result. Why it happens, I can't say, it's too late
here for me to strain my brain :-[

Victor

 

[CrayzeeWulf]

could someone PLEASE tell me why this doesn't work...

-----------------------------------------
#include <complex>
using namespace std;
typedef complex<long double> cld;

void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}

Hi Blair,

There are four different versions of std:ow for complex<T> numbers:

1. template<class T>
complex<T> pow( const complex<T>& x, int y) ;

2. template<class T>
complex<T> pow( const complex<T>& x, const complex<T>& y) ;

3. template<class T>
complex<T> pow( const complex<T>& x, const T& y) ;

4. template<class T>
complext<T> pow( const T& x, const complex<T>& y ) ;

In your case "T" is "long double". Hence, the call "pow(cmplx, 0.5)" does
not directly match any of the above without requiring some implicit
conversion. Overload resolution rules end up matching the call to prototype
#1 in the above list when you really want to use #3. In order to match #3,
you have to help the compiler a little bit by telling it that "0.5" in your
function call is indeed a "long double":

temp = pow( cmplx, 0.5L ) ;

Another option is to use prototype #2 by converting "0.5" into a
complex<long double>:

temp = pow( cmplx, cld(0.5, 0.0) ) ;

if I do this...
void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}

I do not see the difference between this and the version above.

Later,

 

[Blair]

Cy and Victor,
Thanks for the replies. They both worked for a real result case - ie. when
cmplx=cld(4,0) then temp becomes 2. But neither answer works for imaginary
results.

PS. I made a mistake in my original post ... I wrote
--------------------------------
....or

if I do this...
void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}
-----------------------------------
but that code is exactly the same as the first snippet of code that posted!
I mean to say this...
--------------------------------
or

if I do this...
void main()
{
cld cmplx, temp;
cmplx = cld(4.0, 0.0);
temp = pow(cmplx, 0.5);
}
-----------------------------------
to make cmplx = 4.0 + 0j
then temp should = 2.0, but I was getting temp=1+0j. Now, both of your
answers correct the case when sqrt(cmplx) is a real number. But it still
fails when cmplx = -1. (which, of course, should make temp=0+j)

Thanks,
blair

 

[Blair]

Hi Blair,

There are four different versions of std:ow for complex<T> numbers:

1. template<class T>
complex<T> pow( const complex<T>& x, int y) ;

2. template<class T>
complex<T> pow( const complex<T>& x, const complex<T>& y) ;

3. template<class T>
complex<T> pow( const complex<T>& x, const T& y) ;

4. template<class T>
complext<T> pow( const T& x, const complex<T>& y ) ;

In your case "T" is "long double". Hence, the call "pow(cmplx, 0.5)" does
not directly match any of the above without requiring some implicit
conversion. Overload resolution rules end up matching the call to prototype
#1 in the above list when you really want to use #3. In order to match #3,
you have to help the compiler a little bit by telling it that "0.5" in your
function call is indeed a "long double":

temp = pow( cmplx, 0.5L ) ;

Another option is to use prototype #2 by converting "0.5" into a
complex<long double>:

temp = pow( cmplx, cld(0.5, 0.0) ) ;

Hi, and thanks for the response. But, ...darn it... both suggestions were
already posted and neither of them seemed to work for me for some reason. In
fact, if I instead do this...

void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 0.0);
temp = pow(cmplx, 0.5);
}

I get the result: temp = (-1.#IND000000000, 0.00000000000000). HOWEVER, if I
make the imaginary part of cmplx NON-ZERO, as in...

typedef complex<long double> cld;
void main()
{
cld cmplx, temp;
cmplx = cld(-1.0, 1.0);
temp = pow(cmplx, 0.5L);
}

I do get the correct answer: temp = (0.45508986056223, 1.0986841134678)!!!

Strange, huh?

Blair

I do not see the difference between this and the version above.

Yeah, sorry about that. I noticed that too. I posted a correction.

 

[Cy Edmunds]

[snip]

Here is my test program.

void yowza()
{
std::complex<long double> cmplx(-1.0, 0.0), temp;
temp = std:ow(cmplx, 0.5L);
std::cout << temp << '\n';
}

Output is:

(6.12303e-017,1)

Compiler is VC++.Nuts

 

[Blair]

Crap, I messed up again... see correction below...

Hi, and thanks for the response. But, ...darn it... both suggestions were
already posted and neither of them seemed to work for me for some reason. In
fact, if I instead do this...

------insert the following line-------

 

[Blair]

So, it is a compiler issue. I get:

Output: (-1.#IND,0)

from your code.

Thanks for the help. I'd better investigate the compiler, I guess.

Blair

 

[Blair]

Compiler is VC++.Nuts

BTW, what do you mean by .Nuts? I am using VC++ 6.0 and it's not working.
What is .Nuts? Is that what I am missing?

Blair

 

[Victor Bazarov]

BTW, what do you mean by .Nuts? I am using VC++ 6.0 and it's not working.
What is .Nuts? Is that what I am missing?

It's .Net, my guess. You're missing quite some improvements
in the Standard compliance if you're still using 6.0...

 

[Jim West]

Hi Blair,

There are four different versions of std:ow for complex<T> numbers:

1. template<class T>
complex<T> pow( const complex<T>& x, int y) ;

2. template<class T>
complex<T> pow( const complex<T>& x, const complex<T>& y) ;

3. template<class T>
complex<T> pow( const complex<T>& x, const T& y) ;

4. template<class T>
complext<T> pow( const T& x, const complex<T>& y ) ;

In your case "T" is "long double". Hence, the call "pow(cmplx, 0.5)" does
not directly match any of the above without requiring some implicit
conversion. Overload resolution rules end up matching the call to prototype
#1 in the above list when you really want to use #3. In order to match #3,

Wow! This is potentially a big problem for someone who uses complex
numbers extensively. I changed the program to

using namespace std;
typedef complex<float> cld;

int main() {
cld cmplx, temp;
cmplx = cld(-1.0f, 0.0f);
temp = pow(cmplx, 0.5);
}

and ended up getting the same behavior. Dropping the f or L does not
strike me as all that unlikely. The GNU compiler at least warned that
it was passing a double to an int argument, but the Intel compiler
said nothing. There may be a lot of programs out there unknowingly
getting the wrong answer! Could this be considered a flaw in the
standard library definition of the complex pow() function? It would
be better if they left version 1 (with int) out altogether. It might
take longer some times, but at least no one would get the wrong answer
for using the wrong precision in the second argument.

Seriously, this is very, very scary.

 

[CrayzeeWulf]

Hi, and thanks for the response. But, ...darn it... both suggestions were
already posted and neither of them seemed to work for me for some reason.
In fact, if I instead do this...

I do get the correct answer: temp = (0.45508986056223, 1.0986841134678)!!!

Strange, huh?

Hi Blair,

It looks like a compiler issue. The code craps out with Microsoft Visual C++
6.0 but works fine with Microsoft Visual Studio .Net 2003. It also works as
expected with g++ 3.2.3. VC6 has always been buggy when it comes to the STL
(in my limited experience). I have not been able to find a workaround for
VC6, unfortunately.

Later,

 

[CrayzeeWulf]

Jim West wrote:

Wow! This is potentially a big problem for someone who uses complex
numbers extensively. I changed the program to

There may be a lot of programs out there unknowingly
getting the wrong answer! Could this be considered a flaw in the
standard library definition of the complex pow() function?

Seriously, this is very, very scary.

Yeah. It does not make sense. However, I have been looking through the
standard and am not yet convinced if VC.Net and G++ are doing the right
thing. For example, N. M. Josuttis's book (The C++ Standard Library, 3rd
Printing) explicitly states that "pow(c, 1.7)" should result in c^(1.7). Of
course, VC6 is definitely wrong.

Furthermore, this works without any warning under G++ 3.2.3 as well as
VC.Net 2003:

#include <complex>
#include <iostream>
using namespace std;
typedef complex<float> cld;

int main()
{
float half = 0.5 ;
cld cmplx, temp ;
cmplx = cld(-1.0, 0.0) ;
temp = pow(cmplx, half) ;
cout << temp << endl ;
return EXIT_SUCCESS ;
}

 

[CrayzeeWulf]

Yeah. It does not make sense. However, I have been looking through the
standard and am not yet convinced if VC.Net and G++ are doing the right
thing.

Aha. I was looking at the wrong place. The answer lies in the rules of
template argument deduction. In particular, the text of Section 14.8.1
Paragraph 4 of ISO/IEC 14882:2003(E) applies here. According to this
paragraph, during template argument deduction for function templates,
implicit conversions are performed only on function parameters whose types
do not contain template parameters that are used in template argument
deduction.

So if we consider the available function templates:

1. template<class T>
   complex<T> pow( const complex<T>& x, int y) ;

2. template<class T>
   complex<T> pow( const complex<T>& x, const complex<T>& y) ;

3. template<class T>
   complex<T> pow( const complex<T>& x, const T& y) ;

4. template<class T>
   complext<T> pow( const T& x, const complex<T>& y ) ;

and the following call:

using namespace std ;
complex<long double> temp, cmplx(-1, 0) ;
temp = pow(cmplx, 0.5) ;

Function templates #2 and #4 are obviously out. Function template #3 is not
considered because the template argument deduction (i.e. deduction of the
type of T in the above call) does not consider the implicit conversion of
0.5 to "long double". Hence, the only candidate is function template #1.

However, the story is quite different if we change the above call to:

using namespace std ;
complex<long double> temp, cmplx(-1, 0) ;
temp = pow<long double>(cmplx, 0.5) ; // ! ERROR

In this case, there is no need to perform template argument deduction and
both #1 and #3 are candidates. Hence, the call is ambiguous and an error.

Finally, the following uses #3 because 0.5 is a double and no implicit type
conversion is required for overload resolution:

using namespace std;
complex<double> temp, cmplx(-1, 0) ;
temp = pow(cmplx, 0.5) ;

In other words, the complex<T> pow() function templates provide "intuitive"
results only when T is double. In other cases, you have to be careful.

Later,

 

[CrayzeeWulf]

Yeah. It does not make sense. However, I have been looking through the
standard and am not yet convinced if VC.Net and G++ are doing the right
thing.

Aha. I was looking at the wrong place. The answer lies in the rules of
template argument deduction. In particular, the text of Section 14.8.1
Paragraph 4 of ISO/IEC 14882:2003(E) applies here. According to this
paragraph, during template argument deduction for function templates,
implicit conversions are performed only on function parameters whose types
do not contain template parameters that are used in template argument
deduction.

So if we consider the available function templates:

1. template<class T>
   complex<T> pow( const complex<T>& x, int y) ;

2. template<class T>
   complex<T> pow( const complex<T>& x, const complex<T>& y) ;

3. template<class T>
   complex<T> pow( const complex<T>& x, const T& y) ;

4. template<class T>
   complext<T> pow( const T& x, const complex<T>& y ) ;

and the following call:

using namespace std ;
complex<long double> temp, cmplx(-1, 0) ;
temp = pow(cmplx, 0.5) ;

Function template #4 is obviously out. Function template #3 is not
considered because the template argument deduction (i.e. deduction of the
type of T in the above call) does not consider the implicit conversion of
0.5 to "long double". Similarly, function template #2 is not considered even
though a conversion of 0.5 to complex<long double> is available. Hence, the
only candidate is function template #1.

However, the story is quite different if we change the above call to:

using namespace std ;
complex<long double> temp, cmplx(-1, 0) ;
temp = pow<long double>(cmplx, 0.5) ; // ! ERROR

In this case, there is no need to perform template argument deduction and
#1, #2, and #3 are all candidates. Hence, the call is ambiguous and an
error.

Finally, the following uses #3 because 0.5 is a double and no implicit type
conversion is required for overload resolution:

using namespace std;
complex<double> temp, cmplx(-1, 0) ;
temp = pow(cmplx, 0.5) ;

In other words, the complex<T> pow() function templates provide "intuitive"
results only when T is double. In other cases, you have to be careful.

Later,

 

[KTC]

So, it is a compiler issue. I get:

Output: (-1.#IND,0)

from your code.

Thanks for the help. I'd better investigate the compiler, I guess.

Blair

Yeah, you ought to be getting some along the line of (something, 1).
There ususally is something in the Real part coz of small errors and
things working with floating point numbers...

KTC

 

[Cy Edmunds]

Compiler is VC++.Nuts

BTW, what do you mean by .Nuts? I am using VC++ 6.0 and it's not working.
What is .Nuts? Is that what I am missing?

Blair

[snip]

Here is my test program.

void yowza()
{
std::complex<long double> cmplx(-1.0, 0.0), temp;
temp = std:ow(cmplx, 0.5L);
std::cout << temp << '\n';
}

Output is:

(6.12303e-017,1)

Compiler is VC++.Nuts

Sorry. My lame idea of a joke. It's really Visual C++.Net. But on my desktop
the compiler is called Miserable C++.Nuts.

Actually, though, I mostly like it.

 

[Jim West]

In other words, the complex<T> pow() function templates provide "intuitive"
results only when T is double. In other cases, you have to be careful.

In reaction to this thread, today I added to following code to the header
file the I use with most of my numerical processing to make sure I don't
get bitten by this:

inline std::complex<float> pow(const std::complex<float>&x, double n) {
std::cerr << "Warning: Called a mixed mode std::complex pow() function."
<< std::endl;
return pow(x, float(n));
}

I also added several others to cover all combinations of float, double,
and long double. Hopefully I'll never seen the warning come up.


I really do believe that this is a serious problem with the standard
library implementation of complex<T>pow() since the people who will usually
use it will very often be coming from a Fortran environment where
0.5 is single precision, so using it (instead of 0.5f) in a floating
point expression looks "right". I understand that one should know
the syntax of the language that they are using, but the fact that the
double would get demoted to an int in this case to yield an incorrect
answer is extremely counter-intuative to someone who doesn't know the
intimate details of the language.

So, is there any procedure an inexpert C++ programmer like me could
(or better yet, should) do to make sure that the powers-that-be are
aware of what I think is a potential problem?