由于一些很蠢的原因，我写了一份完全用不着的用于计算 Lambert W function 的 C# 代码. 具体原理是很粗暴的牛顿法求解，但有几个特别处理的地方：
首先是估算. W0比较好处理，其实随便给个初始值就好，我这里选择了偏移后的 ex. 而 W-1 则比较特殊，首先在 x < -2 的部分变成了 concave 的，而且越往左斜率越小，所以一律从 x = -2 开始尝试；而在 -2 < x < -1 的区间内其实也是随便给个初始值就好，我这里选择了偏移后的 cos 函数.
其次是精度处理. 通过 log2x 做差可以求出二进制下的有效数位，然而当 y 值较大的时候其实并不能把完整的 52 位双精度浮点数的小数部分求完整，大概只能求到 47 位左右，所以这里我保守选择了 42 位有效数字.
代码：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 |
using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks; using SysMath = System.Math; namespace MartinCl2.Math { /// <summary> /// Lambert W function is the inverse function of x*e^x. /// It has two branches in [-1/e, +inf). /// </summary> class LambertW { /// <summary> /// Use Newton's method to find the numerical solution. /// </summary> /// <param name="branch"></param> /// <param name="input"></param> /// <returns></returns> public static double Evaluate(int branch, double input) { // Check argument range if (input < -1 / SysMath.E) { throw new ArgumentOutOfRangeException ("input", "input should be larger than -1/e."); } if (branch != 0 && branch != -1) { throw new ArgumentOutOfRangeException ("branch", "branch should be 0 or -1."); } if (input > 0 && branch == -1) { throw new ArgumentOutOfRangeException ("Only branch 0 has positive input."); } // Start with an estimation double result = Estimate(branch, input); // Iterate using Newton's method int loopsRemaining = 100; while (TestResult(result, input) && loopsRemaining > 0) { double x0 = result; double y0 = result * SysMath.Exp(result); double k = SysMath.Exp(x0) * (1 + x0); double y1 = input; double x1 = (y1 - y0) / k + x0; result = x1; loopsRemaining--; } if (loopsRemaining == 0) { throw new Exception("Cannot find solution."); } return result; } /// <summary> /// Give a rough estimation of Lambert W function /// </summary> /// <param name="branch"></param> /// <param name="input"></param> /// <returns></returns> private static double Estimate(int branch, double input) { if (branch == -1) { // return SysMath.Max(-2, -0.909091 + (8.97204 * // SysMath.Sqrt(0.00138945 - 0.0102667 * input // * input)) / input); if (input > -0.2706706) { // Cannot go further less otherwise Newton's method // won't work return -2; } else { return -1 - SysMath.Acos(-4.11727 * input - 0.514659); } } else if (branch == 0) { return -1 + SysMath.Log(1 / SysMath.E + SysMath.E + input); } else { throw new ArgumentException(); } } /// <summary> /// Test whether f(x)=x*e^x and y is close enough. /// </summary> /// <param name="x"></param> /// <param name="y"></param> /// <returns></returns> private static bool TestResult(double x, double y) { double fx = x * SysMath.Exp(x); double delta = fx - y; double deltaBits = SysMath.Log(SysMath.Abs(y / delta), 2); // Double-precision floating-point has 52 bits fraction. // We allow 10 bit tolarance return deltaBits < 42; } } } |
0 条评论。