Solution

Answer: 191541795

Let

$$\frac1x=\left(\frac{k}{x}\right)^2(k+x^2)-kx.$$

We are told that $a_k,b_k,c_k$ are the three roots of this equation, and we define

$$S(n)=\sum_{p=1}^n\sum_{k=1}^n (a_k+b_k)^p(b_k+c_k)^p(c_k+a_k)^p.$$

We must compute $S(10^6)\bmod 1,000,000,007$.

Step 1: Convert the equation into a polynomial

Starting from

$$\frac1x=\frac{k^2}{x^2}(k+x^2)-kx,$$

expand:

$$\frac1x=\frac{k^3}{x^2}+k^2-kx.$$

Multiply through by $x^2$:

$$x=k^3+k^2x^2-kx^3.$$

Rearrange:

$$kx^3-k^2x^2+x-k^3=0.$$

Thus $a_k,b_k,c_k$ are the roots of

$$kx^3-k^2x^2+x-k^3=0.$$

Divide by $k$:

$$x^3-kx^2+\frac1k x-k^2=0.$$

So for the monic cubic

$$x^3-s_1x^2+s_2x-s_3,$$

Vieta gives

$$s_1=a+b+c=k,$$

$$s_2=ab+ac+bc=\frac1k,$$

$$s_3=abc=k^2.$$

Step 2: Simplify the expression

We need

$$(a+b)(b+c)(c+a).$$

A standard symmetric identity is

$$(a+b)(b+c)(c+a) =(a+b+c)(ab+ac+bc)-abc.$$

Therefore

$$(a+b)(b+c)(c+a) =s_1s_2-s_3.$$

Substitute the Vieta values:

$$=k\cdot \frac1k-k^2 =1-k^2.$$

Hence

$$(a_k+b_k)^p(b_k+c_k)^p(c_k+a_k)^p =(1-k^2)^p.$$

$$S(n)=\sum_{k=1}^n\sum_{p=1}^n(1-k^2)^p.$$

Step 3: Evaluate the geometric series

For fixed $k$,

$$\sum_{p=1}^n r^p =\frac{r(r^n-1)}{r-1}, \quad r=1-k^2.$$

Thus

$$S(n)=\sum_{k=1}^n \frac{(1-k^2)\big((1-k^2)^n-1\big)}{-k^2}.$$

We only need this modulo

$$M=1,000,000,007.$$

Since $M$ is prime, division by $k^2$ is multiplication by the modular inverse.

Complexity is $O(n\log M)$, easily fast enough for $n=10^6$.

Step 4: Python implementation

MOD = 1_000_000_007
n = 10**6

ans = 0

for k in range(1, n + 1):
    r = (1 - k * k) % MOD

    # Sum_{p=1}^n r^p
    # = r * (r^n - 1) / (r - 1)

    numerator = r * (pow(r, n, MOD) - 1)
    denominator = (r - 1) % MOD

    geom = numerator * pow(denominator, MOD - 2, MOD)
    ans = (ans + geom) % MOD

print(ans % MOD)

Step 5: Code walkthrough

Constants

MOD = 1_000_000_007
n = 10**6

We work modulo the required prime.

Loop over $k$

for k in range(1, n + 1):

Each $k$ contributes one geometric series.

Compute $r=1-k^2$

r = (1 - k * k) % MOD

This equals

$$(a_k+b_k)(b_k+c_k)(c_k+a_k).$$

Geometric sum

numerator = r * (pow(r, n, MOD) - 1)
denominator = (r - 1) % MOD
geom = numerator * pow(denominator, MOD - 2, MOD)

Using

$$\sum_{p=1}^n r^p = \frac{r(r^n-1)}{r-1}.$$

The modular inverse is computed via Fermat’s theorem:

$$a^{-1}\equiv a^{M-2}\pmod M.$$

Accumulate

ans = (ans + geom) % MOD

Step 6: Check against the example

For $n=4$,

$$S(4)=51160,$$

which matches the problem statement.

Step 7: Final computation

Running the program for $n=10^6$ gives

$$S(10^6)\equiv 191541795 \pmod{1,000,000,007}.$$

Answer: 191541795

Project Euler Problem 479