QOJ.ac

QOJ

ID题目提交者结果用时内存语言文件大小提交时间测评时间
#110530#6542. Optimal Quadratic Functionhjroh0315RE 0ms3824kbC++2030.1kb2023-06-02 18:29:102023-06-02 18:29:29

Judging History

你现在查看的是最新测评结果

  • [2023-08-10 23:21:45]
  • System Update: QOJ starts to keep a history of the judgings of all the submissions.
  • [2023-06-02 18:29:29]
  • 评测
  • 测评结果:RE
  • 用时:0ms
  • 内存:3824kb
  • [2023-06-02 18:29:10]
  • 提交

answer

#include<bits/stdc++.h>
using namespace std;

#define assert_msg(x,y) assert((x && y))

class Rng{
private:
    static std::mt19937 engine;

public:
    static std::mt19937& get_engine(){
        return engine;
    }
    template<typename T>
    static void set_seed(T const& seed){
        engine = std::mt19937(seed);
    }
    static void timebased_seed(){
        engine = std::mt19937(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());
    }
    template<typename T>
    static typename std::enable_if<std::is_integral<T>::value, T>::type uniform(T l, T r){
        return std::uniform_int_distribution<T>(l, r)(engine);
    }
    template<typename T>
    static typename std::enable_if<std::is_floating_point<T>::value, T>::type uniform(T l, T r){
        return std::uniform_real_distribution<T>(l, r)(engine);
    }

};
std::mt19937 Rng::engine(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now().time_since_epoch()).count());


class Timer{
public:
    template<typename S, typename T, typename... U>
    static S execute_timed(T func, std::string const& name, U&&... u){
        auto time_begin = std::chrono::high_resolution_clock::now();
        S ret = func(std::forward<U>(u)...);
        auto time_end = std::chrono::high_resolution_clock::now();
        auto timespan = std::chrono::duration_cast<std::chrono::nanoseconds>(time_end - time_begin);
        std::cerr << "Execution of '" << name << "' took " << std::fixed << std::setprecision(6) << timespan.count()*1e-9 << "\n";
        return ret;
    }
};

/**
 *  signed Bigint in base 10^18, used for Input / Output
 *  don't use for computations, doesn't support most operations
 *
 */
class Bigint_base10{
public:
    static constexpr int64_t BASE = 1e18;
    static constexpr int DIGITS = 18;
private:
    bool is_neg;
    vector<int64_t> data;

public:
    Bigint_base10():is_neg(false), data(1, 0){}
    Bigint_base10(int64_t const&val):is_neg(val<0){
        int64_t abs_val = abs(val);
        if(abs_val < BASE){
            data = {abs_val};
        } else {
            data = {abs_val%BASE, abs_val/BASE};
        }
    }

    Bigint_base10 operator+(Bigint_base10 const&o)const{
        assert_msg(is_neg == o.is_neg, "Addition operands need to have equal sign");
        Bigint_base10 ret;
        ret.is_neg = is_neg;
        ret.data.assign(1+max(data.size(), o.data.size()), 0);
        copy(data.begin(), data.end(), ret.data.begin());
        int64_t carry = 0;
        for(unsigned int i=0;i<o.data.size();++i){
            ret.data[i]+=o.data[i] + carry;
            carry = 0;
            if(ret.data[i] >= BASE){
                carry = 1;
                ret.data[i]-=BASE;
            }
        }
        for(unsigned int i=o.data.size();carry;++i){
            ret.data[i]+=carry;
            carry = 0;
            if(ret.data[i] >= BASE){
                carry = 1;
                ret.data[i]-=BASE;
            }
        }
        return ret.trim();
    }

    Bigint_base10 operator*(int64_t const&o)const{
        if(o == 0){
            return Bigint_base10(0);
        }
        if(o<0){
            return operator*(-o).negate();
        }
        if(o&1){
            return operator+(operator*(o-1));
        }
        return operator+(*this)*(o/2);
    }
    Bigint_base10& operator+=(Bigint_base10 const&o){
        *this = operator+(o);
        return *this;
    }
    Bigint_base10& operator*=(int64_t const&o){
        *this = operator*(o);
        return *this;
    }

    Bigint_base10& trim(){
        while(data.size()>1 && data.back() == 0){
            data.pop_back();
        }
        return *this;
    }

