«

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

On Jan 20, 4:32*pm, Out of the China Blue <chine.b...@yahoo.com>
wrote:
> In article <5cd72d26-ce0b-4503-b5f6-861f916a5...@q8g2000prm.googlegroups.com>,
>
> *steve <kar...@comcast.net> wrote:
> > 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);
> > }
>
> 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.

--
steve

Out of the China Blue <> writes:
> In article <5cd72d26-ce0b-4503-b5f6->,
> steve <> wrote:
>
>> 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);
>> }
>
> 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.)

--
Keith Thompson (The_Other_Keith) kst- <http://www.ghoti.net/~kst>
Nokia
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

steve <> writes:

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

--
Ben.

On Jan 20, 5:38*pm, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
> steve <kar...@comcast.net> writes:
> > 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.

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

> > * *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 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)
--
steve

steve <> writes:

> On Jan 20, 5:38Â*pm, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
>> steve <kar...@comcast.net> writes:
>> > Given
>>
>> > #include <stdio.h>
>> > #include <complex.h>
>>
>> > int
>> > main (void)
>> > {
>> > Â* Â*complex z;
<snip replace complex with double complex>
>> > Â* Â*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 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.

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

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.

--
Ben.

On Jan 20, 7:55*pm, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
> steve <kar...@comcast.net> writes:
> > On Jan 20, 5:38*pm, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
> >> steve <kar...@comcast.net> writes:
> >> > Given
>
> >> > #include <stdio.h>
> >> > #include <complex.h>
>
> >> > int
> >> > main (void)
> >> > {
> >> > * *complex z;
>
> <snip replace complex with double complex>
>
> >> > * *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 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.
>
>
>
> >> 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.
>
> 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.

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

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.

--
steve

On 2011-01-21, steve <> wrote:
> 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.

steve <> writes:

> On Jan 20, 7:55Â*pm, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
<snip>
>> Â*In IEEE arithmetic inf times zero is NaN, not 0.
<snip>
>> steve <kar...@comcast.net> writes:
<snip>
>> > 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...]
<snip>

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

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

"io_x" <> writes:

> "Ben Bacarisse" <> ha scritto nel messaggio
> news:0.6d02c8542229707911d7.20110121143620GMT.878v ...
>> steve <> writes:
>>
>>> On Jan 20, 7:55 pm, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
>> <snip>
>>>> In IEEE arithmetic inf times zero is NaN, not 0.
>> <snip>
>>>> steve <kar...@comcast.net> writes:
>> <snip>
>>>> > Note, there already is a long list of special cases in n1256.pdf.
>
>
> my notes on complex numbers say
> def
> x+iy Is in C => e^(x+iy) == e^x(cos(y) + i sin(y))
>
>>>> > The ones involving zero:
>>>> > -- cexp(+-0 + i0) returns 1 + i0.
> e^0(cos(0)+i sin(0)) = 1 + i0.
>>>> > -- cexp(+inf + i0) returns +inf + i0.
> e^(+inf) (cos(0)+i sin(0)) = +inf(1+i 0)
> what is "inf*0"?

NaN in IEEE floating-point arithmetic.

>>>> > -- cexp(+inf + iy) returns +inf * cis(y) for finite nonzero y.
>
> 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...]
>> <snip>
>>
>>> 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.

>> 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
>
> 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"
>
>> 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).
>>
>>>> > -- cexp(NaN + i0) returns NaN + i0.
>
> 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.