#include <bits/stdc++.h>
#define M_PI 3.14159265358979323846
using namespace std;
#define ld long double
ld Weierstrass(ld x) {
    ld y = 0;
    for (int i = 0; i < 30; i++) y += pow(0.14, i) * cos(pow(51, i) * M_PI * x);
    return y;
}
ld Classical_ternary(ld f(ld), ld l = -1, ld r = 1, ld eps = 1e-7) {
    while (r - l > eps) {
        ld midl = l + (r - l) / 3;
        ld midr = r - (r - l) / 3;
        if (f(midl) > f(midr)) r = midr;
        else l = midl;
    }
    return f(l);
}
int main() {
    cout << Classical_ternary(Weierstrass) << endl;
}

ld Blocks(ld f(ld), ld l = -1, ld r = 1, ld eps = 1e-7) {
    int iter = 10000;
    while (r - l > eps) {
        ld midl = l + (r - l) / 3;
        ld midr = r - (r - l) / 3;
        ld vall = -2e9, valr = -2e9;
        for (ld x = midl, i = 0; i < iter; x += eps, i++) vall = max(vall, f(x));
        for (ld x = midr, i = 0; i < iter; x += eps, i++) valr = max(valr, f(x));
        if (vall > valr) r = midr;
        else l = midl;
    }
    return f(l);
}

ld Blocks(ld f(ld), ld l = -1, ld r = 1, ld eps = 1e-7) {
    int iter = 10000;
    ld ans = -2e9;
    while (r - l > eps) {
        ld midl = l + (r - l) / 3;
        ld midr = r - (r - l) / 3;
        ld vall = -2e9, valr = -2e9;
        for (ld x = midl, i = 0; i < iter; x += eps, i++) vall = max(vall, f(x));
        for (ld x = midr, i = 0; i < iter; x += eps, i++) valr = max(valr, f(x));
        ans = max({ans, vall, valr});
        if (vall > valr) r = midr;
        else l = midl;
    }
    return max(ans, f(l));
}

ld find_maximum(ld f(ld), int C = 3, ld l = -1, ld r = 1, ld eps = 1e-7) {
    ld ans = -2e9;
    while (r - l > eps) {
        ld midl = l + (r - l) / C;
        ld midr = r - (r - l) / C;
        ld vall = f(midl), valr = f(midr);
        ans = max({ans, vall, valr});
        if (vall > valr) r = midr;
        else l = midl;
    }
    return max(ans, f(l));
}

Blocks res: 6.6809
Blocks time: 3.645
100-ternary res: 6.27581
100-ternary time: 0.753

Blocks res: 6.6809
Blocks time: 3.365
200-ternary res: 6.48924
200-ternary time: 2.725