    bool is_zero()const{
        for(auto const&e:data) if(e) return false;
        return true;
    }

    Bigint_base10& negate(){
        is_neg = !is_neg;
        if(is_zero()) is_neg = false;
        return *this;
    }

    friend ostream& operator<<(ostream&o, Bigint_base10 const&b){
        if(b.is_neg) o << '-';
        o << b.data.back();
        o << setfill('0');
        for(auto it = next(b.data.rbegin());it != b.data.rend();++it){
            o << setw(9) << *it;
        }
        o << setw(0) << setfill(' ');
        return o;
    }
    friend istream& operator>>(istream&in, Bigint_base10 &b){
        static string tmp;
        in >> tmp;
        assert_msg(in, "input should be readable as a string");
        if(tmp[0] == '-'){
            b.is_neg = true;
            tmp = tmp.substr(1, -1);
        } else {
            b.is_neg = false;
        }
        assert_msg(all_of(tmp.begin(), tmp.end(), [](char const&c){return '0'<=c && c<='9';}), "Input should consist of digits and possibly a '-'");
        assert_msg(!tmp.empty(), "Input should contain at least one digit");

        b.data.resize((tmp.size()+DIGITS-1)/DIGITS);
        unsigned int i, j;
        for(i=tmp.size()-DIGITS, j=0;i>0;i-=DIGITS, ++j){
            b.data[j] = stoll(tmp.substr(i, DIGITS));
        }
        b.data[j] = stoll(tmp.substr(0, i+DIGITS));
        return in;
    }
};


/**
 *  Biginteger with fixed precision
 *  Has 31*len bits, first bit is sign
 *
 *  Is quite fast
 */
template<size_t len>
struct Bigint_Fixedsize{
    unsigned int data[len];
    static constexpr unsigned int bits = 31;
    Bigint_Fixedsize(){memset(data, 0, sizeof(data));}
    Bigint_Fixedsize(long long const&_a){
        memset(data, 0, sizeof(data));
        unsigned long long a = _a;
        data[0] = a&((1u<<bits)-1);
        data[1] = a>>bits;
        data[1]&=~(1u<<bits);
        if(a>~a){ // negative number, use complement
            for(size_t i=2;i<len;++i){
                data[i] = (1u<<bits)-1;
            }
        }
    }
    //__attribute__((optimize("unroll-loops")))
    Bigint_Fixedsize& operator+=(Bigint_Fixedsize const&o){
        unsigned int carry = 0;
        for(size_t i=0;i<len;++i){
            data[i]+=o.data[i] + carry;
            carry = data[i]>>bits;
            data[i]&=~(1u<<bits);
        }
        return *this;
    }
    //__attribute__((optimize("unroll-loops")))
    Bigint_Fixedsize& operator-=(Bigint_Fixedsize const&o){
        unsigned int carry = 0;
        for(size_t i=0;i<len;++i){
            data[i]-=o.data[i] + carry;
            carry = data[i]>>bits;
            data[i]&=~(1u<<bits);
        }
        return *this;
    }
    Bigint_Fixedsize operator+(Bigint_Fixedsize const&o)const{
        Bigint_Fixedsize ret(*this);
        ret+=o;
        return ret;
    }
    Bigint_Fixedsize operator-(Bigint_Fixedsize const&o)const{
        Bigint_Fixedsize ret(*this);
        ret-=o;
        return ret;
    }
    //__attribute__((optimize("unroll-loops")))
    void multo(Bigint_Fixedsize const&o, Bigint_Fixedsize &ret)const{
        static unsigned int tmp[len+1];
        memset(tmp, 0, sizeof(tmp));
        for(size_t i=0;i<len;++i){
            unsigned long long val = 0;
            for(size_t j=0;j<len-i;++j){
                val+= data[i]*(unsigned long long)o.data[j]+tmp[i+j];
                tmp[i+j] = val&((1u<<bits)-1);
                val>>=bits;
            }
        }
        memcpy(ret.data, tmp, sizeof(ret.data));
    }
    Bigint_Fixedsize& operator*=(Bigint_Fixedsize const&o){
        multo(o, *this);
        return *this;
    }
    Bigint_Fixedsize operator*(Bigint_Fixedsize const&o)const{
        Bigint_Fixedsize ret;
        multo(o, ret);
        return ret;
    }
    Bigint_Fixedsize& negate(){
        unsigned int carry = 0;
        for(size_t i=0;i<len;++i){
            data[i] = ~data[i] + !carry;
            carry = (data[i]>>bits);
            data[i]&=~(1u<<bits);
        }
        return *this;
    }

