tests/rtest_lambert_w.mac

   1 /******************************************************************************
   2   rtest_lambert_w.mac
   3   Test for Lambert W function
   4 ******************************************************************************/
   5
   6 kill(all);
   7 done;
   8
   9 (closeto(value,compare,tol):=
  10   block(
  11     [abse],
  12     abse:abs(value-compare),if(abse<tol) then true else abse),
  13     done);
  14 done;
  15
  16
  17 /* Exact values */
  18 lambert_w(-%pi/2)$
  19 %i*%pi/2;
  20 lambert_w(-log(2)/2)$
  21 -log(2);
  22 lambert_w(-1/%e)$
  23 -1;
  24 lambert_w(0)$
  25 0;
  26 lambert_w(%e)$
  27 1;
  28
  29 /* Numerical values.  Results from Mathematica 6 */
  30
  31 /* This point -47/128 is slightly greater than the branch cut at x=-1/%e.
  32    Known difficult region. */
  33 closeto(lambert_w(-0.3671875),-0.93988639805243454196,5.0e-15);
  34 true;
  35
  36 closeto(lambert_w(-0.25),-0.35740295618138890306,1.0e-16);
  37 true;
  38
  39 closeto(lambert_w(0.0),0.0,1.0e-16);
  40 true;
  41
  42 closeto(lambert_w(0.5),0.3517337112491958260,1.0e-16);
  43 true;
  44
  45 closeto(lambert_w(2.0),0.85260550201372549134,2.0e-16);
  46 true;
  47
  48 closeto(lambert_w(16.0),2.05319271746264858727,1.0e-16);
  49 true;
  50
  51 closeto(lambert_w(128.0),3.57739529855165334288,1.0e-16);
  52 true;
  53
  54 closeto(lambert_w(2048.0),5.85698416299949388816,1.0e-15);
  55 true;
  56
  57 closeto(float(lambert_w(-%pi/2)),%i*%pi/2,1.0e-16);
  58 true;
  59
  60 closeto(lambert_w(3.0+2*%i),
  61   1.12605402258467335157 + 0.31512993765865848403*%i,1.0e-16);
  62 true;
  63
  64 /* bfloat evaluation.  Results from Mathematica 6 to 82 sf */
  65 (oldfpprec:fpprec,fpprec:80,done);
  66 done;
  67
  68 closeto(lambert_w(-0.3671875b0),
  69  -0.9398863980524345419695464132107234809629689590715558376016258594821501408118726902b0,5.0b-79);
  70 true;
  71
  72 closeto(lambert_w(-0.25b0),
  73   -0.3574029561813889030688111040559047533165905550760120436276204485896714025961457963b0,5.0b-80);
  74 true;
  75
  76 closeto(lambert_w(0.0b0),0.0b0,1.0e-80);
  77 true;
  78
  79 closeto(lambert_w(0.5b0),
  80   0.3517337112491958260249093009299510651714642155171118040466438461099606107203387109b0,1.0b-80);
  81 true;
  82
  83 closeto(lambert_w(2.0b0),
  84   0.8526055020137254913464724146953174668984533001514035087721073946525150656742630449b0,1.0b-80);
  85 true;
  86
  87 closeto(lambert_w(16.0b0),
  88   2.053192717462648587277573057065703715549879054683959196524461288893671477132320243b0,5.0b-80);
  89 true;
  90
  91 closeto(lambert_w(128.0b0),
  92   3.577395298551653342882928583858277466679549941663670990487212206804013161129017297b0,1.0b-80);
  93 true;
  94
  95 closeto(lambert_w(2048.0b0),
  96   5.856984162999493888162560618072936017287078565867036752738725003619707275074129900b0,5.0b-80);
  97 true;
  98
  99 closeto(lambert_w(3.0b0+2.0b0*%i),
 100 1.1260540225846733515720001139613054043549635214873248102633125218051215201783683422b0
 101 +
 102 0.3151299376586584840326916101687325819105911412081543422333577478256866198890713045b0
 103 *%i ,1.0b-80);
 104 true;
 105
 106 /* Check an argument > most-positive-double-float ~ 1.0e308 */
 107 closeto(lambert_w(1.0b1000),
 108 2294.846671683506869652792785993616789973426699478802684415757740164730983851156757b0,
 109 1.0b-77);
 110 true;
 111
 112 /* Arguments that are exact numbers are not evaluated numerically */
 113 lambert_w(2);
 114 lambert_w(2);
 115
 116 lambert_w(3+2*%i);
 117 lambert_w(3+2*%i);
 118
 119 /* ... unless we set numer evflag*/
 120 closeto(lambert_w(2),0.85260550201372549134,2.0e-16), numer;
 121 true;
 122
 123 diff(lambert_w(x),x);
 124 %e^-lambert_w(x)/(lambert_w(x)+1);
 125
 126 integrate(lambert_w(x),x);
 127 x*(lambert_w(x)^2-lambert_w(x)+1)/lambert_w(x);
 128
 129 /* SF bug report 2468610 the integrator loops endlessly
 130    errcatch() doesn't catch the error for clisp-2.46 on cygwin */
 131 integrate(lambert_w(1/x),x);
 132 'integrate(lambert_w(1/x),x);
 133
 134 taylor(lambert_w(x),x,0,6);
 135 x-x^2+3*x^3/2-(8*x^4/3)+125*x^5/24-(54*x^6/5);
 136
 137 /******************************************************************************
 138   Tests for Generalized Lambert W function
 139 ******************************************************************************/
 140
 141 closeto(generalized_lambert_w(0,-0.25),-0.35740295618138890306,1.0e-16);
 142 true;
 143
 144 closeto(generalized_lambert_w(0,3.0+2*%i),
 145   1.12605402258467335157 + 0.31512993765865848403*%i,1.0e-16);
 146 true;
 147
 148 closeto(generalized_lambert_w(0,-0.25),-0.35740295618138890306,1.0e-16);
 149 true;
 150
 151 /* This point -753/2048 is slightly greater than the branch cut at x=-1/%e.
 152    Known difficult region. */
 153 closeto(generalized_lambert_w(-1,-753.0/2048),-1.033649565301979,5.0e-15);
 154 true;
 155
 156 closeto(generalized_lambert_w(-1,-754.0/2048),
 157    -0.9994844032397146 - 0.0393272270347577*%i,5.32e-15);
 158 true;
 159
 160 closeto(generalized_lambert_w(1,-753.0/2048),
 161   -3.089416043057986+7.461420464927938*%i,5.0e-15);
 162 true;
 163
 164 closeto(generalized_lambert_w(1,-754.0/2048),
 165   -3.088042730569348+7.461585404728464*%i,5.0e-15);
 166 true;
 167
 168 /* Reduced precision when very near branch point z = -1/%e */
 169 closeto(generalized_lambert_w(-1,z:float(-1/%e)-1.0e-12),
 170  -0.999999999998187812 - 2.3316439815951875e-6*%i,1e-9);
 171 true;
 172
 173 /* The branch test was wrong for this point due to roundoff */
 174 closeto(generalized_lambert_w(-1,-1/float(%e)),-1.0,1e-15);
 175 true;
 176
 177 /* bigfloat evaluation near the branch point failed to obtain a starting guess
 178    There is significant (~ 14 decimal digits) loss of precision due to ill-conditioning
 179    The expected answer below checked with 200 and 300 digit bfloats */
 180 closeto(generalized_lambert_w(-1,bfloat(-1/%e)+1b-30),
 181 - 1.000000000000002331643981597126015551421701533827747686376229738168073845444881159b0,
 182  1.0b-64);
 183 true;
 184
 185 closeto(generalized_lambert_w(1,-754.0b0/2048),
 186  -3.088042730569348561456014556706168490039090862315553536854637271843252771228714625b0
 187  +7.461585404728464111653134458076750679263104524925586270112477441975032123188737473b0*%i,
 188  1.0b-77);
 189 true;
 190
 191 /* Check an argument > most-positive-double-float ~ 1.0e308 */
 192 closeto(generalized_lambert_w(0,1.0b1000),
 193 2294.846671683506869652792785993616789973426699478802684415757740164730983851156757b0,
 194 1.0b-77);
 195 true;
 196
 197 closeto(generalized_lambert_w(0,1.0b1000*%i),
 198   2.294846671449550599071991520727747213463470638397366444770394873636319358115746758b3
 199  +1.570112136469240112764525899485121673623387222840770033671975542758566442609326195b0*%i
 200  ,1.0b-77);
 201 true;
 202
 203
 204 /* Arguments that are not evaluated numerically */
 205 generalized_lambert_w(3,1/2);
 206 generalized_lambert_w(3,1/2);
 207
 208 generalized_lambert_w(1/2,3.0);
 209 generalized_lambert_w(1/2,3.0);
 210
 211 generalized_lambert_w(0.5,3.0);
 212 generalized_lambert_w(0.5,3.0);
 213
 214 generalized_lambert_w(1.0,3.0);
 215 generalized_lambert_w(1.0,3.0);
 216
 217 generalized_lambert_w(%pi,3.0);
 218 generalized_lambert_w(%pi,3.0);
 219
 220 /* ... unless we set numer evflag */
 221 closeto(generalized_lambert_w(0,2),0.85260550201372549134,2.0e-16), numer;
 222 true;
 223
 224
 225 diff(generalized_lambert_w(n,x),x);
 226 %e^-generalized_lambert_w(n,x)/(generalized_lambert_w(n,x)+1);
 227
 228 integrate(generalized_lambert_w(n,x),x);
 229 (generalized_lambert_w(n,x)^2-generalized_lambert_w(n,x)+1)*x
 230  /generalized_lambert_w(n,x);
 231
 232 (fpprec:oldfpprec, kill(oldfpprec), done);
 233 done;
 234
 235 /*
 236  * Bug 4156: float of generalized_lambert_w doesn't work
 237  */
 238 float(generalized_lambert_w(-1, 1/2));
 239 generalized_lambert_w(-1, 0.5);
 240
 241 bfloat(generalized_lambert_w(-1, 1/2));
 242 generalized_lambert_w(-1, 0.5b0);
 243
 244 float(generalized_lambert_w(-3/2, 1/2));
 245 generalized_lambert_w(-1.5, 0.5);
 246
 247 bfloat(generalized_lambert_w(-3/2, 1/2));
 248 generalized_lambert_w(-1.5b0, 0.5b0);
 249