Vào thẳng vấn đề luôn :
**Problem :** [AGENT](https://www.codechef.com/problems/AGENTS?tab=statement)
Tóm tắt đề bài : Cho trước hai đa thức $Q(x)$ và $T(x)$ . Tìm đa thức $P(x)$ thỏa mãn :
$$P(x) = Q(x) + \int_0^1 T(x.s)P(s)ds$$
Ban đầu $\deg Q(x) = n, \deg T(x) = m$. Để cho dễ tính toán, ta đặt $k = \max(n, m)$ và thêm các hệ số $0$ vào đa thức $Q(x), T(x)$
- $Q(x) = q_0 + q_1 x + q_2x^2 + \cdots + q_kx^k$ $(q_i = 0$ với $i > n)$
- $T(x) = t_0 + t_1x + t_2x^2 + \cdots + t_kx^k$ $(t_i = 0$ với $i > m)$
Đặt $P(x) = p_0 + p_1x + p_2x^2 + \cdots + p_kx^k$. Khi đó ta có :
- $T(x.s) = t_0 + t_1xs + t_2x^2s^2 + \cdots + t_kx^ks^k$
- $P(s) = p_0 + p_1s + p_2s^2 + \cdots + p_ks^k$
$$\Rightarrow T(x.s)P(s) = \sum_{i = 0}^{k}t_ix^i \sum_{j = 0}^{k} p_js^{i + j}$$
$$\Rightarrow \int_0^1 T(x.s) P(s)ds = \sum_{i = 0}^{k}t_ix^i \sum_{j = 0}^k \frac{p_j}{i + j + 1}$$
Do $P(x) - \int_0^1T(x.s)P(s)ds = Q(x)$ nên ta có $p_i - \sum_{j = 0}^k \frac{t_ip_j}{i + j + 1} = q_i \ \forall \ i = \overline{0, k}$
Từ đó ta có hệ phương trình :
$$
\begin{bmatrix}
\tiny{1-t_0} & -\frac{t_0}{2} & -\frac{t_0}{3} & \dots & -\frac{t_0}{k + 1} \\
-\frac{t_1}{2} & \frac{3-t_1}{3} & -\frac{t_1}{4} & \dots & -\frac{t_1}{k + 2} \\
-\frac{t_2}{3} & -\frac{t_2}{4} & \frac{5-t_2}{5} & \dots & -\frac{t_2}{k + 3} \\
\vdots & \vdots & \vdots &\ddots & \vdots \\
-\frac{t_k}{k + 1} & -\frac{t_k}{k + 2} & -\frac{t_k}{k + 3} & \dots & \frac{2k+ 1 - t_k}{2k + 1} \\
\end{bmatrix}
\times
\begin{bmatrix}
p_0 \\
p_1 \\
p_2 \\
\vdots \\
p_k
\end{bmatrix} = \begin{bmatrix}
q_0 \\
q_1 \\
q_2 \\
\vdots \\
q_k
\end{bmatrix}
$$
Để giải phương trình trên, ta sử dụng [khử Gauss - Jordan](https://en.wikipedia.org/wiki/Gaussian_elimination)
Và khi đó ta tìm được đa thức $P(x)$
code mẫu
```c++
#include <bits/stdc++.h>
using namespace std;
// source : https://cp-algorithms.com/linear_algebra/linear-system-gauss.html
const double EPS = 1e-9;
const int INF = 2; // it doesn't actually have to be infinity or a big number
int gauss (vector <vector<double>> a, vector<double> & ans) {
int n = (int) a.size();
int m = (int) a[0].size() - 1;
vector<int> where (m, -1);
for (int col=0, row=0; col<m && row<n; ++col) {
int sel = row;
for (int i=row; i<n; ++i)
if (abs (a[i][col]) > abs (a[sel][col]))
sel = i;
if (abs (a[sel][col]) < EPS)
continue;
for (int i=col; i<=m; ++i)
swap (a[sel][i], a[row][i]);
where[col] = row;
for (int i=0; i<n; ++i)
if (i != row) {
double c = a[i][col] / a[row][col];
for (int j=col; j<=m; ++j)
a[i][j] -= a[row][j] * c;
}
++row;
}
ans.assign (m, 0);
for (int i=0; i<m; ++i)
if (where[i] != -1)
ans[i] = a[where[i]][m] / a[where[i]][i];
for (int i=0; i<n; ++i) {
double sum = 0;
for (int j=0; j<m; ++j)
sum += ans[j] * a[i][j];
if (abs (sum - a[i][m]) > EPS)
return 0;
}
for (int i=0; i<m; ++i)
if (where[i] == -1)
return INF;
return 1;
}
void solve() {
int n; cin >> n; vector <double> q(n + 1);
for (auto &x : q) cin >> x;
int m; cin >> m; vector <double> t(m + 1);
for (auto &x : t) cin >> x;
int k = max(n, m);
while(q.size() <= k + 1) q.emplace_back(0);
while(t.size() <= k + 1) t.emplace_back(0);
vector <vector <double>> a(k + 1, vector <double> (k + 2));
for (int i = 0; i <= k; ++i) {
for (int j = 0; j <= k + 1; ++j) {
if(j == k + 1) {
a[i][j] = q[i];
} else {
if(i == j) {
a[i][j] = (2 * i + 1 - t[i]) / (i + j + 1);
} else {
a[i][j] = - t[i] / (i + j + 1);
}
}
}
}
vector <double> ans(k + 1);
assert(gauss(a, ans) == 1);
while(ans.back() == 0) {
ans.pop_back();
}
cout << (int) ans.size() - 1 << endl;
cout << fixed << setprecision(6);
for (auto &x : ans) cout << x << " ";
cout << endl;
}
signed main() {
#ifdef LOCAL
freopen("TASK.inp", "r", stdin);
freopen("TASK.out", "w", stdout);
#endif // LOCAL
ios_base::sync_with_stdio(0); cin.tie(0); cout.tie(0);
int T; cin >> T; while(T--) {
solve();
}
}
```