public class Euler677 {

    static final int MOD = 1000000007;

    static int modPow(long a, long e) {
        long r = 1;
        a %= MOD;
        while (e > 0) {
            if ((e & 1L) != 0)
                r = (r * a) % MOD;
            a = (a * a) % MOD;
            e >>= 1L;
        }
        return (int) r;
    }

    static int modNorm(long x) {
        x %= MOD;
        if (x < 0)
            x += MOD;
        return (int) x;
    }

    static void updateSquare(int[] s, int[] sSq, int n, int limitN) {
        int sn = s[n];
        if (sn == 0)
            return;
        int maxK = Math.min(n - 1, limitN - n);
        long add = (2L * sn) % MOD;
        for (int k = 1; k <= maxK; ++k) {
            if (s[k] != 0) {
                sSq[n + k] = (int) ((sSq[n + k] + add * s[k]) % MOD);
            }
        }
        if (2 * n <= limitN) {
            sSq[2 * n] = (int) ((sSq[2 * n] + (long) sn * sn) % MOD);
        }
    }

    static int convCubeAt(int[] s, int[] sSq, int t) {
        if (t <= 1)
            return 0;
        long sum = 0;
        for (int k = 1; k < t; ++k) {
            if (s[k] == 0 || sSq[t - k] == 0)
                continue;
            sum += (long) s[k] * sSq[t - k];
            if (sum >= 8000000000000000000L)
                sum %= MOD; // protect from long overflow
        }
        return (int) (sum % MOD);
    }

    static int convSc2At(int[] s, int t) {
        if (t <= 1)
            return 0;
        long sum = 0;
        for (int u = 2; u < t; u += 2) {
            int k = t - u;
            if (s[k] == 0 || s[u / 2] == 0)
                continue;
            sum += (long) s[k] * s[u / 2];
            if (sum >= 8000000000000000000L)
                sum %= MOD;
        }
        return (int) (sum % MOD);
    }

    static int[] buildShift(int[] s, int k, int nmax) {
        int[] out = new int[nmax + 1];
        for (int i = 1; i * k <= nmax; ++i) {
            out[i * k] = s[i];
        }
        return out;
    }

    static int[] convolutionSingle(int[] a, int[] b, int nmax) {
        int[] out = new int[nmax + 1];
        for (int i = 1; i <= nmax; ++i) {
            if (a[i] == 0)
                continue;
            long ai = a[i];
            int limit = nmax - i;
            for (int j = 1; j <= limit; ++j) {
                if (b[j] == 0)
                    continue;
                long val = out[i + j] + ai * b[j];
                out[i + j] = (int) (val % MOD);
            }
        }
        return out;
    }

    static int[] convolutionStride(int[] a, int[] base, int stride, int nmax) {
        int[] out = new int[nmax + 1];
        int maxK = nmax / stride;
        for (int k = 1; k <= maxK; ++k) {
            long coeff = base[k];
            if (coeff == 0)
                continue;
            int offset = stride * k;
            int limit = nmax - offset;
            for (int j = 1; j <= limit; ++j) {
                if (a[j] == 0)
                    continue;
                long val = out[j + offset] + a[j] * coeff;
                out[j + offset] = (int) (val % MOD);
            }
        }
        return out;
    }

    static class PlantedData {
        int[] p_r, p_b, p_y, s_all, s_no_y, s_all_sq, s_no_y_sq;
    }

    static PlantedData buildPlanted(int nmax) {
        int inv2 = modPow(2, MOD - 2);
        int inv6 = modPow(6, MOD - 2);

        int[] p_r = new int[nmax + 1];
        int[] p_b = new int[nmax + 1];
        int[] p_y = new int[nmax + 1];
        int[] s_all = new int[nmax + 1];
        int[] s_no_y = new int[nmax + 1];
        int[] s_all_sq = new int[nmax + 1];
        int[] s_no_y_sq = new int[nmax + 1];

        for (int n = 1; n <= nmax; ++n) {
            int t = n - 1;
            int base = (t == 0) ? 1 : 0;

            int s1_all = s_all[t];
            int s2_all = s_all_sq[t];
            int s3_all = convCubeAt(s_all, s_all_sq, t);
            int sc2_all = convSc2At(s_all, t);
            int c2_all = (t % 2 == 0) ? s_all[t / 2] : 0;
            int c3_all = (t % 3 == 0) ? s_all[t / 3] : 0;

            int m1_all = s1_all;
            int m2_all = (int) ((1L * (s2_all + c2_all) % MOD) * inv2 % MOD);
            int m3_all = (int) ((1L * (s3_all + 3L * sc2_all + 2L * c3_all) % MOD) * inv6 % MOD);

            int s1_no = s_no_y[t];
            int s2_no = s_no_y_sq[t];
            int s3_no = convCubeAt(s_no_y, s_no_y_sq, t);
            int sc2_no = convSc2At(s_no_y, t);
            int c2_no = (t % 2 == 0) ? s_no_y[t / 2] : 0;
            int c3_no = (t % 3 == 0) ? s_no_y[t / 3] : 0;

            int m1_no = s1_no;
            int m2_no = (int) ((1L * (s2_no + c2_no) % MOD) * inv2 % MOD);
            int m3_no = (int) ((1L * (s3_no + 3L * sc2_no + 2L * c3_no) % MOD) * inv6 % MOD);

            p_r[n] = modNorm((long) base + m1_all + m2_all + m3_all);
            p_b[n] = modNorm((long) base + m1_all + m2_all);
            p_y[n] = modNorm((long) base + m1_no + m2_no);

            s_all[n] = modNorm((long) p_r[n] + p_b[n] + p_y[n]);
            s_no_y[n] = modNorm((long) p_r[n] + p_b[n]);

            updateSquare(s_all, s_all_sq, n, nmax);
            updateSquare(s_no_y, s_no_y_sq, n, nmax);
        }

        PlantedData pd = new PlantedData();
        pd.p_r = p_r;
        pd.p_b = p_b;
        pd.p_y = p_y;
        pd.s_all = s_all;
        pd.s_no_y = s_no_y;
        pd.s_all_sq = s_all_sq;
        pd.s_no_y_sq = s_no_y_sq;
        return pd;
    }

