`Decimal.lambertw` occasionally spikes with inputs around the ee31 range

[James103]

When running the following:

[
Decimal.lambertw("ee31.666733435898778").toString(),
Decimal.lambertw("ee31.666733435898788").toString(),
Decimal.lambertw("ee31.666733435898798").toString()
]

the output is

[
"1.0689296521861698e32",
"2.575091926685474e3133",
"1.0689296521862222e32"
]

The second value is many orders of magnitude higher that of the first and third values, effectively a large spike.
The correct answer for the second expression is approximately midway between the first value and the third value.

Some more test cases, all of which return incorrect values

W(ee31.965173139544905) == 3.590892194141079e6229
W(ee31.752610949621566) == 3.296534240252085e3818
W(ee31.913848146747807) == 1.5268311171822413e5535
W(ee31.640295850532457) == 2.2637793318601247e2948
W(ee31.764948764454818) == 3.5762246079696114e3928
W(ee31.652926931263902) == 2.3139757541974464e3035
W(ee31.929958343490664) == 2.332630066079301e5744
W(ee31.98205373630235) == 2.878867493568092e6476
W(ee31.98283664180519) == 1.3960688507821357e6488
W(ee31.787391560044053) == 8.01963757821285e4136
W(ee31.823874021678865) == 2.7187091896742213e4499
W(ee31.825232704314676) == 3.4104034609250986e4513
W(ee31.863878434083848) == 3.8176698551518697e4933
W(ee31.918276882779058) == 8.276372977284206e5591
W(ee31.669172925011136) == 1.151060360788823e3151
W(ee31.909157710719125) == 5.305882880703211e5475
W(ee31.954588884894793) == 3.7058668964725516e6079
W(ee31.690265691504234) == 7.565092881039769e3307
W(ee31.696486703711212) == 4.0107213474029315e3355
W(ee31.907517771777723) == 1.2214311038222574e5455

Code to generate the above test cases

let count = 0;
while (count < 100) {
	const n = Math.random() * 3 + 30
	const Wn = Decimal.lambertw("ee" + n)
	if (Wn.gt(1e40)) {
		console.log("W(ee"+n+") == "+Wn.toString())
		count++
	}
	else
		count+=(1/1024)
}

[James103]

I've added some debug lines to the d_lambertw function to try to trace down the cause of this issue.

The cause appears to be that ew (derived from w.neg().exp() where w is Decimal.ln(z) assuming the principal branch) is not completely cancelled out by z when it normally is at these large numbers, causing z.mul(ew) to be some large number (≈ ee17) instead of 1, which then causes later calculations during that same iteration to break and the final result to be many thousands of orders of magnitude higher than what it should be.

Although z.mul(Decimal.ln(z).neg().exp()) $= z e^{-\ln(z)}$ simplifies to 1 for all z, it can in practice be very far away from 1 (like ee17 units off) due to incomplete cancellation of the number z and its reciprocal calculated with Decimal.ln(z).neg().exp().

d_lambertw output when running code in OP using the principal branch

17:40:17.881 lambertw trace: z = ee31.666733435898777
17:40:17.881   Iteration 0:
17:40:17.881     w = 1.0689296521861698e32, w.neg() = -1.0689296521861698e32, ew = w.neg().exp() = ee-31.666733435898777
17:40:17.881     z.mul(ew) = 1, wewz = w.sub(z.mul(ew)) = 1.0689296521861698e32
17:40:17.881     Decimal.mul(2, w).add(2) = 2.13785930437234e32
17:40:17.882     w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)) = 5.344648260930848e31
17:40:17.882     w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))) = 5.344648260930848e31
17:40:17.882     wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)))) = 2.0000000000000004
17:40:17.882     wn = w.sub(wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))))) = 1.0689296521861698e32
17:40:17.882   Returning: 1.0689296521861698e32
17:40:17.882 lambertw trace: z = ee31.666733435898788
17:40:17.882   Iteration 0:
17:40:17.882     w = 1.0689296521861873e32, w.neg() = -1.0689296521861873e32, ew = w.neg().exp() = ee-31.666733435898784
17:40:17.882     z.mul(ew) = ee17.581375385438758, wewz = w.sub(z.mul(ew)) = -ee17.581375385438758
17:40:17.883     Decimal.mul(2, w).add(2) = 2.137859304372375e32
17:40:17.883     w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)) = -ee17.581375385438754
17:40:17.883     w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))) = ee17.581375385438754
17:40:17.883     wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)))) = -2.575091926685474e3133
17:40:17.883     wn = w.sub(wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))))) = 2.575091926685474e3133
17:40:17.883   Iteration 1:
17:40:17.883     w = 2.575091926685474e3133, w.neg() = -2.575091926685474e3133, ew = w.neg().exp() = ee-3133.0485770485766
17:40:17.883     z.mul(ew) = ee-3133.0485770485766, wewz = w.sub(z.mul(ew)) = 2.575091926685474e3133
17:40:17.883     Decimal.mul(2, w).add(2) = 5.15018385337196e3133
17:40:17.883     w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)) = 1.287545963342484e3133
17:40:17.883     w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))) = 1.287545963342484e3133
17:40:17.883     wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)))) = 2.000000000000393
17:40:17.883     wn = w.sub(wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))))) = 2.575091926685474e3133
17:40:17.884   Returning: 2.575091926685474e3133
17:40:17.884 lambertw trace: z = ee31.6667334358988
17:40:17.884   Iteration 0:
17:40:17.884     w = 1.0689296521862222e32, w.neg() = -1.0689296521862222e32, ew = w.neg().exp() = ee-31.6667334358988
17:40:17.884     z.mul(ew) = 1, wewz = w.sub(z.mul(ew)) = 1.0689296521862222e32
17:40:17.884     Decimal.mul(2, w).add(2) = 2.1378593043724448e32
17:40:17.884     w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)) = 5.344648260931111e31
17:40:17.884     w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))) = 5.344648260931111e31
17:40:17.884     wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2)))) = 2.0000000000000004
17:40:17.884     wn = w.sub(wewz.div(w.add(1).sub(w.add(2).mul(wewz).div(Decimal.mul(2, w).add(2))))) = 1.0689296521862222e32
17:40:17.884   Returning: 1.0689296521862222e32

---

EDIT: Found it! Change this line in the d_lambertw function:

break_eternity.js/src/index.ts

Line 346 in 9690197

wewz = w.sub(z.mul(ew));

By replacing z.mul(ew) in the .sub() call of wewz = w.sub(z.mul(ew)) with 1, which is what it simplifies to after undoing the substitutions. Note: Only tested with certain numbers in the principal branch; needs testing for the real non-principal branch and over a wider range of numbers.
EDIT 2: The above fix does not appear to converge for smaller numbers (layer 1). A better solution may be needed.