cexp implementation questions

Discussion in 'C Programming' started by steve, Jan 20, 2011.

  1. steve

    steve Guest

    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
    steve, Jan 20, 2011
    #1
    1. Advertising

  2. steve

    steve Guest

    On Jan 20, 4:32 pm, Out of the China Blue <>
    wrote:
    > In article <>,
    >
    >  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.
    >


    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
    steve, Jan 21, 2011
    #2
    1. Advertising

  3. Out of the China Blue <> writes:
    > In article <>,
    > 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) <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"
    Keith Thompson, Jan 21, 2011
    #3
  4. 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.
    Ben Bacarisse, Jan 21, 2011
    #4
  5. steve

    steve Guest

    On Jan 20, 5:38 pm, Ben Bacarisse <> wrote:
    > 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.


    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, Jan 21, 2011
    #5
  6. steve <> writes:

    > On Jan 20, 5:38 pm, Ben Bacarisse <> wrote:
    >> steve <> 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.
    Ben Bacarisse, Jan 21, 2011
    #6
  7. steve

    steve Guest

    On Jan 20, 7:55 pm, Ben Bacarisse <> wrote:
    > steve <> writes:
    > > On Jan 20, 5:38 pm, Ben Bacarisse <> wrote:
    > >> steve <> 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
    steve, Jan 21, 2011
    #7
  8. steve

    Ike Naar Guest

    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.
    Ike Naar, Jan 21, 2011
    #8
  9. steve <> writes:

    > On Jan 20, 7:55 pm, Ben Bacarisse <> wrote:

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

    <snip>
    >> steve <> 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.
    Ben Bacarisse, Jan 21, 2011
    #9
  10. "io_x" <> writes:

    > "Ben Bacarisse" <> ha scritto nel messaggio
    > news:...
    >> steve <> writes:
    >>
    >>> On Jan 20, 7:55 pm, Ben Bacarisse <> wrote:

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

    >> <snip>
    >>>> steve <> 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.
    Ben Bacarisse, Jan 21, 2011
    #10
  11. [I've removed the cross post to sci.maths. That's just going to
    complicate matters because C's infinities are not mathematical
    infinities.]

    "io_x" <> writes:

    > "io_x" <> ha scritto nel messaggio
    > news:4d39c3db$0$19264$...
    >> "Ben Bacarisse" <> ha scritto nel messaggio
    >> news:...
    >>> steve <> writes:
    >>> 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
    >>
    >> the problem is exp(inf+i0)

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

    --
    Ben.
    Ben Bacarisse, Jan 21, 2011
    #11
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. John Rowell
    Replies:
    0
    Views:
    549
    John Rowell
    Apr 22, 2004
  2. Replies:
    2
    Views:
    352
    Ian Collins
    Nov 13, 2007
  3. Michael Tsang
    Replies:
    32
    Views:
    1,103
    Richard Bos
    Mar 1, 2010
  4. Michael Tsang
    Replies:
    54
    Views:
    1,189
    Phil Carmody
    Mar 30, 2010
  5. sanket
    Replies:
    7
    Views:
    996
    Tsung
    Nov 3, 2011
Loading...

Share This Page