Mercurial > octave-nkf

comparison scripts/specfun/nthroot.m @ 10415:976e76b77fa0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
rewrite nth_root, move to specfun
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 16 Mar 2010 15:33:35 +0100
parents scripts/general/nthroot.m@1bf0ce0930be
children fbd7843974fa
comparison
equal deleted inserted replaced
10414:2a8b1db1e2ca 10415:976e76b77fa0
1 ## Copyright (C) 2004, 2006, 2007, 2009 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{n} must be a scalar. If @var{n} is not an even integer and @var{X} has
39 ## negative entries, an error is produced.
40 ##
41 ## @end deftypefn
42
43 function y = nthroot (x, n)
44
45 if (nargin != 2)
46 print_usage ();
47 endif
48
49 if (! isscalar (n))
50 error ("nthroot: n must be a nonzero scalar");
51 endif
52
53 if (n == 3)
54 y = cbrt (x);
55 elseif (n == -3)
56 y = 1 ./ cbrt (x);
57 elseif (n < 0)
58 y = 1 ./ nthroot (x, -n);
59 else
60 ## Compute using power.
61 if (n == round (n) && mod (n, 2) == 1)
62 y = abs (x) .^ (1/n) .* sign (x);
63 elseif (any (x(:) < 0))
64 error ("nthroot: if x contains negative values, n must be an odd integer");
65 else
66 y = x .^ (1/n);
67 endif
68
69 if (finite (n) && n > 0 && n == round (n))
70 ## Correction.
71 y = ((n-1)*y + x ./ (y.^(n-1))) / n;
72 y = merge (finite (y), y, x);
73 endif
74
75 endif
76
77 endfunction
78
79 %!assert(nthroot(-32,5), -2);
80 %!assert(nthroot(81,4), 3);
81 %!assert(nthroot(Inf,4), Inf);
82 %!assert(nthroot(-Inf,7), -Inf);
83 %!assert(nthroot(-Inf,-7), 0);