# SIMARRAY - Editorial

Author: Sundar
Tester: Riley Borgard

MEDIUM

# PROBLEM:

You are given two arrays of positive integers, A and B, both of length N. Find the minimum possible value of \displaystyle \sum_{i=1}^{N} (A_i - x_i B_i)^2, over all possible N-tuples of real numbers (x_1, x_2, \ldots x_N) satisfying x_1 \ge x_2 \ge \cdots \ge x_N

# EXPLANATION:

For each 1 \le L \le R \le N, define F_{[L,R]}(x) = \sum_{i=L}^{R} (A_i - B_ix)^2 = Q_{[L,R]} x^2 - 2 P_{[L,R]} x + C_{[L,R]}, where Q_{[L,R]} is the sum of B_i^2, P_{[L, R]} is the sum of A_iB_i and C_{[L,R]} is the sum of A_i^2 in the range [L,R] respectively.

It is minimized at t_{[L,R]} = \dfrac{P_{[L,R]}}{Q_{[L,R]}}. For all x < t_{[L,R]}, it is decreasing and for all x > t_{[L,R]}, it is increasing. Let the minimum value be V_{[L,R]}

Define x_0 = \infty, x_{N+1} = -\infty. Lets prove the following statement:

Lemma 1 : Lets partition the optimal solution into maximal segments of equal values. Consider any such segment [L,R]. Then, x_i = t_{[L,R]} for all i in [L,R].

Proof : Assume that all x_i for i in range [L,R] were equal to u \ne t_{[L,R]}. Let r = \min(|u - t_{[L,R]}|, \min(x_{L-1} - u, u - x_{R+1})). Clearly, r > 0. If u < t_{[L,R]} we can change all values to u+r, else we can change all values to u-r to improve the objective, while maintaining the non increasing property.

So, any optimal solution can be thought of as segments [1=L_1,R_1], [L_2=R_1+1,R_2] \ldots [L_{k} = R_{k-1}+1,R_k = N] with the segment [L_i,R_i] having a value x_i = t_{[L_i,R_i]}, such that x_1\geq x_2 \ldots \geq x_k.

There are exactly 2^{N-1} ways of partitioning the array into segments (for each index i > 1, choose whether a new segment starts at i). So, we can brute-force over these partitions and find the minimum objective over all the valid partitions. This passes subtasks 1 and 2.

For subtask 3, let’s define dp[L][R] to be the minimum value of the sums of deviation over first R indices, if the last segment was [L,R]. Then, let’s iterate over the left end L' of the second last segment, such that t_{[L',L-1]} \geq t_{[L,R]}. Over all such L', we have to minimize dp[L'][L-1] + V_{[L,R]}. The total time complexity would be O(N^3).

