From bug-octave-request at bevo dot che dot wisc dot edu  Thu Sep  7 15:51:38 2000
Subject: log(-Inf) + exp(Inf+pi*i) problems
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Rolf Fabian <fabian at tu-cottbus dot de>
Cc: "'bug-octave UWISC'" <bug-octave at bevo dot che dot wisc dot edu>
Date: Thu, 7 Sep 2000 15:50:10 -0500 (CDT)

On  6-Sep-2000, Rolf Fabian <fabian at tu-cottbus dot de> wrote:

| again, I don't know, whether the the following bugs have disappeared
| in the current development version 2.1.31 ?
| 
| --  A ----
| :> version	|- ans = 2.1.14		#OS/2 (ported by Klaus Gebhardt)
| :> log(-Inf)	|- ans = Inf + 3i	<<
| :> imag(ans)-pi	|- ans = 0		<<	conclusion 3==pi?
| This bug is not severe, but may be irritating.
| Looks like a formatting problem.

What values do you have set for output_precision and
output_max_field_width?  On my system, I see

  octave:4> log(-Inf)
  ans = Inf + 3.14i

| --  B ----
| :1> exp(Inf-pi*i) |- ans = -Inf - Infi	<< BUG
| :2> exp(Inf+pi*i) |- ans = -Inf + Infi	<< BUG
| 
| Both results are incorrect, because complex exponential function
| has a period of P=2*pi*i, so that
| 	exp(x-pi*i)   and   exp(x+pi*i)
| are expected to return the SAME result, because both arguments differ by P,
| even if x is Inf .

Hmm.  On my system I see

  octave:6> exp (Inf + i*pi)
  ans = NaN + NaNi

which is also wrong.

Octave uses the exp function from the standard C++ library (typically,
this is the GNU implementation).  In the sources for the complex class
that is distributed with gcc 2.95, the result is computed using (with
extraneous template declarations removed):

  complex
  exp (const complex& x)
  {
    return polar (exp (real (x)), imag (x));
  }

and

  complex
  polar (double r, double t)
  {
    return complex (r * cos (t), r * sin (t));
  }

so, the reason that you see a result of -Inf +/- Inf*i from this is
that sin (pi) does not return exactly zero.  Even if it did, the
result of r * 0 would be NaN, not 0.  This is yet another case of
"what should happen when the imaginary part of a complex number
becomes zero and that complex number is multiplied by Infinity?"

Now, why is the result even worse (NaN + NaN*i) on my system?  It
seems that there are inline math functions that are available to the
compiler which are simply wrong, and the inline version of the exp
function returns NaN for exp(Inf).  Octave tries to deal with this by
defining __NO_MATH_INLINES (which I have defined when building Octave
on my system), but because the code for the double complex class is
compiled into the standard C++ library (apparently using the inline
functions), it is still computed incorrectly here.

I'm afraid that the tools just aren't up to the task yet.

I believe that the fact that the gcc/libc inline math functions are
buggy has been reported to the gcc/libc maintainers, but I don't think
the problem has been fixed yet.  Please correct me if I'm wrong, and
tell me which version I should be using.  (Sometimes I think that
compiler developers never do and don't care much about numerical
computations.)

I think that fixing all of this on all systems would be hard.

Even if it is fixed for Octave by using a special complex class, the
results would still be wrong in compiled C, C++, or Fortan code (for
example) if that code used the buggy system libraries.

jwe



-----------------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.che.wisc.edu/octave/octave.html
How to fund new projects:  http://www.che.wisc.edu/octave/funding.html
Subscription information:  http://www.che.wisc.edu/octave/archive.html
-----------------------------------------------------------------------