My code

… was written in C++11 and can be compiled with G++. You can download it, too. Or just jump to my GitHub repository.

       
#include <iostream>
       
#include <vector>
       
#include <unordered_map>
       
 
       
// I had overflows with 64 bit arithmetic, so I switched to GCC's 128 bit extension

       
typedef __int128 Number;
       
 
       
// compute the exponent of prime in the factorization of f!

       
unsigned long long countFactors(Number f, unsigned int prime)
       
{
       
  // Legendre's formula, https://en.wikipedia.org/wiki/Legendre%27s_formula
       
  unsigned long long result = 0;
       
  Number power  = prime;
       
  while (power <= f)
       
  {
       
    result += f / power;
       
    power  *= prime;
       
  }
       
 
       
  return result;
       
}
       
 
       
// variables of the problem statement:
       
// p = prime      = 61
       
// q = iterations = 10^7
       
// moduloExponent = 10 (because of 61^10)

       
unsigned long long solve(unsigned int prime, unsigned int iterations, unsigned int exponent)
       
{
       
  // compute prime^moduloExponent => 61^10
       
  unsigned long long modulo = 1;
       
  for (unsigned int i = 1; i <= exponent; i++)
       
    modulo *= prime;
       
 
       
  // will be 61^i => 61^0, 61^1, 61^2, ..., 61^10 and keep it at 61^10
       
  Number maxPower = 1;
       
 
       
  // seed of pseudo-random number generator
       
  unsigned long long s = 290797;
       
 
       
  unsigned long long result = 0;
       
  for (unsigned int i = 0; i <= iterations; i++)
       
  {
       
    // T(i, 61) * 61^i
       
    auto t = s % prime;
       
    // => but 61^i limited to 61^10
       
    auto product = t * maxPower;
       
 
       
    // I observed only 71 different products (in 1000000 iterations !), let's cache them
       
    static std::unordered_map<Number, unsigned long long> cache;
       
    auto lookup = cache.find(product);
       
    if (lookup == cache.end())
       
    {
       
      // new, unknown parameters, must call countFactors()
       
      auto current   = countFactors(product, prime);
       
      // cache the value
       
      cache[product] = current;
       
      // add it to the result, too
       
      result        += current;
       
    }
       
    else
       
    {
       
      // add looked up value
       
      result += lookup->second;
       
    }
       
 
       
    // keep under 61^10
       
    result %= modulo;
       
 
       
    // 61^i => 61^(i+1), but limit to 61^10
       
    if (maxPower < modulo)
       
      maxPower *= prime;
       
 
       
    // next random number
       
    s *= s;
       
    s %= 50515093;
       
  }
       
 
       
  return result;
       
}
       
 
       
int main()
       
{
       
  unsigned int prime      = 61;
       
  unsigned int iterations = 10000000;
       
  unsigned int exponent   = 10;
       
  std::cin >> prime >> iterations >> exponent;
       
  std::cout << solve(prime, iterations, exponent) << std::endl;
       
  return 0;
       
}

This solution contains 13 empty lines, 20 comments and 3 preprocessor commands.