# HG changeset patch # User Jaroslav Hajek # Date 1268750015 -3600 # Node ID 976e76b77fa061b169524990018c810e56eee893 # Parent 2a8b1db1e2ca81b011bf59dbb3e362689a0041e3 rewrite nth_root, move to specfun diff -r 2a8b1db1e2ca -r 976e76b77fa0 scripts/ChangeLog --- a/scripts/ChangeLog Tue Mar 16 15:16:32 2010 +0100 +++ b/scripts/ChangeLog Tue Mar 16 15:33:35 2010 +0100 @@ -1,3 +1,10 @@ +2010-03-16 Jaroslav Hajek + + * general/nthroot.m: Remove. + * general/module.mk: Update. + * specfun/nthroot.m: New source. + * specfun/module.mk: Update. + 2010-03-16 Jaroslav Hajek * miscellaneous/intwarning.m: Deprecate. diff -r 2a8b1db1e2ca -r 976e76b77fa0 scripts/general/module.mk --- a/scripts/general/module.mk Tue Mar 16 15:16:32 2010 +0100 +++ b/scripts/general/module.mk Tue Mar 16 15:33:35 2010 +0100 @@ -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 \ diff -r 2a8b1db1e2ca -r 976e76b77fa0 scripts/general/nthroot.m --- a/scripts/general/nthroot.m Tue Mar 16 15:16:32 2010 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -## Copyright (C) 2004, 2006, 2007, 2009 Paul Kienzle -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . -## -## Original version by Paul Kienzle distributed as free software in the -## public domain. - -## -*- texinfo -*- -## @deftypefn {Function File} {} nthroot (@var{x}, @var{n}) -## -## Compute the n-th root of @var{x}, returning real results for real -## components of @var{x}. For example -## -## @example -## @group -## nthroot (-1, 3) -## @result{} -1 -## (-1) ^ (1 / 3) -## @result{} 0.50000 - 0.86603i -## @end group -## @end example -## -## @end deftypefn - -function y = nthroot (x, m) - - if (nargin != 2) - print_usage (); - endif - - y = x.^(1./m); - - if (isscalar (x)) - x *= ones (size (m)); - endif - - if (isscalar (m)) - m *= ones (size (x)); - endif - - idx = (mod (m, 2) == 1 & imag (x) == 0 & x < 0); - - if (any (idx(:))) - y(idx) = -(-x(idx)).^(1./m(idx)); - 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)); diff -r 2a8b1db1e2ca -r 976e76b77fa0 scripts/specfun/module.mk --- a/scripts/specfun/module.mk Tue Mar 16 15:16:32 2010 +0100 +++ b/scripts/specfun/module.mk Tue Mar 16 15:33:35 2010 +0100 @@ -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 \ diff -r 2a8b1db1e2ca -r 976e76b77fa0 scripts/specfun/nthroot.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/specfun/nthroot.m Tue Mar 16 15:33:35 2010 +0100 @@ -0,0 +1,83 @@ +## Copyright (C) 2004, 2006, 2007, 2009 Paul Kienzle +## Copyright (C) 2010 VZLU Prague +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . +## +## Original version by Paul Kienzle distributed as free software in the +## public domain. + +## -*- texinfo -*- +## @deftypefn {Function File} {} nthroot (@var{x}, @var{n}) +## +## Compute the n-th root of @var{x}, returning real results for real +## components of @var{x}. For example +## +## @example +## @group +## nthroot (-1, 3) +## @result{} -1 +## (-1) ^ (1 / 3) +## @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, n) + + if (nargin != 2) + print_usage (); + endif + + if (! isscalar (n)) + error ("nthroot: n must be a nonzero scalar"); + endif + + 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 (finite (n) && n > 0 && n == round (n)) + ## Correction. + y = ((n-1)*y + x ./ (y.^(n-1))) / n; + y = merge (finite (y), y, x); + endif + + endif + +endfunction + +%!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);