c_cpp/etc/calc/cal/lambertw.cal

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