Strange NSolve failure [duplicate]

Question

In[22]:= NSolve[-LogIntegral[2] + LogIntegral[x] == 119002, x]

During evaluation of In[22]:= NSolve::ifun: Inverse functions are being used by NSolve,
so some solutions may not be found; use Reduce for complete solution information.

Out[22]= {{x -> 1.56751*10^6}}

There are some natural numbers where this fails:

In[20]:= NSolve[-LogIntegral[2] + LogIntegral[x] == 119003, x]

During evaluation of In[20]:= NSolve::nsmet: This system cannot be solved
with the methods available to NSolve.

Out[20]= NSolve[-LogIntegral[2] + LogIntegral[x] == 119003, x]

What is so special about 119003?

In[23]:= NSolve[-LogIntegral[2] + LogIntegral[x] == 119004, x]

During evaluation of In[23]:= NSolve::ifun: Inverse functions are being used by NSolve,
so some solutions may not be found; use Reduce for complete solution information.

Out[23]= {{x -> 1.56754*10^6}}

Some other such integers are 128005, 133003, 137004. What could be the reason?

In case it is version dependent, mine is 11.0.1.0 (on Windows 10)

Answer 1

As Bob says, NSolve has an issue with using machine numbers, so why not use Solve instead? In the following, I also use the InverseFunctions option to avoid the Solve::ifun message, since I do want Solve to use inverse functions, and I include N to convert the Root object output of Solve into a real number:

Map[
    N @ Solve[-LogIntegral[2]+LogIntegral[x]==#, x, InverseFunctions->True]&,
    {119003, 128005, 133003, 137004}
]

{{{x -> 1.56753*10^6}}, {{x -> 1.6963*10^6}}, {{x -> 1.76809*10^6}}, {{x -> 1.82572*10^6}}}

Answer 2

You need to use arbitrary-precision rather than machine precision.

With machine precision:

Quiet[NSolve[-LogIntegral[2] + LogIntegral[x] == #, x] & /@ {119003, 128005, 
   133003, 137004}]

(* {NSolve[-LogIntegral[2] + LogIntegral[x] == 119003, 
  x], {{x -> 1.6963*10^6 - 1.62507*10^-15 I}}, {{x -> 
    1.76809*10^6 + 4.29213*10^-15 I}}, 
 NSolve[-LogIntegral[2] + LogIntegral[x] == 137004, x]} *)

With arbitrary-precision

Quiet[NSolve[-LogIntegral[2] + LogIntegral[x] == #, x,
    WorkingPrecision -> 20] & /@ {119003, 128005, 133003, 137004}]

(* {{{x -> 1.5675261048618991866*10^6}}, {{x -> 
    1.6962994222144938141*10^6}}, {{x -> 1.7680947906149612601*10^6}}, {{x -> 
    1.8257153022466093236*10^6}}} *)