1/*
2 * lambertw - Lambert's W-function
3 *
4 * Copyright (C) 2013,2021 Christoph Zurnieden
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	2013/08/11 01:31:28
21 * File existed as early as:	2013
22 */
23
24
25static resource_debug_level;
26resource_debug_level = config("resource_debug", 0);
27
28
29/*
30
31     R. M. Corless and G. H. Gonnet and D. E. G. Hare and D. J. Jeffrey and
32     D. E. Knuth, "On the Lambert W Function", Advances n Computational
33     Mathematics, 329--359, (1996)
34     http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.112.6117
35
36     D. J. Jeffrey, D. E. G. Hare, R. M. Corless, "Unwinding the branches of the
37     Lambert W function", The Mathematical Scientist, 21, pp 1-7, (1996)
38     http://www.apmaths.uwo.ca/~djeffrey/Offprints/wbranch.pdf
39
40     Darko Verebic, "Having Fun with Lambert W(x) Function"
41     arXiv:1003.1628v1, March 2010, http://arxiv.org/abs/1003.1628
42
43     Winitzki, S. "Uniform Approximations for Transcendental Functions",
44     In Part 1 of Computational Science and its Applications - ICCSA 2003,
45     Lecture Notes in Computer Science, Vol. 2667, Springer-Verlag,
46     Berlin, 2003, 780-789. DOI 10.1007/3-540-44839-X_82
47     A copy may be found by Google.
48
49
50*/
51static true = 1;
52static false = 0;
53
54/* Branch 0, Winitzki (2003) , the well known Taylor series*/
55define __CZ__lambertw_0(z,eps){
56  local a=2.344e0, b=0.8842e0, c=0.9294e0, d=0.5106e0, e=-1.213e0;
57  local y=sqrt(2*exp(1)*z+2);
58  return (2*ln(1+b*y)-ln(1+c*ln(1+d*y))+e)/(1+1/(2*ln(1+b*y)+2*a));
59}
60/* branch -1 */
61define __CZ__lambertw_m1(z,eps){
62  local wn k;
63  /* Cut-off found in Maxima */
64  if(z < 0.3) return __CZ__lambertw_app(z,eps);
65  wn = z;
66  /*  Verebic (2010) eqs. 16-18*/
67  for(k=0;k<10;k++){
68    wn = ln(-z)-ln(-wn);
69  }
70  return wn;
71}
72
73/*
74  generic approximation
75
76  series for 1+W((z-2)/(2 e))
77
78  Corless et al (1996) (4.22)
79  Verebic (2010) eqs. 35-37; more coefficients given at the end of sect. 3.1
80  or online
81     http://www.wolframalpha.com/input/?
82       i=taylor+%28+1%2Bproductlog%28+%28z-2%29%2F%282*e%29+%29+%29
83  or by using the function lambertw_series_print() after running
84  lambertw_series(z,eps,branch,terms) at least once with the wanted number of
85  terms and z = 1 (which might throw an error because the series will not
86  converge in anybodies lifetime for something that far from the branch point).
87
88
89*/
90define __CZ__lambertw_app(z,eps){
91  local  b0=-1, b1=1, b2=-1/3, b3=11/72;
92  local  y=sqrt(2*exp(1)*z+2);
93  return b0 + ( y * (b1 + (y * (b2 + (b3 * y)))));
94}
95
96static __CZ__Ws_a;
97static __CZ__Ws_c;
98static __CZ__Ws_len=0;
99
100define lambertw_series_print(){
101  local k;
102  for(k=0;k<__CZ__Ws_len;k++){
103    print num(__CZ__Ws_c[k]):"/":den(__CZ__Ws_c[k]):"*p^":k;
104  }
105}
106
107/*
108   The series is fast but only if _very_ close to the branch point
109   The exact branch must be given explicitly, e.g.:
110
111; lambertw(-exp(-1)+.001)-lambertw_series(-exp(-1)+.001,epsilon()*1e-10,0)
112	-0.14758879113205794065490184399030194122136720202792-
113                  0.00000000000000000000000000000000000000000000000000i
114; lambertw(-exp(-1)+.001)-lambertw_series(-exp(-1)+.001,epsilon()*1e-10,1)
115	0.00000000000000000000000000000000000000000000000000-
116                  0.00000000000000000000000000000000000000000000000000i
117 */
118define lambertw_series(z,eps,branch,terms){
119  local k l limit tmp sum A C P PP epslocal;
120  if(!isnull(terms))
121    limit = terms;
122  else
123    limit = 100;
124
125  if(isnull(eps))
126    eps = epsilon(epsilon()*1e-10);
127  epslocal = epsilon(eps);
128
129  P = sqrt(2*(exp(1)*z+1));
130  if(branch != 0) P = -P;
131  tmp=0;sum=0;PP=P;
132
133    __CZ__Ws_a   = mat[limit+1];
134    __CZ__Ws_c   = mat[limit+1];
135    __CZ__Ws_len = limit;
136   /*
137      c0 = -1; c1 = 1
138      a0 =  2; a1 =-1
139    */
140    __CZ__Ws_c[0] = -1; __CZ__Ws_c[1] =  1;
141    __CZ__Ws_a[0] =  2; __CZ__Ws_a[1] = -1;
142    sum += __CZ__Ws_c[0];
143    sum += __CZ__Ws_c[1] * P;
144    PP *= P;
145    for(k=2;k<limit;k++){
146      for(l=2;l<k;l++){
147         __CZ__Ws_a[k] += __CZ__Ws_c[l]*__CZ__Ws_c[k+1-l];
148      }
149
150      __CZ__Ws_c[k] = (k-1) * (  __CZ__Ws_c[k-2]/2
151                            +__CZ__Ws_a[k-2]/4)/
152                      (k+1)-__CZ__Ws_a[k]/2-__CZ__Ws_c[k-1]/(k+1);
153      tmp = __CZ__Ws_c[k] * PP;
154      sum += tmp;
155      if(abs(tmp) <= eps){
156        epsilon(epslocal);
157        return sum;
158      }
159      PP *= P;
160    }
161    epsilon(epslocal);
162    return
163      newerror(strcat("lambertw_series: does not converge in ",
164                        str(limit)," terms" ));
165}
166
167/* */
168define lambertw(z,branch){
169  local eps epslarge ret branchpoint bparea w we ew w1e wn k places m1e;
170  local closeness;
171
172  eps = epsilon();
173  if(branch == 0){
174    if(!im(z)){
175      if(abs(z) <= eps) return 0;
176      if(abs(z-exp(1)) <= eps) return 1;
177      if(abs(z - (-ln(2)/2)) <= eps ) return -ln(2);
178      if(abs(z - (-pi()/2)) <= eps ) return 1i*pi()/2;
179    }
180  }
181
182  branchpoint = -exp(-1);
183  bparea = .2;
184  if(branch == 0){
185    if(!im(z) && abs(z-branchpoint) == 0) return -1;
186    ret = __CZ__lambertw_0(z,eps);
187   /* Yeah, C&P, I know, sorry */
188   ##ret = ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
189  }
190  else if(branch == 1){
191    if(im(z)<0 && abs(z-branchpoint) <= bparea)
192      ret = __CZ__lambertw_app(z,eps);
193   /* Does calc have a goto? Oh, it does!  */
194    ret =ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
195  }
196  else if(branch == -1){##print "-1";
197    if(!im(z) && abs(z-branchpoint) == 0) return -1;
198    if(!im(z) && z>branchpoint && z < 0){##print "0";
199      ret = __CZ__lambertw_m1(z,eps);}
200    if(im(z)>=0 && abs(z-branchpoint) <= bparea){##print "1";
201      ret = __CZ__lambertw_app(z,eps);}
202    ret =ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
203  }
204  else
205    ret = ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
206
207  /*
208     Such a high precision is only needed _very_ close to the branchpoint
209     and might even be insufficient if z has not been computed with
210     sufficient precision itself (M below was calculated by Mathematica and also
211     with the series above with epsilon(1e-200)):
212    ; epsilon(1e-50)
213	0.00000000000000000001
214    ; display(50)
215	20
216    ; M=-0.9999999999999999999999997668356018402875796636464119050387
217    ; lambertw(-exp(-1)+1e-50,0)-M
218	-0.00000000000000000000000002678416515423276355643684
219    ; epsilon(1e-60)
220	0.0000000000000000000000000000000000000000000000000
221    ; A=-exp(-1)+1e-50
222    ; epsilon(1e-50)
223	0.00000000000000000000000000000000000000000000000000
224    ; lambertw(A,0)-M
225	-0.00000000000000000000000000000000000231185460220585
226    ; lambertw_series(A,epsilon(),0)-M
227	-0.00000000000000000000000000000000000132145133161626
228    ; epsilon(1e-100)
229	0.00000000000000000000000000000000000000000000000001
230    ; A=-exp(-1)+1e-50
231    ; epsilon(1e-65)
232	0.00000000000000000000000000000000000000000000000000
233    ; lambertw_series(A,epsilon(),0)-M
234	0.00000000000000000000000000000000000000000000000000
235    ; lambertw_series(-exp(-1)+1e-50,epsilon(),0)-M
236	-0.00000000000000000000000000000000000000002959444084
237    ; epsilon(1e-74)
238	0.00000000000000000000000000000000000000000000000000
239    ; lambertw_series(-exp(-1)+1e-50,epsilon(),0)-M
240	-0.00000000000000000000000000000000000000000000000006
241   */
242  closeness = abs(z-branchpoint);
243  if( closeness< 1){
244    if(closeness != 0)
245      eps = epsilon(epsilon()*( closeness));
246    else
247      eps = epsilon(epsilon()^2);
248  }
249  else
250    eps = epsilon(epsilon()*1e-2);
251
252
253  epslarge =epsilon();
254
255  places = highbit(1 + int(1/epslarge)) + 1;
256  w = ret;
257  for(k=0;k<100;k++){
258    ew = exp(w);
259    we = w*ew;
260    if(abs(we-z)<= 4*epslarge*abs(z))break;
261    w1e = (1+w)*ew;
262    wn = bround(w- ((we - z) / ( w1e - ( (w+2)*(we-z) )/(2*w+2)  ) ),places++) ;
263    if( abs(wn - w) <= epslarge*abs(wn)) break;
264    else  w = wn;
265  }
266
267  if(k==100){
268    epsilon(eps);
269    return newerror("lambertw: Halley iteration does not converge");
270  }
271  /* The Maxima coders added a check if the iteration converged to the correct
272     branch. This coder deems it superfluous. */
273
274  epsilon(eps);
275  return wn;
276}
277
278
279config("resource_debug", resource_debug_level),;
280if (config("resource_debug") & 3) {
281    print "lambertw(z,branch)";
282    print "lambertw_series(z,eps,branch,terms)";
283    print "lambertw_series_print()";
284}
285