lambert_w 2.0.2

Fast and accurate evaluation of the Lambert W function by the method of T. Fukushima.
Documentation
// Copyright 2024-2026 Johanna Sörngård
// SPDX-License-Identifier: MIT OR Apache-2.0

//! This module contains elementary and rational functions used in the Lambert W function approximations.
//! They are generic over all types that implement the [`Float`] trait.

use num_traits::Float;

// The `inline` annotations are motivated by benchmarks on an AMD 5800X3D processor.
// The `lambert_w0` function is sped up by around 50% when the inline annotations are present.
// The reason it's not `inline(always)` is to let the compiler have more of a say on platforms where
// it might decide that inlining is not a good idea.

/// Evaluate a rational function at `x` using [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
///
/// The coefficients are assumed to be sorted in ascending order by the degree of the associated x-term.
#[inline]
pub(crate) fn rational_function<T: Float, const N: usize, const D: usize>(
    x: T,
    numerator_coefficients: [T; N],
    denominator_coefficients: [T; D],
) -> T {
    // According to both `cargo-show-asm` and Compiler Explorer these two polynomials are evaluated in parallel using vector instructions.
    // See the Compiler Explorer link below for how the `lambert_w0` function in this crate gets compiled.
    // https://godbolt.org/#z:OYLghAFBqd5QCxAYwPYBMCmBRdBLAF1QCcAaPECAMzwBtMA7AQwFtMQByARg9KtQYEAysib0QXACx8BBAKoBnTAAUAHpwAMvAFYTStJg1DEArgoKkl9ZATwDKjdAGFUtEywYgATKUcAZPAZMADl3ACNMYhAAVi8ATlIAB1QFQjsGFzcPbySUtIEAoNCWCKjYhKtMG3ShAiZiAkz3Tx9K6oFa%2BoJCkPDImPjLOoam7Nbh7sDekv7ygEpLVBNiZHYOAHp1gGohBCXadC2IrYB2DS8AOmiAZkkADjOANhP7k%2BeuOK2mZDRifCMtkQAKQaACCmy2CAIBESChAmwA7kiLgjXFRiKwxIkEEwLmgWOtAokTAQgdcAGJ4MkAET8rFKBAA6kCvNFogAhDQs2JOFns7leO5cAC0Xm50WwwpuAriACoxayJVxzpIBeSFWzqSDwdtDIdCFs8ApTucrrcHhpnq8rdriWEtlQGFs0AxzAB9TAANzEHtUrES9Agcy2wrJ2Adj0kWyBJ35YK2Ca2BhmDTdCI01EjIBAAFlQQANLbbS4aObamNasG2kz2x1J%2BmRAhpjMALxAEckwdD13DVEj0dj2sTW3Q6eQEBb0a87K2XAuGiLzpI7BAfck2Zd5jhIGwZarJ0roO1m22LeQ0eu1K2k75s/WmG1ddHGnHZ/ba67YY7A7joOHeCoa9zzJJwtg0ecfyHYdEzXbNglBYIoITCstkwWglENQCzwvXlLy2S5okkOI3TuOJrjdOJLTdaINHIyQTnIx5IPjaCEwhRkcJpMCLlI65SC2fM3S4csWNYjFbAEH0qBMBh2gYCAkNYxMzwuBQAEcGiDUhFKU6MOR03TE2FcC4lMiizNM4jLPMqztNEwzh1FC4GKsx5HhItySK8DQ3WuNy7L/BzoPA853K4SM3W8k43ROLx3NMnzpAMhzHgua4zgo8KSI0aLSMkCiNC4AKgugk4LmVLhfMFN1HmiSqTjuaI3UkWjipKxNLjiO4fI0O4vBqhrfI0HzHka3y2va7j7juIS4iauJRqE9KSK6mKJvakKvHIk4bjdWimoa6LrmOmqiuSwyQuGvaTiO2qYseHqyJi2Lmr487hyBaJqXW3TPt/SaEznH6guiecTkqx5lRI%2B5mq4aLJEFYGHLiC4uvcu5rhmoUZuuQjqMR96lJRtyfJ2pqPn6qKSIeJHDMkZy4ny9L%2BtI4iuBOyGfEJ1i53C/q4gy25opuxLYnG7ngvKrgZqI/LHlZPact85a2Vp3TLro65%2Btq/q6vIrwWsirxdbVxNPu%2B8690C5CD1Q9DMEwoDOLwyRrguSHyLiA3mvlmrHmI9LmOt6D2Odq9Lh26QBMikTg%2BHcT0ikmS5IU%2BzdJU9TNIWCW/olxyTIswbCpm4aDtIoTTdY4znI0fKuGiaK6vm2uSL8yvoMuLgtqN9zev17rIrFu52%2BHS5MeI55e4y%2BJyIbrm06C8Cbui2ryK4aW9poyK4kq4e88TELy9G9z/Zm0zGMkfKdpHg/51oob3INmaG6a9K5aSheHMunqFeVdy3kSicNmWsSKxHbubcB%2BlP6GSBvvBMbt9o1XODVPq2Vn6SDOtA3SbspB0V2oVcm1Vnh7ywUpOczNfL1xqtEE%2B0RiIGyvjfBM4FWQy0kPFF4QkG7NQRhXOB3FCrETZJVPmg9d4fD2pIaITCBFXVeDFOhyt%2BqX1foVSKjwZFjyynVCGdw65K2uFQ6ImBhQaJzl9WmVsPq2zQhhACTsQJcToXcC4UjZ4nXBjNBqLMg66VDo4l2PEvBFWjtcWOukE6SVoG6aSskJLyQlhnDSBAtLmP%2BgDauXUJ6cO8uRChDFpH8PAnQpqDx0ZUUtJ7V2JEZFzgRrrXqQ1xHk0ZuLUhrFwKXyxtVUiuj9aYz4e0yWtd8o0XyozPJIDriUT2m9IZw5LqVViLvIBzUMbUQwRRbhJC46Lzvj1LKOVKaPFnjQ6iyDjYyLdhQ46R1kG0RGnVExOySoQLSfwucGhanuzmhRQ2W1iLsO1hIk4MjwLdQnmNBidc2SKwSEU8qhtzg%2BVuFZVacMfJbTBXfAO3Dpa5SIr5VZhUzHzNvsNTFIC6r5QePNE6pFhbz12Q5Oc3Ve5w2uuTLe8Q4gmNJcywyBEdp7U8gxbaai6H9WuCY%2BFQy3nQKsWbGx9tHbYQCVeEuRFXE7womNHKB1IbUV8Upfx1xcJXjKrXEJglVRVmgZE5g0TYkp0ScgVSyTUlyqgQK3S1cbprwijvKyL1LRSu%2BezYBeSt4I2ySi7FrsRGGrhrrMWGKYofx9UpEKuNnrkRyjjCRtFG4aEKWS5h%2ByKIRQwe5FqpSS3URmbKzNrEXECylWNKRxE7jPGog8Plmiri%2B1xiiqQsN6ovWMcKF5AM6m3GppVdKBKIYnKEiY9ekCLZevSZNT52KXjuTniRThNDyaX0itiruj8BbPRRXOj4s9sXDQ4TCyq8R5pnMMdOya38hKG14r5L2PDTmnO%2BecA61TBQFqavEJm/b%2BGXHiJixa4NiKFUql7aVwom0Ay2CjLwqyEZrzUVrGtk72YbssYpFCtiHb2LVWari68REG1RnjLur6OV9WNaxU15rcNpW6vxQS0RwlKQdUnOJ6RU7NqAu6rOkDt0bXKt5a6lUSUxSMX5am2KMazzis1RpPKhLVTmTJhZYMfJSH6jdeaWVGY4wzThxZay5Y5qAegxuXdmqPquoVQBXjuF%2BTClDeupazMdSuODXt9DVpnGisqbazzvlkRlsEnh0UFojXos1NdmDwsJkuHoq%2BLdKKe3Zn8rga6wuvIse8sts55zYtZJ7HtmMZZnI%2BK%2B0FCLCos2mZFAZlEjrXBHV%2BpTFLIrIJ2ttD9v8NBWTGyVXmMM31IKOmenekg4P1a0VjZalofI0PhvprbWHvkI3i2Zat1EclnNO13b50Q7gy3Wf7fKALmqrMeGu/lNXN0CsVTbK8NHVXAQY3hOcT3oq135kKSbM1jbEROJgJiMZFPDl44x1l82hM1VE2JJg8SJMuqGUk%2BTdX8sNdiv1deiVqnn2eiI7F0y2a/MZizCKeiWY%2BYy896ieNCr5UtM/VqPXfO8Oe0j3hzXB4ZaxfwsqQo5aoOxpvWHyOp0yPpl1dDK6ymvSauwu4uWZEozeE1XqZ8osIwPZh%2BuVya7QbMrFGayahIvW%2B8KD4FGKc4d3T1zhw14YDPXiXFTqsxdqZOsOrhI1/bPTrlr1xWnep1tni9HalXTEyNSt47er7DE1WWsqU72GAZu1GprJqzxSYrsZryr3TKcOpT7qdd9vd9MPXu79mdPFR1owoidZxxnrMmO8j7hVVHlV2KwmDvjY9Nn0R8noyzDwNnfZONxkO2wOLqtnPTEbw9o4nHx9BcTTrk7xOkyVMnKTs5bv4W7Ha/NC8tSldUqRTML14x3nRfTpWbsS5ut6sf0doT5X8qIdpUMlYy8d15x8MhIHoKJu46F4pr1jds9%2BF6YGpX5vY3hG4non010vkFceIQFuohcbo9oEDCINcpAHcJkho8lhUMY2Za411Fsgo3Y%2Bp%2B5H5w9PJzgZpPd5cvV/sat0cgpYEQD4DLMvMupX4BlGp8pL5HMAYf0cor4vNBQwo8YMZLNVDJoXF7gepFFuoUV49MZMNOCHI3YMEEdeEME81lp/ZKp7s8scN6ZnsT4FZf9K1UNfkyNqsAZO5jZlZbptZ1l2FZ411YD2o5w3J3tVojtjNe4ntfIx8DDDJ5UAcp9gcVU6M58uI3Y2QV5stXZl8rpTIa1MA7gt8Mcd8w5ZxUpu1pFo47hT945CdE4L9JMBBr8gpb9PUZNc56sW8JFap3t60aFHke9v15xdoiIT4Xo9F6p0DZjxtfMIpbhl9CV8NX5tEB0Ph4p5s9ooZhohEFoYokt4MrhAN3MVpxU8kqI7cm9e8coMsJEhQwoB565R8vd1ils0oi4HCTiKoYpy4G9gkHd0p0Zr1DELdbtJ18MJ9hjvU/cLgfNoMB4FocYMpCMwjvN%2BFQZlR%2BoS1LM0sUstlPcATJCrg1EMZX14gjZod1ks9YjATkiRY/YL4F0N4MDyMiSGZ6EW4u5yJekhJnsOCc8rhytDFocqDaJSSGox9TN0StZDtGlCpPY0ttiaox9rDoJsjdJAcBw8iZ8HFwcrxQY2RLMot5YJ4so/JpVPg0dzpMcIcUYjdcc4gOjExz8YlL8pNXU5M78FN%2BEjDqpQpQFHordH0OVQsFF0FH5VlwZVSAZQZW83JSTPIuoX8vBttKcD8c0iJzC2Z1lWEcsvcSD6s5xBR9F6pR0aFRk0i7dMjATbh0MoZJ4KJj0%2BdPdvcsC0oQt489EmpCJ6Etox8gjJpLhIxKZqlbhKZGkhQG9DEUTxCPkMTI9fJuFvIRphVO1qIRtBlKdrldpYgmYYYrMap8zhRgDKcCIVN/Z4pvY5yUjKz11%2BFc9yt5ZtYECRtUi7caSUo0oQF6k/YS4O1FphD3De8hQWYsoXgepiMqDhDgLDIUYgtJs8160hRGIMDjp1ylITTqN8jZ9GiXFu14Y8T1kFoDplTlQ6jEx3Tw5O4Hto5lRfSEx/TnUr9gzM5QzfcZ0TJC9KkBtAFlEupMSjZlFu5P9qDmyH1H8hTTpDsotpY1N2YTF7ycMXEHkbtcU8lNlaJ%2BzXjJpQZfCBZn5fYzhUN%2BydLe8bgEd7kc05onjMNLlNzjYsZWkME/z0Z159Tvlu0D1GkvZocBCvNMNcYiKHIRjKcpDKcf1bh%2B5ocFYdEezSlDj%2BsbguVd5dpntYNNdByS1lEOVi5uSuEA4100yZznIN5ryE0jz08s9pZpS3EVohC9cc1xyMik9IxAVDV0pn4%2BdSIhCTERtpSRlJtlEtNYorJpkJqHK4ratJ97JSLzT6M%2BM6lssboyyD0EDEN0iu4mK2IGi98WMgFrUhJhI7UBUeLAy%2Bj%2BKPV79USJCWVijkFGDepl8noHhYK5jhollqoojqERo515YT4B1wZq9LiCkvJH451J0DTdI6l15ZlAV60u53sspbyBSdtB01LNZmN0MGIYjvk3yu4A4XpL5dZ6UpyHctpJLlD49cFe1TstZvlvJDdPI4o65f5Bpor2SlIjSgp4r0SfMeooUJUvMHhUNy52FvljpmNlDqkX51F6EZVpT8MfLlEz1aJn5cYqtQMYZhqSIqFjZoMVNTtx8vzypWlUVzJvDxlwox9lq4ipZLMzkGpe4515td4Jrar2pijy4GJXdZbjzSJqjhQMFYroISLp9aNyKLrXF0apBSkV1WQRFGpMB2ZTqixzrLT8IyoaIvkOKxQ7qIkuiokAzeiElSc3UBKhiNyazUou4gCwoe1XZUNwpJbXpGJ75Gp6oIoa81pNzu0J5DZMZ4tUF5tbzhaHIx5WbVkG4T4ZkpAvFjbNyyTKEyztZLj4TKyoTNzIwwo7jAskDuyPKUalIypfC%2BoV5GkhZmpYgJr0LdIyppZRly5aImYMpwoNdL447fo0Te8%2B7zgRE8ZFR1EiN%2BY2ygoXFlQRp0bCJ3splBQTd%2BFPSjFVoHpKY7MG84YHcHhdYvMgFMUmSbgVSZFQYTkWbyYUHuES0NdJqbj7he5fllRUs1NaIlrktBo9DPttpeEbpN6Y7pysjVqcj1rE7QdGj6Zeokd75Ix3skDeplDc6oxXToEWKtgRsLgtoNEOKwlK6xNq7HVa6ScZNBjXqW6EqXEgF3IRsREpi512G%2B714u1fLo9Lj2E5YF6YEDH2ECTn8PIV5u10iqziGKSpEUE60o86D3aSp6ZlQxyLDCCFZvIMCDZHtFFYg15Ymu9oZo61zx6qFwYjpYURthZTIJqAnsFUYB4gEkc3hXpgVAGb7kIpG7GJbBTwo2C15VoMZ5ovYCyPCyCU0Z6pVC90pWqzKPb0pSZRrOEXCapAN%2BTOnoJyEDyzlU9FY9bsmJHkmeI4dpYEdr0aFvjA6aGDG6oPI64VdyrgkTFY7t6JFgl0NFFDE6JhVPciJgHEwE6zSk6LS%2BNLVYU6F9CV5BpTJkd6587dGtZUYFpccpAuKtgHq67%2BiHIbGwydtIcrons%2Bs49ht4Gv4K1FSeF0MINlp4b8MkmEGwYa0D7NlK8FEyN5mSoLKYYit%2BdkHKZSJybNybga0uyN5Yollu0x936yFXF74tY14mSzhKUxHbgGCfC198MC81LAGg6SoSxfZKIhCBlLbc09SpROWRbumxbQG4Cua7mbhh73sFyc7MCCaGTZlXC0t7h6ocoarkteEsLIwIYorIwx9qzCz6qOcIY5Fym/IgqSqFy50XIJTPZfK37QMHXR0pAmZWS%2BSXmyW6ZUZC8Bq9onpCJGJANPdCIAWEwgW7ZNrCiIcDGtMSZ2WYpr1ngWZc7N9tGBVEXQY9rUWRNTGCciceirGb9G6Xq8WHz6YtpSkvMbgyqOEZWW13Ysp65DszIG5wmxn0zyp5Xql65vj0Z8N/XwzUYIoPssKntxkw3RQAb2prSFZXZRXXcTiXhbzBgvKW4aJop8NlF60DtzX0oZE9LlojMKnqEIjTt6JQMStR1eoON35qqpQ129IxCbX3rAmHdho5ZWmv7Ptxl7htKKbjyd5EoEC9nCJo76CbjRp3sECEtqCvGPLLWWVB1dcAPM7YnntXcJqI2cM3ZgkA0z575aIkdcG6mB1hp%2BkJ5vYbha8ZYXnNnR5bnSSECpEeaV5joTEnta2th62QcCjGjO4LJlpzQO215c7ai%2B2/FC6%2BNMYrg%2BpUXHh0XMXJ2Bjp3ycH8dtUpTJXc51n6up4ZC3dIEM7nTWbp/5dyQb93JocFU1u59VdU17/411wvZXDFXD74G5Sa7L2PksqFiEDM8lLjDHqIJqn2SoUYo6YoroRtX4z1SOY6avaSmzJF20u0z09E6J9OhOAZ6Z5YQ2A4W5LQa0%2B1TF2vWJRaVqcPUatz6tP7dye0TlocZkdoMC1PEw8Mt4d6Ho2YUGdofsHd5ZPY4c5oEdBodpG4%2BrBSKqbpX1CVMYeoEZrm2Gzknsz7uld4/WY6OPDIK81a7o%2BuAM2Y3J9PAfdIzcINCVLQhdNNTsTlDPjOyLQWiir3Db6FPi3G1FvsPgEXHOuIME5Wo5BI4YPPzHic%2BKG6Qzm7sP4Mzc9F18Gu%2BtSabmgF4tWnu7Ip4y9O7zksqJQvIpVk4pLM5ot6ayeJPILQ9o189n9cPKsvWJUoEYTC1Nflbky2OaVfSoeIlYsyfYLc4cGUC2B1KJ%2BY9d49/Zq9C8Nca2bjBcJToc%2BOzlIYMDapDOPpbWPaKbANL4nGN5a4OFDtrj6tiTBYnomqu5DdCI10GWHJsC19659YXoT3FYNc4pHtPJL5Dtyy8Ym4JqjnaShReTtpVotogv9YXmS%2BQLbhDcqFBC/khF3v0Ok8g9t4a1kELyeE%2BVE/DTrXWI0fG3Gi0brMlZqV8Z4Y7p8yS8iethd8i6toD8DOtg5AJSqfx3LHafrG3VaB5JbHGfpCXg6451TI5Y1E4pcp40kOwKtNpYTC2Y9fzNxfjNq9WWe1oVEDH1vXo0mG2cSBGDH3TqheJu4N0bWNlhLRswk8uMGfrinRiF5AuCXD2pVwApIJSkgGNNBrh24JgUYgiWGFfFaYnIfIVRNdNDzNhD9JGC3WVszjnQPBw6%2B5AjpZl/4zZEBTMVpA3BGjIkxch2aXB%2BmPKKNroR0JPK7FyirRS4suStLeQw7Dhay1SWqOgiRyy8s8uA/jOcF1ioJvu90NeFK0bwyIY6BjetA1EYjUcoYbaIVqIUowyNgWC/BMBvw%2BBb9uiO/IMnTwP4M95u/CauJDFPT81Xch1TttimeBSciMC6LTCNgDiPoj4Y0XGOKm/qlIKBHSO%2BEsirStIN6AGM%2BJjBkRlRv%2BCUF3ubV07Ukua5ccKFfDxhuQrIatcPpTktQ/0qUQzOEnDlMqm4eImyAWKhhOgXYeyeNDDnN2oGblH0cJK6HFG2hQwWorhF/uSnx5Qw/IAHThLy28iyDEwiuPApwiQZDRhYKAwEjBmQJhRI616dmHLGqE4ZLUuDHVOxkHzvYMY5AloVzhvIi4mOE/bvAYMhyXFnesUaDJfXrTCE9efQ4cEC3LAHhY4LIMJByECC0ApgEAMQAiCYAABPBQGWC%2BiPgnQnneJCBE3AEAtgwQdsGYDwAthMA/EDEVsGpA4jUg%2BIsMBAHUAdh%2BIDAdwJEEJwkA3QaATAFQBoDIA8AjAAgHCD0jsg1wZIGcIhAsQjhGAqAFgIEAZHEAmRqAFkWyI5GCBuRf0PkdcBnBagvon4HsN%2BHs6Jh6AmI2kWwHEgkALwV4ZILQFhEMBRReAMQJSJpF0iDRko5kayLwDsjORCI/kafh1HCjzRYo5gEQGIBGitgJos0RaKtGqB%2BIWAL0eKN9FSiZRTouUVyLLDKjT8eo%2Bkb6MXDhiLRPokgICMPDlgtof0cEVMHNjIiAxrgIMWKLEDoiBA5gLEaSLxEPgewlI98JGEJHSjHRzo%2BUe2EVGRh%2BRWI82OqN7D9gtRCYB0bKJdHnQLggQIgG6EICRAgw444gF6DnHQILg/AA4BABCj8QYwTgb4MgEJFbivgPwLYLKC2CqApwM4ZAFbArAcAFgtATgNEF4CeAOAWgUgKgE4AAAlMwJiIUBLAVgDsEETwFIAEBNA14hYAAGsBg9VQqpGBLRaxpAt4jgJIAfHASXxnAXgHCC%2BRASnx140gHAFgBIB8QiQOgJEHICUACJREqIMQA%2BAtQ%2BAdAAgJEDhAQAwgyEsIOKOICwjOAAElicwDYkAB5MINoCqCYSAJ%2BINgIIB4kMBTRyErAGEBMDABtxtAdCBxN4BYAWAhgYAOICwmkB8AC46oJ6EwBwhNJmAVQFUBJBrAAJk4tCMhIhFhAMQbElwFgGQkEBiAeAFgEpNIB6TiAYQFIJgGpCYBVJRgCEUYGAkLAqABgYAAoAABqHIhEDxMSCMA3J/AQQCIDEDsB6CiU%2BQEoDUDITdARUAwEFNMDmB9AeAMIHCEgALBUAiQeJAZOFA8TeAqADyc5KwClSoAzANgCAEwD4B4k7ksQCYDWDeQDYIEywGhEEnpAHADAZwK4GaB6B/AUwYoKUD0DJBUg8SUYJ4CKiLT8gDAHoHNP6BFQ2g8SToCMEmnZBdpw0uSAdMmBFA%2BgUQXaRMBWl6BzAXQLaVdIkALBvxywVYC9P0B3ikJmk18RwAxafjzwc4RmBBAgC4BCAho/8XMF4CYStAcwBYAgEwBMAsAUQLSOBPogGMy%2BuVRmJPC%2BkISfpz4v6WhJAAYTgpN4zgF4AJl1TUJgEsme5PomjTJAQAA%3D%3D
    // Note the vector instructions with names of the form "v*pd", the "pd" at the end means "packed doubles".
    // That they are operating on the 128-bit xmm registers means that two 64-bit floats are operated on at once,
    // and thus this is a way of formulating this operation that makes the compiler likely to vectorize it.

    let numerator = polynomial(x, numerator_coefficients);

    let denominator = polynomial(x, denominator_coefficients);

    numerator / denominator
}

/// Evaluate a polynomial at `x` using [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
///
/// The coefficients are assumed to be sorted in ascending order by the degree of the associated x-term.
#[inline]
fn polynomial<T: Float, const N: usize>(x: T, coefficients: [T; N]) -> T {
    coefficients
        .into_iter()
        .rev()
        .fold(T::zero(), |acc, c| acc * x + c)
}

// The functions below are wrappers around the [`num-traits`] crate,
// mainly to ensure that we can always call them regardless of the presence of
// the standard library. I do not just import the trait in the files where the
// lambert w functions are defined because the standard library is available during testing,
// which means that the crate would produce warnings about the unused imports.

/// Compute the square root of `x`.
///
/// Just wraps the [`sqrt`](Float::sqrt) function from [`num_traits`].
#[inline]
pub(crate) fn sqrt<T: Float>(x: T) -> T {
    Float::sqrt(x)
}

/// Compute the natural logarithm of `x`.
///
/// Just wraps the [`ln`](Float::ln) function from [`num_traits`].
#[inline]
pub(crate) fn ln<T: Float>(x: T) -> T {
    Float::ln(x)
}