cexp implementation questions

[steve]

Given

#include <stdio.h>
#include <complex.h>

int
main (void)
{
complex z;
z = 1000. + I * 0;
z = cexp(z);
printf("%lf %lf\n", creal(z), cimag(z));
return (0);
}

What is the expected output and what is correct output?

 

[steve]

In general, exp(x+yi) = exp(x)(cos y + i sin y)

exp (1000+0i) = exp(1000)(cos 0 + i sin 0) = exp(1000) (1+0i) = exp(1000)

I think exp(1000) is out of range for IEEE reals.

Yes, it is out of range, and many (most? all?) compilers will
return exp(1000) = +Inf (unless one uses an option to trap
FP exceptions). So, your mathematics may break down.

exp(x) * cos(y) + I exp(x) * sin(y)

become

+Inf * 1 + I * (+Inf) * 0 = +Inf + I NaN

Now, mathematically, exp(1000) is a finite number, and a finite number
times zero is zero. So, is the correct answer +Inf + I Nan or +Inf +
I 0.

 

[Keith Thompson]

In general, exp(x+yi) = exp(x)(cos y + i sin y)

exp (1000+0i) = exp(1000)(cos 0 + i sin 0) = exp(1000) (1+0i) = exp(1000)


I think exp(1000) is out of range for IEEE reals.

In addition, "complex z;" is the wrong syntax; it should be
"double complex z;". The use of "complex" to mean "double complex"
is a gcc extension. (It's also an incorrect example in the original
C99 standard, which may be why gcc implemented it that way.)

 

[Ben Bacarisse]

Given

#include <stdio.h>
#include <complex.h>

int
main (void)
{
complex z;

Technically, this is not C. The three complex types are float complex,
double complex and long double complex. Since you use cexp, I presume
you intended to use double complex here.

z = 1000. + I * 0;
z = cexp(z);
printf("%lf %lf\n", creal(z), cimag(z));
return (0);
}

What is the expected output and what is correct output?

That's a very negative question -- one should always expect the correct
output!

I'd expect "inf nan" because exp(1000) is inf (in most implementations)
and since

cexp(a + ib) = exp(a)(cos(b) + i sin(b))

we get inf * (1 + i*0) and inf * 0 = nan in IEEE arithmetic.

The output I get from my implementation ("inf inf") seems curious. I
can't say if this is correct or not. I can't find anything in the C
standard that pins down what cexp should return in this case, so maybe
"inf inf" is as correct as any other result.

Maybe you expected "inf 0"? This is in some sense "closest" to the
mathematical result, but if this were the required result, cexp
implementations would probably have to special-case complex values with
a zero imaginary part. It may well be a permitted result, but I doubt
it would be required.

 

[steve]

Technically, this is not C.  The three complex types are float complex,
double complex and long double complex.  Since you use cexp, I presume
you intended to use double complex here.

Yes. Sorry about the abridged source. I was typing from
memory and rush to get out to door before rush hour.

That's a very negative question -- one should always expect the correct
output!

I don't see it as negative. What's the result of sqrt(-0.)?
The correct output is -0; while one might expect NaN.

I'd expect "inf nan" because exp(1000) is inf (in most implementations)
and since

  cexp(a + ib) = exp(a)(cos(b) + i sin(b))

we get inf * (1 + i*0) and inf * 0 = nan in IEEE arithmetic.

The output I get from my implementation ("inf inf") seems curious.  I
can't say if this is correct or not.  I can't find anything in the C
standard that pins down what cexp should return in this case, so maybe
"inf inf" is as correct as any other result.

Maybe you expected "inf 0"?  This is in some sense "closest" to the
mathematical result, but if this were the required result, cexp
implementations would probably have to special-case complex values with
a zero imaginary part.  It may well be a permitted result, but I doubt
it would be required.

That's the crux the problem. exp(a) is finite for all finite values
of a. Zero times a finite number is zero.

Note, there already is a long list of special cases in n1256.pdf.
The ones involving zero:

-- cexp(+-0 + i0) returns 1 + i0.

-- cexp(+inf + i0) returns +inf + i0.

-- cexp(+inf + iy) returns +inf * cis(y) for finite nonzero y.

