HP49-50G: MLSV OK for -1E-495, not OK for -1E-496, ROOT OK for both

[Gil]

I have a program and may get a complex number, for which I look for a solution.

The function is x*EXP(x).

I want to solve x*EXP(x)=-1E-495.

With ROOT
'x*EXP(x)+1E-495' (in stack level 3)
'x' (in stack level 2)
-1000 (an approximate guess in stack 1)
ROOT
and the output: -1146.82437196

With MLSV
['x*EXP(x)+1E-495'] (in stack level 3)
['x'] (in stack level 2)
[-1000] (an approximate guess in stack 1)
MLSV
and the outputs:
[ 'X*EXP(X)+1.E-495' ]
[ 'X' ]
[ -1146.82437276 ]

Fine, the final result found with MLSV is almost the same as with ROOT.

Now repeat the two exercices for x*EXP(x)=-1E-496

With ROOT
'x*EXP(x)+1E-496' (in stack level 3)
'x' (in stack level 2)
-1000 (an approximate guess in stack 1)
ROOT
and the output: --1148.9899614

With MLSV
['x*EXP(x)+1E-496'] (in stack level 3)
['x'] (in stack level 2)
[-1000] (an approximate guess, even --1148.98 or -1148.9899614, in stack 1)
MLSV
and then an endless "search" will occur, with no output.

Why?
And is there a way, in this case, to get the wished value of about -1148.99 with the MLSV (as already mentioned, I am not supposed to know, beforehand, that x is a real and that I can use instead ROOT command without the brackets) ?

[Albert Chan]

(01-25-2024 11:08 AM)Gil Wrote:  I have a program and may get a complex number, for which I look for a solution.

The function is x*EXP(x).

Newton's method is not good for solving x*e^x = a, if |x| is big (*)

Instead, we may solve for the tempered version, by taking ln

(x*e^x = a) --> (ln(x) + x = ln(a)) --> (x + ln(x/a) = 0)

HP71B
> a = -1e-496
> x = ln(-a) ! branch -1
> for i=1 to 3 @ x = (1-ln(x/a))/(1+1/x) @ x @ next i
-1149.12898462
-1149.12896563
-1149.12896564

(*) OTTH, if we already have an extremely good guess, good for finishing touch.

[Gil]

I don't think that the built-in MSLV solver (multi solver for several variables and equation ) of the HP50G calculator uses the Newton's algorithm.

For complex numbers, we must use it instead of the ROOT command reserved for a real unknown.

What is strange is that both methods find the appropriate solution for negative numbers very close to zero,
ie up to -1E-495.

For negative numbers still nearer to zero,
for instance -1E-496, only the command ROOT gives an answer.

I wished to know the reason — in a word, what is so particular with -1E-496 in comparison to -1E-495.

[Gil]

I check your result with Wolfram.

And, as usual, your result is correct.

Though, the ROOT solver seems to find all solutions for negative values close to a — up to -1E-495 —, the values it finds for a= -1E-496,-1E-497, -1E-498 and -1E-496 are all the same (-1148.98996141)... and incorrect!

Thanks, Albert, for your insight : I will soon correct my program according to your modified formulae when taking the log.

But strange anyway that the results are correct up to -1E-494, almost correct for -1E-495 and wrong for negative numbers still nearer to zero.

[Gil]

I used your method, Albert.

It gives correct digits for -1E-492, -1-E493, -1E-494, -1E-495 & - -1E-496 — and better results than normal method regarding last digits for - 1E-492, -1-E493, -1E-494 and -1E-495.

But your method on my HP50G gives

{ { for -1.E-497: -1151.2925465 } { for -1.E-498: -1151.2925465 } {for -1.E-499: -1151.2925465 } }.

Three same results, and incorrect.

Could you check on your HP calculator what you get?

x/a gives the max of the calculator capacity: 9.99999999999E499,
and we will take always the Ln of that maximum value, to get always the same result... 1151.2925465.

We should try and alternative solution.

Just changing your formula Ln (x/a) by Ln (x) - Ln(a):
'(1-(LN(x)-LN(a)))/(1+1/x)'

[Albert Chan]

(01-28-2024 12:59 AM)Gil Wrote:  Though, the ROOT solver seems to find all solutions for negative values close to a — up to -1E-495 —, the values it finds for a= -1E-496,-1E-497, -1E-498 and -1E-496 are all the same (-1148.98996141)... and incorrect!