Blocks res: 6.6809
Blocks time: 3.442
400-ternary res: 7.05544
400-ternary time: 10.751

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ld long double
#define fi first
#define se second
#define pb push_back
#pragma GCC optimize("unroll-loops,O3")
#pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#define cok cout << (ok ? "YES\n" : "NO\n");
#define dbg(x) cout << (#x) << ": " << x << endl;
#define dbga(x,l,r) cout << (#x) << ": "; for (int ii=l;ii<r;ii++) cout << x[ii] << " "; cout << endl;
#define int long long
#define pi pair<int, int>
const int N = 1e6+9;
int n, q;
struct pt {
    int x, y;
    pt operator-(pt b) {return {x - b.x, y - b.y};}
    int operator%(pt b) {return x * b.y - y * b.x;}
    int len() {return x * x + y * y;}
};
pt a[N];
bool cmp(pt a, pt b) {
    int sign = a % b;
    if (sign > 0) return 1;
    if (sign < 0) return 0;
    return a.len() < b.len() ? 1 : 0;
}
vector<pt> hull;
int mx[N];
ld find(int lx, int ly, int rx, int ry, int px, int py) {
    ld k = (ld)py / px;
    return (lx * ry - ly * rx) / (k * (lx - rx) + ry - ly);
}
ld findgood(ld lx, ld ly, ld rx, ld ry, ld px, ld py) {
    ld l = 0, r = 1e9;
    for (int o = 0; o < 200; o++) {
        ld mid = (l + r) / 2;
        if ((lx - px / mid) * (ry - py / mid) - (ly - py / mid) * (rx - px / mid) >= 0) r = mid;
        else l = mid;
    }
    return px / l;
}
signed main() {
    cin.tie(0); ios_base::sync_with_stdio(0);
    cout.precision(20);
    cin >> n;
    int longest = 0;
    for (int i = 0; i < n; i++) cin >> a[i].x >> a[i].y, longest = max(longest, a[i].y);
    sort(a, a + n, cmp);
    hull.pb({0, 0});
    hull.pb(a[0]);
    for (int i = 1; i < n; i++) {
        while (hull.size() >= 2 && (a[i] - hull[hull.size() - 2]) % (hull.back() - hull[hull.size() - 2]) >= 0) hull.pop_back();
        hull.pb(a[i]);
    }
    mx[hull.size()] = 0;
    for (int i = hull.size() - 1; i > 0; i--) {
        mx[i] = max(mx[i + 1], hull[i].y);
    }
    cin >> q;
    for (int o = 0; o < q; o++) {
        pt P;
        cin >> P.x >> P.y;
        if (P % hull[1] >= 0 && 1) {
            cout << (long double)P.y / longest << "\n";
            continue;
        }
        if (P % hull.back() < 0 && 1) {
            cout << "-1\n";
            continue;
        }
        int l = 1, r = hull.size() - 1;
        while (r - l > 1) {
            int mid = l + r >> 1;
            pt current = hull[mid];
            if (current % P >= 0) l = mid;
            else r = mid;
        }
        ld X = find(hull[l].x, hull[l].y, hull[r].x, hull[r].y, P.x, P.y);
        cout << (long double)min(P.x / X, P.y / (ld)mx[r]) << "\n";
    }
}

В тетраэдре проведены внешние векторы нормалей ко всем граням. Длина каждого вектора численно равна площади соответствующей грани. Доказать, что сумма векторов равна 0

Закачаем в тетраэдр газ давления p. Тогда сила, с которой газ давит на грань, перпендикулярна грани и пропорциональна ее площади: F = pS. Из условия равновесия, сумма сил равна 0, и тогда сумма наших веркторов тоже равна 0. Это доказательство легко обобщается на произвольный многогранник.

Общеизвестно

Ясно, что прямоугольный треугольник определяется гипотенузой c и углом при вершине α. Тогда, из соображений размерности, его площадь S равна f(α) * c^2, где f(x) — какая-то функция от угла (на самом деле, f(x) = 0.5 * sin(x)). Проведём высоту к гипотенузе, разделив треугольник на два маленьких треугольника с площадями S1 и S2. В обоих треугольниках есть угол α. Суммируя их площади по формуле, сокращаем на f(α) и получаем требуемое равенство.

Внутри остроугольного треугольника ABC отмечена точка P, минимизирующая сумму расстояний PA+PB+PC. Доказать, что стороны треугольника видны из точки P под равными углами.