    Bigint_Fixedsize operator-()const{
        Bigint_Fixedsize ret(*this);
        ret.negate();
        return ret;
    }
    bool operator<(Bigint_Fixedsize const&o)const{
        // treat sign bit
        if(data[len-1]>>(bits-1)!=o.data[len-1]>>(bits-1)){
            return data[len-1]>>(bits-1);
        }
        for(size_t i=len-1;~i;--i){
            if(data[i]!=o.data[i]) return data[i]<o.data[i];
        }
        return false;
    }
    bool operator>(Bigint_Fixedsize const&o)const{
        return o<*this;
    }
    bool operator<=(Bigint_Fixedsize const&o)const{
        return !(operator>(o));
    }
    bool operator>=(Bigint_Fixedsize const&o)const{
        return !(operator<(o));
    }
    bool operator==(Bigint_Fixedsize const&o)const{
        for(size_t i=0;i<len;++i){
            if(data[i] !=o.data[i]) return false;
        }
        return true;
    }
    bool operator!=(Bigint_Fixedsize const&o)const{
        return !operator==(o);
    }
    bool operator!()const{
        for(size_t i=0;i<len;++i){
            if(data[i]) return false;
        }
        return true;
    }
    bool is_negative()const{
        return data[len-1]>>(bits-1);
    }
    void print_binary(ostream&o, Bigint_Fixedsize const&b){
        o << "[";
        for(size_t i=len;i>0;--i){
            o << bitset<bits>(b.data[i-1]);
        }
        o << "]";
    }
    friend ostream& operator<<(ostream&o, Bigint_Fixedsize const&b){
        if(b.is_negative()){
            return o << '-' << -b << "\n";
        }
        Bigint_base10 ret(0);
        int64_t base = 1u<<bits;
        for(int i = len-1;i>=0;--i){
            ret*=base;
            ret+=Bigint_base10(b.data[i]);
        }
        o << ret;

        return o;
    }
    explicit operator long double()const{
        if(is_negative()){
            return (long double)operator-();
        }
        long double ret = 0.0;
        long double base = 1u<<bits;
        for(int i = len-1;i>=0;--i){
            ret = ret * base + data[i];
        }
        return ret;
    }
    /// TODO: implement for larger inputs
    friend istream& operator>>(istream&i, Bigint_Fixedsize &b){
        int64_t tmp;
        i >> tmp;
        b = Bigint_Fixedsize(tmp);
        return i;
    }
};

/**
 *  Biginteger with fixed precision
 *  Has 32*len bits, is signed
 *
 *  Trying out more optimizations.
 */