    public static String solve() {
        int nmax = 10000;
        int inv2 = modPow(2, MOD - 2);
        int inv6 = modPow(6, MOD - 2);
        int inv24 = modPow(24, MOD - 2);

        PlantedData pd = buildPlanted(nmax);
        int[] p_all = pd.s_all;
        int[] p_no_y = pd.s_no_y;
        int[] p_y = pd.p_y;

        int[] c2_all = buildShift(p_all, 2, nmax);
        int[] c3_all = buildShift(p_all, 3, nmax);
        int[] c4_all = buildShift(p_all, 4, nmax);

        int[] c2_no = buildShift(p_no_y, 2, nmax);
        int[] c3_no = buildShift(p_no_y, 3, nmax);

        int[] all_cube = convolutionSingle(pd.s_all_sq, p_all, nmax);
        int[] all_four = convolutionSingle(pd.s_all_sq, pd.s_all_sq, nmax);
        int[] all_c2 = convolutionStride(p_all, p_all, 2, nmax);
        int[] all_sq_c2 = convolutionStride(pd.s_all_sq, p_all, 2, nmax);
        int[] all_c3 = convolutionStride(p_all, p_all, 3, nmax);

        int[] no_cube = convolutionSingle(pd.s_no_y_sq, p_no_y, nmax);
        int[] no_c2 = convolutionStride(p_no_y, p_no_y, 2, nmax);

        int[] c2_sq_all = new int[nmax + 1];
        for (int i = 1; 2 * i <= nmax; ++i) {
            c2_sq_all[2 * i] = pd.s_all_sq[i];
        }

        int[] m1_all = p_all;
        int[] m2_all = new int[nmax + 1];
        int[] m3_all = new int[nmax + 1];
        int[] m4_all = new int[nmax + 1];

        int[] m1_no = p_no_y;
        int[] m2_no = new int[nmax + 1];
        int[] m3_no = new int[nmax + 1];

        for (int t = 0; t <= nmax; ++t) {
            m2_all[t] = (int) ((1L * (pd.s_all_sq[t] + c2_all[t]) % MOD) * inv2 % MOD);
            m2_no[t] = (int) ((1L * (pd.s_no_y_sq[t] + c2_no[t]) % MOD) * inv2 % MOD);

            m3_all[t] = (int) ((1L * (all_cube[t] + 3L * all_c2[t] + 2L * c3_all[t]) % MOD) * inv6 % MOD);
            m3_no[t] = (int) ((1L * (no_cube[t] + 3L * no_c2[t] + 2L * c3_no[t]) % MOD) * inv6 % MOD);

            m4_all[t] = (int) ((1L
                    * (all_four[t] + 6L * all_sq_c2[t] + 3L * c2_sq_all[t] + 8L * all_c3[t] + 6L * c4_all[t]) % MOD)
                    * inv24 % MOD);
        }

        int[] a_all = new int[nmax + 1];
        int[] a_r = new int[nmax + 1];
        int[] a_b = new int[nmax + 1];
        int[] a_y = new int[nmax + 1];

        for (int n = 1; n <= nmax; ++n) {
            int t = n - 1;
            int base = (t == 0) ? 1 : 0;
            a_r[n] = modNorm((long) base + m1_all[t] + m2_all[t] + m3_all[t] + m4_all[t]);
            a_b[n] = modNorm((long) base + m1_all[t] + m2_all[t] + m3_all[t]);
            a_y[n] = modNorm((long) base + m1_no[t] + m2_no[t] + m3_no[t]);
            a_all[n] = modNorm((long) a_r[n] + a_b[n] + a_y[n]);
        }

        int[] p_y_sq = convolutionSingle(p_y, p_y, nmax);
        int[] p_all_x2 = buildShift(p_all, 2, nmax);
        int[] p_y_x2 = buildShift(p_y, 2, nmax);

        int[] g = new int[nmax + 1];
        for (int n = 1; n <= nmax; ++n) {
            int d = modNorm((long) pd.s_all_sq[n] - p_y_sq[n]);
            int e = modNorm((long) pd.s_all_sq[n] + p_all_x2[n] - p_y_sq[n] - p_y_x2[n]);
            int e_half = (int) ((1L * e) * inv2 % MOD);
            g[n] = modNorm((long) a_all[n] + e_half - d);
        }

        return Integer.toString(g[nmax]);
    }

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