Скрепим три нити в одной точке (точке P), которая может свободно скользить по плоскости. Сделаем отверстия в вершинах треугольника и пропустим через них нити. Подвесим на нижние концы нитей грузы равной массы. Минимум потенциальной энергии соответствует максимальной суммарной длине вертикальных отрезков нитей, то есть минимальной суммарной длине отрезков PA, PB, PC. Точка P находится в покое, и на неё действуют три одинаковые силы. Тогда из симметрии углы между каждой парой нитей равны по 120 градусов.

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ld long double
#define fi first
#define se second
#define pb push_back
#define cok cout << (ok ? "YES\n" : "NO\n");
#define dbg(x) cout << (#x) << ": " << (x) << endl;
#define dbga(x,l,r) cout << (#x) << ": "; for (int ii=l;ii<r;ii++) cout << x[ii] << " "; cout << endl;
// #define int long long
#define pi pair<int, int>
const int N = 7, C = 1e7, MAX_DEG = 4, MAX_PRODUCT = 1e5;
const ld EPS = 1e-9, EPS_CHECK = 1e-9;
const string SEP = "  (", END = ")\n";
const bool APPROXIMATION = true;
array <string, N> names;
array <int, N> max_exp, powers, current_converted, cur_exp;
array<vector<ll>, N> POINTS;
ll DIV[N][MAX_DEG + 1][MAX_DEG + 1], PW[N][MAX_DEG + 1][MAX_DEG + 1];
ld SUM[MAX_PRODUCT];
ld F_CACHE[MAX_PRODUCT];
ll pow(ll a, int b)
{
	if (b == 0) return 1;
	if (b == 1) return a;
	ll s = pow(a, b / 2);
	s *= s;
	if (b & 1) s *= a;
	return s;
}
ld approximate(ld k)
{
	int k_ = k;
	int k__ = k_ + abs(k) / k;
	if (abs(k - k_) < EPS) return k_;
	else if (abs(k - k__) < EPS) return k__;
	else
	{
		int i = 1, j = 1;
		ld ka = abs(k);
		while (i < C && j < C)
		{
			ld p = ka * j;
			if (abs(p - i) < EPS) break;
			if (p < i) j++;
			else i++;
		}
		if (i >= C || j >= C) return k;
		if (k < 0) i = -i;
		return (ld)i / j;
	}	
}
void normalize(ld k)
{
    if (!APPROXIMATION)
    {
        cout << k << SEP;
        return;
    }
	int k_ = k;
	int k__ = k_ + abs(k) / k;
	if (abs(k - k_) < EPS) cout << k_ << SEP;
	else if (abs(k - k__) < EPS) cout << k__ << SEP;
	else
	{
		int i = 1, j = 1;
		ld ka = abs(k);
		while (i < C && j < C)
		{
			ld p = ka * j;
			if (abs(p - i) < EPS) break;
			if (p < i) j++;
			else i++;
		}
		if (i >= C || j >= C)
		{
			cout << k << SEP;
			return;
		}
		if (k < 0) i = -i;
		cout << i << "/" << j << SEP;
	}
}
struct monom
{
	array<int, N> exp;
	ld k;
	int deg;
	monom(array<int, N> v, ld k_)
	{
		k = k_;
		exp = v;
		deg = 0;
		for (int i=0;i<N;i++) deg += exp[i];
	}
	void display()
	{
		normalize(k);
		if (deg == 0) { cout << "1" << END; return;}
		bool go = 0;
		for (int i=0;i<N;i++)
		{
			if (go && exp[i]) cout << " * ";
			if (exp[i]) go = 1, cout << names[i] + "^" + to_string(exp[i]);
		}
		cout << END;
	}
	ld operator()(array<int, N> v)
	{
		ll res = 1;
		for (int i=0;i<N;i++) res *= PW[i][v[i]][exp[i]];
		return k * res;
	}
	ld getRandom(array<ll, N> v)
	{
		ld res = 1;
		for (int i=0;i<N;i++) res *= pow(v[i], exp[i]);
		return k * res;
	}
};
bool operator<(monom a, monom b)
{
	if (a.deg > b.deg) return 1;
	if (a.deg < b.deg) return 0;
	if (a.exp > b.exp) return 1;
	if (a.exp < b.exp) return 0;
	return a.k > b.k;
}
struct polynom
{
	vector<monom> st;
	void add(monom m)
	{
		if (abs(m.k) < EPS) return;
		st.pb(m);
	}
	void print() { if(st.size() == 0) {cout << "Polynom is 0\n"; return;} sort(st.begin(), st.end()); for (monom &m: st) m.display();}
	ld operator()(array<ll, N> v)
	{
		ld res = 0;
		for (auto &m: st) res += m.getRandom(v);
		return res;
	}
};
ld gen(int index=0, int current_hash=0)
{
	if (index == N)
	{
		ll div = 1;
		for (int i=0;i<N;i++) div *= DIV[i][current_converted[i]][cur_exp[i]];
		return (ld)(F_CACHE[current_hash] - SUM[current_hash]) / div;
	}
	ld res = 0;
	for (int i=0;i<=cur_exp[index];i++)
	{
		current_converted[index] = i;
		res += gen(index + 1, current_hash + i * powers[index]);
	}
	return res;
}
array<int, N> convert(int h)
{
	array<int, N> res;
	for (int i=0;i<N;i++) res[i] = h / powers[i], h -= res[i] * powers[i];
	return res;
}
array<ll, N> convert_points(int h)
{
	array<ll, N> res;
	for (int i=0;i<N;i++) res[i] = POINTS[i][h / powers[i]], h %= powers[i];
	return res;
}
polynom interpolate(ld f(array<ll, N>))
{
    int max_pow = -2e9, sum = 0, h_max = 0;
    set<int> remaining_points, st;
	polynom res;
    for (int x: max_exp) max_pow = max(max_pow, x), sum += x, h_max = h_max * (x + 1) + x;

    powers[N - 1] = 1;
    for (int i=N-2;i>-1;i--) powers[i] = powers[i + 1] * (max_exp[i + 1] + 1);

    for (int i=0;i<max_exp.size();i++) for (int j=0;j<=max_exp[i];j++) POINTS[i].pb(j);
	
    for (int i=0;i<N;i++) for (int j=0;j<=max_exp[i];j++) for (int u=0;u<=max_exp[i];u++) DIV[i][j][u] = (u ? DIV[i][j][u - 1] : 1) * (u == j ? 1 : (POINTS[i][j] - POINTS[i][u]));

    for (int i=0;i<N;i++) for (int j=0;j<=max_exp[i];j++) for (int u=0;u<=max_pow;u++) PW[i][j][u] = u ? PW[i][j][u - 1] * POINTS[i][j] : 1;

    for (int i=0;i<=h_max;i++) F_CACHE[i] = f(convert_points(i)), remaining_points.insert(i);
    st.insert(h_max);
	
    while (st.size())
	{
		int v = *st.rbegin();
		st.erase(v);
		remaining_points.erase(v);
		cur_exp = convert(v);
		ld k = gen();
		if (abs(k) > EPS)
		{
			if (APPROXIMATION) k = approximate(k);
			monom mn = monom(cur_exp, k);
			res.add(mn);
			for (int i: remaining_points) SUM[i] += mn(convert(i));
		}
		for (int i=0;i<N;i++) if (cur_exp[i]) st.insert(v - powers[i]);
	}
	return res;
}
ld f(array<ll, N> v)
{
	auto [a, b, c, d, e, f, g] = v;
	ld res = 0;
	for (int i=0;i<a;i++)
		for (int j=0;j<b;j++)
			for (int u=0;u<c;u++)
				for (int x=0;x<d;x++)
					for (int y=0;y<e;y++)
						for (int z=0;z<f;z++)
							for (int k=0;k<g;k++)
								res += 13ll * i * j * u * i * i * u - 49ll * k * k * z * z * y + 90ll * c * u * k * x * x * x;
	return res;
}
void check(polynom p, ld(array<ll, N> f))
{
	mt19937 rnd(228);
	for (int i=0;i<10000;i++)
	{
		int t = clock();
		array<ll, N> ex;
		for (int j=0;j<N;j++) ex[j] = rnd() % 20 + 2;
		ld F = f(ex);
		ld P = p(ex);
		if (abs(F - P) > max(EPS_CHECK, EPS_CHECK * abs(F)))
		{
			cout << "Polynom is wrong, test " << i << endl;
			cout << F << endl << P << endl;
			for (int x: ex) cout << x << " ";
			cout << endl;
			return;
		}
		cout << "Test " << i << " has been passed, time = " << (ld)(clock() - t) / CLOCKS_PER_SEC << "s" << endl;
	}
	cout << "Polynom is OK" << endl;
}
signed main()
{
    cin.tie(0); ios_base::sync_with_stdio(0);
    cout << setprecision(20) << fixed;

    names = {"a", "b", "c", "d", "e", "f", "g"};
    max_exp = {4, 2, 3, 4, 2, 3, 3};
    
    polynom P = interpolate(f);
    P.print();
    //cout << "Checking polynom..." << endl;
    //check(P, f);
}