template<size_t len>
struct Bigint_Fixedsize_Fast{
    unsigned int data[len];
    uint16_t siz;
    bool sign;
    static constexpr unsigned int bits = 32;
    Bigint_Fixedsize_Fast(){
        data[0] = 0;
        siz = 1;
        sign = false;
    }
    Bigint_Fixedsize_Fast(long long a){
        sign = false;
        if(a<0){
            sign = true;
            a=-a;
        }
        siz = 0;
        do{
            long long b = a>>bits;
            data[siz] = a - (b<<bits);
            a = b;
            ++siz;
        } while(a);
    }
    void trim(){
        while(siz>1 && !data[siz-1]) --siz;
        if(siz == 1 && data[0] == 0) sign=false;
    }
    int comp_unsigned(Bigint_Fixedsize_Fast const&o)const{
        uint16_t lim = min(siz, o.siz);
        for(unsigned int i=lim;i<siz;++i){
            if(data[i]){
                return 1;
            }
        }
        for(unsigned int i=lim;i<o.siz;++i){
            if(o.data[i]){
                return -1;
            }
        }
        for(unsigned int i=lim-1;i+1;--i){
            if(data[i]!=o.data[i]){
                return data[i] < o.data[i]?-1:1;
            }
        }
        return 0;
    }
    int comp(Bigint_Fixedsize_Fast const&o)const{
        int sign_mul = 1-2*sign;
        if(sign != o.sign){
            return sign_mul;
        }
        return sign_mul * comp_unsigned(o);
    }
    bool operator<(Bigint_Fixedsize_Fast const&o)const{
        return comp(o)<0;
    }
    bool operator>(Bigint_Fixedsize_Fast const&o)const{
        return comp(o)>0;
    }
    bool operator<=(Bigint_Fixedsize_Fast const&o)const{
        return comp(o)<=0;
    }
    bool operator>=(Bigint_Fixedsize_Fast const&o)const{
        return comp(o)>=0;
    }
    bool operator==(Bigint_Fixedsize_Fast const&o)const{
        return comp(o)==0;
    }
    bool operator!=(Bigint_Fixedsize_Fast const&o)const{
        return comp(o)!=0;
    }
    bool operator!()const{
        return operator==(ZERO);
    }
    Bigint_Fixedsize_Fast operator-()const{
        Bigint_Fixedsize_Fast ret(*this);
        if(!!ret){
            ret.sign ^=1;
        }
        return ret;
    }
    Bigint_Fixedsize_Fast operator*(Bigint_Fixedsize_Fast const&o)const{
        Bigint_Fixedsize_Fast ret;
        ret.siz = min(siz+o.siz, (int)len);
        ret.sign = (sign!=o.sign);
        fill(ret.data, ret.data+ret.siz, 0);
        for(unsigned int i=0;i<siz;++i){
            unsigned long long carry = 0, carry_2;
            for(unsigned int j=0;j<o.siz;++j){
                carry+= data[i]*(unsigned long long)o.data[j] + ret.data[i+j];
                carry_2 = carry >> bits;
                ret.data[i+j] = carry - (carry_2<<bits);
                carry = carry_2;
            }
            for(unsigned int j=i+o.siz;carry;++j){
                carry+= ret.data[j];
                carry_2 = carry >> bits;
                ret.data[j] = carry - (carry_2<<bits);
                carry = carry_2;
            }
        }
        ret.trim();
        return ret;
    }
    Bigint_Fixedsize_Fast& operator*=(Bigint_Fixedsize_Fast const&o){
        *this = operator*(o);
        return *this;
    }
    static void unsigned_add(Bigint_Fixedsize_Fast &ret, Bigint_Fixedsize_Fast const&A, Bigint_Fixedsize_Fast const&B){
        const Bigint_Fixedsize_Fast *a = &A, *b = &B;
        if(a->siz < b->siz) swap(a, b);
        ret.sign = A.sign;
        unsigned long long carry = 0, carry_2;
        unsigned int j;
        for(j=0;j<b->siz;++j){
            carry+=(unsigned long long)a->data[j] + (unsigned long long)b->data[j];
            carry_2 = carry>>bits;
            ret.data[j] = carry - (carry_2<<bits);
            carry = carry_2;
        }
        for(;j<a->siz;++j){
            carry+=a->data[j];
            carry_2 = carry>>bits;
            ret.data[j] = carry - (carry_2<<bits);
            carry = carry_2;
        }
        if(carry){
            ret.data[j++] = carry;
        }
        ret.siz = j;
        ret.trim();
    }
    static void unsigned_subtract(Bigint_Fixedsize_Fast &ret, Bigint_Fixedsize_Fast const&A, Bigint_Fixedsize_Fast const&B){
        int com = A.comp_unsigned(B);
        if(com == 0){
            ret.sign = false;
            ret.siz = 1;
            ret.data[0] = 0;
            return;
        }
        ret.sign = A.sign;
        const Bigint_Fixedsize_Fast *a = &A, *b = &B;
        if(com < 0){
            ret.sign ^= true;
            swap(a, b);
        }
        // deal with case then o is not trimed.
        unsigned int min_siz = min(A.siz, B.siz);
        unsigned long long carry = 0, carry_2;
        unsigned int j;
        for(j=0;j<min_siz;++j){
            carry+=(unsigned long long)a->data[j] - (unsigned long long)b->data[j];
            carry_2 = carry>>bits;
            ret.data[j] = carry - (carry_2<<bits);
            carry = -(carry_2 & 1u);
        }
        for(;carry;++j){
            assert(j < a->siz);
            carry+=a->data[j];
            carry_2 = carry>>bits;
            ret.data[j] = carry - (carry_2<<bits);
            carry = -(carry_2 & 1u);
        }
        copy(a->data+j, a->data+a->siz, ret.data+j);
        ret.siz = a->siz;
        ret.trim();
    }
    static void add(Bigint_Fixedsize_Fast &ret, Bigint_Fixedsize_Fast const&A, Bigint_Fixedsize_Fast const&B){
        if(A.sign == B.sign){
            unsigned_add(ret, A, B);
        } else {
            unsigned_subtract(ret, A, B);
        }
    }
    static void sub(Bigint_Fixedsize_Fast &ret, Bigint_Fixedsize_Fast const&A, Bigint_Fixedsize_Fast const&B){
        if(A.sign != B.sign){
            unsigned_add(ret, A, B);
        } else {
            unsigned_subtract(ret, A, B);
        }
    }
    Bigint_Fixedsize_Fast operator+(Bigint_Fixedsize_Fast const&o)const{
        Bigint_Fixedsize_Fast ret;
        add(ret, *this, o);
        return ret;
    }
    Bigint_Fixedsize_Fast& operator+=(Bigint_Fixedsize_Fast const&o){
        add(*this, *this, o);
        return *this;
    }
    Bigint_Fixedsize_Fast operator-(Bigint_Fixedsize_Fast const&o)const{
        Bigint_Fixedsize_Fast ret;
        sub(ret, *this, o);
        return ret;
    }
    Bigint_Fixedsize_Fast operator-=(Bigint_Fixedsize_Fast const&o){
        sub(*this, *this, o);
        return *this;
    }
    /// TODO: more operators
    void print_binary(ostream&o, Bigint_Fixedsize_Fast const&b){
        o << "[";
        o << sign << "/" << len << "/";
        for(size_t i=siz;i>0;--i){
            o << bitset<bits>(b.data[i-1]);
        }
        o << "]";
    }
    friend ostream& operator<<(ostream&o, Bigint_Fixedsize_Fast const&b){
        if(b.sign){
            return o << '-' << -b;
        }
        Bigint_base10 ret(0);
        int64_t base = 1ll<<bits;
        for(int i = b.siz-1;i>=0;--i){
            ret*=base;
            ret+=Bigint_base10(b.data[i]);
        }
        o << ret;

        return o;
    }
    explicit operator long double()const{
        if(sign){
            return (long double)operator-();
        }
        long double ret = 0.0;
        long double base = 1ll<<bits;
        for(int i = siz-1;i>=0;--i){
            ret = ret * base + data[i];
        }
        return ret;
    }
    /// TODO: implement for larger inputs
    friend istream& operator>>(istream&i, Bigint_Fixedsize_Fast &b){
        int64_t tmp;
        i >> tmp;
        b = Bigint_Fixedsize_Fast(tmp);
        return i;
    }
    static const Bigint_Fixedsize_Fast ZERO;
};
template<size_t len>
const Bigint_Fixedsize_Fast<len> Bigint_Fixedsize_Fast<len>::ZERO(0);

