Re: plotting even function

[Mike Miller]

On Sun, 20 Mar 2005 address@hidden wrote:

> OK, if I understood correctly the very interesting discussion, in
> principle (a^b)^c==(a^c)^b (a is a real scalar) if, and only if
> b*c is a rational number.

I think, mathematically speaking, (a^b)^c==(a^c)^b when b and c are real, 
not just for rational b*c.  Computationally, however, the two may differ 
greatly depending on choice of algorithms.  As noted earlier, 1/3 cannot 
be represented as a finite binary expansion without some imprecision. 
This means that the order matters:  ((-1)^2) = 1 and ((-1)^2)^(1/3) = 1, 
approximately, but (-1)^(1/3) is not -1, and it is complex, so 
((-1)^(1/3))^2 is also complex, and not equal to 1.


> As octave treats evry inputed value as a double defined to some finite 
> precision, and as within this uncertainty there are both numbers 
> commesurate with each other and non-commensurate ones, it has to make 
> some arbitrary decision.

I don't know that the decision should be called "arbitrary."  It seems 
like a good set of choices to me.  It probably uses De Moivre's theorem:

http://scholar.hw.ac.uk/site/maths/topic17.asp

This also provides consistency with things like...

log(-1)
exp(pi*i)


> Now consider the folowing simple example:
> 588~> x=[-10:10];
> 588~> y=x.^(1/3);
> 589~> plot(x,y) % the plot is asymetric with respect to the
> origin
> 590~> floor(1)
> ans = 1
> 591~> floor(3)
> ans = 3
> 592~> y=x.^(floor(1)/floor(3));
> 589~> plot(x,y) % I get the same plot, although floor(1) and
> floor(3) are integers.
>
> It's a pity, because that mighty have been a way to obtain from the 
> software the expected behaviour.

When you're working with x^(a/b), and a,b are scalar integers, I think you 
are expecting something like this:

when a,b both odd:
sign(x).*(abs(x).^(a/b))

when a is even:
abs(x).^(a/b)

otherwise (a odd, b even):
x.^(a/b)   # the usual behavior

altogether using indicators...

abs(rem(a,2))*abs(rem(b,2))*sign(x).*(abs(x).^(a/b)) + \
(1-abs(rem(a,2)))*abs(x).^(a/b) + \
abs(rem(a,2))*(1-abs(rem(b,2)))*x.^(a/b)

...on a single line:

abs(rem(a,2))*abs(rem(b,2))*sign(x).*(abs(x).^(a/b)) + 
(1-abs(rem(a,2)))*abs(x).^(a/b) + abs(rem(a,2))*(1-abs(rem(b,2)))*x.^(a/b)

That might not be what everyone expects for things like (-8)^(2/3), and it 
isn't what Octave currently does, nor is it what MATLAB does, so I guess 
you should just code it like that by hand when needed.  I think it works 
for integers a,b up to about 16 digits long.  Does that seem like a good 
strategy?


> Cheers and thanks to all hte participants for the pleasant hour I spent 
> recalling long forgotten basic math, Avraham

Same here, believe me.  I don't use much of this every day.

Mike

--
Michael B. Miller, Ph.D.
Assistant Professor
Division of Epidemiology and Community Health
and Institute of Human Genetics
University of Minnesota
http://taxa.epi.umn.edu/~mbmiller/



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