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(); } } ```