Inria has published a paper factoring 50-200 bits integers.
The four methods of choice are: SQUFOF, ECM, SIQS and CFRAC.

Within the paper a short pseudocode is given for each method, so you may be able to judge whether they will suit your computation restrictions.

From their experiments, SQUFOF is indeed the fastest up to $2^{60}$ range, where it is tied against SIQS. Note that these experiments are done on integers known to have few prime factors (Something like 2-4).

For your case, it appears that you are factoring random integers up to $2^{60}$.
Hence, except for ECM, it is more efficient to trial factor small primes up to a certain bound.
(See comments below on the case for ECM)

One additional point: ECM (specifically GMP-ECM, also from Inria) might be a good idea since the software is highly optimized. (With a research team coding it for a few years)
On the other hand, the arithmetic is done modulo $N$.
This might mean that there will be overflow during the multiplication steps.

If you choose ECM, I can elaborate a little more.
A common problem is the fact that ECM is probabilistic.
This can be resolved if you know which curves to choose to cover all the prime factors.
Inria has also published its work on this area, showing that 124 carefully chosen curves finds all prime factors up to $2^{32}$.
The numbers are quite suitable for your case.
Statistically speaking though, you find the factors with probability around 0.15 per curve, so usually 20 curves will do the job.
As mentioned earlier, using these parameters you can skip the trial factorization part.
During the search for bigger primes, you will also find the smaller ones in the process.

If you want to push the efficiency even higher, the optimal approach is to first do a filtering based on $P-1$ and $P+1$ factorization method, followed by the actual ECM.
This is a little tougher since you need to calibrate the parameters correctly.

EDIT: Did not notice that all numbers are of the form $x^3-y^3$ with $ 0 < y < x < 50000 $.
There are about $50000^2/2=2^{36.86}$ possibilities.
Out of all these, a portion will be primes.
A portion will be smooth with factors $< B$ for some $B$ that can be easily trial divided.
The remain cases are composites with primes $> B$.

Suppose you have $n$ remaining entries.
Sort them in the following format, by $N_i$:
$N_i L_i$, where $N_i=x^3-y^3$ is the entry's value and $L_i$ is the offset to find the factorization.
At offset $L_i$, store the distinct prime factors (may choose not to store multiplicity by doing trial division).

Each time you get a number:

1. Do primality test
2. Do trial division
3. Do another primality test
4. If residue remains, do a binary search on the factorization

I think this should be able to handle your 4 million factorization.
Clearly there are more ways to tweak the method, but the idea here is you can do a one time computation and store the difficult computations as a look up table.

Inria has published a paper factoring 50-200 bits integers.
The four methods of choice are: SQUFOF, ECM, SIQS and CFRAC.

Within the paper a short pseudocode is given for each method, so you may be able to judge whether they will suit your computation restrictions.

From their experiments, SQUFOF is indeed the fastest up to $2^{60}$ range, where it is tied against SIQS. Note that these experiments are done on integers known to have few prime factors (Something like 2-4).

For your case, it appears that you are factoring random integers up to $2^{60}$.
Hence, except for ECM, it is more efficient to trial factor small primes up to a certain bound.
(See comments below on the case for ECM)

One additional point: ECM (specifically GMP-ECM, also from Inria) might be a good idea since the software is highly optimized. (With a research team coding it for a few years)
On the other hand, the arithmetic is done modulo $N$.
This might mean that there will be overflow during the multiplication steps.

If you choose ECM, I can elaborate a little more.
A common problem is the fact that ECM is probabilistic.
This can be resolved if you know which curves to choose to cover all the prime factors.
Inria has also published its work on this area, showing that 124 carefully chosen curves finds all prime factors up to $2^{32}$.
The numbers are quite suitable for your case.
Statistically speaking though, you find the factors with probability around 0.15 per curve, so usually 20 curves will do the job.
As mentioned earlier, using these parameters you can skip the trial factorization part.
During the search for bigger primes, you will also find the smaller ones in the process.

If you want to push the efficiency even higher, the optimal approach is to first do a filtering based on $P-1$ and $P+1$ factorization method, followed by the actual ECM.
This is a little tougher since you need to calibrate the parameters correctly.

EDIT: Did not notice that all numbers are of the form $x^3-y^3$ with $ 0 < y < x < 50000 $.
There are about $50000^2/2=2^{36.86}$ possibilities.
Out of all these, a portion will be primes.
A portion will be smooth with factors $< B$ for some $B$ that can be easily trial divided.
The remain cases are composites with primes $> B$.

