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 5 files changed, 91 insertions(+), 71 deletions(-) [+]
line wrap: on
line diff
--- 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  <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	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 \
--- 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
-## <http://www.gnu.org/licenses/>.
-##
-## 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));
--- 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 \
--- /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
+## <http://www.gnu.org/licenses/>.
+##
+## 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);