int stupid(int a, int b)
{
    int ans = 0;
    for (int x1=0;x1<a;x1++)
        for (int x2=0;x2<a;x2++)
            for (int x3=0;x3<a;x3++)
                for (int y1=0;y1<b;y1++)
                    for (int y2=0;y2<b;y2++)
                        for (int y3=0;y3<b;y3++)
                        {
                            int a = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
                            int b = (x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3);
                            int c = (x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2);
                            if ((a + b == c) && min(a, min(b, c))) ans++;
                        } 
    return ans / 2;
}
signed main()
{
    int C1 = 5, C2 = 5;
    ld res = 0;
    for (int a1=4;a1<=C1+3;a1++)
    {
    	for (int a2=4;a2<=C2+3;a2++)
    	{
    		int d = 1;
    		for (int ai=4;ai<=C1+3;ai++) if (ai != a1) d *= a1 - ai;
    		for (int ai=4;ai<=C2+3;ai++) if (ai != a2) d *= a2 - ai;
    		res += (ld)(stupid(a1, a2)) / d;
    	}
    }
    cout << setprecision(20) << fixed << res;
}

int fast(int a, int b)
{
    int ans = 0;
    for (int a1=-a+1;a1<=a-1;a1++)
    {
	for (int b1=-b+1;b1<=b-1;b1++)
	{
       	    if (a1 == 0 && b1 == 0) continue; 
       	    for (int a2=-a+1;a2<=a-1;a2++)
	    {
	        for (int b2=-b+1;b2<=b-1;b2++)
		{
		    if (a2 == 0 && b2 == 0) continue;

		    // первое условие - скалярное произведение равно 0
		    // второе условие - векторное произведение < 0 - чтобы не считать одну пару дважды
		    if (b1 * b2 + a1 * a2 == 0 && a1 * b2 - b1 * a2 < 0)
		    {
		        int cnta = a - max(max(abs(a1), abs(a2)), abs(a1 - a2));
		        int cntb = b - max(max(abs(b1), abs(b2)), abs(b1 - b2));
		        if (cnta > 0 && cntb > 0) ans += cnta * cntb;
		    }
	        }
	    }
        }
    }
    return ans;
}

