Mercurial > hg > octave-lyh
changeset 10415:976e76b77fa0
rewrite nth_root, move to specfun
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 16 Mar 2010 15:33:35 +0100 |
parents | 2a8b1db1e2ca |
children | f06ede00fd8f |
files | scripts/ChangeLog scripts/general/module.mk scripts/general/nthroot.m scripts/specfun/module.mk scripts/specfun/nthroot.m |
diffstat | 4 files changed, 40 insertions(+), 20 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,10 @@ +2010-03-16 Jaroslav Hajek <highegg@gmail.com> + + * general/nthroot.m: Remove. + * general/module.mk: Update. + * specfun/nthroot.m: New source. + * specfun/module.mk: Update. + 2010-03-16 Jaroslav Hajek <highegg@gmail.com> * miscellaneous/intwarning.m: Deprecate.
--- a/scripts/general/module.mk +++ b/scripts/general/module.mk @@ -54,7 +54,6 @@ general/nargchk.m \ general/nargoutchk.m \ general/nextpow2.m \ - general/nthroot.m \ general/num2str.m \ general/perror.m \ general/pol2cart.m \
--- a/scripts/specfun/module.mk +++ b/scripts/specfun/module.mk @@ -11,6 +11,7 @@ specfun/isprime.m \ specfun/legendre.m \ specfun/nchoosek.m \ + specfun/nthroot.m \ specfun/perms.m \ specfun/pow2.m \ specfun/primes.m \
rename from scripts/general/nthroot.m rename to scripts/specfun/nthroot.m --- a/scripts/general/nthroot.m +++ b/scripts/specfun/nthroot.m @@ -1,4 +1,5 @@ ## Copyright (C) 2004, 2006, 2007, 2009 Paul Kienzle +## Copyright (C) 2010 VZLU Prague ## ## This file is part of Octave. ## @@ -33,38 +34,50 @@ ## @result{} 0.50000 - 0.86603i ## @end group ## @end example +## +## @var{n} must be a scalar. If @var{n} is not an even integer and @var{X} has +## negative entries, an error is produced. ## ## @end deftypefn -function y = nthroot (x, m) +function y = nthroot (x, n) if (nargin != 2) print_usage (); endif - - y = x.^(1./m); - if (isscalar (x)) - x *= ones (size (m)); + if (! isscalar (n)) + error ("nthroot: n must be a nonzero scalar"); endif - if (isscalar (m)) - m *= ones (size (x)); - endif - - idx = (mod (m, 2) == 1 & imag (x) == 0 & x < 0); + if (n == 3) + y = cbrt (x); + elseif (n == -3) + y = 1 ./ cbrt (x); + elseif (n < 0) + y = 1 ./ nthroot (x, -n); + else + ## Compute using power. + if (n == round (n) && mod (n, 2) == 1) + y = abs (x) .^ (1/n) .* sign (x); + elseif (any (x(:) < 0)) + error ("nthroot: if x contains negative values, n must be an odd integer"); + else + y = x .^ (1/n); + endif - if (any (idx(:))) - y(idx) = -(-x(idx)).^(1./m(idx)); - endif + if (finite (n) && n > 0 && n == round (n)) + ## Correction. + y = ((n-1)*y + x ./ (y.^(n-1))) / n; + y = merge (finite (y), y, x); + endif - ## If result is all real, make sure it looks real - if (all (imag (y) == 0)) - y = real (y); endif endfunction -%!assert(nthroot(-1,[3,-3]), [-1,-1],eps); -%!assert(nthroot([-1,1],[3.1,-3]), [-1,1].^(1./[3.1,-3])); -%!assert(nthroot([-1+1i,-1-1i],3), [-1+1i,-1-1i].^(1/3)); +%!assert(nthroot(-32,5), -2); +%!assert(nthroot(81,4), 3); +%!assert(nthroot(Inf,4), Inf); +%!assert(nthroot(-Inf,7), -Inf); +%!assert(nthroot(-Inf,-7), 0);