search

Patashu / break_eternity.js

A Javascript numerical library to represent numbers as large as 10^^1e308 and as small as 10^-10^^1e308. Sequel to break_infinity.js, designed for incremental games.
MIT License
120 stars 43 forks source link

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

Open James103 opened 5 months ago

James103 commented 5 months ago

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 ```javascript 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 commented 1 month ago

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 ```js 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: https://github.com/Patashu/break_eternity.js/blob/96901974c175cb28f66c7164a5a205cdda783872/src/index.ts#L346 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.