Calcul précis de la branche principale de la fonction W de Lambert avec la bibliothèque mathématique C standard

Question

En Lambert W est la fonction inverse de f(w) = we w . Il s'agit d'une fonction multivaluée qui possède une infinité de branches sur les nombres complexes, mais qui n'a que deux branches sur les nombres réels, désignées par W ₀ et W _-1 . W ₀ est considérée comme la branche principale, avec une entrée domaine de [ -1/e, ] et W _-1 a un domaine d'entrée de ( -1/e, 0 ). Les implémentations correspondantes sont souvent appelées lambert_w0() y lambert_wm1() .

Un proche parent de la fonction a été identifié pour la première fois par Leonard Euler [1] lors du suivi des travaux de Johann Heinrich Lambert [2] . Euler a examiné la solution d'équations transcendantes x -x \= ( - ) v x + et dans le processus, nous avons considéré un cas simplifié ln x = v x . Il a ainsi introduit une fonction d'aide avec l'expansion en série suivante autour de zéro :

y = 1 + (2 1 /(1-2))u + (3 2 /(1-2-3))u 2 + (4 3 /(1-2-3-4))u 3 + (5 4 /((1-2 -5))u 4 +

En termes modernes, cette fonction (non nommée par Gauss) représente W(-x)/x et la solution de ln x = v x es x=(-W(-v)/-v) (1/) .

Bien que la fonction W de Lambert ait fait une apparition occasionnelle dans la littérature, par ex. [3] elle n'a pas été nommée et reconnue comme un élément constitutif important avant les travaux fondamentaux de Robert Corless dans les années 1990, par exemple [4] . Par la suite, l'applicabilité de la fonction W de Lambert à la fois aux mathématiques et aux sciences physiques a été élargie par des recherches continues, dont certains exemples sont donnés dans [5] .

La fonction W de Lambert ne fait pas actuellement partie de la bibliothèque mathématique standard de l'ISO-C, et il ne semble pas y avoir de projet immédiat pour l'ajouter. *Comment la branche principale de la fonction W de Lambert, W ₀ être mis en œuvre avec précision en utilisant la bibliothèque mathématique standard ISO-C ?

Une implémentation fidèlement arrondie est probablement trop ambitieuse, mais maintenir une limite d'erreur de 4 ulp (telle que choisie par le profil LA de la bibliothèque mathématique Intel) semble réalisable et souhaitable. La prise en charge de l'arithmétique binaire à virgule flottante IEEE-754 (2008) et la prise en charge des opérations de multiplication-addition fusionnées (FMA) accessibles par l'intermédiaire de la bibliothèque mathématique Intel. fma() y fmaf() les fonctions de la bibliothèque standard peuvent être supposées.

---

[1] Leonard Euler, "De serie Lambertina plurimisque eius insignibus proprietatibus," (Sur la série de Lambert et ses nombreuses propriétés distinctives) Acta Academiae Scientiarum Imperialis Petropolitanae pro Anno MDCCLXXIX Tomus III, Pars II, (Actes de l'Académie impériale des sciences de Saint-Pétersbourg pour l'année 1779, volume 3, partie 2, juil. - déc.), Saint-Pétersbourg : Académie des sciences 1783, pp. 29-51 ( numérisation en ligne à la Bibliothèque d'État de Bavière, Munich)

[2] Johann Heinrich Lambert, "Observationes variae in mathesin puram" (Observations diverses sur les mathématiques pures). Acta Helveticae physico-mathematico-anatomico-botanico-medica Vol. 3, Bâle : J. R. Imhof 1758, pp. 128-168 ( numérisation en ligne à la bibliothèque du patrimoine de la biodiversité)

[3] F.N. Fritsch, R.E. Shafer et W.P. Crowley, "Algorithm 443 : Solution of the Transcendental Equation". nous x \=x ", Communications de l'ACM , vol. 16, n° 2, février 1973, p. 123-124.

[4] R.M. Corless, et al, "On the Lambert W function," Progrès dans le domaine des mathématiques computationnelles , vol. 5, n° 1, décembre 1996, pp. 329-359

[5] Iordanis Kesisoglou, Garima Singh et Michael Nikolaou, "The Lambert Function Should Be in the Engineering Mathematical Toolbox", Génie informatique et chimique , 148, mai 2021

Answer 1

D'après la littérature, il est clair que la méthode la plus courante de calcul de la fonction Lambert W sur les nombres réels est l'itération fonctionnelle. La méthode de Newton du second ordre est le plus simple de ces schémas :

w _i+1 \= w _i - (w _i exp(w _i ) - x) / (exp (w _i ) + w _i exp(w _i )).

Une grande partie de la littérature préfère les méthodes d'ordre supérieur, telles que celles de Halley, Fritsch et Schroeder. Au cours de travaux exploratoires, j'ai constaté que, lorsqu'elles sont effectuées en virgule flottante à précision finie, les propriétés numériques de ces itérations d'ordre supérieur ne sont pas aussi favorables que l'ordre peut le suggérer. À titre d'exemple, trois itérations de Newton étaient systématiquement plus performantes que deux itérations de Halley en termes de précision. Pour cette raison, j'ai choisi l'itération de Newton comme bloc de construction principal, et j'ai utilisé des implémentations personnalisées de exp() qui n'ont qu'à fournir des résultats qui sont représentables comme des normales positives dans la terminologie IEEE-754 pour regagner un peu de performance.

J'ai également constaté que l'itération de Newton ne converge que lentement pour les grands opérandes, en particulier lorsque l'approximation de départ n'est pas très précise. Comme un nombre élevé d'itérations n'est pas propice à de bonnes performances, j'ai cherché une alternative et j'ai trouvé un candidat supérieur dans le schéma d'itération basé sur le logarithme de Iacono et Boyd. [*] qui a également une convergence du second ordre :

w _i+1 \= (w _i / (1 + w _i )) * (1 + log (x / w _i )

De nombreuses implémentations de la fonction W de Lambert semblent utiliser différentes approximations de départ pour différentes parties du domaine d'entrée. Ma préférence allait à une seule approximation de départ pour l'ensemble du domaine d'entrée, dans l'optique de futures implémentations vectorisées.

Heureusement, Iacono et Boyd fournissent également une approximation de départ universelle qui fonctionne dans tout le domaine d'entrée de W_. 0 qui, bien que ne tenant pas entièrement ses promesses, fonctionne très bien. Je l'ai affiné pour l'implémentation en simple précision qui traite un domaine d'entrée beaucoup plus étroit, en utilisant une heuristique de recherche optimisante pour obtenir la meilleure précision possible. J'ai également utilisé des implémentations personnalisées de log() qui ne doivent traiter que les entrées qui sont des normales positives.

Il faut faire attention à la fois à l'approximation de départ et à l'itération de Newton pour éviter les débordements et les sous-débordements dans les calculs intermédiaires. Ceci est facilement et économiquement réalisable par une mise à l'échelle par des puissances de deux appropriées.

Bien que les schémas d'itération qui en résultent donnent des résultats précis en général, des erreurs de plusieurs ulp se produisent pour les arguments proches de zéro et pour les arguments proches de -1/e -0.367879. J'ai résolu le premier problème en utilisant les premiers termes de l'expansion de la série de Taylor autour de zéro : x - x 2 + (3/2)x 3 . Le fait que W ₀ (1+ex)-1 sur [ -1/e, 0 ] suggère l'utilisation d'une approximation polynomiale minimaxe p(t) con t=x+1/e) qui s'avère fonctionner raisonnablement bien près de -1/e . J'ai généré cette approximation avec l'algorithme de Remez.

La précision atteinte pour les deux IEEE-754 binary32 correspond à float et IEEE-754 binary64 correspond à double est bien dans la limite de l'erreur spécifiée. L'erreur maximale dans le demi-plan positif est inférieure à 1,5 ulp, et l'erreur maximale dans le demi-plan négatif est inférieure à 2,7 ulp. Le code en simple précision a été testé de manière exhaustive, tandis que le code en double précision a été testé avec un milliard d'arguments aléatoires.

---

[*] Roberto Iacono et John P. Boyd. "New approximations to the principal real-valued branch of the Lambert W-function". Avancées en mathématiques computationnelles , vol. 43, n° 6, décembre 2017, p. 1403-1436.

L'implémentation en simple précision de la méthode Lambert W ₀ est la suivante :

float expf_scale_pos_normal (float a, int scale);
float logf_pos_normal (float a);

/*
  Compute the principal branch of the Lambert W function, W_0. The maximum
  error in the positive half-plane is 1.49874 ulps and the maximum error in 
  the negative half-plane is 2.56002 ulps
*/
float lambert_w0f (float z) 
{
    const float em1_fact_0 = 0.625529587f; // exp(-1)_factor_0
    const float em1_fact_1 = 0.588108778f; // exp(-1)_factor_1
    const float qe1 = 2.71828183f / 4.0f; //  0x1.5bf0a8p-01 // exp(1)/4
    float e, w, num, den, rden, redz, y, r;

    if (isnan (z) || (z == INFINITY) || (z == 0.0f)) return z + z;
    if (fabsf (z) < 1.220703125e-4f) return fmaf (-z, z, z); // 0x1.0p-13
    redz = fmaf (em1_fact_0, em1_fact_1, z); // z + exp(-1)
    if (redz < 0.0244140625f) { // expansion at -(exp(-1))
        r = sqrtf (redz);
        w =             -7.88867188f;  // -0x1.f8e000p+2
        w = fmaf (w, r,  5.89442062f); //  0x1.793e30p+2
        w = fmaf (w, r, -4.16159487f); // -0x1.0a5792p+2
        w = fmaf (w, r,  3.06752682f); //  0x1.88a4b8p+1
        w = fmaf (w, r, -2.35340595f); // -0x1.2d3c68p+1
        w = fmaf (w, r,  1.93658125f); //  0x1.efc3cap+0
        w = fmaf (w, r, -1.81217659f); // -0x1.cfeacep+0
        w = fmaf (w, r,  2.33164287f); //  0x1.2a7346p+1
        w = fmaf (w, r, -1.00000000f); // -0x1.000000p+0
    } else {
        /* Compute initial approximation. Based on: Roberto Iacono and John 
           Philip Boyd, "New approximations to the principal real-valued branch
           of the Lambert W function", Advances in Computational Mathematics, 
           Vol. 43, No. 6, December 2017, pp. 1403-1436
        */
        y = fmaf (2.0f, sqrtf (fmaf (qe1, z, 0.25f)), 1.0f);
        y = logf_pos_normal (fmaf (1.15262585f, y, -0.15262585f) / 
                             fmaf (0.45906518f, logf_pos_normal (y), 1.0f));
        w = fmaf (2.0390625f, y, -1.0f);

        /* perform Newton iterations to refine approximation to full accuracy */
        for (int i = 0; i < 3; i++) {
            e = expf_scale_pos_normal (w, -3); // 0.125f * expf (w);
            num = fmaf (w, e, -0.125f * z);
            den = fmaf (w, e, e);
            rden = 1.0f / den;
            w = fmaf (-num, rden, w);
        }
    }
    return w;
}

float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/* exp(a) * 2**scale; positive normal results only! Maximum error 0.86565 ulp */
float expf_scale_pos_normal (float a, int scale)
{
    const float flt_int_cvt = 12582912.0f; // 0x1.8p23 
    float f, r, j, t;
    uint32_t i;

    /* exp(a) = 2**i * exp(f); i = rintf (a / log(2)) */
    j = fmaf (1.442695f, a, flt_int_cvt); // // 0x1.715476p0 // log2(e)
    t = j - flt_int_cvt; 
    f = fmaf (t, -6.93145752e-1f, a); // -0x1.62e400p-1  // log_2_hi 
    f = fmaf (t, -1.42860677e-6f, f); // -0x1.7f7d1cp-20 // log_2_lo 
    i = float_as_uint32 (j);

    /* approximate r = exp(f) on interval [-log(2)/2, +log(2)/2] */
    r =             1.37805939e-3f;  // 0x1.694000p-10
    r = fmaf (r, f, 8.37312452e-3f); // 0x1.125edcp-7
    r = fmaf (r, f, 4.16695364e-2f); // 0x1.555b5ap-5
    r = fmaf (r, f, 1.66664720e-1f); // 0x1.555450p-3
    r = fmaf (r, f, 4.99999851e-1f); // 0x1.fffff6p-2
    r = fmaf (r, f, 1.00000000e+0f); // 0x1.000000p+0
    r = fmaf (r, f, 1.00000000e+0f); // 0x1.000000p+0

    /* exp(a) = 2**(i+scale) * r; */
    r = uint32_as_float (((i + scale) << 23) + float_as_uint32 (r));
    return r;
}

/* compute natural logarithm of positive normals; maximum error: 0.85089 ulp */
float logf_pos_normal (float a)
{
    const float ln2 = 0.693147182f; // 0x1.62e430p-1 // log(2)
    const float two_to_m23 = 1.19209290e-7f; // 0x1.0p-23
    float m, r, s, t, i, f;
    int32_t e;

    /* log(a) = log(m * 2**i) = i * log(2) + log(m) */
    e = (float_as_uint32 (a) - float_as_uint32 (0.666666667f)) & 0xff800000;
    m = uint32_as_float (float_as_uint32 (a) - e);
    i = (float)e * wo_to_m23;

    /* log(m) = log1p(f) */
    f = m - 1.0f;
    s = f * f;

    /* compute log1p(f) for f in [-1/3, 1/3] */
    r =             -0.130310059f;  // -0x1.0ae000p-3
    t =              0.140869141f;  //  0x1.208000p-3
    r = fmaf (r, s, -0.121483363f); // -0x1.f1988ap-4
    t = fmaf (t, s,  0.139814854f); //  0x1.1e5740p-3
    r = fmaf (r, s, -0.166846141f); // -0x1.55b36ep-3
    t = fmaf (t, s,  0.200120345f); //  0x1.99d8b2p-3
    r = fmaf (r, s, -0.249996200f); // -0x1.fffe02p-3
    r = fmaf (t, f, r);
    r = fmaf (r, f,  0.333331972f); //  0x1.5554fap-2
    r = fmaf (r, f, -0.500000000f); // -0x1.000000p-1
    r = fmaf (r, s, f);

    /* log(a) = i * log(2) + log(m) */
    r = fmaf (i, ln2, r);
    return r;
}

L'implémentation en double précision est structurellement équivalente à l'implémentation en simple précision, sauf qu'elle utilise le schéma d'itération Iacono-Boyd :

double exp_scale_pos_normal (double a, int scale);
double log_pos_normal (double a);

/* Compute the principal branch of the Lambert W function, W_0. Maximum error:
   positive half-plane: 1.48025 ulp
   negative half-plane: 2.67268 ulp
*/
double lambert_w0 (double z) 
{
    const double qe1 = 2.7182818284590452 * 0.25; //  0x1.5bf0a8b145769p-01 // exp(1)/4
    const double em1_hi =  3.6787944117144233e-1; //  0x1.78b56362cef38p-02 // exp(-1)_hi
    const double em1_lo = -1.2428753672788363e-17;// -0x1.ca8a4270fadf5p-57 // exp(-1)_lo
    double e, r, t, w, y, num, den, rden, redz;
    int i;

    if (isnan (z) || (z == INFINITY) || (z == 0.0)) return z + z;
    if (fabs (z) < 0x1.0p-19) return fma (fma (1.5, z, -1.0) * z, z, z);
    redz = (z + em1_hi) + em1_lo;
    if (redz < 0.01025390625) { // expansion at -(exp(-1))
        r = sqrt (redz);
        w =             12.2549150361339980;  //  0x1.88284393ee88ap+3
        w = fma (w, r, -15.0789952403151450); // -0x1.e2872106b62ecp+3
        w = fma (w, r,  11.8725001236137670); //  0x1.7beb856115aabp+3
        w = fma (w, r,  -8.3706388286143945); // -0x1.0bdc45f5f0d9bp+3
        w = fma (w, r,   5.8564285634808346); //  0x1.76cfb9bfe0ab2p+2
        w = fma (w, r,  -4.1752813239140965); // -0x1.0b37cf2873e15p+2
        w = fma (w, r,   3.0668577430022212); //  0x1.888ecb65d6e6ap+1
        w = fma (w, r,  -2.3535511874120569); // -0x1.2d412a51b2c8cp+1
        w = fma (w, r,   1.9366311143993924); //  0x1.efc70e84c2eccp+0
        w = fma (w, r,  -1.8121878856391300); // -0x1.cfeb8b9707071p+0
        w = fma (w, r,   2.3316439815971242); //  0x1.2a734f5b6ffbep+1
        w = fma (w, r,  -1.0000000000000000); // -0x1.0000000000000p+0
        return w;
    }
    /* Roberto Iacono and John Philip Boyd, "New approximations to the 
       principal real-valued branch of the Lambert W function", Advances 
       in Computational Mathematics, Vol. 43, No. 6, December 2017, 
       pp. 1403-1436
     */
    y = fma (2.0, sqrt (fma (qe1, z, 0.25)), 1.0);
    y = log_pos_normal (fma (1.14956131, y, -0.14956131) / 
                        fma (0.4549574, log_pos_normal (y), 1.0));
    w = fma (2.036, y, -1.0);

    /* Use iteration scheme w = (w / (1 + w)) * (1 + log (z / w) from
       Roberto Iacono and John Philip Boyd, "New approximations to the 
       principal real-valued branch of the Lambert W function", Advances
       in Computational Mathematics, Vol. 43, No. 6, December 2017, pp. 
       1403-1436
    */
    for (i = 0; i < 3; i++) {
        t = w / (1.0 + w);
        w = fma (log_pos_normal (z / w), t, t);
    }

    /* Fine tune approximation with a single Newton iteration */
    e = exp_scale_pos_normal (w, -3); // 0.125 * exp (w)
    num = fma (w, e, -0.125 *z);
    den = fma (w, e, e);
    rden = 1.0 / den;
    w = fma (-num, rden, w);

    return w;
}

int double2hiint (double a)
{
    unsigned long long int t;
    memcpy (&t, &a, sizeof t);
    return (int)(t >> 32);
}

int double2loint (double a)
{
    unsigned long long int t;
    memcpy (&t, &a, sizeof t);
    return (int)(unsigned int)t;
}

double hiloint2double (int hi, int lo) 
{
    double r;
    unsigned long long int t;
    t = ((unsigned long long int)(unsigned int)hi << 32) | (unsigned int)lo;
    memcpy (&r, &t, sizeof r);
    return r;
}

/* exp(a) * 2**scale; pos. normal results only! Max. err. found: 0.89028 ulp */
double exp_scale_pos_normal (double a, int scale)
{
    const double ln2_hi = 6.9314718055829871e-01; // 0x1.62e42fefa00000p-01
    const double ln2_lo = 1.6465949582897082e-12; // 0x1.cf79abc9e3b3a0p-40
    const double l2e = 1.4426950408889634; // 0x1.71547652b82fe0p+00 // log2(e)
    const double cvt = 6755399441055744.0; // 0x1.80000000000000p+52 // 3*2**51
    double f, r;
    int i;

    /* exp(a) = exp(i + f); i = rint (a / log(2)) */
    r = fma (l2e, a, cvt);
    i = double2loint (r);
    r = r - cvt;
    f = fma (r, -ln2_hi, a);
    f = fma (r, -ln2_lo, f);

    /* approximate r = exp(f) on interval [-log(2)/2,+log(2)/2] */
    r =            2.5022018235176802e-8;  // 0x1.ade0000000000p-26
    r = fma (r, f, 2.7630903497145818e-7); // 0x1.28af3fcbbf09bp-22
    r = fma (r, f, 2.7557514543490574e-6); // 0x1.71dee623774fap-19
    r = fma (r, f, 2.4801491039409158e-5); // 0x1.a01997c8b50d7p-16
    r = fma (r, f, 1.9841269589068419e-4); // 0x1.a01a01475db8cp-13
    r = fma (r, f, 1.3888888945916566e-3); // 0x1.6c16c1852b805p-10
    r = fma (r, f, 8.3333333334557735e-3); // 0x1.11111111224c7p-7
    r = fma (r, f, 4.1666666666519782e-2); // 0x1.55555555502a5p-5
    r = fma (r, f, 1.6666666666666477e-1); // 0x1.5555555555511p-3
    r = fma (r, f, 5.0000000000000122e-1); // 0x1.000000000000bp-1
    r = fma (r, f, 1.0000000000000000e+0); // 0x1.0000000000000p+0
    r = fma (r, f, 1.0000000000000000e+0); // 0x1.0000000000000p+0

    /* exp(a) = 2**(i+scale) * r */
    r = hiloint2double (double2hiint (r) + ((i + scale) << 20), 
                        double2loint (r));
    return r;
}

/* compute natural logarithm of positive normals; max. err. found: 0.86902 ulp*/
double log_pos_normal (double a)
{
    const double ln2_hi = 6.9314718055994529e-01; // 0x1.62e42fefa39efp-01
    const double ln2_lo = 2.3190468138462996e-17; // 0x1.abc9e3b39803fp-56
    double m, r, i, s, t, p, q;
    int e;

    /* log(a) = log(m * 2**i) = i * log(2) + log(m) */
    e = (double2hiint (a) - double2hiint (0.70703125)) & 0xfff00000;
    m = hiloint2double (double2hiint (a) - e, double2loint (a));
    t = hiloint2double (0x41f00000, 0x80000000 ^ e);
    i = t - (hiloint2double (0x41f00000, 0x80000000));

    /* m now in [181/256, 362/256]. Compute q = (m-1) / (m+1) */
    p = m + 1.0;
    r = 1.0 / p;
    q = fma (m, r, -r);
    m = m - 1.0;

    /* compute (2*atanh(q)/q-2*q) as p(q**2), q in [-75/437, 53/309] */
    s = q * q;
    r =            1.4794533702196025e-1;  // 0x1.2efdf700d7135p-3
    r = fma (r, s, 1.5314187748152339e-1); // 0x1.39a272db730f7p-3
    r = fma (r, s, 1.8183559141306990e-1); // 0x1.746637f2f191bp-3
    r = fma (r, s, 2.2222198669309609e-1); // 0x1.c71c522a64577p-3
    r = fma (r, s, 2.8571428741489319e-1); // 0x1.24924941c9a2fp-2
    r = fma (r, s, 3.9999999999418523e-1); // 0x1.999999998006cp-2
    r = fma (r, s, 6.6666666666667340e-1); // 0x1.5555555555592p-1
    r = r * s;

    /* log(a) = 2*atanh(q) + i*log(2) = ln2_lo*i + p(q**2)*q + 2q + ln2_hi * i.
       Use K.C. Ng's trick to improve the accuracy of the computation, like so:
       p(q**2)*q + 2q = p(q**2)*q + q*t - t + m, where t = m**2/2.
    */
    t = m * m * 0.5;
    r = fma (q, t, fma (q, r, ln2_lo * i)) - t + m;
    r = fma (ln2_hi, i, r);

    return r;
}