/*
Tool chi dung trong truong hop ham f thoa man 2 dieu kien:
- Ham f bieu dien duoc duoi dang 1/2 * xT * A * x + bT * x + c
- A la ma tran xac dinh duong
Input:
- Dong dau la 'n': So chieu cua khong gian
- Dong thu hai gom (n+1)*n/2 so: Cac he so cua ham trong bieu dien phan bac 2
- Dong thu ba gom n so: Cac he so cua ham trong bieu dien phan bac 1
- Dong thu tu gom 1 so: He so tu do
- Dong thu nam gom n so: Diem xuat phat
Output:
- Bieu dien qua trinh thuc hien quy tac Exact Line Search/Armijo duoi dang bang
Vi du:
Input:
4
5 1 0 1 4 1 0 3 1 2
2 4 -7 8
6
4 2 3 7
-> f(x,y,z,t) = 5x^2 + 1xy + 0xz + 1xt
+ 4y^2 + 1yz + 0yt
+ 3z^2 + 1zt
+ 2t^2
+ 2x + 4y -7z + 8t
+ 6
Xuat phat: (x0; y0; z0; t0) = (4; 2; 3; 7);
Output:
-> Nghiem toi uu cho bai toan min f(x,y,z,t) tren R4:
(0.118055555589, -0.726851851816, 1.696759259221, -2.453703703590)
Gia tri nho nhat dat duoc la -11.089120370370
*/
//------- Khoi tao ------------
#include <bits/stdc++.h>
#define up(i,a,b) for (int i = (int)a; i <= (int)b; i++)
#define PB push_back
#define ld long double
using namespace std;
const int maxn = 1e5 + 10;
ld epsilon;
int precision_digit;
template <typename T>
string to_string_precs(const T a_value, const int n = precision_digit){
std::ostringstream out;
out.precision(n);
out << std::fixed << a_value;
return std::move(out).str();
}
vector<ld> x[maxn];
ld t[maxn], res[maxn];
vector<string> fin[maxn];
string description[maxn];
string indent[maxn];
//--------- Cac ham quan trong ---------------
//Tich vo huong cua hai vector
ld DOT(const vector<ld>& A, const vector<ld>& B, int n){
ld sum = 0;
up(i, 0, n-1) sum += A[i] * B[i];
return sum;
}
//Gia tri ham f doi voi ma tran A, B, C tai diem X, khong gian n chieu.
//f = 1/2 * xT * A * x + bT * x + c
ld f(const vector<vector<ld> >& A, const vector<ld>& B, const ld& C, const vector<ld>& X, const int& n){
vector<ld> temp(n, 0);
up(j, 0, n-1){
for (int i = 0; i <= j; i++){
ld H = A[i][j];
if (i == j) H /= 2;
temp[j] += X[i] * H;
}
}
ld sum = DOT(temp, X, n) + DOT(B, X, n) + C;
return sum;
}
//Gradient tai diem X
vector<ld> GRAD(const vector<vector<ld> >& A, const vector<ld>& B, const vector<ld>& X, const int& n){
vector<ld> temp(n, 0);
up(j, 0, n-1) up(i, 0, n-1){
temp[j] += X[i] * A[i][j];
}
up(i, 0, n-1) temp[i] += B[i];
return temp;
}
//Am Gradient tai diem X
vector<ld> MINUS_GRAD(const vector<vector<ld> >& A, const vector<ld>& B, const vector<ld>& X, const int& n){
vector<ld> temp(n, 0);
temp = GRAD(A, B, X, n);
up(i, 0, n-1) temp[i] = -temp[i];
return temp;
}
//Tich cua mot vector [1 x n] voi mot ma tran [n x n] -> vector [1 x n]
vector<ld> MXV(const vector<vector<ld> >& A, const vector<ld>& X, const int& n){
vector<ld> temp(n, 0);
up(j, 0, n-1) up(i, 0, n-1){
temp[j] += X[i] * A[i][j];
}
return temp;
}
//X_(k+1) = X_(k) + t_(k)*d_(k)
vector<ld> FIND_NEXT(const vector<ld>& X, const ld& T, const vector<ld>& D, const int& n){
vector<ld> temp(n, 0);
up(i, 0, n-1) temp[i] = X[i] + T*D[i];
return temp;
}
//------- Exact line search starts ---------
ld Exact_line_search (const vector<vector<ld> >& A, const vector<ld>& B, const ld& C, const vector<ld>& X, const vector<ld>& d, const int& n){
return (DOT(d, d, n))/DOT(MXV(A, d, n), d, n);
}
//----- Armijo starts ---------
int cnt = 0;
ld t0, m, alpha, beta;
ld Armijo_step_size(const vector<vector<ld> >& A, const vector<ld>& B, const ld& C, const vector<ld>& X, const vector<ld>& d, const int& n) {
ld t = t0; // Reset t ve t0 cap nhat duoc gan nhat trong moi lan tim co buoc phu hop
ld fx = f(A, B, C, X, n);
ld grad_dot_d = DOT(GRAD(A, B, X, n), d, n);
while (true) {
vector<ld> X_new = FIND_NEXT(X, t, d, n);
ld fx_new = f(A, B, C, X_new, n);
if (fx_new - fx <= m * t * grad_dot_d) {
t0 = t;
return t;
}
t *= alpha;
if (t < beta) break;
++cnt;
}
return t;
}
//------ GRADIENT DESCENT CONCEPT --------
void Gradient_Descent(vector<vector<ld> >& A, vector<ld>& B, ld& C, vector<ld>& X, int n, int& lim,
ld (*kind)(const vector<vector<ld> >&, const vector<ld>&, const ld&, const vector<ld>&, const vector<ld>&, const int&)){
x[0] = X;
res[0] = f(A, B, C, x[0], n);
for (int i = 1; i <= lim; i++){
vector<ld> d = MINUS_GRAD(A, B, x[i-1], n);
if (sqrt(DOT(d, d, n)) < epsilon){
cout << "Desired stationary point found at iteration: " << i-1 << "\n\n";
lim = i-1;
break;
}
t[i-1] = kind(A, B, C, X, d, n);
x[i] = FIND_NEXT(x[i-1], t[i-1], d, n);
res[i] = f(A, B, C, x[i], n);
}
}
//------ INPUT + PREPROCESSING ----------
void input(vector<vector<ld> >& A, vector<ld>& B, ld& C, vector<ld>& X, int& n){
up(i, 0, n-1) up(j, i, n-1) cin >> A[i][j];
//xx xy xz yy yz zz
up(i, 0, n-1) cin >> B[i];
//x y z
cin >> C;
//C
up(i, 0, n-1) cin >> X[i];
//x1 x2 x3
up(i, 0, n-1){
up(j, i, n-1){
if (i == j) A[i][j] *= 2;
else A[j][i] = A[i][j];
}
}
// A in 1/2(xT*A*x)
}
// ------ OUTPUT ---------
void parse_indent(const vector<string>& fin, const int& max_length){
int k = (int)fin.size();
description[0] = "i";
description[k-1] = "delta f";
description[k-1 - 1] = "ti";
description[k-1 - 2] = "f";
for (int i = 1; i < k-1 - 2; i++){
string cur_description = "x" + to_string(i);
description[i] = cur_description;
}
up(i, 0, k-1) indent[i] = string(max_length+2, ' ');
}
void print_columns(const vector<string>& fin, const int& max_length){
up(i, 0, fin.size()-1){
cout << description[i];
cout << indent[i];
}
cout << "\n";
}
void print_result_table(const int& lim, const int& max_length){
up(i,0,lim){
for (int j = 0; j < (int)fin[i].size(); j++){
cout << fin[i][j] << " ";
if ((int)fin[i][j].size() < max_length){
int dotdot = max_length - (int)fin[i][j].size();
if (j == 0) dotdot -= (max_length - 7);
cout << string(dotdot, ' ');
}
cout << "| ";
}
cout << "\n";
}
cout << "\n";
}
//--------- MAIN -------------
signed main(){
ios_base::sync_with_stdio(false);
cin.tie(0);
#define Task "A"
if (fopen(Task".inp", "r")){
freopen(Task".inp", "r", stdin);
freopen(Task".out", "w", stdout);
}
cout << fixed << setprecision(12);
// --------- INPUT ------------
int n;
cin >> n;
vector<vector<ld> > A(n, vector<ld>(n, 0));
vector<ld> B(n, 0);
ld C;
vector<ld> X(n, 0);
input(A, B, C, X, n);
// --------- ALGORITHM ----------
/**
int lim = 10000; // So lan lap toi da
ld t0 = 1; // t0 ban dau trong thuat toan Armijo
ld m = 0.5; // He so dieu kien: nen dat trong khoang [0.5 -> 0.7]
ld alpha = 0.9; // Learning rate: nen dat trong khoang [0.5 -> 0.95]
ld beta = 1e-6; // Co buoc toi thieu: cham qua thi thoi
ld epsilon = 1e-9; // Gia tri toi thieu can dat duoc cua NORM(GRAD(diem dung)): cang chuan -> chay cang cham
int precision_digit = 12; // So chu so sau dau phay khi hien thi ket qua
**/
int lim = 50000;
t0 = 1;
m = 0.5;
alpha = 0.8;
beta = 1e-6;
epsilon = 1e-5;
precision_digit = 5;
Gradient_Descent(A, B, C, X, n, lim, Armijo_step_size);
// Gradient_Descent(A, B, C, X, n, lim, Exact_line_search);
//-------- Parsing Output -----------
up(i, 0, lim){
fin[i].PB(to_string_precs(i));
up(j, 0, n-1){
fin[i].PB(to_string_precs(x[i][j]));
}
fin[i].PB(to_string_precs(res[i]));
fin[i].PB(to_string_precs(t[i]));
fin[i].PB(to_string_precs(res[i] - res[i-1]));
}
int max_length = 0;
up(i, 0, lim){
for (int j = 0; j < (int)fin[i].size(); j++){
max_length = max(max_length, (int)fin[i][j].size());
}
}
//---- Printing the result table ----
parse_indent(fin[0], max_length);
print_columns(fin[0], max_length);
print_result_table(lim, max_length);
cout << "Tong so buoc t = alpha*t trong thuat toan Armijo: " << cnt;
}
LyoKVG9vbCBjaGkgZHVuZyB0cm9uZyB0cnVvbmcgaG9wIGhhbSBmIHRob2EgbWFuIDIgZGlldSBraWVuOgotIEhhbSBmIGJpZXUgZGllbiBkdW9jIGR1b2kgZGFuZyAxLzIgKiB4VCAqIEEgKiB4ICsgYlQgKiB4ICsgYwotIEEgbGEgbWEgdHJhbiB4YWMgZGluaCBkdW9uZwoKCgoKSW5wdXQ6Ci0gRG9uZyBkYXUgbGEgJ24nOiBTbyBjaGlldSBjdWEga2hvbmcgZ2lhbgotIERvbmcgdGh1IGhhaSBnb20gKG4rMSkqbi8yIHNvOiBDYWMgaGUgc28gY3VhIGhhbSB0cm9uZyBiaWV1IGRpZW4gcGhhbiBiYWMgMgotIERvbmcgdGh1IGJhIGdvbSBuIHNvOiBDYWMgaGUgc28gY3VhIGhhbSB0cm9uZyBiaWV1IGRpZW4gcGhhbiBiYWMgMQotIERvbmcgdGh1IHR1IGdvbSAxIHNvOiBIZSBzbyB0dSBkbwotIERvbmcgdGh1IG5hbSBnb20gbiBzbzogRGllbSB4dWF0IHBoYXQKCk91dHB1dDoKLSBCaWV1IGRpZW4gcXVhIHRyaW5oIHRodWMgaGllbiBxdXkgdGFjIEV4YWN0IExpbmUgU2VhcmNoL0FybWlqbyBkdW9pIGRhbmcgYmFuZwoKCgoKClZpIGR1OgpJbnB1dDoKNAo1IDEgMCAxIDQgMSAwIDMgMSAyCjIgNCAtNyA4CjYKNCAyIDMgNwoKLT4gZih4LHkseix0KSA9IDV4XjIgKyAxeHkgICsgMHh6ICArIDF4dAoJCQkgICAgICAgICArIDR5XjIgKyAxeXogICsgMHl0CgkJCSAgICAgICAgCSAgICArIDN6XjIgKyAxenQKCQkJICAgICAgICAJCQkgICArIDJ0XjIKCQkJICArIDJ4ICAgKyA0eSAgIC03eiAgICArIDh0CgkJCSAgKyA2CgpYdWF0IHBoYXQ6ICh4MDsgeTA7IHowOyB0MCkgPSAoNDsgMjsgMzsgNyk7CgoKT3V0cHV0OgotPiBOZ2hpZW0gdG9pIHV1IGNobyBiYWkgdG9hbiBtaW4gZih4LHkseix0KSB0cmVuIFI0OgooMC4xMTgwNTU1NTU1ODksICAtMC43MjY4NTE4NTE4MTYsICAxLjY5Njc1OTI1OTIyMSwgIC0yLjQ1MzcwMzcwMzU5MCkKR2lhIHRyaSBuaG8gbmhhdCBkYXQgZHVvYyBsYSAtMTEuMDg5MTIwMzcwMzcwCiovCgoKCi8vLS0tLS0tLSBLaG9pIHRhbyAtLS0tLS0tLS0tLS0KI2luY2x1ZGUgPGJpdHMvc3RkYysrLmg+CiNkZWZpbmUgdXAoaSxhLGIpIGZvciAoaW50IGkgPSAoaW50KWE7IGkgPD0gKGludCliOyBpKyspCiNkZWZpbmUgUEIgcHVzaF9iYWNrCiNkZWZpbmUgbGQgbG9uZyBkb3VibGUKdXNpbmcgbmFtZXNwYWNlIHN0ZDsKCmNvbnN0IGludCBtYXhuID0gMWU1ICsgMTA7CmxkIGVwc2lsb247CmludCBwcmVjaXNpb25fZGlnaXQ7Cgp0ZW1wbGF0ZSA8dHlwZW5hbWUgVD4Kc3RyaW5nIHRvX3N0cmluZ19wcmVjcyhjb25zdCBUIGFfdmFsdWUsIGNvbnN0IGludCBuID0gcHJlY2lzaW9uX2RpZ2l0KXsKICAgIHN0ZDo6b3N0cmluZ3N0cmVhbSBvdXQ7CiAgICBvdXQucHJlY2lzaW9uKG4pOwogICAgb3V0IDw8IHN0ZDo6Zml4ZWQgPDwgYV92YWx1ZTsKICAgIHJldHVybiBzdGQ6Om1vdmUob3V0KS5zdHIoKTsKfQoKdmVjdG9yPGxkPiB4W21heG5dOwpsZCB0W21heG5dLCByZXNbbWF4bl07CnZlY3RvcjxzdHJpbmc+IGZpblttYXhuXTsKc3RyaW5nIGRlc2NyaXB0aW9uW21heG5dOwpzdHJpbmcgaW5kZW50W21heG5dOwoKCgoKCgovLy0tLS0tLS0tLSBDYWMgaGFtIHF1YW4gdHJvbmcgLS0tLS0tLS0tLS0tLS0tCi8vVGljaCB2byBodW9uZyBjdWEgaGFpIHZlY3RvcgpsZCBET1QoY29uc3QgdmVjdG9yPGxkPiYgQSwgY29uc3QgdmVjdG9yPGxkPiYgQiwgaW50IG4pewogICAgbGQgc3VtID0gMDsKICAgIHVwKGksIDAsIG4tMSkgc3VtICs9IEFbaV0gKiBCW2ldOwogICAgcmV0dXJuIHN1bTsKfQoKLy9HaWEgdHJpIGhhbSBmIGRvaSB2b2kgbWEgdHJhbiBBLCBCLCBDIHRhaSBkaWVtIFgsIGtob25nIGdpYW4gbiBjaGlldS4KLy9mID0gMS8yICogeFQgKiBBICogeCArIGJUICogeCArIGMKbGQgZihjb25zdCB2ZWN0b3I8dmVjdG9yPGxkPiA+JiBBLCBjb25zdCB2ZWN0b3I8bGQ+JiBCLCBjb25zdCBsZCYgQywgY29uc3QgdmVjdG9yPGxkPiYgWCwgY29uc3QgaW50JiBuKXsKICAgIHZlY3RvcjxsZD4gdGVtcChuLCAwKTsKICAgIHVwKGosIDAsIG4tMSl7CiAgICAgICAgZm9yIChpbnQgaSA9IDA7IGkgPD0gajsgaSsrKXsKICAgICAgICAgICAgbGQgSCA9IEFbaV1bal07CiAgICAgICAgICAgIGlmIChpID09IGopIEggLz0gMjsKICAgICAgICAgICAgdGVtcFtqXSArPSBYW2ldICogSDsKICAgICAgICB9CiAgICB9CgogICAgbGQgc3VtID0gRE9UKHRlbXAsIFgsIG4pICsgRE9UKEIsIFgsIG4pICsgQzsKICAgIHJldHVybiBzdW07Cn0KCi8vR3JhZGllbnQgdGFpIGRpZW0gWAp2ZWN0b3I8bGQ+IEdSQUQoY29uc3QgdmVjdG9yPHZlY3RvcjxsZD4gPiYgQSwgY29uc3QgdmVjdG9yPGxkPiYgQiwgY29uc3QgdmVjdG9yPGxkPiYgWCwgY29uc3QgaW50JiBuKXsKICAgIHZlY3RvcjxsZD4gdGVtcChuLCAwKTsKICAgIHVwKGosIDAsIG4tMSkgdXAoaSwgMCwgbi0xKXsKICAgICAgICB0ZW1wW2pdICs9IFhbaV0gKiBBW2ldW2pdOwogICAgfQogICAgdXAoaSwgMCwgbi0xKSB0ZW1wW2ldICs9IEJbaV07CgogICAgcmV0dXJuIHRlbXA7Cn0KCi8vQW0gR3JhZGllbnQgdGFpIGRpZW0gWAp2ZWN0b3I8bGQ+IE1JTlVTX0dSQUQoY29uc3QgdmVjdG9yPHZlY3RvcjxsZD4gPiYgQSwgY29uc3QgdmVjdG9yPGxkPiYgQiwgY29uc3QgdmVjdG9yPGxkPiYgWCwgY29uc3QgaW50JiBuKXsKICAgIHZlY3RvcjxsZD4gdGVtcChuLCAwKTsKICAgIHRlbXAgPSBHUkFEKEEsIEIsIFgsIG4pOwogICAgdXAoaSwgMCwgbi0xKSB0ZW1wW2ldID0gLXRlbXBbaV07CiAgICByZXR1cm4gdGVtcDsKfQoKLy9UaWNoIGN1YSBtb3QgdmVjdG9yIFsxIHggbl0gdm9pIG1vdCBtYSB0cmFuIFtuIHggbl0gLT4gdmVjdG9yIFsxIHggbl0KdmVjdG9yPGxkPiBNWFYoY29uc3QgdmVjdG9yPHZlY3RvcjxsZD4gPiYgQSwgY29uc3QgdmVjdG9yPGxkPiYgWCwgY29uc3QgaW50JiBuKXsKICAgIHZlY3RvcjxsZD4gdGVtcChuLCAwKTsKICAgIHVwKGosIDAsIG4tMSkgdXAoaSwgMCwgbi0xKXsKICAgICAgICB0ZW1wW2pdICs9IFhbaV0gKiBBW2ldW2pdOwogICAgfQogICAgcmV0dXJuIHRlbXA7Cn0KCi8vWF8oaysxKSA9IFhfKGspICsgdF8oaykqZF8oaykKdmVjdG9yPGxkPiBGSU5EX05FWFQoY29uc3QgdmVjdG9yPGxkPiYgWCwgY29uc3QgbGQmIFQsIGNvbnN0IHZlY3RvcjxsZD4mIEQsIGNvbnN0IGludCYgbil7CiAgICB2ZWN0b3I8bGQ+IHRlbXAobiwgMCk7CiAgICB1cChpLCAwLCBuLTEpIHRlbXBbaV0gPSBYW2ldICsgVCpEW2ldOwogICAgcmV0dXJuIHRlbXA7Cn0KCgoKCgoKCgoKLy8tLS0tLS0tIEV4YWN0IGxpbmUgc2VhcmNoIHN0YXJ0cyAtLS0tLS0tLS0KbGQgRXhhY3RfbGluZV9zZWFyY2ggKGNvbnN0IHZlY3Rvcjx2ZWN0b3I8bGQ+ID4mIEEsIGNvbnN0IHZlY3RvcjxsZD4mIEIsIGNvbnN0IGxkJiBDLCBjb25zdCB2ZWN0b3I8bGQ+JiBYLCBjb25zdCB2ZWN0b3I8bGQ+JiBkLCBjb25zdCBpbnQmIG4pewogICAgcmV0dXJuIChET1QoZCwgZCwgbikpL0RPVChNWFYoQSwgZCwgbiksIGQsIG4pOwp9CgovLy0tLS0tIEFybWlqbyBzdGFydHMgLS0tLS0tLS0tCmludCBjbnQgPSAwOwpsZCB0MCwgbSwgYWxwaGEsIGJldGE7CmxkIEFybWlqb19zdGVwX3NpemUoY29uc3QgdmVjdG9yPHZlY3RvcjxsZD4gPiYgQSwgY29uc3QgdmVjdG9yPGxkPiYgQiwgY29uc3QgbGQmIEMsIGNvbnN0IHZlY3RvcjxsZD4mIFgsIGNvbnN0IHZlY3RvcjxsZD4mIGQsIGNvbnN0IGludCYgbikgewogICAgbGQgdCA9IHQwOyAgLy8gUmVzZXQgdCB2ZSB0MCBjYXAgbmhhdCBkdW9jIGdhbiBuaGF0IHRyb25nIG1vaSBsYW4gdGltIGNvIGJ1b2MgcGh1IGhvcAoKICAgIGxkIGZ4ID0gZihBLCBCLCBDLCBYLCBuKTsKICAgIGxkIGdyYWRfZG90X2QgPSBET1QoR1JBRChBLCBCLCBYLCBuKSwgZCwgbik7CiAgICB3aGlsZSAodHJ1ZSkgewogICAgICAgIHZlY3RvcjxsZD4gWF9uZXcgPSBGSU5EX05FWFQoWCwgdCwgZCwgbik7CiAgICAgICAgbGQgZnhfbmV3ID0gZihBLCBCLCBDLCBYX25ldywgbik7CiAgICAgICAgaWYgKGZ4X25ldyAtIGZ4IDw9IG0gKiB0ICogZ3JhZF9kb3RfZCkgewogICAgICAgICAgICB0MCA9IHQ7CiAgICAgICAgICAgIHJldHVybiB0OwogICAgICAgIH0KICAgICAgICB0ICo9IGFscGhhOwogICAgICAgIGlmICh0IDwgYmV0YSkgYnJlYWs7CiAgICAgICAgKytjbnQ7CiAgICB9CiAgICByZXR1cm4gdDsKfQoKCgoKCgoKLy8tLS0tLS0gR1JBRElFTlQgREVTQ0VOVCBDT05DRVBUIC0tLS0tLS0tCnZvaWQgR3JhZGllbnRfRGVzY2VudCh2ZWN0b3I8dmVjdG9yPGxkPiA+JiBBLCB2ZWN0b3I8bGQ+JiBCLCBsZCYgQywgdmVjdG9yPGxkPiYgWCwgaW50IG4sIGludCYgbGltLAogICAgICAgICAgIGxkICgqa2luZCkoY29uc3QgdmVjdG9yPHZlY3RvcjxsZD4gPiYsIGNvbnN0IHZlY3RvcjxsZD4mLCBjb25zdCBsZCYsIGNvbnN0IHZlY3RvcjxsZD4mLCBjb25zdCB2ZWN0b3I8bGQ+JiwgY29uc3QgaW50JikpewogICAgeFswXSA9IFg7CiAgICByZXNbMF0gPSBmKEEsIEIsIEMsIHhbMF0sIG4pOwoKICAgIGZvciAoaW50IGkgPSAxOyBpIDw9IGxpbTsgaSsrKXsKICAgICAgICB2ZWN0b3I8bGQ+IGQgPSBNSU5VU19HUkFEKEEsIEIsIHhbaS0xXSwgbik7CiAgICAgICAgaWYgKHNxcnQoRE9UKGQsIGQsIG4pKSA8IGVwc2lsb24pewogICAgICAgICAgICBjb3V0IDw8ICJEZXNpcmVkIHN0YXRpb25hcnkgcG9pbnQgZm91bmQgYXQgaXRlcmF0aW9uOiAiIDw8IGktMSA8PCAiXG5cbiI7CiAgICAgICAgICAgIGxpbSA9IGktMTsKICAgICAgICAgICAgYnJlYWs7CiAgICAgICAgfQoKICAgICAgICB0W2ktMV0gPSBraW5kKEEsIEIsIEMsIFgsIGQsIG4pOwogICAgICAgIHhbaV0gPSBGSU5EX05FWFQoeFtpLTFdLCB0W2ktMV0sIGQsIG4pOwogICAgICAgIHJlc1tpXSA9IGYoQSwgQiwgQywgeFtpXSwgbik7CiAgICB9Cn0KCgoKCgoKCgovLy0tLS0tLSBJTlBVVCArIFBSRVBST0NFU1NJTkcgLS0tLS0tLS0tLQp2b2lkIGlucHV0KHZlY3Rvcjx2ZWN0b3I8bGQ+ID4mIEEsIHZlY3RvcjxsZD4mIEIsIGxkJiBDLCB2ZWN0b3I8bGQ+JiBYLCBpbnQmIG4pewogICAgdXAoaSwgMCwgbi0xKSB1cChqLCBpLCBuLTEpIGNpbiA+PiBBW2ldW2pdOwogICAgLy94eCB4eSB4eiB5eSB5eiB6egoKICAgIHVwKGksIDAsIG4tMSkgY2luID4+IEJbaV07CiAgICAvL3ggeSB6CgogICAgY2luID4+IEM7CiAgICAvL0MKCiAgICB1cChpLCAwLCBuLTEpIGNpbiA+PiBYW2ldOwogICAgLy94MSB4MiB4MwoKICAgIHVwKGksIDAsIG4tMSl7CiAgICAgICAgdXAoaiwgaSwgbi0xKXsKICAgICAgICAgICAgaWYgKGkgPT0gaikgQVtpXVtqXSAqPSAyOwogICAgICAgICAgICBlbHNlIEFbal1baV0gPSBBW2ldW2pdOwogICAgICAgIH0KICAgIH0KICAgIC8vIEEgaW4gMS8yKHhUKkEqeCkKfQoKCgoKCgoKLy8gLS0tLS0tIE9VVFBVVCAtLS0tLS0tLS0Kdm9pZCBwYXJzZV9pbmRlbnQoY29uc3QgdmVjdG9yPHN0cmluZz4mIGZpbiwgY29uc3QgaW50JiBtYXhfbGVuZ3RoKXsKICAgIGludCBrID0gKGludClmaW4uc2l6ZSgpOwogICAgZGVzY3JpcHRpb25bMF0gPSAiaSI7CiAgICBkZXNjcmlwdGlvbltrLTFdID0gImRlbHRhIGYiOwogICAgZGVzY3JpcHRpb25bay0xIC0gMV0gPSAidGkiOwogICAgZGVzY3JpcHRpb25bay0xIC0gMl0gPSAiZiI7CgogICAgZm9yIChpbnQgaSA9IDE7IGkgPCBrLTEgLSAyOyBpKyspewogICAgICAgIHN0cmluZyBjdXJfZGVzY3JpcHRpb24gPSAieCIgKyB0b19zdHJpbmcoaSk7CiAgICAgICAgZGVzY3JpcHRpb25baV0gPSBjdXJfZGVzY3JpcHRpb247CiAgICB9CgogICAgdXAoaSwgMCwgay0xKSBpbmRlbnRbaV0gPSBzdHJpbmcobWF4X2xlbmd0aCsyLCAnICcpOwp9Cgp2b2lkIHByaW50X2NvbHVtbnMoY29uc3QgdmVjdG9yPHN0cmluZz4mIGZpbiwgY29uc3QgaW50JiBtYXhfbGVuZ3RoKXsKICAgIHVwKGksIDAsIGZpbi5zaXplKCktMSl7CiAgICAgICAgY291dCA8PCBkZXNjcmlwdGlvbltpXTsKICAgICAgICBjb3V0IDw8IGluZGVudFtpXTsKICAgIH0KICAgIGNvdXQgPDwgIlxuIjsKfQoKdm9pZCBwcmludF9yZXN1bHRfdGFibGUoY29uc3QgaW50JiBsaW0sIGNvbnN0IGludCYgbWF4X2xlbmd0aCl7CiAgICB1cChpLDAsbGltKXsKICAgICAgICBmb3IgKGludCBqID0gMDsgaiA8IChpbnQpZmluW2ldLnNpemUoKTsgaisrKXsKICAgICAgICAgICAgY291dCA8PCBmaW5baV1bal0gPDwgIiAiOwogICAgICAgICAgICBpZiAoKGludClmaW5baV1bal0uc2l6ZSgpIDwgbWF4X2xlbmd0aCl7CiAgICAgICAgICAgICAgICBpbnQgZG90ZG90ID0gbWF4X2xlbmd0aCAtIChpbnQpZmluW2ldW2pdLnNpemUoKTsKICAgICAgICAgICAgICAgIGlmIChqID09IDApIGRvdGRvdCAtPSAobWF4X2xlbmd0aCAtIDcpOwogICAgICAgICAgICAgICAgY291dCA8PCBzdHJpbmcoZG90ZG90LCAnICcpOwogICAgICAgICAgICB9CiAgICAgICAgICAgIGNvdXQgPDwgInwgICI7CiAgICAgICAgfQogICAgICAgIGNvdXQgPDwgIlxuIjsKICAgIH0KICAgIGNvdXQgPDwgIlxuIjsKfQoKCgoKCgoKCi8vLS0tLS0tLS0tIE1BSU4gLS0tLS0tLS0tLS0tLQpzaWduZWQgbWFpbigpewogICAgaW9zX2Jhc2U6OnN5bmNfd2l0aF9zdGRpbyhmYWxzZSk7CiAgICBjaW4udGllKDApOwogICAgI2RlZmluZSBUYXNrICJBIgogICAgaWYgKGZvcGVuKFRhc2siLmlucCIsICJyIikpewogICAgICAgIGZyZW9wZW4oVGFzayIuaW5wIiwgInIiLCBzdGRpbik7CiAgICAgICAgZnJlb3BlbihUYXNrIi5vdXQiLCAidyIsIHN0ZG91dCk7CiAgICB9CgogICAgY291dCA8PCBmaXhlZCA8PCBzZXRwcmVjaXNpb24oMTIpOwoKLy8gLS0tLS0tLS0tIElOUFVUIC0tLS0tLS0tLS0tLQogICAgaW50IG47CiAgICBjaW4gPj4gbjsKCiAgICB2ZWN0b3I8dmVjdG9yPGxkPiA+IEEobiwgdmVjdG9yPGxkPihuLCAwKSk7CiAgICB2ZWN0b3I8bGQ+IEIobiwgMCk7CiAgICBsZCBDOwogICAgdmVjdG9yPGxkPiBYKG4sIDApOwogICAgaW5wdXQoQSwgQiwgQywgWCwgbik7CgoKCgoKCgoKCgoKLy8gLS0tLS0tLS0tIEFMR09SSVRITSAtLS0tLS0tLS0tCiAgICAvKioKICAgICAgICBpbnQgbGltID0gMTAwMDA7ICAgICAgICAgICAgLy8gU28gbGFuIGxhcCB0b2kgZGEKICAgICAgICBsZCB0MCA9IDE7ICAgICAgICAgICAgICAgICAgLy8gdDAgYmFuIGRhdSB0cm9uZyB0aHVhdCB0b2FuIEFybWlqbwogICAgICAgIGxkIG0gPSAwLjU7ICAgICAgICAgICAgICAgICAvLyBIZSBzbyBkaWV1IGtpZW46IG5lbiBkYXQgdHJvbmcga2hvYW5nIFswLjUgLT4gMC43XQogICAgICAgIGxkIGFscGhhID0gMC45OyAgICAgICAgICAgICAvLyBMZWFybmluZyByYXRlOiBuZW4gZGF0IHRyb25nIGtob2FuZyBbMC41IC0+IDAuOTVdCiAgICAgICAgbGQgYmV0YSA9IDFlLTY7ICAgICAgICAgICAgIC8vIENvIGJ1b2MgdG9pIHRoaWV1OiBjaGFtIHF1YSB0aGkgdGhvaQogICAgICAgIGxkIGVwc2lsb24gPSAxZS05OyAgICAgICAgICAvLyBHaWEgdHJpIHRvaSB0aGlldSBjYW4gZGF0IGR1b2MgY3VhIE5PUk0oR1JBRChkaWVtIGR1bmcpKTogY2FuZyBjaHVhbiAtPiBjaGF5IGNhbmcgY2hhbQogICAgICAgIGludCBwcmVjaXNpb25fZGlnaXQgPSAxMjsgICAvLyBTbyBjaHUgc28gc2F1IGRhdSBwaGF5IGtoaSBoaWVuIHRoaSBrZXQgcXVhCiAgICAqKi8KICAgIGludCBsaW0gPSA1MDAwMDsKICAgIHQwID0gMTsKICAgIG0gPSAwLjU7CiAgICBhbHBoYSA9IDAuODsKICAgIGJldGEgPSAxZS02OwogICAgZXBzaWxvbiA9IDFlLTU7CiAgICBwcmVjaXNpb25fZGlnaXQgPSA1OwogICAgR3JhZGllbnRfRGVzY2VudChBLCBCLCBDLCBYLCBuLCBsaW0sIEFybWlqb19zdGVwX3NpemUpOwovLyAgICAgR3JhZGllbnRfRGVzY2VudChBLCBCLCBDLCBYLCBuLCBsaW0sIEV4YWN0X2xpbmVfc2VhcmNoKTsKCgoKCgoKCgoKCgoKCi8vLS0tLS0tLS0gUGFyc2luZyBPdXRwdXQgLS0tLS0tLS0tLS0KICAgIHVwKGksIDAsIGxpbSl7CiAgICAgICAgZmluW2ldLlBCKHRvX3N0cmluZ19wcmVjcyhpKSk7CiAgICAgICAgdXAoaiwgMCwgbi0xKXsKICAgICAgICAgICAgZmluW2ldLlBCKHRvX3N0cmluZ19wcmVjcyh4W2ldW2pdKSk7CiAgICAgICAgfQogICAgICAgIGZpbltpXS5QQih0b19zdHJpbmdfcHJlY3MocmVzW2ldKSk7CiAgICAgICAgZmluW2ldLlBCKHRvX3N0cmluZ19wcmVjcyh0W2ldKSk7CiAgICAgICAgZmluW2ldLlBCKHRvX3N0cmluZ19wcmVjcyhyZXNbaV0gLSByZXNbaS0xXSkpOwogICAgfQoKICAgIGludCBtYXhfbGVuZ3RoID0gMDsKICAgIHVwKGksIDAsIGxpbSl7CiAgICAgICAgZm9yIChpbnQgaiA9IDA7IGogPCAoaW50KWZpbltpXS5zaXplKCk7IGorKyl7CiAgICAgICAgICAgIG1heF9sZW5ndGggPSBtYXgobWF4X2xlbmd0aCwgKGludClmaW5baV1bal0uc2l6ZSgpKTsKICAgICAgICB9CiAgICB9CgovLy0tLS0gUHJpbnRpbmcgdGhlIHJlc3VsdCB0YWJsZSAtLS0tCiAgICBwYXJzZV9pbmRlbnQoZmluWzBdLCBtYXhfbGVuZ3RoKTsKICAgIHByaW50X2NvbHVtbnMoZmluWzBdLCBtYXhfbGVuZ3RoKTsKICAgIHByaW50X3Jlc3VsdF90YWJsZShsaW0sIG1heF9sZW5ndGgpOwogICAgY291dCA8PCAiVG9uZyBzbyBidW9jIHQgPSBhbHBoYSp0IHRyb25nIHRodWF0IHRvYW4gQXJtaWpvOiAiIDw8IGNudDsKfQo=