This is a late reply after leaving comments. I was looking back at my 'stackoverflow snippets' directory, and realised I had implemented a solution for this question (67689345), and forgot about it...
The key is to find an analytic expression for the gcd sum, which requires the factorization of (n). see: Pillai's function, or: The gcd-sum Function
Tabulating all prime factors for (n <= 1000000) at runtime is so fast that it's almost negligible to the total running time. However, I used a utility of mine to generate a table of all primes (p) <= (1000). It actually finds all primes up to (1024), but that's fine. sp_factor yields a prime factorization with multiplicity if a prime factor occurs more than once, so any repeated factors must be handled in the gcdsum function. These are required for the prime power exponents in any case.
Note: each (n) is factored to find the summation of gcd terms. We exploit the transform of the index of summation to yield:
Note, that the summation for (n <= 1000000) may exceed 32-bits. So a 64-bit accumulator is used. On an old laptop I'm using just now - a 2Ghz dual-core i7 - this takes less than (0.25) seconds, which should be well within the contest requirements.
The results for: n = {1, .., 5} match the contest values: {0, 1, 3, 7, 11}, while: n = 1000000 yields: 4071628673912
/******************************************************************************/
#include <inttypes.h>
/******************************************************************************/
static const uint16_t sp_lut[] =
{
0x0002, 0x0003, 0x0005, 0x0007, 0x000b, 0x000d, 0x0011, 0x0013,
0x0017, 0x001d, 0x001f, 0x0025, 0x0029, 0x002b, 0x002f, 0x0035,
0x003b, 0x003d, 0x0043, 0x0047, 0x0049, 0x004f, 0x0053, 0x0059,
0x0061, 0x0065, 0x0067, 0x006b, 0x006d, 0x0071, 0x007f, 0x0083,
0x0089, 0x008b, 0x0095, 0x0097, 0x009d, 0x00a3, 0x00a7, 0x00ad,
0x00b3, 0x00b5, 0x00bf, 0x00c1, 0x00c5, 0x00c7, 0x00d3, 0x00df,
0x00e3, 0x00e5, 0x00e9, 0x00ef, 0x00f1, 0x00fb, 0x0101, 0x0107,
0x010d, 0x010f, 0x0115, 0x0119, 0x011b, 0x0125, 0x0133, 0x0137,
0x0139, 0x013d, 0x014b, 0x0151, 0x015b, 0x015d, 0x0161, 0x0167,
0x016f, 0x0175, 0x017b, 0x017f, 0x0185, 0x018d, 0x0191, 0x0199,
0x01a3, 0x01a5, 0x01af, 0x01b1, 0x01b7, 0x01bb, 0x01c1, 0x01c9,
0x01cd, 0x01cf, 0x01d3, 0x01df, 0x01e7, 0x01eb, 0x01f3, 0x01f7,
0x01fd, 0x0209, 0x020b, 0x021d, 0x0223, 0x022d, 0x0233, 0x0239,
0x023b, 0x0241, 0x024b, 0x0251, 0x0257, 0x0259, 0x025f, 0x0265,
0x0269, 0x026b, 0x0277, 0x0281, 0x0283, 0x0287, 0x028d, 0x0293,
0x0295, 0x02a1, 0x02a5, 0x02ab, 0x02b3, 0x02bd, 0x02c5, 0x02cf,
0x02d7, 0x02dd, 0x02e3, 0x02e7, 0x02ef, 0x02f5, 0x02f9, 0x0301,
0x0305, 0x0313, 0x031d, 0x0329, 0x032b, 0x0335, 0x0337, 0x033b,
0x033d, 0x0347, 0x0355, 0x0359, 0x035b, 0x035f, 0x036d, 0x0371,
0x0373, 0x0377, 0x038b, 0x038f, 0x0397, 0x03a1, 0x03a9, 0x03ad,
0x03b3, 0x03b9, 0x03c7, 0x03cb, 0x03d1, 0x03d7, 0x03df, 0x03e5,
0x03f1, 0x03f5, 0x03fb, 0x03fd, 0x0000
};
/* 172 primes < 2^10 factor 83.87% of all odd integers. */
/******************************************************************************/
static inline unsigned int sp_factor (uint32_t p[], uint32_t n)
{
uint32_t sp = sp_lut[0], q;
unsigned int np = 1;
/* assert(n > 1 && n < (UINT32_C(1) << (20))); */
for (unsigned int i = 1; (q = n / sp) >= sp; )
{
if (q * sp == n)
np++, n = q, *p++ = sp;
else if ((sp = sp_lut[i++]) == 0) /* EOT entry: */
break;
}
*p++ = n;
return np; /* the number of prime factors (with multiplicity). */
}
/******************************************************************************/
static uint64_t gcdsum (uint32_t n_max)
{
uint64_t sum = 0;
/* assert(n < (UINT32_C(1) << (20))); */
for (uint32_t n = 2; n <= n_max; n++)
{
uint32_t ptab[(20)], pn, gn = 2 * n - 1;
if ((pn = sp_factor(ptab, n)) != 1) /* composite: */
{
uint32_t etab[(20)], un, p, i;
etab[0] = 1, un = 1;
for (p = ptab[0], i = 1; i < pn; i++)
{
uint32_t pi = ptab[i];
if (pi == p)
etab[un - 1]++;
else
{
ptab[un] = (p = pi), etab[un] = 1;
un++;
}
}
for (gn = 1, i = 0; i < un; i++)
{
uint32_t pi = ptab[i], ei = etab[i];
gn *= ((ei + 1) * pi - ei);
while (--ei) gn *= pi; /* pi ^ (ei - 1) */
}
}
sum += gn - n; /* exclude (n, n) term from summation. */
}
return sum;
}
/******************************************************************************/
