/usr/share/octave/packages/specfun-1.1.0/lambertw.m is in octave-specfun 1.1.0-1build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 | %% Copyright (C) 1998 by Nicol N. Schraudolph
%%
%% This program is free software; you can redistribute and/or
%% modify it under the terms of the GNU General Public
%% License as published by the Free Software Foundation;
%% either version 3, or (at your option) any later version.
%%
%% This program 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 this software; see the file COPYING. If not,
%% see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{x} = } lambertw (@var{z})
## @deftypefnx {Function File} {@var{x} = } lambertw (@var{z}, @var{n})
## Compute the Lambert W function of @var{z}.
##
## This function satisfies W(z).*exp(W(z)) = z, and can thus be used to express
## solutions of transcendental equations involving exponentials or logarithms.
##
## @var{n} must be integer, and specifies the branch of W to be computed;
## W(z) is a shorthand for W(0,z), the principal branch. Branches
## 0 and -1 are the only ones that can take on non-complex values.
##
## If either @var{n} or @var{z} are non-scalar, the function is mapped to each
## element; both may be non-scalar provided their dimensions agree.
##
## This implementation should return values within 2.5*eps of its
## counterpart in Maple V, release 3 or later. Please report any
## discrepancies to the author, Nici Schraudolph <schraudo@@inf.ethz.ch>.
##
## For further details, see:
##
## Corless, Gonnet, Hare, Jeffrey, and Knuth (1996), `On the Lambert
## W Function', Advances in Computational Mathematics 5(4):329-359.
## @end deftypefn
%% Author: Nicol N. Schraudolph <schraudo@inf.ethz.ch>
%% Version: 1.0
%% Created: 07 Aug 1998
%% Keywords: Lambert W Omega special transcendental function
function w = lambertw(b,z)
if (nargin == 1)
z = b;
b = 0;
else
%% some error checking
if (nargin != 2)
print_usage;
else
if (any(round(real(b)) != b))
usage('branch number for lambertw must be integer')
end
end
end
%% series expansion about -1/e
%
% p = (1 - 2*abs(b)).*sqrt(2*e*z + 2);
% w = (11/72)*p;
% w = (w - 1/3).*p;
% w = (w + 1).*p - 1
%
% first-order version suffices:
%
w = (1 - 2*abs(b)).*sqrt(2*e*z + 2) - 1;
%% asymptotic expansion at 0 and Inf
%
v = log(z + ~(z | b)) + 2*pi*I*b;
v = v - log(v + ~v);
%% choose strategy for initial guess
%
c = abs(z + 1/e);
c = (c > 1.45 - 1.1*abs(b));
c = c | (b.*imag(z) > 0) | (~imag(z) & (b == 1));
w = (1 - c).*w + c.*v;
%% Halley iteration
%
for n = 1:10
p = exp(w);
t = w.*p - z;
f = (w != -1);
t = f.*t./(p.*(w + f) - 0.5*(w + 2.0).*t./(w + f));
w = w - t;
if (abs(real(t)) < (2.48*eps)*(1.0 + abs(real(w)))
&& abs(imag(t)) < (2.48*eps)*(1.0 + abs(imag(w))))
return
end
end
warning('iteration limit reached, result of lambertw may be inaccurate');
end
|