int fast(int a, int b)
{
    int ans = 0;
    for (int a1=-a+1;a1<=a-1;a1++)
	for (int b1=-b+1;b1<=b-1;b1++)
	{
       	    if (a1 == 0 && b1 == 0) continue; 
	    for (int a2=-a+1;a2<=a-1;a2++)
	        for (int b2=-b+1;b2<=b-1;b2++)
		{
		    if (a2 == 0 && b2 == 0) continue;
		    if (b1 * b2 + a1 * a2 == 0 && a1 * b2 - b1 * a2 < 0)
		    {
		        int cnta = a - max(max(abs(a1), abs(a2)), abs(a1 - a2));
			int cntb = b - max(max(abs(b1), abs(b2)), abs(b1 - b2));
			if (cnta > 0 && cntb > 0) ans += cnta * cntb;
		    }
	        }
	}
    return ans;
}
struct polynom
{
    int x0 = 0, x1 = 0, x2 = 0;
    void print()
    {
	cout << x2 << "x^2 + " << x1 << "x - " << -x0 << endl;
    }
    ld get(ld x)
    {
	return x2 * x * x + x1 * x + x0;
    }
};
void interpolate(int b)
{
    polynom P = {0, 0, 0};
    vector<int> A = {b * b + 1, b * b + 2, b * b + 3};
    ld res = 0;
    vector<ld> B = {fast(A[0], b), fast(A[1], b), fast(A[2], b)};
    for (int i=0;i<3;i++)
    {
	int d = 1;
	for (int j=0;j<3;j++) if (A[i] != A[j]) d *= A[i] - A[j];
	res += (B[i] - P.get(A[i])) / d;
    }
    P.x2 = res;
    res = 0;
    A.pop_back();
    for (int i=0;i<2;i++)
    {
	int d = 1;
	for (int j=0;j<2;j++) if (A[i] != A[j]) d *= A[i] - A[j];
	res += (B[i] - P.get(A[i])) / d;
    }
    P.x1 = res;
    P.x0 = B[0] - P.get(A[0]);
    P.print();
}

