Calculate the Leyland number index of a given Leyland (x,y) pair

Question

A Leyland number L(x,y) is x^y + y^x where x>=y>1. OEIS A076980 prepends L(2,1) as its first term. The first ten Leyland number (x,y) pairs ordered by magnitude are (2,1), (2,2), (3,2), (4,2), (3,3), (5,2), (6,2), (4,3), (7,2), (8,2). We wish to efficiently calculate the A076980 index of a given (x,y) pair, understanding that the procedure may be somewhat onerous for large L(x,y). For example, the recently discovered L(101311,90816) probable prime has, by my own admittedly inefficient Mathematica procedure, an index of 6137957218.

It fails for the initial (2,1), but here is the code that I've been using since 2015:

L[{x_, y_}] := x^y + y^x;
LN[{x_, y_}] := (
  q = L[{x, y}]; 
  nn = q*.001; 
  n = 1.0;
  count = 2; 
  While[
    n = n + 1.0;
    num = 2*n^n;
    num < nn || 2*Round[n]^Round[n] < q,
    count++;
    k = n; 
    While[
      k = k + 1.0; 
      num = L[{n, k}]; 
      num < nn || L[{Round[n], Round[k]}] < q,
      count++
    ]
  ];
  count
  )

Here's a neat application. Assuming the conjecture that (for d > 11) the smallest d-digit Leyland number is L(d-1,10), then the number of d-digit Leyland numbers is LN[{d,10}] - LN[{d-1,10}]. For d = 1000000, this works out to 39542.

Answer 1

@Hans: Sorry for causing a nuisance yesterday with my poor code. I deleted it to prevent clutter in your thread. Your code was much, much faster even when I corrected it.
Here is a new suggestion: Afraid your code was a bit difficult for me to follow so I re-wrote it with the intention of more easily compiling it. Here is your code for $L(1000,1000)$ and the timing at 4.5 GHz:

L[{x_, y_}] := x^y + y^x; 
LN[{x_, y_}] := (q = L[{x, y}]; nn = q*.001; n = 1.0; count = 2; 
  While[n = n + 1.0; num = 2*n^n; num < nn || 2*Round[n]^Round[n] < q,
    count++; k = n; 
   While[k = k + 1.0; num = L[{n, k}]; 
    num < nn || L[{Round[n], Round[k]}] < q, count++]]; count);

AbsoluteTiming[
 index = LN[{1000, 1000}];
 ]
index

Out[18]= {8.91903, Null}

Out[19]= 725836

And here is your code with some very minor improvements but modularized:

myLN[x_, y_] := Module[{q, nn, n, count, nR, num, k, kR},
   q = x^y + y^x;
   nn = 0.001 q;
   n = 1.;
   count = 2;
   While[
    n = n + 1.0;
    nR = Round[n];
    num = 2 n^n;
    num < nn || (2 nR^nR < q),
    count++;
    k = n;
    While[
     k = k + 1.0;
     kR = Round[k];
     num = n^k + k^n;
     num < nn || (nR^kR + kR^nR) < q,
     count++;
     ];
    ];
   count
   ];

AbsoluteTiming[
 myLN[1000, 1000]
 ]
(* {8.84537, 725836} *)

A modicum of improvements for sure but in a form I think is easier to compile but not sure that's doable. I'm a bit weak on compiling code. Maybe someone else here can do so. I'll try a bit later as well. $L(2500,2500)$ results are $57.6$ and $55.6$ respectively.