#define lp_debug(x) do{}while(0)

/**
 *  Randomized LP in expected
 *  O(d! 4^d n)
 *  Does exact calculations.
 */
template<typename FLOAT>
class Lp_Seidel{
private:


    // orthogonal projection of 'vec' into 'plane'
    vector<FLOAT> proj_down(vector<FLOAT> const&vec, vector<FLOAT> const&plane, size_t j){
        assert(vec.size() <= plane.size() && plane.size()<=vec.size()+1);
        assert(j+1 < plane.size());
        assert(!!plane[j]);
        vector<FLOAT> ret (vec.size()-1);
        //FLOAT tmp;
        if(plane[j] < FLOAT(0)){
            for(size_t i=0;i+1<vec.size();++i){
                ret[i] = vec[j] * plane[i+(i>=j)] - vec[i+(i>=j)]* plane[j];
            }
        } else {
            for(size_t i=0;i+1<vec.size();++i){
                ret[i] = vec[i+(i>=j)]*plane[j] - vec[j]*plane[i+(i>=j)];
            }
        }
        return ret;
    }

    // orthogonal projection of 'vec' out of 'plane'
    vector<FLOAT> proj_up(vector<FLOAT> const&vec, vector<FLOAT> const&plane, size_t j){
        assert(vec.size()+1 == plane.size());
        assert(j+1 < plane.size());
        assert(!!plane[j]);
        vector<FLOAT> ret(vec.size()+1);
        copy(vec.begin(), vec.begin()+j, ret.begin());
        copy(vec.begin()+j, vec.end(), ret.begin()+j+1);
        for(size_t i=0;i<vec.size();++i){
            ret[j]+=vec[i]*plane[i+(i>=j)];
        }
        FLOAT denom = plane[j];
        if(denom < FLOAT(0)){
            denom = -denom;
        }
        for(size_t i=0;i<vec.size();++i){
            ret[i+(i>=j)]*=denom;
        }
        if(plane[j] >= FLOAT(0)){
            ret[j] = -ret[j];
        }
        return ret;
    }
    FLOAT planescal(vector<FLOAT> const&x, vector<FLOAT> const&a){
        assert(x.size() == a.size());
        FLOAT ret=0;
        for(size_t i=0;i<x.size();++i){
            ret+=x[i]*a[i];
        }
        return ret;
    }

