(28/48/50) Lambert W Function

03202023, 08:43 PM
(This post was last modified: 06282023 06:23 PM by John Keith.)
Post: #1




(28/48/50) Lambert W Function
NOTE: The first and last programs have been obsoleted by the program in post #22, which should be used in most cases. That program gives accurate results over the entire complex plane for any branch, positive or negative, but accuracy is limited for values very close to the branch point 1/e.
For very accurate results near the branch point in branches 0 and 1, use the program in post #28, which is a translation of Albert Chan's HP71 program in post # 19. The following text refers to the (mostly obsolete) programs in this post only.  This program computes the Lambert W function W(z) for branches 0 and 1. For branch 1, the program covers the realvalued range 1/e < z < 0. For branch 0, the program covers the entire complex plane. For most values, the results are accurate to within 2 ULPs. For values very close to 1/e, the results are approximate. The program expects the branch (0 or 1) on level 2 and z on level 1. Some notes and background information: The gold standard for Lambert W programs on the HP 50g is Egan Ford's SPECIAL.LIB. It is very accurate for almost any input value, and can return results for any branch. It does have some problems, however. It is an HPGCC program and only works on the 49g+ and 50g. It requires two large programs on the calculators SD card. Also, it is rather slow. I believe that the slowness is caused by the need to fetch the executable from the SD card, as one would expect the HPGCC code to be quite fast. The program described here is comparatively limited since no branches other than 0 and 1 are covered, and then only the realvalued range of branch 1. However, it is about 5 times as fast as SPECIAL.LIB and about 9 times as fast as programs using the calculator's builtin solvers. It also runs on all RPL calculators. If anyone here has an idea of how to get an accurate estimate for other branches, or to improve accuracy around the branch point at 1/e, please chime in! The code that provides the initial estimate for branch 0 comes from page 8 of Iacono and Boyd. The estimation code for branch 1 comes from formula 27 on page 12 of https://www.sciencedirect.com/science/ar...via%3Dihub. The recurrence used in the main loop occurs in both papers and is mentioned in Wikipedia. Many additional references here. Code:
Also, a streamlined version for the 49 and 50, 25 bytes smaller: Code:
Now, here are a couple of experimental programs for those interested in further exploration. The first one is a variation of the recurrence from the main program using the LongFloat Library. The program expects z on level 2 and an estimate (such as the result of the main program) on level 1. It is rather slow, about 22 seconds on my 50g but reasonable on emulators. The value of c assumes the default value of DIGITS which is 50. Removing the R\<\>F at the end will leave a LongFloat number (type 27) on the stack. Code:
The last program implements the Halley recurrence shown at the Wikipedia link and in the Lóczi paper. The reason for including it is that the recurrence used in the main program works poorly for branches other than branch 0 and the realvalued part of branch 1. It tends to "drag" the values back to branch 0 even if given an accurate estimate for the value in another branch. The Halley recurrence converges rapidly even if the estimate is not very close, although it can be fooled if the estimate is vague. As above, the program expects z on level 2 and an estimate (usersupplied) on level 1. For HP28 and 48, change PICK3 to 3 PICK. Code:


« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)