From octave-sources-request at bevo dot che dot wisc dot edu Fri Nov 22 09:27:01 2002 Subject: submission: Lambert W function From: schraudo at inf dot ethz dot ch To: octave-sources at bevo dot che dot wisc dot edu Date: Fri, 22 Nov 2002 06:14:25 -0600 This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. Send mail to mime at docserver dot cac dot washington dot edu for more info. --8323328-1616375311-1037967260=:1186 Content-Type: TEXT/PLAIN; charset=US-ASCII Dear Octavers, a few years ago I needed the Lambert W function, used to solve certain transcendental equations, which available as LambertW() in Maple and ProductLog[] in Mathematica. However, I needed it in a Matlab/octave context, where to my knowledge it doesn't exist. So I went back to the seminal paper describing the algorithm, and implemented it. The result is the attached W.m function, which I think would fit nicely in octave's "specfun" library. I hereby submit it for inclusion in future releases of octave. For more information on this function, see http://mathworld.wolfram.com/LambertsW-Function.html With best regards, -- Dr. Nicol N. Schraudolph http://www.inf.ethz.ch/~schraudo/ Institute of Computational Science mobile: +41-76-585-3877 ETH Zentrum, HRS H5 (Hirschengraben 84) office: -1-632-7942 CH-8092 Zuerich, Switzerland fax: -1374 --8323328-1616375311-1037967260=:1186 Content-Type: TEXT/PLAIN; charset=US-ASCII; name="W.m" Content-Transfer-Encoding: BASE64 Content-ID: Content-Description: Lambert W function Content-Disposition: attachment; filename="W.m" JSUgdXNhZ2U6IFcoeikgb3IgVyhuLHopDQolJQ0KJSUgQ29tcHV0ZSB0aGUg TGFtYmVydCBXIGZ1bmN0aW9uIG9mIHouICBUaGlzIGZ1bmN0aW9uIHNhdGlz Zmllcw0KJSUgVyh6KS4qZXhwKFcoeikpID0geiwgYW5kIGNhbiB0aHVzIGJl IHVzZWQgdG8gZXhwcmVzcyBzb2x1dGlvbnMNCiUlIG9mIHRyYW5zY2VuZGVu dGFsIGVxdWF0aW9ucyBpbnZvbHZpbmcgZXhwb25lbnRpYWxzIG9yIGxvZ2Fy aXRobXMuDQolJQ0KJSUgbiBtdXN0IGJlIGludGVnZXIsIGFuZCBzcGVjaWZp ZXMgdGhlIGJyYW5jaCBvZiBXIHRvIGJlIGNvbXB1dGVkOw0KJSUgVyh6KSBp cyBhIHNob3J0aGFuZCBmb3IgVygwLHopLCB0aGUgcHJpbmNpcGFsIGJyYW5j aC4gIEJyYW5jaGVzDQolJSAwIGFuZCAtMSBhcmUgdGhlIG9ubHkgb25lcyB0 aGF0IGNhbiB0YWtlIG9uIG5vbi1jb21wbGV4IHZhbHVlcy4NCiUlDQolJSBJ ZiBlaXRoZXIgbiBvciB6IGFyZSBub24tc2NhbGFyLCB0aGUgZnVuY3Rpb24g aXMgbWFwcGVkIHRvIGVhY2gNCiUlIGVsZW1lbnQ7IGJvdGggbWF5IGJlIG5v bi1zY2FsYXIgcHJvdmlkZWQgdGhlaXIgZGltZW5zaW9ucyBhZ3JlZS4NCiUl DQolJSBUaGlzIGltcGxlbWVudGF0aW9uIHNob3VsZCByZXR1cm4gdmFsdWVz IHdpdGhpbiAyLjUqZXBzIG9mIGl0cw0KJSUgY291bnRlcnBhcnQgaW4gTWFw bGUgViwgcmVsZWFzZSAzIG9yIGxhdGVyLiAgUGxlYXNlIHJlcG9ydCBhbnkN CiUlIGRpc2NyZXBhbmNpZXMgdG8gdGhlIGF1dGhvciwgTmljaSBTY2hyYXVk b2xwaCA8c2NocmF1ZG9AaW5mLmV0aHouY2g+Lg0KJSUNCiUlIEZvciBmdXJ0 aGVyIGRldGFpbHMsIHNlZToNCiUlDQolJSBDb3JsZXNzLCBHb25uZXQsIEhh cmUsIEplZmZyZXksIGFuZCBLbnV0aCAoMTk5NiksICJPbiB0aGUgTGFtYmVy dA0KJSUgVyBGdW5jdGlvbiIsIEFkdmFuY2VzIGluIENvbXB1dGF0aW9uYWwg TWF0aGVtYXRpY3MgNSg0KTozMjktMzU5Lg0KDQolJSBBdXRob3I6ICAgTmlj b2wgTi4gU2NocmF1ZG9scGggPHNjaHJhdWRvQGluZi5ldGh6LmNoPg0KJSUg VmVyc2lvbjogIDEuMA0KJSUgQ3JlYXRlZDogIDA3IEF1ZyAxOTk4DQolJSBL ZXl3b3JkczogTGFtYmVydCBXIE9tZWdhIHNwZWNpYWwgdHJhbnNjZW5kZW50 YWwgZnVuY3Rpb24NCg0KJSUgQ29weXJpZ2h0IChDKSAxOTk4IGJ5IE5pY29s IE4uIFNjaHJhdWRvbHBoDQolJQ0KJSUgVGhpcyBwcm9ncmFtIGlzIGZyZWUg c29mdHdhcmU7IHlvdSBjYW4gcmVkaXN0cmlidXRlIGFuZC9vcg0KJSUgbW9k aWZ5IGl0IHVuZGVyIHRoZSB0ZXJtcyBvZiB0aGUgR05VIEdlbmVyYWwgUHVi bGljDQolJSBMaWNlbnNlIGFzIHB1Ymxpc2hlZCBieSB0aGUgRnJlZSBTb2Z0 d2FyZSBGb3VuZGF0aW9uOw0KJSUgZWl0aGVyIHZlcnNpb24gMiwgb3IgKGF0 IHlvdXIgb3B0aW9uKSBhbnkgbGF0ZXIgdmVyc2lvbi4NCiUlDQolJSBUaGlz IHByb2dyYW0gaXMgZGlzdHJpYnV0ZWQgaW4gdGhlIGhvcGUgdGhhdCBpdCB3 aWxsIGJlDQolJSB1c2VmdWwsIGJ1dCBXSVRIT1VUIEFOWSBXQVJSQU5UWTsg d2l0aG91dCBldmVuIHRoZSBpbXBsaWVkDQolJSB3YXJyYW50eSBvZiBNRVJD SEFOVEFCSUxJVFkgb3IgRklUTkVTUyBGT1IgQSBQQVJUSUNVTEFSDQolJSBQ VVJQT1NFLiAgU2VlIHRoZSBHTlUgR2VuZXJhbCBQdWJsaWMgTGljZW5zZSBm b3IgbW9yZQ0KJSUgZGV0YWlscy4NCiUlDQolJSBZb3Ugc2hvdWxkIGhhdmUg cmVjZWl2ZWQgYSBjb3B5IG9mIHRoZSBHTlUgR2VuZXJhbCBQdWJsaWMNCiUl IExpY2Vuc2UgYWxvbmcgd2l0aCBPY3RhdmU7IHNlZSB0aGUgZmlsZSBDT1BZ SU5HLiAgSWYgbm90LA0KJSUgd3JpdGUgdG8gdGhlIEZyZWUgU29mdHdhcmUg Rm91bmRhdGlvbiwgNTkgVGVtcGxlIFBsYWNlIC0NCiUlIFN1aXRlIDMzMCwg Qm9zdG9uLCBNQSAwMjExMS0xMzA3LCBVU0EuDQoNCmZ1bmN0aW9uIHcgPSBX KGIseikNCiAgICBpZiAobmFyZ2luID09IDEpDQogICAgICAgIHogPSBiOw0K ICAgICAgICBiID0gMDsNCiAgICBlbHNlDQogICAgICAgICUlIHNvbWUgZXJy b3IgY2hlY2tpbmcNCiAgICAgICAgJQ0KICAgICAgICBpZiAobmFyZ2luICE9 IDIpDQogICAgICAgICAgICB1c2FnZSgncmVzdWx0ID0gVyhicmFuY2gsIGFy Z3VtZW50KScpDQogICAgICAgIGVsc2UNCiAgICAgICAgICAgIGlmIChhbnko cm91bmQocmVhbChiKSkgIT0gYikpDQogICAgICAgICAgICAgICAgdXNhZ2Uo J2JyYW5jaCBudW1iZXIgZm9yIFcgbXVzdCBiZSBpbnRlZ2VyJykNCiAgICAg ICAgICAgIGVuZA0KICAgICAgICBlbmQNCiAgICBlbmQNCg0KICAgICUlIHNl cmllcyBleHBhbnNpb24gYWJvdXQgLTEvZQ0KICAgICUNCiAgICAlIHAgPSAo MSAtIDIqYWJzKGIpKS4qc3FydCgyKmUqeiArIDIpOw0KICAgICUgdyA9ICgx MS83MikqcDsNCiAgICAlIHcgPSAodyAtIDEvMykuKnA7DQogICAgJSB3ID0g KHcgKyAxKS4qcCAtIDENCiAgICAlDQogICAgJSBmaXJzdC1vcmRlciB2ZXJz aW9uIHN1ZmZpY2VzOg0KICAgICUNCiAgICB3ID0gKDEgLSAyKmFicyhiKSku KnNxcnQoMiplKnogKyAyKSAtIDE7DQoNCiAgICAlJSBhc3ltcHRvdGljIGV4 cGFuc2lvbiBhdCAwIGFuZCBJbmYNCiAgICAlDQogICAgdiA9IGxvZyh6ICsg fih6IHwgYikpICsgMipwaSpJKmI7DQogICAgdiA9IHYgLSBsb2codiArIH52 KTsNCg0KICAgICUlIGNob29zZSBzdHJhdGVneSBmb3IgaW5pdGlhbCBndWVz cw0KICAgICUNCiAgICBjID0gYWJzKHogKyAxL2UpOw0KICAgIGMgPSAoYyA+ IDEuNDUgLSAxLjEqYWJzKGIpKTsNCiAgICBjID0gYyB8IChiLippbWFnKHop ID4gMCkgfCAofmltYWcoeikgJiAoYiA9PSAxKSk7DQogICAgdyA9ICgxIC0g YykuKncgKyBjLip2Ow0KDQogICAgJSUgSGFsbGV5IGl0ZXJhdGlvbg0KICAg ICUNCiAgICBmb3IgbiA9IDE6MTANCiAgICAgICAgcCA9IGV4cCh3KTsNCiAg ICAgICAgdCA9IHcuKnAgLSB6Ow0KICAgICAgICBmID0gKHcgIT0gLTEpOw0K ICAgICAgICB0ID0gZi4qdC4vKHAuKih3ICsgZikgLSAwLjUqKHcgKyAyLjAp Lip0Li8odyArIGYpKTsNCiAgICAgICAgdyA9IHcgLSB0Ow0KICAgICAgICBp ZiAoYWJzKHJlYWwodCkpIDwgKDIuNDgqZXBzKSooMS4wICsgYWJzKHJlYWwo dykpKQ0KICAgICAgICAgICAgJiYgYWJzKGltYWcodCkpIDwgKDIuNDgqZXBz KSooMS4wICsgYWJzKGltYWcodykpKSkNCiAgICAgICAgICAgIHJldHVybg0K ICAgICAgICBlbmQNCiAgICBlbmQNCiAgICB3YXJuaW5nKCdpdGVyYXRpb24g bGltaXQgcmVhY2hlZCwgcmVzdWx0IG9mIFcgbWF5IGJlIGluYWNjdXJhdGUn KTsNCmVuZA0KDQo= --8323328-1616375311-1037967260=:1186--