    // solve lp recursively
    vector<FLOAT> solve(vector<vector<FLOAT> > const &A, vector<FLOAT> const&c, int d, FLOAT const& barier_loc){
        int n=A.size();
        if(d==1){ // base case: single dimension
            vector<FLOAT> ret(2);
            ret[0] = (c[0]<FLOAT(0) ? -barier_loc : barier_loc);
            ret[1] = 1ull;
            for(int i=0;i<n;++i){
                if(ret[0]*A[i][0]+ret[1]*A[i].back()>FLOAT(0)){
                    if(!A[i][0]){
                        lp_debug("infeasible single\n");
                        return vector<FLOAT>();
                    }
                    ret[0] = -A[i].back();
                    ret[1] = A[i][0];
                    /*if(ret[1].is_negative()){
                        ret[1].negate();
                        ret[0].negate();
                    }*/
                    if(ret[1] < FLOAT(0)){
                        ret[1] = -ret[1];
                        ret[0] = -ret[0];
                    }
                    lp_debug(" -> " << i << "\n");
                }
            }
            for(int i=0;i<n;++i){
                if(ret[0]*A[i][0]+ret[1]*A[i].back()>FLOAT(0)){
                    lp_debug("infeasible\n");
                    return vector<FLOAT>();
                }
            }
            return ret;
        }
        FLOAT barier_next = barier_loc * barier_loc;
        // initial solution
        vector<FLOAT> x(d+1);
        for(int i=0;i<d;++i){
            x[i] = (c[i]<FLOAT(0)?-barier_loc:barier_loc);
        }
        x.back() = FLOAT(1);
        for(size_t i=0;i<A.size();++i){
            if(planescal(x, A[i])>FLOAT(0)){
                int k = 0;
                while(k<d && !A[i][k]) ++k;
                // recurse
                if(k==d) {lp_debug("what?\n"); return vector<FLOAT>();} // degenerate failing plane??????
                vector<vector<FLOAT> > A2(i);
                for(size_t j=0;j<A2.size();++j){
                    A2[j] = proj_down(A[j], A[i], k);
                }
                shuffle(A2.begin(), A2.end(), Rng::get_engine());
                lp_debug(string(2*d, ' ') << i << "\n");
                vector<FLOAT> c2 = proj_down(c, A[i], k);
                vector<FLOAT> x2 = solve(A2, c2, d-1, barier_next);
                if(x2.empty()) return x2; // infeasible
                x = proj_up(x2, A[i], k);
                lp_debug(string(2*d, ' ') << ":");
                lp_debug("";for(auto const&e:x) lp_debug(" " << e));
                lp_debug("\n");
            }
        }
        return x;
    }

public:
    vector<FLOAT> solve(vector<vector<FLOAT> > const &A, vector<FLOAT> const&c, const FLOAT& barier = FLOAT((long long)1e9)){
        assert(A.empty() || A[0].size() == c.size()+1);
        return solve(A, c, c.size(), barier);
    }
    /**
     *  Maximize c^T x
     *  subject to Ax <= b
     *
     *  Returns empty vector if infeasible
     */
    vector<FLOAT> solve(vector<vector<FLOAT> > A, vector<FLOAT> const&b, vector<FLOAT> const&c){
        assert(A.size() == b.size());
        for(unsigned int i=0;i<A.size();++i){
            A[i].push_back(-b[i]);
        }
        return solve(A, c);
    }
};