-- cexp(NaN + i0) returns NaN + i0.

cis(y) = cos(y) + i sin(y)

Note the fun with y = n * pi / 2 for all n. Note also the
2nd special case is the mathematically correct answer in
the limit that exp(a) goes to exp(inf).

PS: the 1st special is also interesting in that if a
compiler supports a signed zero, then this rule tells
us that exp(+-0 + i0) = 1 + i0, but it does not address
the question of a signed zero for the imaginary part
(which may be important if you want to be on the right
Riemann sheet)

 

[Ben Bacarisse]

I don't see it as negative. What's the result of sqrt(-0.)?
The correct output is -0; while one might expect NaN.

I think we may have been using "correct" in two different ways. I took
your question to be "what is the expected output and what should a
conforming C implementation output?". With that meaning it is natural
to expect the correct output.

That's the crux the problem. exp(a) is finite for all finite values
of a. Zero times a finite number is zero.

Yes, but the problem is about unrepresentable values. exp(1000) is not
representable as an IEEE double and the result is inf. In IEEE
arithmetic inf times zero is NaN, not 0.

This rule is slightly odd. In fact it contradicts Goldberg's[1] "simple
rule":

"The rule for determining the result of an operation that has infinity
as an operand is simple: Replace infinity with a finite number x and
take the limit as x -> oo."

The rule that inf * 0 == NaN is given elsewhere in that paper but I
would have liked to see an explanation of why this seemingly obvious
rule does not apply in all cases.

Note, there already is a long list of special cases in n1256.pdf.
The ones involving zero:

-- cexp(+-0 + i0) returns 1 + i0.

-- cexp(+inf + i0) returns +inf + i0.

-- cexp(+inf + iy) returns +inf * cis(y) for finite nonzero y.

[Aside: this is not a case involving zero. If it were not for the
exclusion of y == 0 this case would have answered my question (though
maybe not yours) about what a conforming implementation should do -- it
boils down to my guess of "inf Nan".]

-- cexp(NaN + i0) returns NaN + i0.

cis(y) = cos(y) + i sin(y)

It still looks as if the C standard does not say what the answer should
be, though "inf 0" and "inf Nan" seem to be more helpful than "inf inf"
which is what my gcc + glibc gives me.

It is at times like this that I miss Dik Winter's posts.

[1] "What Every Computer Scientist Should Know About Floating-Point
Arithmetic", ACM Computing Surveys, Vol 23, No 1, March 1991.

 

[steve]

I don't see it as negative.  What's the result of sqrt(-0.)?
The correct output is -0; while one might expect NaN.

I think we may have been using "correct" in two different ways.  I took
your question to be "what is the expected output and what should a
conforming C implementation output?".  With that meaning it is natural
to expect the correct output.

That's the crux the problem.  exp(a) is finite for all finite values
of a.  Zero times a finite number is zero.

Yes, but the problem is about unrepresentable values.  exp(1000) is not
representable as an IEEE double and the result is inf.  In IEEE
arithmetic inf times zero is NaN, not 0.

This rule is slightly odd.  In fact it contradicts Goldberg's[1] "simple
rule":

  "The rule for determining the result of an operation that has infinity
  as an operand is simple: Replace infinity with a finite number x and
  take the limit as x -> oo."

This is similar to language in an old working draft of the IEEE 754
revision that I have on hand. In Sec. 6.1, it contains

"The behavior of infinity in floating-point arithmetic is
derived from the limiting cases of real arithmetic with
operands of arbitrarily large magnitude, when such a limit
exists."

However, later in 7.2, one finds

For operations producing results in floating-point format, the
default result of an invalid exception operation shall be a quiet
NaN that should provide some diagnostic information (see 6.2).
These invalid 5 exception operations are:

b) multiplication: multiplication(0, inf) or multiplication(inf, 0)

The rule that inf * 0 == NaN is given elsewhere in that paper but I
would have liked to see an explanation of why this seemingly obvious
rule does not apply in all cases.

Note, there already is a long list of special cases in n1256.pdf.
The ones involving zero:

-- cexp(+-0 + i0) returns 1 + i0.

-- cexp(+inf + i0) returns +inf + i0.

