200_success's solution cannot be much improved on, except by computing the handful of primes up to 1000 (i.e. sqrt(1,000,000)) up front to reduce the number of trial divisions that need to be performed.
I find it surprising though that this solution should TLE. On my aging Lynnfield a C# rendering of it (with pre-sieved primes) takes 4 ms for 20,000 random numbers in the target range, and 10 ms for 20,000 random primes from that range...
In a compiled language like C# there are additional optimisations that could be applied, like treating the prime 2 separately (elementary bit ops instead of division), or computing the modular inverses of above-mentioned small primes so that a divisibility test reduces to a multiplication followed by a comparison. However, in an interpreted language like Python there is little to be gained from that, because the interpreter overhead dwarfs the cost of these elementary operations.
Be that as it may, the high number of queries makes up-front computation of a phi table a viable alternative here. A suitable algorithm is discussed in Linear time Euler's Totient Function calculation. It uses the multiplicative properties of phi() to compute all values at little more than one division per result, by enumerating the target range in the form of least factor decompositions. It is worth keeping in mind because it can easily be adapted to handle other multiplicative functions.
Since my code in that topic is an early draft at early levels of understanding, here is a more mature rendering that I used my SPOJ submission (0.03 s):
static uint compute_phi_table (uint limit = 1000000)
var phi = new uint[limit + 1];
phi = 1;
var small_primes = new System.Collections.Generic.List<uint>();
// For each number in the target range, the outer loop has to supply what
// remains after dividing by the smallest factor; hence we have to go up to N/2.
uint n = 2, small_prime_limit = (uint)Math.Sqrt(limit);
for (uint e = limit / 2; n <= e; ++n)
uint phi_n = phi[n];
if (phi_n == 0)
phi_n = phi[n] = n - 1;
if (n <= small_prime_limit)
foreach (uint current_prime in small_primes)
uint nth_multiple = n * current_prime;
if (nth_multiple > limit)
if (n % current_prime != 0) // current_prime is not a factor of n
phi[nth_multiple] = phi_n * (current_prime - 1);
else // current_prime is a factor of n (the smallest, to be precise)
phi[nth_multiple] = phi_n * current_prime;
// process the primes in the rest of the range up to N (stepping mod 6)
n += 1 - (n & 1);
n += n % 3 == 0 ? 2u : 0u;
for (uint d = (3 - n % 3) * 2; n <= limit; n += d, d ^= 6)
if (phi[n] == 0)
phi[n] = n - 1;
The last step - the loop that replaces null entries with
n - 1 - can be skipped for SPOJ ETF, if null values retrieved from the table are treated accordingly.
Note, however, that this solution is slower than the earlier one. On my aged PC it takes 17 ms for computing the table up to 10^6, whereas factorisation with pre-computed primes can answer 20,000 random point queries in 4 ms - a mere fraction.
In other words, it seems that 200_success's solution combined with pre-computing a list of primes up to 1,000 is already the best that can be done in Python...
The table-based solution isn't very competitive in this particular instance (especially in Python, with its interpreter overhead), but in similar challenges it might well be. So it is worth keeping it in mind, especially as other multiplicative functions can be handled in the same fashion.