template<typename Big_Int, bool use_two_phase = true>
class Lp_Clarkson{
private:
    /**
     *  Returns a sub-multiset of size siz uniformly at random
     *  out of the set where i is present weight[i] times.
     *
     *  Runs in O(|weight| + siz^2) expected time.
     *  Could be optimized
     */
    vector<int> sample_subset(vector<int64_t> const&weight, unsigned int siz){
        int64_t total_weight = accumulate(weight.begin(), weight.end(), 0ll);
        vector<int64_t> samples;
        while(samples.size() < siz){
            int64_t new_sample = Rng::uniform<int64_t>(0, total_weight-1);
            if(find(samples.begin(), samples.end(), new_sample) == samples.end()){
                samples.push_back(new_sample);
            }
        }
        sort(samples.begin(), samples.end());
        vector<int> ret;
        int64_t left_weight = 0;
        for(unsigned int i=0, j=0;i<weight.size() && j<samples.size();){
            if(samples[j] < left_weight + weight[i]){
                ret.push_back(i);
                ++j;
            } else {
                left_weight+=weight[i];
                ++i;
            }
        }
        return ret;
    }
    /// violation check
    bool is_satisfied(vector<Big_Int> const&x, vector<Big_Int> const&a){
        assert(x.size() == a.size());
        Big_Int ret=0;
        for(size_t i=0;i<x.size();++i){
            ret+=x[i]*a[i];
        }
        return ret <= Big_Int(0);
    }
    vector<Big_Int> solve_two(vector<vector<Big_Int> > const&A, vector<Big_Int> const&c){
        const unsigned int sample_size = c.size()*c.size()*4;
        Lp_Seidel<Big_Int> sub_lp;
        // to few constrains -> use other solver
        if(A.size() < sample_size){
            return sub_lp.solve(A, c);
        } else {
            int constraints = A.size();
            int variables = c.size();
            vector<int64_t> weight(constraints, 1);
            vector<Big_Int> x;
            vector<vector<Big_Int> > subproblem_A;
            vector<char> is_violated(constraints, 0);
            for(unsigned int iteration=1;;++iteration){
                subproblem_A.clear();
                vector<int> subspace = sample_subset(weight, sample_size);
                for(int const&e:subspace){
                    subproblem_A.push_back(A[e]);
                }

                x = sub_lp.solve(subproblem_A, c);
                // infeasible case
                if(x.empty()){
                    return x;
                }

                int64_t total_violated = 0;
                for(int i=0;i<constraints;++i){
                    is_violated[i] = !is_satisfied(x, A[i]);
                    if(is_violated[i]){
                        total_violated+=weight[i];
                    }
                }
                if(total_violated == 0){
                    cerr << "Iterations: " << iteration;
                    cerr << ", max weight: " << *max_element(weight.begin(), weight.end()) << "\n";
                    break;
                }
                if(total_violated*3*variables <= accumulate(weight.begin(), weight.end(), 0ll)){
                    for(int i=0;i<constraints;++i){
                        if(is_violated[i]){
                            weight[i]*=2;
                        }
                    }
                    assert_msg(accumulate(weight.begin(), weight.end(), 0ll) < (1ll<<62), "Weight overflow");
                }
            }
            return x;
        }
    }
    vector<Big_Int> solve_one(vector<vector<Big_Int> > const&A, vector<Big_Int> const&c){
        const unsigned int constraints = A.size(), variables = c.size();
        if(constraints <= variables*variables*6){
            return solve_two(A, c);
        } else {
            const unsigned int sqrt_constraints = sqrt(constraints);
            const unsigned int sample_size = variables * sqrt(constraints);
            vector<Big_Int> x;
            vector<vector<Big_Int> > subproblem_A;
            vector<int> violations;
            for(unsigned int iteration=1;;++iteration){
                //function<vector<int>()> samp = [&, this](){return sample_subset(vector<int64_t>(constraints, 1), sample_size);};
                //vector<int> subspace = Timer::execute_timed<vector<int>>(samp, "Sampling");
                vector<int> subspace = sample_subset(vector<int64_t>(constraints, 1), sample_size);
                for(int const&e:subspace){
                    subproblem_A.push_back(A[e]);
                }

                x = solve_two(subproblem_A, c);
                // infeasible case
                if(x.empty()){
                    return x;
                }
                violations.clear();
                for(unsigned int i=0;i<constraints;++i){
                    if(!is_satisfied(x, A[i])){
                        violations.push_back(i);
                    }
                }
                cerr << "Violations: " << violations.size() << " / " << 2*sqrt_constraints << "\n";
                if(violations.empty()){
                    cerr << "Iterations: " << iteration;
                    cerr << ", used constraints:" << subproblem_A.size() << "\n";
                    break;
                }
                subproblem_A.erase(subproblem_A.end()-sample_size, subproblem_A.end());
                if(violations.size() <= 2*sqrt_constraints){
                    for(int const&e : violations){
                        subproblem_A.push_back(A[e]);
                    }
                }
            }
            return x;
        }
    }



public:
    vector<Big_Int> solve(vector<vector<Big_Int> > const&A, vector<Big_Int> const&c){
        if(use_two_phase){
            return solve_one(A, c);
        } else {
            return solve_two(A, c);
        }
    }