-- cexp(+inf + iy) returns +inf * cis(y) for finite nonzero y.

[Aside: this is not a case involving zero.  If it were not for the
exclusion of y == 0 this case would have answered my question (though
maybe not yours) about what a conforming implementation should do -- it
boils down to my guess of "inf Nan".]

I disgree. The case involves a zero at every y = n * pi / 2 for
integral n.

cexp(+inf + i0) = exp(+inf) * (cos(0) + i sin(0))
= exp(+inf) * (1 + i 0)
= +inf * (1 + i0)
= +inf + i0 ?

cexp(+inf + iy) = exp(+inf) * (cos(y) + i sin(y))
= exp(+inf) * (1 + i 0) for y = n * pi
= +inf * (1 + i0)
= +inf + i0 ?

I have a newer draft of IEEE 754 on some computer. But
given the above quotes, a result of inf + i NAN seems
correct to the y = n * pi case.

It still looks as if the C standard does not say what the answer should
be, though "inf 0" and "inf Nan" seem to be more helpful than "inf inf"
which is what my gcc + glibc gives me.

Can't help with gcc+glibc. Now, if you use FreeBSD, I may
be of help eventually (depends on whether my cexp[fl] make
it into libm).

Thanks for the help.

 

[Ike Naar]

cexp(+inf + iy) = exp(+inf) * (cos(y) + i sin(y))
= exp(+inf) * (1 + i 0) for y = n * pi
= +inf * (1 + i0)
= +inf + i0 ?

I have a newer draft of IEEE 754 on some computer. But
given the above quotes, a result of inf + i NAN seems
correct to the y = n * pi case.

For nonzero x one should expect sin(x) to be nonzero because
no nonzero floatingpoint number is an exact multiple of pi.

 

[Ben Bacarisse]

On Jan 20, 7:55 pm, Ben Bacarisse <[email protected]> wrote:

 In IEEE arithmetic inf times zero is NaN, not 0.

steve <[email protected]> writes:

Note, there already is a long list of special cases in n1256.pdf.
The ones involving zero:
-- cexp(+-0 + i0) returns 1 + i0.
-- cexp(+inf + i0) returns +inf + i0.
-- cexp(+inf + iy) returns +inf * cis(y) for finite nonzero y.

[Aside: this is not a case involving zero...]

I disgree. The case involves a zero at every y = n * pi / 2 for
integral n.

I don't want to argue about what "involving a zero means". The point is
that the third special case does not tell us what cexp(inf + i0) should
be -- only the previous case does that. The fact that, mathematically,
cexp(x + i*0) = cexp(x + i*2*pi) does mean that C's cexp must treat
these two in the same way. In fact, the rules you quoted show that they
must be treated differently.

[For completeness, G.6.3.1 lists another case involving zero:

-- cexp(-inf + iy) returns +0 * cis(y), for finite y]

cexp(+inf + i0) = exp(+inf) * (cos(0) + i sin(0))
= exp(+inf) * (1 + i 0)
= +inf * (1 + i0)
= +inf + i0 ?

cexp(inf + i0) = inf + i0 because Annex G says so (by the second rule
you quoted). You can't use a derivation to see what you should get. If
this derivation really did apply to that case, the last line above would
read "+inf + iNan".

cexp(+inf + iy) = exp(+inf) * (cos(y) + i sin(y))
= exp(+inf) * (1 + i 0) for y = n * pi

I think it's exp(+inf) * (-1^n + i0) but I don't think that alters your
point.

= +inf * (1 + i0)
= +inf + i0 ?

That last line should be

-1^n * inf + i*Nan

which is what the third rule mandates when y is not zero. The
derivation helps to explain the rule.

However, as Ian has pointed out, since pi can't be represented what
should be i0 will really be a small non-zero imaginary part. The result
is therefore likely to be

-1^n * inf + i*inf

with an arbitrary sign for the imaginary part. (In fact, I think you
can predict the sign of the imaginary infinity if you know which side of
pi your approximation uses.)

I have a newer draft of IEEE 754 on some computer. But
given the above quotes, a result of inf + i NAN seems
correct to the y = n * pi case.

Except that y is never = n * pi so sin(y) is unlikely to be zero. I
don't think there is anything stopping an implementation returning zero
for sin(y) when y is very close to n * pi, but I don't think it's
required and it seems unlikely in practise.

In this case, glibc's return of "-inf inf" for cexp(inf + i*pi) looks
correct (for a close approximation to pi).

It still looks as if the C standard does not say what the answer should
be, though "inf 0" and "inf Nan" seem to be more helpful than "inf inf"
which is what my gcc + glibc gives me.

Can't help with gcc+glibc. Now, if you use FreeBSD, I may
be of help eventually (depends on whether my cexp[fl] make
it into libm).

What result have you decided on for cexp(1000 + i0)?

<snip>

 

[Ben Bacarisse]

my notes on complex numbers say
def
x+iy Is in C => e^(x+iy) == e^x(cos(y) + i sin(y))

e^(+inf) (cos(0)+i sin(0)) = +inf(1+i 0)
what is "inf*0"?

NaN in IEEE floating-point arithmetic.

e^(+inf) (cos(y)+i sin(y)) = +inf*(v+i n) could be
inf+i inf, inf-i inf , -inf + i inf, -inf - i inf
depend from y from the sign of cos(y) and sin(y)

[Aside: this is not a case involving zero...]

I disgree. The case involves a zero at every y = n * pi / 2 for
integral n.

I don't want to argue about what "involving a zero means". The point is
that the third special case does not tell us what cexp(inf + i0) should

be -- only the previous case does that. The fact that, mathematically,
cexp(x + i*0) = cexp(x + i*2*pi) does mean that C's cexp must treat
these two in the same way. In fact, the rules you quoted show that they
must be treated differently.

for the above defintion they are the same
but possible i not understand well

[For completeness, G.6.3.1 lists another case involving zero:

-- cexp(-inf + iy) returns +0 * cis(y), for finite y]

cexp(+inf + i0) = exp(+inf) * (cos(0) + i sin(0))
= exp(+inf) * (1 + i 0)
= +inf * (1 + i0)
= +inf + i0 ?

inf * 0 i don't know who is

See above.

for me the last line should be
+inf + i NaN

When y = n*pi, cos(y) is not always 1. Hence the -1^n.

but it is many time i not use complex numbers
nor i'm very good on that: i see complex numbers in math in very short
"corso universitario"


Nan+i NaN
0*NaN == NaN or not?

Yes, 0 * NaN is NaN, but G.6.3.1 recommends that cexp(NaN + i0) be NaN
+ i0. There may be reasons to dispute that choice, but it's certainly a
reasonable one.

cis(y) = cos(y) + i sin(y)

It still looks as if the C standard does not say what the answer should
be, though "inf 0" and "inf Nan" seem to be more helpful than "inf inf"
which is what my gcc + glibc gives me.

Can't help with gcc+glibc. Now, if you use FreeBSD, I may
be of help eventually (depends on whether my cexp[fl] make
it into libm).

What result have you decided on for cexp(1000 + i0)?

easy
exp(1000+i 0)= exp(1000)(cos(0) + i sin(0))=exp(1000)(1 + i 0)=
=exp(1000) + i 0

But exp(1000) overflows (on most systems). cexp(1000 + i0) could be inf
+ i0 but that is not the only choice. inf * (1 + i0) could also, quite
reasonably, be inf + iNaN.

the problem is exp(inf+i0)

Annex G recommends inf + i0 as the result.

 

[Ben Bacarisse]

[I've removed the cross post to sci.maths. That's just going to
complicate matters because C's infinities are not mathematical
infinities.]

but where is the problem in define
DeF
exp(inf+i0)=inf+i 0
?

There is no problem with that definition. In fact, annex G recommends
that that be result of cexp(inf + i0). It was you that said this case
was a problem.

so in all case
if x!=NaN => exp(x+i0)=exp(x)+i 0
exp(NaN+i 0)=NaN + i NaN

Annex G recommends NaN + i0 here. In other words, it suggests that
cimag(cexp(x + i0)) == 0 for all x (finite, infinite and NaN).

None the less, your suggestion (that cexp(NaN + i0) be NaN + i*NaN) is
quite reasonable. In fact it is what glibc does.

have u done some course in complex analisys?

Yes, but I don't see how that matters.