public class Euler1008 {
    private static final long MOD = 1_000_000_007L;
    private static final int TARGET = 10_000_000;
    // P_N has a factor x, so its x^10 coefficient comes from degree nine.
    private static final int WANTED_DEGREE = 9;
    private static final int BATCH_SIZE = 32_768;

    private static long add(long lhs, long rhs) {
        long sum = lhs + rhs;
        return sum >= MOD ? sum - MOD : sum;
    }

    private static long subtract(long lhs, long rhs) {
        return lhs >= rhs ? lhs - rhs : lhs + MOD - rhs;
    }

    private static long multiply(long lhs, long rhs) {
        // Both operands are below MOD, so their product fits in signed long.
        return lhs * rhs % MOD;
    }

    private static long modPow(long base, long exponent) {
        long result = 1;
        while (exponent != 0) {
            if ((exponent & 1L) != 0) {
                result = multiply(result, base);
            }
            base = multiply(base, base);
            exponent >>= 1;
        }
        return result;
    }

    private static long targetCoefficient(int n) {
        return targetCoefficient(n, BATCH_SIZE);
    }

    private static long targetCoefficient(int n, int batchSize) {
        if (n < 1 || 2L * n - 1 >= MOD) {
            throw new IllegalArgumentException(
                    "this implementation requires 1 <= n and 2n-1 < MOD");
        }
        if (batchSize < 1) {
            throw new IllegalArgumentException("batchSize must be positive");
        }

        // coefficients[d] is [x^d] A_k, where
        // A_k(x) = product(1 - x/j^2, j=1..k).
        long[] coefficients = new long[WANTED_DEGREE + 1];
        coefficients[0] = 1;
        long interpolationCoefficient = 0;
        long factorial = 1;
        long[] prefix = new long[Math.min(batchSize, n) + 1];

        for (int first = 1; first <= n; first += batchSize) {
            int count = Math.min(batchSize, n - first + 1);

            // Prefix products of v_k = k(2k-1).
            prefix[0] = 1;
            for (int index = 1; index <= count; ++index) {
                long k = (long) first + index - 1;
                long value = multiply(k % MOD, (2 * k - 1) % MOD);
                prefix[index] = multiply(prefix[index - 1], value);
            }

            // One Fermat inverse followed by a backward pass supplies every
            // 1/v_k in the batch.  Store those weights back into prefix[].
            long suffixInverse = modPow(prefix[count], MOD - 2);
            for (int index = count; index >= 1; --index) {
                long k = (long) first + index - 1;
                long value = multiply(k % MOD, (2 * k - 1) % MOD);
                long weight = multiply(suffixInverse, prefix[index - 1]);
                suffixInverse = multiply(suffixInverse, value);
                prefix[index] = weight;
            }

            for (int index = 1; index <= count; ++index) {
                long k = (long) first + index - 1;
                long weight = prefix[index]; // 1/(k(2k-1))
                long inverseK = multiply((2 * k - 1) % MOD, weight);
                long inverseSquare = multiply(inverseK, inverseK);

                // Newton interpolation gives
                // P_N/x = sum A_(k-1)/(k(2k-1)).  Accumulate the summand
                // before changing A_(k-1) into A_k.
                interpolationCoefficient = add(
                        interpolationCoefficient,
                        multiply(weight, coefficients[WANTED_DEGREE]));
                for (int degree = WANTED_DEGREE; degree >= 1; --degree) {
                    coefficients[degree] = subtract(
                            coefficients[degree],
                            multiply(inverseSquare, coefficients[degree - 1]));
                }

                factorial = multiply(factorial, k % MOD);
            }
        }

        long factorialSquare = multiply(factorial, factorial);
        long leading = multiply(
                n,
                modPow(multiply((2L * n - 1) % MOD, factorialSquare), MOD - 2));
        if ((n & 1) == 0) {
            leading = subtract(0, leading);
        }

        // If P_N is not monic, the lowest-degree monic solution is obtained
        // by adding x*product(x-k^2, k=1..N).
        if (leading == 1) {
            return interpolationCoefficient;
        }

        long vanishingScale = (n & 1) == 0
                ? factorialSquare
                : subtract(0, factorialSquare);
        return add(
                interpolationCoefficient,
                multiply(vanishingScale, coefficients[WANTED_DEGREE]));
    }

    private static long[] appendRoot(long[] polynomial, long root) {
        long[] result = new long[polynomial.length + 1];
        for (int degree = 0; degree < polynomial.length; ++degree) {
            result[degree] = subtract(
                    result[degree], multiply(root, polynomial[degree]));
            result[degree + 1] = add(result[degree + 1], polynomial[degree]);
        }
        return result;
    }

    // Independent cubic-time Lagrange interpolation for small checkpoints.
    private static long interpolateDirectly(int n) {
        long[] polynomial = new long[n + 2];

        for (int i = 1; i <= n; ++i) {
            long[] basis = {1};
            long denominator = 1;
            long nodeI = (long) i * i % MOD;

            for (int j = 0; j <= n; ++j) {
                if (j == i) {
                    continue;
                }
                long nodeJ = (long) j * j % MOD;
                basis = appendRoot(basis, nodeJ);
                denominator = multiply(denominator, subtract(nodeI, nodeJ));
            }

            long scale = multiply(i, modPow(denominator, MOD - 2));
            for (int degree = 0; degree < basis.length; ++degree) {
                polynomial[degree] = add(
                        polynomial[degree], multiply(scale, basis[degree]));
            }
        }

        if (polynomial[n] != 1) {
            long[] vanishing = {1};
            for (int i = 0; i <= n; ++i) {
                vanishing = appendRoot(vanishing, (long) i * i % MOD);
            }
            for (int degree = 0; degree < vanishing.length; ++degree) {
                polynomial[degree] = add(polynomial[degree], vanishing[degree]);
            }
        }

        return polynomial.length > 10 ? polynomial[10] : 0;
    }

    private static void requireCheckpoint(boolean condition, String description) {
        if (!condition) {
            throw new AssertionError("checkpoint failed: " + description);
        }
    }

    private static void runCheckpoints() {
        for (int n = 1; n <= 40; ++n) {
            long expected = interpolateDirectly(n);
            requireCheckpoint(
                    targetCoefficient(n) == expected,
                    "interpolation check for n=" + n);
            requireCheckpoint(
                    targetCoefficient(n, 7) == expected,
                    "batch check for n=" + n);
        }
    }

    public static void main(String[] args) {
        runCheckpoints();
        System.out.println(targetCoefficient(TARGET));
    }
}