    /**
     *  Maximize c^T x
     *  Subject to Ax <= b
     *
     *  Returns empty vector if infeasible
     */
    vector<Big_Int> solve(vector<vector<Big_Int> > A, vector<Big_Int> const&b, vector<Big_Int> const&c){
        assert(A.size() == b.size());
        for(unsigned int i=0;i<A.size();++i){
            A[i].push_back(-b[i]);
        }
        return solve(A, c);
    }
};
using ll=long long;
using lf=long double;
using elf=Bigint_Fixedsize_Fast<10>;

vector<vector<elf>>A;
vector<elf>B;
vector<elf>x;

int main()
{
    cin.tie(0)->sync_with_stdio(0);
	// simplex with 7 dimensions, each being r, a1, a2, b1, b2, c1, c2 (1 and 2 exist due to bypassing nonnegativity constraints)
    // max -r
    // s.t. -r - a1*x^2 - b1*x - c1 <= -y
    //      -r + a1*x^2 + b1*x + c1 <= y
    vector<elf>c{-1,0,0,0};
    int q;cin>>q;
    while(q--)
    {
        int n;cin>>n;
        while(n--)
        {
            ll x,y;
            cin>>x>>y;
            B.push_back(-y);
            A.push_back({-1,-x*x,-x,-1});
            B.push_back(y);
            A.push_back({-1,x*x,x,1});
        }
        Lp_Clarkson<elf>solver;
        auto res=solver.solve(A,B,c);
        lf ans=(lf)res[0]/(lf)res[4];
        cout<<setprecision(12)<<fixed<<ans*ans<<"\n";
    }
}

詳細信息

Test #1:

score: 100
Accepted
time: 0ms
memory: 3824kb

input:

1
4
0 0
1 3
2 9
3 0

output:

5.062500000000

result:

ok found '5.0625000', expected '5.0625000', error '0.0000000'

Test #2:

score: -100
Runtime Error

input:

60
1
1000 -990
2
171 -638
949 -99
2
633 227
-257 -602
3
634 -994
633 999
-374 995
3
445 -110
586 -121
462 29
9
-995 -224
-458 -833
691 -670
456 -259
-376 55
-563 -12
834 827
-826 -220
299 744
17
997 991
997 976
997 988
998 -986
999 -982
999 -980
999 -996
998 -988
998 -991
997 987
1000 996
999 -1000
...

output:

0.000000000000
0.000000000000
159849.477731154145
991594.870178950218
991594.870178950218
992269.005360952082
995976.553062163949
995976.553062163949
995979.578552127102
999000.250000000000
999000.250000000000
999000.250000000000
999000.250000000000
1000000.000000000000
1000000.000000000000
1000000....

result: