octave_packages/m/specfun/nthroot.m

   1 ## Copyright (C) 2004-2012 Paul Kienzle
   2 ## Copyright (C) 2010 VZLU Prague
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19 ##
  20 ## Original version by Paul Kienzle distributed as free software in the
  21 ## public domain.
  22
  23 ## -*- texinfo -*-
  24 ## @deftypefn {Function File} {} nthroot (@var{x}, @var{n})
  25 ##
  26 ## Compute the n-th root of @var{x}, returning real results for real
  27 ## components of @var{x}.  For example:
  28 ##
  29 ## @example
  30 ## @group
  31 ## nthroot (-1, 3)
  32 ## @result{} -1
  33 ## (-1) ^ (1 / 3)
  34 ## @result{} 0.50000 - 0.86603i
  35 ## @end group
  36 ## @end example
  37 ##
  38 ## @var{x} must have all real entries.  @var{n} must be a scalar.
  39 ## If @var{n} is an even integer and @var{X} has negative entries, an
  40 ## error is produced.
  41 ## @seealso{realsqrt, sqrt, cbrt}
  42 ## @end deftypefn
  43
  44 function y = nthroot (x, n)
  45
  46   if (nargin != 2)
  47     print_usage ();
  48   endif
  49
  50   if (any (iscomplex (x(:))))
  51     error ("nthroot: X must not contain complex values");
  52   endif
  53
  54   if (! isscalar (n) || n == 0)
  55     error ("nthroot: N must be a nonzero scalar");
  56   endif
  57
  58   if (n == 3)
  59     y = cbrt (x);
  60   elseif (n == -3)
  61     y = 1 ./ cbrt (x);
  62   elseif (n < 0)
  63     y = 1 ./ nthroot (x, -n);
  64   else
  65     ## Compute using power.
  66     if (n == fix (n) && mod (n, 2) == 1)
  67       y = abs (x) .^ (1/n) .* sign (x);
  68     elseif (any (x(:) < 0))
  69       error ("nthroot: if X contains negative values, N must be an odd integer");
  70     else
  71       y = x .^ (1/n);
  72     endif
  73
  74     if (finite (n) && n > 0 && n == fix (n))
  75       ## Correction.
  76       y = ((n-1)*y + x ./ (y.^(n-1))) / n;
  77       y = merge (finite (y), y, x);
  78     endif
  79
  80   endif
  81
  82 endfunction
  83
  84 %!assert (nthroot(-32,5), -2);
  85 %!assert (nthroot(81,4), 3);
  86 %!assert (nthroot(Inf,4), Inf);
  87 %!assert (nthroot(-Inf,7), -Inf);
  88 %!assert (nthroot(-Inf,-7), 0);
  89
  90 %% Test input validation
  91 %!error (nthroot ())
  92 %!error (nthroot (1))
  93 %!error (nthroot (1,2,3))
  94 %!error (nthroot (1+j,2))
  95 %!error (nthroot (1,[1 2]))
  96 %!error (nthroot (1,0))
  97 %!error (nthroot (-1,2))
  98