Re: plotting even function

[Geraint Paul Bevan]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Mike Miller wrote:

> I am not a mathematician, but I believe that you are missing something.
> The equation y^3 + 1 = 0 has three solutions, as noted above.  The
> equation y^3 - 1 = 0 also has three solutions:
>
> y = 1
> y = -0.5-i*sqrt(3)/2
> y = -0.5+i*sqrt(3)/2
>
> How does Octave "know what we want" (1) in the latter case, but it does
> not know what we want (-1) in the former case?  I think the answer has
> to do with numerical analysis.  The binary representation of 1/3 is never
> exactly 1/3, but it is only at the exact value of 1/3 that (-1)^(1/3)
> equals -1.  This is not a problem for (+1)^(1/3).
>
> For (-1)^(333/1000), Octave gives a very similar answer to what it gives
> for (-1)^(1/3), which seems appropriate.  Octave doesn't seem to know or
> care that there are 1000 possible solutions in the complex plane!
>
> Mike


I believe that Octave uses the standard pow function within the system's
math library to perform this operation. The cause of the difference in
behaviour can be found in the Octave sources, xpow.cc:

octave_value
xpow (double a, double b)
{
  if (a < 0.0 && static_cast<int> (b) != b)
    {
      // XXX FIXME XXX -- avoid apparent GNU libm bug by converting
      // A and B to complex instead of just A.

      Complex atmp (a);
      Complex btmp (b);

      return pow (atmp, btmp);
    }
  else
    return pow (a, b);
}

If the base is not negative, or the exponent is an integer, Octave uses
the standard double version of pow. This routine never finds a complex
solution -- it either finds a real one or fails:

a.cc:

#include <cmath>
#include <iostream>

int main(void) {
  std::cout << "+1^2.0\t" << pow(+1.0,2.0) << std::endl;
  std::cout << "-1^2.0\t" << pow(-1.0,2.0) << std::endl;
  std::cout << "+1^0.5\t" << pow(+1.0,0.5) << std::endl;
  std::cout << "-1^0.5\t" << pow(-1.0,0.5) << std::endl;
  return 0;
}

$ g++ a.cc && ./a.out
+1^2.0  1
- -1^2.0        1
+1^0.5  1
- -1^0.5        nan

As you can see from the results above, the double version of pow cannot
cope with a negative base and non-integer exponent. Therefore in this
case Octave uses the Complex (i.e. complex<double>) version of pow. It
is this routine which is finding complex solutions for negative x
instead of the real root that would be required for symmetry about the
y-axis.

- --
Geraint Bevan
http://homepage.ntlworld.com/geraint.bevan

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.4 (GNU/Linux)

iEYEARECAAYFAkI9xogACgkQcXV3N50QmNPprACfZA5fdTDy2IHLef4XmUeyd+yW
4fkAniXoa1L3VlyCxC5JyiWNsM1UBjua
=eDGK
-----END PGP SIGNATURE-----



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

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