We need to minimize dp[L'][L-1] over all L' with t_{[L',L-1]} \geq t_{[L,R]}. Let’s sort the intervals [L',L-1] according to t_{[L,L-1]}. Then, we need the minimum of a suffix. The suffix minimums can be precomputed. We don’t have to sort these segments ending at L-1 for every R. We can do this once we have processed all segments with right end L-1. The time complexity would be O(N^2 \log(N))

Let’s consider a more general problem. Say we have m quadratic functions H_1(x), H_2(x), \ldots H_m(x) each having a positive coefficient of x^2. Let the respective minimizing points be v_1, v_2, \ldots, v_m. We want to minimize \sum H_i (x_i) over all tuples x with x_1 \geq x_2 \ldots \geq x_m .

Lemma 2: Consider some i with v_i < v_{i+1}. Then, in the optimal solution, x_i = x_{i+1}

Proof : Assume that in the optimal solution x_i > x_{i+1}. If x_i > v_i, we can decrement x_i by a small enough value \epsilon improving the objective without violating the non increasing constraint. Similarly, if x_{i+1} < v_{i+1}, we can increment x_{i+1} by a small enough \epsilon to improve the objective. Therefore, x_i \le v_i and x_{i+1} \geq v_{i+1}. But then, v_{i+1} \le x_{i+1} < x_i \leq v_i, a contradiction.

Using this, we can get a simple O(N^2) algorithm. If there exists an index with v_i < v_{i+1}, remove H_i(x), H_{i+1}(x) from the sequence and add H_i(x) + H_{i+1}(x).

Once v_i \geq v_{i+1} for all i, the solution x_i = v_i for all i is valid and clearly optimal.

To improve the above algorithm, we can just keep a stack of quadratic functions. For each i, we first add H_i(x) to the stack and then keep on merging the top two elements of the stack while their minimizing points are out of order. This leads to a simple O(N) algorithm.

# SOLUTIONS:

Setter's Solution
#include <bits/stdc++.h>
using namespace std;

#define ll long long
#define pii pair<int, int>
#define all(c) ((c).begin()), ((c).end())
#define sz(x) ((int)(x).size())

#ifdef LOCAL
#include <print.h>
#else
#define trace(...)
#endif

long long readInt(long long l,long long r,char endd){
long long x=0;
int cnt=0;
int fi=-1;
bool is_neg=false;
while(true){
char g=getchar();
if(g=='-'){
assert(fi==-1);
is_neg=true;
continue;
}
if('0'<=g && g<='9'){
x*=10;
x+=g-'0';
if(cnt==0){
fi=g-'0';
}
cnt++;
assert(fi!=0 || cnt==1);
assert(fi!=0 || is_neg==false);

assert(!(cnt>19 || ( cnt==19 && fi>1) ));
} else if(g==endd){
if(is_neg){
x= -x;
}
if(!(l<=x && x<=r))cerr<<l<<"<="<<x<<"<="<<r<<endl;
assert(l<=x && x<=r);
return x;
} else {
assert(false);
}
}
}
string readString(int l,int r,char endd, char minc = 'a', char maxc = 'z'){
string ret="";
int cnt=0;
while(true){
char g=getchar();
assert(g!=-1);
if(g==endd){
break;
}
assert(g >= minc && g <= maxc);
cnt++;
ret+=g;
}
assert(l<=cnt && cnt<=r);
return ret;
}
long long readIntSp(long long l,long long r){
}
long long readIntLn(long long l,long long r){
}
string readStringLn(int l,int r, char minc = 'a', char maxc = 'z'){
}
string readStringSp(int l,int r, char minc = 'a', char maxc = 'z'){
}
template<class T>
vector<T> readVector(int n, long long l, long long r){
vector<T> ret(n);
for(int i = 0; i < n; i++){
ret[i] = i == n - 1 ? readIntLn(l, r) : readIntSp(l, r);
}
return ret;
}

struct Rational
{
ll num, den;
int cnt;
Rational(ll n, ll d, int cnt) : num(n), den(d), cnt(cnt){}
bool operator<(const Rational &r) const{
return num * __int128_t(r.den) < den * __int128_t(r.num);
};
Rational operator + (const Rational & r) const{
return Rational(num + r.num, den + r.den, cnt + r.cnt);
}
};

const int MAXN = 500000;
int main(){
int sn = 0;
while(t--){
sn += n;
assert(sn <= MAXN);
double ans = 0;

for(ll & x : a){
ans += x * x;
}

vector<Rational> stk;

for(int i = 0; i < n; i++){
a[i] *= b[i];
b[i] *= b[i];
// b * x^2 - 2 * a * x
stk.push_back(Rational(a[i], b[i], 1));
while(stk.size() >= 2){
Rational x = stk.back(); stk.pop_back();
Rational y = stk.back(); stk.pop_back();
if(y < x){
stk.push_back(x + y);
} else{
stk.push_back(y);
stk.push_back(x);
break;
}
}
}
int j = 0;
for(int i = 0; i < stk.size(); i++){
double x = stk[i].num / (double) stk[i].den;
for(int k = 0; k < stk[i].cnt; k++, j++){
ans += x * (b[j] * x - 2 * a[j]);
}
}
trace(stk.size());
cout << setprecision(20) << fixed << ans << endl;
}
}

Tester's Solution
#include <bits/stdc++.h>

#define ll long long
#define sz(x) ((int) (x).size())
#define all(x) (x).begin(), (x).end()
#define vi vector<int>
#define pii pair<int, int>
#define rep(i, a, b) for(int i = (a); i < (b); i++)
using namespace std;
template<typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;

using ftype = long double;

struct node {
ftype ab, bb;
int n;
node(ftype a, ftype b) {
n = 1;
ab = a * b;
bb = b * b;
}
node(ftype ab, ftype bb, int n) : ab(ab), bb(bb), n(n) {}
node operator+(const node &o) const {
return node(ab + o.ab, bb + o.bb, n + o.n);
}
ftype getx() const {
return ab / bb;
}
};

void solve() {
int n;
cin >> n;
vector<ftype> a(n), b(n);
rep(i, 0, n) cin >> a[i];
rep(i, 0, n) cin >> b[i];
vector<node> st;
rep(i, 0, n) {
node x(a[i], b[i]);
while(!st.empty() && st.back().getx() < x.getx()) {
x = x + st.back();
st.pop_back();
}
st.push_back(x);
}
vector<ftype> r;
for(node &x : st) {
ftype val = x.getx();
rep(i, 0, x.n) r.push_back(val);
}
ftype ans = 0;
rep(i, 0, n) ans += (a[i] - r[i] * b[i]) * (a[i] - r[i] * b[i]);
cout << ans << '\n';
}

int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout << fixed << setprecision(10);
int te;
cin >> te;
while(te--) solve();
}