Radiation

gistfile1.txt

#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <limits>
#include <stdexcept>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
#include <CGAL/Gmpz.h>

using namespace std;

typedef CGAL::Gmpz IT;
typedef CGAL::Gmpz ET;
typedef CGAL::Quadratic_program<IT> Program;
typedef CGAL::Quadratic_program_solution<ET> Solution;

struct Cell
{
    long x = 0;
    long y = 0;
    long z = 0;
};

typedef std::vector<IT> VIT;
typedef std::vector<VIT> VVIT;

VVIT pow_memo;

bool run(int D, const std::vector<Cell> &h_cell, const std::vector<Cell> &t_cell) {
    int h = h_cell.size();
    int t = t_cell.size();

    Program lp(CGAL::SMALLER, false, 0, false, 0);

    for (int i = 0; i < h; i++)
    {
        lp.set_b(i, -1);
    }
    
    for (int i = 0; i < t; i++)
    {
        lp.set_b(i + h, -1);
    }

    int var_idx = 0;

    for (int d = 0; d <= D; d++)
    {
        for (int i = d; i >= 0; i--)
        {
            int jk = d - i;
            for (int j = 0; j <= jk; j++)
            {
                int k = jk - j;

                for (int cell_idx = 0; cell_idx < t; cell_idx++)
                {
                    Cell c = t_cell[cell_idx];

                    if (pow_memo[c.x+1024][i] == -1)
                    {
                        pow_memo[c.x+1024][i] = std::pow(c.x, i);
                    }

                    if (pow_memo[c.y+1024][j] == -1)
                    {
                        pow_memo[c.y+1024][j] = std::pow(c.y, j);
                    }

                    if (pow_memo[c.z+1024][k] == -1)
                    {
                        pow_memo[c.z+1024][k] = std::pow(c.z, k);
                    }

                    IT coef = pow_memo[c.x+1024][i] * pow_memo[c.y+1024][j] * pow_memo[c.z+1024][k];

                    if (i == 0 && j == 0 && k == 0)
                    {
                        coef = 1;
                    }

                    lp.set_a(var_idx, cell_idx + h, -coef);
                }

                for (int cell_idx = 0; cell_idx < h; cell_idx++)
                {
                    Cell c = h_cell[cell_idx];
                    
                    if (pow_memo[c.x+1024][i] == -1)
                    {
                        pow_memo[c.x+1024][i] = std::pow(c.x, i);
                    }

                    if (pow_memo[c.y+1024][j] == -1)
                    {
                        pow_memo[c.y+1024][j] = std::pow(c.y, j);
                    }

                    if (pow_memo[c.z+1024][k] == -1)
                    {
                        pow_memo[c.z+1024][k] = std::pow(c.z, k);
                    }

                    IT coef = pow_memo[c.x+1024][i] * pow_memo[c.y+1024][j] * pow_memo[c.z+1024][k];

                    if (i == 0 && j == 0 && k == 0)
                    {
                        coef = 1;
                    }

                    lp.set_a(var_idx, cell_idx, coef);
                }
                
                var_idx++;
            }
        }

    }

    CGAL::Quadratic_program_options options;
    options.set_pricing_strategy(CGAL::QP_BLAND);
    Solution s = CGAL::solve_linear_program(lp, ET(), options);
    assert(s.solves_linear_program(lp));

    return s.is_infeasible();
}

void T()
{
    // The exercise description asks for a polynomial p(x,y,z) that depends on three variables (hence 3D)
    // and that has degree d, where d should be as small as possible. A polynomial of degree d is a
    // polynomial where all monomials have degree at most d. And a monomial of degree at most d is a term
    // of the form x^i y^j z^k, where i, j, k are natural numbers such that i+j+k <= d.
    int h, t;
    std::cin >> h >> t;

    std::vector<Cell> h_cell;
    h_cell.reserve(h);
    std::vector<Cell> t_cell;
    t_cell.reserve(t);

    for (int i = 0; i < h; i++)
    {
        Cell c;
        std::cin >> c.x >> c.y >> c.z;
        h_cell.push_back(c);
    }

    for (int i = 0; i < t; i++)
    {
        Cell c;
        std::cin >> c.x >> c.y >> c.z;
        t_cell.push_back(c);
    }

    pow_memo = VVIT(2049, VIT(31, -1));

    for (int i = 0; i <= 10; i++)
    {
        if (!run(i, h_cell, t_cell))
        {
            std::cout << i << "\n";
            return;
        }
    }

    int l = 10;
    int r = 30;

    int new_deg;
    bool infeasible;
    while (l <= r)
    {
        new_deg = (l + r) / 2;

        infeasible = run(new_deg, h_cell, t_cell);

        if (infeasible && new_deg == 30)
        {
            std::cout << "Impossible!\n";
            return;
        }

        if (infeasible)
        {
            l = new_deg + 1;
        }
        else
        {
            r = new_deg - 1;
        }
    }

    if (infeasible)
    {
        std::cout << new_deg + 1 << "\n";
    }
    else
    {
        std::cout << new_deg << "\n";
    }
}

int main(int argc, char const *argv[])
{
    std::ios_base::sync_with_stdio(false);
    int t;
    std::cin >> t;

    while (t--)
        T();
}