x*e^x = a

min(e^x) = 1e-499      // if e^x = 0 then x*e^x = -inf * 0 = nan
min(x) = ln(1e-499) = -1148.9899614

(01-28-2024 02:01 AM)Gil Wrote:  But your method on my HP50G gives

{ { for -1.E-497: -1151.2925465 } { for -1.E-498: -1151.2925465 } {for -1.E-499: -1151.2925465 } }.

Three same results, and incorrect.

x + ln(x/a) = 0

max(x/a) = 9.99999999999E499
min(x) = -ln(9.99999999999E499) = -1151.2925465

Quote:We should try and alternative solution.

Just changing your formula Ln (x/a) by Ln (x) - Ln(a):
'(1-(LN(x)-LN(a)))/(1+1/x)'

Since x and a are negative real, try (ln(-x) - ln(-a))
Or, with (1+1/x) ≈ 1, just iterate x = ln(-a) - ln(-x)

Cas> c := ln(10) * -100000.; /* = ln(-a) */
−230258.509299

Cas> c - ln(-Ans)
−230270.856257
−230270.856311
−230270.856311 /* = W_-1(-1e-100000) */

[Gil]

Nicely seen that ln(x/a) = ln(-x/-a) =ln(-x) - ln(-a).

So we get rid of complexes, instead of having a final pseudo complex =a+0×i that would oblige to take only the real part.

With ln(-x) & ln(-a), we remain in R.

Much better!

Thanks.

[Gil]

Hi, Albert.

One more "challenge" for you.

How would it possible to get, on your calculator, about 10 or more significant correct digits for W0(x) and W-1(x), 'x being a real number very near of -1/e, but > -1/e'?

Thanks in advance for your solutions.

[Albert Chan]

(01-28-2024 01:51 PM)Gil Wrote:  How would it possible to get, on your calculator, about 10 or more significant correct digits for W0(x) and W-1(x),
'x being a real number very near of -1/e, but > -1/e'?

It is a solved problem, incorporated in both Lua I.W(z,k) and Python W(z,k)

I also made HP-71B version (real z, branch 0,-1) (post#18 solved with y=e^x, post#19 wth accurate ln(1+x)-x)
John Keith had kindly translated post#19 HP71B version to RPL (post #28)

All versions relied on precise 1/e = float(1/e) + eps, both constants stored in advance.

[Gil]

As you have noted, I am no math guy.

I tried this "journey", in which frankly I got somewhat lost.

I should start the whole program process from new basis.

But I am afraid I haven't got the courage for it... nor the capacity.

Sometimes, I even doesn't know why we (you!) took a guess value, and not another one, despite your long explanations.

I am like a parrot, with nice colours — or a beautiful toys in the hands, the great EMU48, but with half a brain.

Of course, tomorrow is a new week.

[Gil]

For x near -1/e, I tried to copy your code, Albert

« a NEG 1 60
START 'y' STO y 'y*LN(y)-a' —>NUM 'y-1/e' —>NUM e * LNP1 —>NUM / -
NEXT LN
»

Code:

\<< a NEG 1 60
  START 'y' STO y 'y*LN(y)-a' \->NUM 'y-1/e' \->NUM e * LNP1 \->NUM / -
  NEXT LN
\>>

May I ask what your program gives for
W-1(a=-.3678794411).

My result with the above code on my HP50G is
-1.000019 63105

But Wolfram differs:
-1.000019 708014416801577108775088535244577659214984302969737

[Albert Chan]

(01-28-2024 11:55 PM)Gil Wrote:  May I ask what your program gives for
W-1(a=-.3678794411).

Keep in mind Lua code has decimal to binary conversion error.
a = -0x1.78b5636194be9p-2 = -0.36787944109999998199...

lua> a = -.3678794411
lua> I.W(a, -1)
(-1.000019708016901-0*I)

Compare apples to apples, with hexfloat of a, it agreed with WA (error = 1 ULP)

---

HP71B code matched WA number exactly.

>run
a, k? -.3678794411, -1
e^W, W = .36787219107      -1.00001970801

BTW, my code does not do (y-1/e).
Instead, it does (y-r-eps), with *doubled* precision of float(1/e)