
Well, I found a some kind of 'magical' solution:
I said that the answer polynomial is (59x^5  11x^3) / 3, but actually it's same modulo to 666666691x^5 + 333333332x^3 (let this the answer polynomial).
It seems like the sample input only provides (x,y)=(1,16),(3,4680),(5,61000), but it's not true. Because: 1. You can also get (x, y) = (0, 0), obviously. 2. Doing "newton's polynomial interpolation" with mod 1000000007, you can add (x, y) = (12345, 598011702), (999993007, 763641672) from sample input (999993007 is equal to 1000000000000 mod 1000000007).
Now you have 6 points! So you can get the answer polynomial!!!!!
Also, the code of newton's polynomial interpolation (my library) is this:
vector<int> interpolation(vector<int> x, int m) {
int n = x.size()  1;
vector<int> inv(n + 1); inv[1] = 1;
for (int i = 2; i <= n; i++) inv[i] = 1LL * inv[m % i] * (m  m / i) % m;
vector<vector<int> > dp(n + 1, vector<int>(n + 2));
for (int i = 0; i <= n; i++) dp[i][i + 1] = x[i];
for (int i = 2; i <= n + 1; i++) {
for (int r = i; r <= n + 1; r++) {
int l = r  i;
dp[l][r] = 1LL * (dp[l + 1][r]  dp[l][r  1] + m) * inv[i  1] % m;
}
}
vector<int> ret(n + 1); ret[0] = x[0];
vector<int> mul({ 1 });
for (int i = 0; i < n; i++) {
vector<int> mul2(i + 2, 0);
for (int j = 0; j <= i; j++) mul2[j + 1] = mul[j];
for (int j = 0; j <= i; j++) mul2[j] = (mul2[j] + 1LL * (m  i) * mul[j]) % m;
mul = mul2;
for (int j = 0; j <= i + 1; j++) ret[j] = (ret[j] + 1LL * mul[j] * dp[0][i + 2]) % m;
}
return ret;
}