Выберем одну вертикальную прямую x = z, z < a - 1 (a - 1 способ), вертикальную прямую x = a - 1 и две горизонтальные прямые (b * (b - 1) / 2 способов). Пересечением наших прямых являются 4 точки. Мы можем выкинуть из них любую, и после этого оставшиеся 3 точки будут образовывать прямоугольный треугольник второго типа. Ясно, что любой треугольник второго типа может быть представлен в виде пересечения таких прямых и отбрасывания одной из точек пересечения.

Лемма. Максимальный модуль разности абсцисс вершин треугольника третьего типа не превосходит (b - 1) ^ 2 + 1

Для доказательства посмотрим на рисунок.

У нашего треугольника есть хотя бы две стороны, непараллельные осям координат. Пусть они AB и BC. Пусть угловой коэффициент прямой AB равен k, тогда у прямой BC он равен -1/k. Тогда разность абсцисс верши A и C не больше, чем (b - 1) / |k| + (b - 1) / |-1/k| = (b - 1) * (|k| + 1 / |k|).

Известно, что если мы возьмём два положительных числа и будем их сближать на числовой прямой с сохранением произведения, то их сумма будет уменьшаться. Тогда для максимизации |k| + 1 / |k| нужно минимизировать или максимизировать |k|. Из ограничения по y-координатам, |k| < b - 1 и 1 / |k| < b - 1. Тогда максимальное значение нашего выражения (b - 1) * ((b - 1) + 1 / (b - 1)) = (b - 1) ^ 2 + 1. Лемма доказана.

Нужно доказать, что при a > (b - 1) ^ 2 + 1 количество треугольников третьего типа не зависит от a. Это очевидно следует из леммы, ведь каждому такому треугольнику для a мы можем сопоставить треугольник для a + 1, полученный из первого сдвигом на 1 вправо. И каждому треугольнику 3 типа для a + 1 аналогично сопоставим треугольник 3 типа для a, сдвинув его на 1 влево.

ll fast(int a, int b)
{
    ll ans = 0;
    for (int a1=-a+1;a1<=a-1;a1++)
	for (int b1=-b+1;b1<=b-1;b1++)
	{
       	    if (a1 == 0 && b1 == 0) continue; 
	    for (int a2=-a+1;a2<=a-1;a2++)
	        for (int b2=-b+1;b2<=b-1;b2++)
		{
		    if (a2 == 0 && b2 == 0) continue;
		    if (b1 * b2 + a1 * a2 == 0 && a1 * b2 - b1 * a2 < 0)
		    {
		        int cnta = a - max(max(abs(a1), abs(a2)), abs(a1 - a2));
			int cntb = b - max(max(abs(b1), abs(b2)), abs(b1 - b2));
			if (cnta > 0 && cntb > 0) ans += cnta * cntb;
		    }
	        }
	}
    return ans;
}
struct polynom
{
    ll x0 = 0, x1 = 0, x2 = 0;
    void print()
    {
	cout << x2 << "x^2 + " << x1 << "x - " << -x0 << endl;
    }
    ld get(ld x)
    {
	return x2 * x * x + x1 * x + x0;
    }
};
ll get(int a, int b)
{
    if (a < b) swap(a, b);
    if (a <= (b - 1) * (b - 1) + 1) return fast(a, b);
    int x = fast((b - 1) * (b - 1) + 1, b);
    int y = fast((b - 1) * (b - 1) + 2, b);
    // y = x + c + (a - 1) * b * (b - 1) / 2 * 4;
    int c = y - x - ((b - 1) * (b - 1) + 1) * (b - 1) * b * 2;
    polynom P;
    P.x2 = b * b - b;
    P.x1 = c - b * b + b;
    P.x0 =-b*b*b*b*b*b + 5*b*b*b*b*b - 11*b*b*b*b + 13*b*b*b + b*b*(-8 - c) + b * (2 * c + 2) - 2 * c + x;
    return P.get(a);
}