Suppose you have $n$ remaining entries.
Sort them in the following format, by $N_i$:
$N_i L_i$, where $N_i$ is the entry's value, each a composite of factors $> B$, and $L_i$ is the offset to find the factorization.
At offset $L_i$, store the distinct prime factors (may choose not to store multiplicity by doing trial division).

Each time you get a number:

1. Do primality test
2. Do trial division
3. Do another primality test
4. If residue remains, do a binary search on the factorization

I think this should be able to handle your 4 million factorization.
Clearly there are more ways to tweak the method, but the idea here is you can do a one time computation and store the difficult computations as a look up table.

Inria has published a paper factoring 50-200 bits integers.
The four methods of choice are: SQUFOF, ECM, SIQS and CFRAC.

Within the paper a short pseudocode is given for each method, so you may be able to judge whether they will suit your computation restrictions.

From their experiments, SQUFOF is indeed the fastest up to $2^{60}$ range, where it is tied against SIQS. Note that these experiments are done on integers known to have few prime factors (Something like 2-4).

For your case, it appears that you are factoring random integers up to $2^{60}$.
Hence, except for ECM, it is more efficient to trial factor small primes up to a certain bound.
(See comments below on the case for ECM)

One additional point: ECM (specifically GMP-ECM, also from Inria) might be a good idea since the software is highly optimized. (With a research team coding it for a few years)
On the other hand, the arithmetic is done modulo $N$.
This might mean that there will be overflow during the multiplication steps.

If you choose ECM, I can elaborate a little more.
A common problem is the fact that ECM is probabilistic.
This can be resolved if you know which curves to choose to cover all the prime factors.
Inria has also published its work on this area, showing that 124 carefully chosen curves finds all prime factors up to $2^{32}$.
The numbers are quite suitable for your case.
Statistically speaking though, you find the factors with probability around 0.15 per curve, so usually 20 curves will do the job.
As mentioned earlier, using these parameters you can skip the trial factorization part.
During the search for bigger primes, you will also find the smaller ones in the process.

If you want to push the efficiency even higher, the optimal approach is to first do a filtering based on $P-1$ and $P+1$ factorization method, followed by the actual ECM.
This is a little tougher since you need to calibrate the parameters correctly.

Inria has published a paper factoring 50-200 bits integers.
The four methods of choice are: SQUFOF, ECM, SIQS and CFRAC.

Within the paper a short pseudocode is given for each method, so you may be able to judge whether they will suit your computation restrictions.

From their experiments, SQUFOF is indeed the fastest up to $2^{60}$ range, where it is tied against SIQS. Note that these experiments are done on integers known to have few prime factors (Something like 2-4).

For your case, it appears that you are factoring random integers up to $2^{60}$.
Hence, except for ECM, it is more efficient to trial factor small primes up to a certain bound.
(See comments below on the case for ECM)

One additional point: ECM (specifically GMP-ECM, also from Inria) might be a good idea since the software is highly optimized. (With a research team coding it for a few years)
On the other hand, the arithmetic is done modulo $N$.
This might mean that there will be overflow during the multiplication steps.

If you choose ECM, I can elaborate a little more.
A common problem is the fact that ECM is probabilistic.
This can be resolved if you know which curves to choose to cover all the prime factors.
Inria has also published its work on this area, showing that 124 carefully chosen curves finds all prime factors up to $2^{32}$.
The numbers are quite suitable for your case.
Statistically speaking though, you find the factors with probability around 0.15 per curve, so usually 20 curves will do the job.
As mentioned earlier, using these parameters you can skip the trial factorization part.
During the search for bigger primes, you will also find the smaller ones in the process.

If you want to push the efficiency even higher, the optimal approach is to first do a filtering based on $P-1$ and $P+1$ factorization method, followed by the actual ECM.
This is a little tougher since you need to calibrate the parameters correctly.

EDIT: Did not notice that all numbers are of the form $x^3-y^3$ with $ 0 < y < x < 50000 $.
There are about $50000^2/2=2^{36.86}$ possibilities.
Out of all these, a portion will be primes.
A portion will be smooth with factors $< B$ for some $B$ that can be easily trial divided.
The remain cases are composites with primes $> B$.

Suppose you have $n$ remaining entries.
Sort them in the following format, by $N_i$:
$N_i L_i$, where $N_i=x^3-y^3$ is the entry's value and $L_i$ is the offset to find the factorization.
At offset $L_i$, store the distinct prime factors (may choose not to store multiplicity by doing trial division).

Each time you get a number:

1. Do primality test
2. Do trial division
3. Do another primality test
4. If residue remains, do a binary search on the factorization

I think this should be able to handle your 4 million factorization.
Clearly there are more ways to tweak the method, but the idea here is you can do a one time computation and store the difficult computations as a look up table.