Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - lambert.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28695-49bb1ac00f) Lines: 237 264 89.8 %
Date: 2023-09-24 07:47:42 Functions: 18 22 81.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2021  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : /***********************************************************************/
      19             : /**                 LAMBERT's W_K FUNCTIONS                           **/
      20             : /***********************************************************************/
      21             : /* roughly follows Veberic, https://arxiv.org/abs/1003.1628 */
      22             : 
      23             : static double
      24          21 : dblL1L2(double L1)
      25             : {
      26          21 :   double L2 = log(-L1), LI = 1 / L1, N2, N3, N4, N5;
      27          21 :   N2 = (L2-2.)/2.; N3 = (6.+L2*(-9.+2.*L2))/6.;
      28          21 :   N4 = (-12.+L2*(36.+L2*(-22.+3*L2)))/12.;
      29          21 :   N5 = (60.+L2*(-300.+L2*(350.+L2*(-125.+12*L2))))/60.;
      30          21 :   return L1-L2+L2*LI*(1+LI*(N2+LI*(N3+LI*(N4+LI*N5))));
      31             : }
      32             : 
      33             : /* rough approximation to W0(a > -1/e), < 1% relative error */
      34             : double
      35      474466 : dbllambertW0(double a)
      36             : {
      37      474466 :   if (a < -0.2583)
      38             :   {
      39           0 :     const double c2 = -1./3, c3 = 11./72, c4 = -43./540, c5 = 769./17280;
      40           0 :     double p = sqrt(2 * (M_E * a + 1));
      41           0 :     if (a < -0.3243) return -1+p*(1+p*(c2+p*c3));
      42           0 :     return -1+p*(1+p*(c2+p*(c3+p*(c4+p*c5))));
      43             :   }
      44             :   else
      45             :   {
      46      474466 :     double Wd = log(1.+a);
      47      474466 :     Wd *= (1.-log(Wd/a))/(1.+Wd);
      48      474466 :     if (a < 0.6482 && a > -0.1838) return Wd;
      49      419401 :     return Wd*(1.-log(Wd/a))/(1.+Wd);
      50             :   }
      51             : }
      52             : /* uniform approximation to W0, at least 15 bits. */
      53             : static double
      54         532 : dbllambertW0init(double a)
      55             : {
      56         532 :   if (a < -0.323581)
      57             :   {
      58          35 :     const double c2 = 1./3., c3 = 11./72., c4 = 43./540., c5 = 769./17280.;
      59          35 :     const double c6 = 221./8505., c7 = 680863./43545600.;
      60          35 :     const double c8 = 1963./204120., c9 = 226287557./37623398400.;
      61          35 :     double p = M_E * a + 1;
      62          35 :     if (p <= 0) return -1;
      63          23 :     p = -sqrt(2 * p);
      64          23 :     return -(1.+p*(1.+p*(c2+p*(c3+p*(c4+p*(c5+p*(c6+p*(c7+p*(c8+p*c9)))))))));
      65             :   }
      66         497 :   if (a < 0.145469)
      67             :   {
      68         112 :     const double a1 = 5.931375, a2 = 11.392205, a3 = 7.338883, a4 = 0.653449;
      69         112 :     const double b1 = 6.931373, b2 = 16.823494, b3 = 16.430723, b4 = 5.115235;
      70         112 :     double n = 1.+a*(a1+a*(a2+a*(a3+a*a4)));
      71         112 :     double d = 1.+a*(b1+a*(b2+a*(b3+a*b4)));
      72         112 :     return a * n / d;
      73             :   }
      74         385 :   if (a < 8.706658)
      75             :   {
      76         378 :     const double a1 = 2.445053, a2 = 1.343664, a3 = 0.148440, a4 = 0.000804;
      77         378 :     const double b1 = 3.444708, b2 = 3.292489, b3 = 0.916460, b4 = 0.053068;
      78         378 :     double n = 1.+a*(a1+a*(a2+a*(a3+a*a4)));
      79         378 :     double d = 1.+a*(b1+a*(b2+a*(b3+a*b4)));
      80         378 :     return a * n / d;
      81             :   }
      82             :   else
      83             :   {
      84           7 :     double w = log(1.+a);
      85           7 :     w *= (1.-log(w/a)) / (1.+w);
      86           7 :     return w * (1.-log(w/a)) / (1.+w);
      87             :   }
      88             : }
      89             : 
      90             : /* rough approximation to W_{-1}(0 > a > -1/e), < 1% relative error */
      91             : double
      92       66962 : dbllambertW_1(double a)
      93             : {
      94       66962 :   if (a < -0.2464)
      95             :   {
      96         280 :     const double c2 = -1./3, c3 = 11./72, c4 = -43./540, c5 = 769./17280;
      97         280 :     double p = -sqrt(2 * (M_E * a + 1));
      98         280 :     if (a < -0.3243) return -1+p*(1+p*(c2+p*c3));
      99         175 :     return -1+p*(1+p*(c2+p*(c3+p*(c4+p*c5))));
     100             :   }
     101             :   else
     102             :   {
     103             :     double Wd;
     104       66682 :     a = -a; Wd = -log(a);
     105       66682 :     Wd *= (1.-log(Wd/a))/(1.-Wd);
     106       66682 :     if (a < 0.0056) return -Wd;
     107         658 :     return -Wd*(1.-log(Wd/a))/(1.-Wd);
     108             :   }
     109             : }
     110             : /* uniform approximation to W_{-1}, at least 15 bits. */
     111             : static double
     112         112 : dbllambertW_1init(double a)
     113             : {
     114         112 :   if (a < -0.302985)
     115             :   {
     116          21 :     const double c2 = 1./3., c3 = 11./72., c4 = 43./540., c5 = 769./17280.;
     117          21 :     const double c6 = 221./8505., c7 = 680863./43545600.;
     118          21 :     const double c8 = 1963./204120., c9 = 226287557./37623398400.;
     119          21 :     double p = M_E * a + 1;
     120          21 :     if (p <= 0) return -1;
     121          21 :     p = sqrt(2 * p);
     122          21 :     return -(1.+p*(1.+p*(c2+p*(c3+p*(c4+p*(c5+p*(c6+p*(c7+p*(c8+p*c9)))))))));
     123             :   }
     124          91 :   if (a <= -0.051012)
     125             :   {
     126          77 :     const double a0 = -7.814176, a1 = 253.888101, a2 = 657.949317;
     127          77 :     const double b1 = -60.439587, b2 = 99.985670, b3 = 682.607399;
     128          77 :     const double b4 = 962.178439, b5 = 1477.934128;
     129          77 :     double n = a0+a*(a1+a*a2);
     130          77 :     double d = 1+a*(b1+a*(b2+a*(b3+a*(b4+a*b5))));
     131          77 :     return n / d;
     132             :   }
     133          14 :   return dblL1L2(log(-a));
     134             : }
     135             : 
     136             : /* uniform approximation to more than 46 bits, 50 bits away from -1/e;
     137             :  * branch = -1 or 0 */
     138             : static double
     139         651 : dbllambertWfritsch(GEN ga, int branch)
     140             : {
     141             :   double a, z, w1, q, w;
     142         651 :   if (expo(ga) >= 0x3fe)
     143             :   { /* branch = 0 */
     144           7 :     double w = dbllog2(ga) * M_LN2; /* ~ log(1+a) ~ log a */
     145           7 :     return w * (1.+w-log(w)) / (1.+w);
     146             :   }
     147         644 :   a = rtodbl(ga);
     148         644 :   w = branch? dbllambertW_1init(a): dbllambertW0init(a);
     149         644 :   if (w == -1.|| w == 0.) return w;
     150         632 :   z = log(a / w) - w; w1 = 1. + w;
     151         632 :   q = 2. * w1 * (w1 + (2./3.) * z);
     152         632 :   return w * (1 + (z / w1) * (q - z) / (q - 2 * z));
     153             : }
     154             : 
     155             : static double
     156           7 : dbllambertWhalleyspec(double loga)
     157             : {
     158           7 :   double w = dblL1L2(loga);
     159             :   for(;;)
     160           0 :   {
     161           7 :     double n = w + log(-w) - loga, d = 1 - w, r = n / (d + n / d);
     162           7 :     w *= 1 - r; if (r < 2.e-15) return w;
     163             :   }
     164             : }
     165             : /* k = 0 or -1. */
     166             : static GEN
     167         658 : lambertW(GEN z, long k, long prec)
     168             : {
     169         658 :   pari_sp av = avma;
     170         658 :   long bit = prec2nbits(prec), L = -(bit / 3 + 10), ct = 0, p, pb;
     171             :   double wd;
     172             :   GEN w, vp;
     173             : 
     174         658 :   if (gequal0(z) && !k) return real_0(prec);
     175         658 :   z = gtofp(z, prec);
     176         658 :   if (k == -1)
     177             :   {
     178         119 :     long e = expo(z);
     179         119 :     if (signe(z) >= 0) pari_err_DOMAIN("lambertw", "z", ">", gen_0, z);
     180         119 :     wd = e < -512? dbllambertWhalleyspec(dbllog2(z) * M_LN2)
     181         119 :                  : dbllambertWfritsch(z, -1);
     182             :   }
     183             :   else
     184         539 :     wd = dbllambertWfritsch(z, 0);
     185         658 :   if (fabs(wd + 1) < 1e-5)
     186             :   {
     187          14 :     long prec2 = prec + EXTRAPREC64;
     188          14 :     GEN Z = rtor(z, prec2);
     189          14 :     GEN t = addrs(mulrr(Z, gexp(gen_1, prec2)), 1);
     190          14 :     if (signe(t) <= 0) { set_avma(av); return real_m1(prec); }
     191           0 :     if (realprec(t) < prec)
     192             :     {
     193           0 :       prec2 += prec - realprec(t);
     194           0 :       Z = rtor(z, prec2);
     195           0 :       t = addrs(mulrr(Z, gexp(gen_1, prec2)), 1);
     196             :     }
     197           0 :     t = sqrtr(shiftr(t, 1));
     198           0 :     w = gprec_w(k == -1? subsr(-1, t) : subrs(t, 1), prec);
     199           0 :     p = prec - 2; vp = NULL;
     200             :   }
     201             :   else
     202             :   { /* away from -1/e: can reduce accuracy and self-correct */
     203         644 :     w = wd == 0.? z: dbltor(wd);
     204         644 :     vp = cgetg(30, t_VECSMALL); ct = 0; pb = bit;
     205        1425 :     while (pb > BITS_IN_LONG * 3/4)
     206         781 :     { vp[++ct] = (pb + BITS_IN_LONG-1) >> TWOPOTBITS_IN_LONG; pb = (pb + 2) / 3; }
     207         644 :     p = vp[ct]; w = gprec_w(w, p + 2);
     208             :   }
     209         644 :   if ((k == -1 && (bit < 192 || bit > 640)) || (k == 0 && bit > 1024))
     210             :   {
     211             :     for(;;)
     212          13 :     {
     213             :       GEN t, ew, n, d;
     214         104 :       ew = mplog(divrr(w, z)); n = addrr(w, ew); d = addrs(w, 1);
     215         104 :       t = divrr(n, shiftr(d, 1));
     216         104 :       w = mulrr(w, subsr(1, divrr(n, addrr(d, t))));
     217         104 :       if (p >= prec-2 && expo(n) - expo(d) - expo(w) <= L) break;
     218          13 :       if (vp) { if (--ct) p = vp[ct]; w = gprec_w(w, ct? p + 2: prec); }
     219             :     }
     220             :   }
     221             :   else
     222             :   {
     223             :     for(;;)
     224         185 :     {
     225             :       GEN t, ew, wew, n, d;
     226         738 :       ew = mpexp(w); wew = mulrr(w, ew); n = subrr(wew, z); d = addrr(ew, wew);
     227         738 :       t = divrr(mulrr(addrs(w, 2), n), shiftr(addrs(w, 1), 1));
     228         738 :       w = subrr(w, divrr(n, subrr(d, t)));
     229         738 :       if (p >= prec-2 && expo(n) - expo(d) - expo(w) <= L) break;
     230         185 :       if (vp) { if (--ct) p = vp[ct]; w = gprec_w(w, ct? p + 2: prec); }
     231             :     }
     232             :   }
     233         644 :   return gerepileupto(av, w);
     234             : }
     235             : 
     236             : /*********************************************************************/
     237             : /*                       Complex branches                            */
     238             : /*********************************************************************/
     239             : 
     240             : /* x *= (1 - (x + log(x) - L) / (x + 1)); L = log(z) + 2IPi * k */
     241             : static GEN
     242      632752 : lamaux(GEN x, GEN L, long *pe, long prec)
     243             : {
     244      632752 :   GEN n = gsub(gadd(x, glog(x, prec)), L);
     245      632752 :   if (pe) *pe = maxss(4, -gexpo(n));
     246      632752 :   if (gequal0(imag_i(n))) n = real_i(n);
     247      632752 :   return gmul(x, gsubsg(1, gdiv(n, gaddsg(1, x))));
     248             : }
     249             : 
     250             : /* Complex branches, experimental */
     251             : static GEN
     252       78428 : lambertWC(GEN z, long branch, long prec)
     253             : {
     254       78428 :   pari_sp av = avma;
     255             :   GEN w, pii2k, zl, lzl, L, Lz;
     256       78428 :   long bit0, si, j, fl = 0, lim = 6, lp = DEFAULTPREC, bit = prec2nbits(prec);
     257             : 
     258       78428 :   si = gsigne(imag_i(z)); if (!si) z = real_i(z);
     259       78428 :   pii2k = gmulsg(branch, PiI2(lp));
     260       78428 :   zl = gtofp(z, lp); lzl = glog(zl, lp);
     261             :   /* From here */
     262       78428 :   if (branch == 0 || branch * si < 0
     263       26579 :       || (si == 0 && gsigne(z) < 0 && branch == -1))
     264             :   {
     265       51989 :     GEN lnzl1 = gaddsg(1, glog(gneg(zl), lp));
     266       51989 :     if (si == 0) si = gsigne(lnzl1);
     267       51989 :     if ((branch == 0 || branch * si < 0) && gexpo(lnzl1) < -1)
     268             :     { /* close to -1/e */
     269        2408 :       w = gaddsg(1, gmul(z, gexp(gen_1, prec)));
     270        2408 :       w = gprec_wtrunc(w, lp);
     271        2408 :       w = gsqrt(gmul2n(w, 1), lp);
     272        2408 :       w = branch * si < 0? gsubsg(-1, w): gaddsg(-1, w);
     273        2408 :       lim = 10; fl = 1;
     274             :     }
     275       51989 :     if (branch == 0 && !fl && gexpo(lzl) < 0) { w = zl; fl = 1; }
     276             :   }
     277       78428 :   if (!fl)
     278             :   {
     279       68152 :     if (branch)
     280             :     {
     281       51212 :       GEN lr = glog(pii2k, lp);
     282       51212 :       w = gadd(gsub(gadd(pii2k, lzl), lr), gdiv(gsub(lr, lzl), pii2k));
     283             :     }
     284             :     else
     285             :     {
     286       16940 :       GEN p = gaddsg(1, gmul(z, gexp(gen_1, lp)));
     287       16940 :       w = gexpo(p) > 0? lzl: gaddgs(gsqrt(p, lp), -1);
     288             :     }
     289             :   }
     290             :   /* to here: heuristic */
     291       78428 :   L = gadd(lzl, pii2k);
     292      480200 :   for (j = 1; j < lim; j++) w = lamaux(w, L, NULL, lp);
     293       78428 :   Lz = NULL;
     294       78428 :   if (branch == 0 || branch == -1)
     295             :   {
     296       52199 :     Lz = glog(z, prec);
     297       52199 :     if (branch == -1)
     298             :     {
     299       25970 :       long flag = 1;
     300       25970 :       if (!si && signe(z) <= 0 && signe(addrs(Lz, 1))) flag = 0;
     301       25970 :       if (flag) Lz = gsub(Lz, PiI2(prec));
     302             :     }
     303             :   }
     304       78428 :   w = lamaux(w, L, &bit0, lp);
     305      230980 :   while (bit0 < bit || (Lz && gexpo(gsub(gadd(w, glog(w, prec)), Lz)) > 16-bit))
     306             :   {
     307      152552 :     long p = nbits2prec(bit0 <<= 1);
     308      152552 :     L = gadd(gmulsg(branch, PiI2(p)), glog(gprec_w(z, p), p));
     309      152552 :     w = lamaux(gprec_w(w, p), L, NULL, p);
     310             :   }
     311       78428 :   return gerepilecopy(av, gprec_w(w, nbits2prec(bit)));
     312             : }
     313             : 
     314             : /* exp(t (1 + O(t^n))), n >= 0 */
     315             : static GEN
     316         154 : serexp0(long v, long n)
     317             : {
     318         154 :   GEN y = cgetg(n+3, t_SER), t;
     319             :   long i;
     320         154 :   y[1] = evalsigne(1) | evalvarn(v) | evalvalser(0);
     321         154 :   gel(y,2) = gen_1; if (!n) return y;
     322         147 :   gel(y,3) = gen_1; if (n == 1) return y;
     323        1295 :   for (i=2, t = gen_2; i < n; i++, t = muliu(t,i)) gel(y,i+2) = mkfrac(gen_1,t);
     324         119 :   gel(y,i+2) = mkfrac(gen_1,t); return y;
     325             : }
     326             : 
     327             : /* series expansion of W at -1/e */
     328             : static GEN
     329           7 : Wbra(long N)
     330             : {
     331           7 :   GEN v = cgetg(N + 2, t_VEC);
     332             :   long n;
     333           7 :   gel(v, 1) = gen_m1;
     334           7 :   gel(v, 2) = gen_1;
     335          56 :   for (n = 2; n <= N; n++)
     336             :   {
     337          49 :     GEN t = gel(v,n), a = gen_0;
     338          49 :     long k, K = (n - 1) >> 1;
     339         133 :     for (k = 1; k <= K; k++) t = gadd(t, gmul2n(gel(v,n-2*k), -k));
     340         196 :     for (k = 2; k < n; k++) a = gadd(a, gmul(gel(v,k+1), gel(v,n+2-k)));
     341          49 :     gel(v,n+1) = gsub(gdivgs(t, -n-1), gmul2n(a, -1));
     342             :   }
     343           7 :   return RgV_to_RgX(v, 0);
     344             : }
     345             : 
     346             : static GEN
     347         154 : reverse(GEN y)
     348             : {
     349         154 :   GEN z = ser2rfrac_i(y);
     350         154 :   long l = lg(z);
     351         154 :   return RgX_to_ser(RgXn_reverse(z, l-2), l-1);
     352             : }
     353             : static GEN
     354         182 : serlambertW(GEN y, long branch, long prec)
     355             : {
     356             :   long n, vy, val, v;
     357         182 :   GEN t = NULL;
     358             : 
     359         182 :   if (!signe(y)) return gcopy(y);
     360         182 :   v = valser(y);
     361         182 :   if (v < 0) pari_err_DOMAIN("lambertw","valuation", "<", gen_0, y);
     362         175 :   if (v > 0 && branch)
     363           0 :     pari_err_DOMAIN("lambertw [k != 0]", "x", "~", gen_0, y);
     364         175 :   vy = varn(y); n = lg(y)-3;
     365         483 :   for (val = 1; val < n; val++)
     366         434 :     if (!gequal0(polcoef_i(y, val, vy))) break;
     367         175 :   if (v)
     368             :   {
     369          70 :     t = serexp0(vy, n / val);
     370          70 :     setvalser(t, 1); t = reverse(t); /* rev(x exp(x)) */
     371             :   }
     372             :   else
     373             :   {
     374         105 :     GEN y0 = gel(y,2), x = glambertW(y0, branch, prec);
     375         105 :     if (val > n) return scalarser(x, vy, n+1);
     376          98 :     y = serchop0(y);
     377          98 :     if (gequalm1(x))
     378             :     { /* y0 ~ -1/e, branch = 0 or -1 */
     379          14 :       GEN p = gmul(shiftr(gexp(gen_1,prec), 1), y);
     380          14 :       if (odd(val)) pari_err(e_MISC, "odd valuation at branch point");
     381           7 :       p = gsqrt(p, prec); if (odd(branch)) p = gneg(p);
     382           7 :       n -= val >> 1;
     383           7 :       t = RgXn_eval(Wbra(n), ser2rfrac_i(p), n);
     384           7 :       return gtoser(t, varn(t), lg(p));
     385             :     }
     386          84 :     t = serexp0(vy, n / val);
     387             :     /* (x + t) exp(x + t) = (y0 + t y0/x) * exp(t) */
     388          84 :     t = gmul(deg1pol_shallow(gdiv(y0,x), y0, vy), t);
     389          84 :     t = gadd(x, reverse(serchop0(t)));
     390             :   }
     391         154 :   return normalizeser(gsubst(t, vy, y));
     392             : }
     393             : 
     394             : static GEN
     395          14 : lambertp(GEN x)
     396             : {
     397          14 :   pari_sp av = avma;
     398             :   long k;
     399             :   GEN y;
     400             : 
     401          14 :   if (gequal0(x)) return gcopy(x);
     402          14 :   if (!valp(x)) { x = leafcopy(x); setvalp(x, 1); }
     403          14 :   k = Qp_exp_prec(x);
     404          14 :   if (k < 0) return NULL;
     405          14 :   y = gpowgs(cvstop2(k, x), k - 1);
     406         266 :   for (k--; k; k--)
     407         252 :     y = gsub(gpowgs(cvstop2(k, x), k - 1), gdivgu(gmul(x, y), k + 1));
     408          14 :   return gerepileupto(av, gmul(x, y));
     409             : }
     410             : 
     411             : /* y a t_REAL */
     412             : static int
     413        2219 : useC(GEN y, long k)
     414             : {
     415        2219 :   if (signe(y) > 0 || (k && k != -1)) return k ? 1: 0;
     416         980 :   return gsigne(addsr(1, logr_abs(y))) > 0;
     417             : }
     418             : static GEN
     419       80808 : glambertW_i(void *E, GEN y, long prec)
     420             : {
     421             :   pari_sp av;
     422       80808 :   long k = (long)E, p;
     423             :   GEN z;
     424       80808 :   if (gequal0(y))
     425             :   {
     426          21 :     if (k) pari_err_DOMAIN("glambertW","argument","",gen_0,y);
     427          14 :     return gcopy(y);
     428             :   }
     429       80787 :   switch(typ(y))
     430             :   {
     431        2219 :     case t_REAL:
     432        2219 :       p = minss(prec, realprec(y));
     433        2219 :       return useC(y, k)? lambertWC(y, k, p): lambertW(y, k, p);
     434          14 :     case t_PADIC: z = lambertp(y);
     435          14 :       if (!z) pari_err_DOMAIN("glambertW(t_PADIC)","argument","",gen_0,y);
     436          14 :       return z;
     437       76867 :     case t_COMPLEX:
     438       76867 :       p = precision(y);
     439       76867 :       return lambertWC(y, k, p? p: prec);
     440        1687 :     default:
     441        1687 :       av = avma; if (!(z = toser_i(y))) break;
     442         182 :       return gerepileupto(av, serlambertW(z, k, prec));
     443             :   }
     444        1505 :   return trans_evalgen("lambert", E, glambertW_i, y, prec);
     445             : }
     446             : 
     447             : GEN
     448       79303 : glambertW(GEN y, long k, long prec)
     449             : {
     450       79303 :   return glambertW_i((void*)k, y, prec);
     451             : }
     452             : GEN
     453           0 : mplambertW(GEN y, long prec) { return lambertW(y, 0, prec); }
     454             : 
     455             : /*********************************************************************/
     456             : /*                        Application                                */
     457             : /*********************************************************************/
     458             : /* Solve x - a * log(x) = b with a > 0 and b >= a * (1 - log(a)). */
     459             : GEN
     460           0 : mplambertx_logx(GEN a, GEN b, long bit)
     461             : {
     462           0 :   pari_sp av = avma;
     463           0 :   GEN e = gexp(gneg(gdiv(b, a)), nbits2prec(bit));
     464           0 :   return gerepileupto(av, gmul(gneg(a), lambertW(gneg(gdiv(e, a)), -1, bit)));
     465             : }
     466             : /* Special case a = 1, b = log(y): solve e^x / x = y with y >= exp(1). */
     467             : GEN
     468           0 : mplambertX(GEN y, long bit)
     469             : {
     470           0 :   pari_sp av = avma;
     471           0 :   return gerepileupto(av, gneg(lambertW(gneg(ginv(y)), -1, bit)));
     472             : }
     473             : 
     474             : /* Solve x * log(x) - a * x = b; if b < 0, assume a >= 1 + log |b|. */
     475             : GEN
     476           0 : mplambertxlogx_x(GEN a, GEN b, long bit)
     477             : {
     478           0 :   pari_sp av = avma;
     479           0 :   long s = gsigne(b);
     480             :   GEN e;
     481           0 :   if (!s) return gen_0;
     482           0 :   e = gexp(gneg(a), nbits2prec(bit));
     483           0 :   return gerepileupto(av, gdiv(b, lambertW(gmul(b, e), s > 0? 0: -1, bit)));
     484             : }

Generated by: LCOV version 1.14