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