水果忍者

题目描述

    你知道水果忍者吗？这是一个智能手机上面很流行的游戏。

    这个出色的游戏迎合了人类最喜爱的运动轨迹：抛物线。就是各种各样物体被扔出去之后都会形成的曲线运动。现在这里的问题也是关于抛物线的。

    现在有N个水果，假设它们就是一个个在平面上的的圆，有着同样的半径。每个水果有着各自的位置和移动速度。由于重力的影响（这里假设重力常数g=10)，所有水果都从一开始就沿着抛物线运动。你的目标就是在给定的时间限制之内，一次切到最多的水果，即是找到在某一时刻的一条直线穿过最多的圆。

输入格式

输入包含多给测试数据。

输入的第一行给出了测试数据个数T。

每个测试数据首先给出三个整数，水果数目N(<=100)，水果的半径R(<=1000)，还有时间限制E(<=1000)。

接下来有N行，每行给出4个实数，分别对应圆心坐标(xi,yi)和速度(vxi,vyi)。

所有数都是采用标准单位，所以你可以放心使用以下著名的物理公式来计算某一时刻水果所处的位置：

x' = x + vx * t

y' = y + vy * t - 0.5 * 10 * t * t

数据保证所有的浮点数的绝对值小于100000。

输出格式

对应每个测试数据，输出一个整数，即在时限内可以同时切到的水果数目的最大值。

样例输入

2

4 1 2

2.0 4.0 1.0 0.0

5.0 4.0 2.0 0.0

2.0 12.0 1.5 -4.0

-2.0 -5.0 10.0 10.0

4 1 1

2.0 4.0 1.0 0.0

5.0 4.0 2.0 0.0

2.0 12.0 1.5 -4.0

-2.0 -5.0 10.0 10.0

样例输出

4

3

提示

请注意一条直线“穿过”一个圆意味着，圆心与直线的距离小于圆的半径加上精度误差(1e-5)。

#include <cstdio>

#include <cstring>

#include <cmath>

#include <cassert>

#include <algorithm>

#include <vector>

using namespace std;

const int maxn = 101;

const double INF = 1e8;

const double EPS = 1e-5;

const double MAXLOOP = 100;

struct seg {

    seg() { l = 1; r = 0; }

    seg(double il, double ir) { l = il; r = ir; }

    bool empty() { return r < l - EPS; }

    double l, r;

};

struct point {

    int mrk;

    double pos;

}points[maxn * 2];

struct vect {

   

    double d[20];

    int n;

    int size() {

        return n;

    }

    vect() {

        n= 0;

    }

    void push_back(double x) {

        d[n++] = x;

    }

    double operator [] (int index) {

        return d[index];

    }

};

int pn;

double x[maxn], y[maxn], vx[maxn],vy[maxn];

int n, r, e; // number, radius, time limit

int best; // maximum number of cutting

void add_seg(seg &s) {

    if (abs(s.l - s.r)<EPS) s.r = s.l;

    points[pn].mrk = 1;

    points[pn].pos = s.l;

    ++pn;

    points[pn].mrk = 0;

    points[pn].pos = s.r;

    ++pn;

}

vectdifferentiate(vect co) {

    vect res;

    for (int i = 1; i < co.size(); ++i) {

        res.push_back(co * i);

    }

    return res;

}

double substitute(vect co, double x) {

    double res = 0, tmp = 1;

    for (int i = 0; i < co.size(); ++i) {

        res += tmp * co;

        tmp *= x;

    }

    return res;

}

vectremove_dup(vect vec) {

    vect res;

    if (vec.size() == 0) return res;

    //sort(vec.begin(), vec.end());

    sort(vec.d, vec.d+vec.n);

    res.push_back(vec[0]);

    for (int i=1;i<vec.size();++i) {

        if (vec - vec[i-1] > EPS * 10)res.push_back(vec);

    }

    return res;

}

bool newton(vect co, vect dco, double x, double &res_x) {

    double last = -INF;

    int count = 0;

    while (count++ < MAXLOOP) {

        double f = substitute(co, x);

        double fp = substitute(dco, x);

        if ( abs(fp) < EPS ) {

            // no idea whether this works

            x += EPS;

            continue;

        }

        x -= f / fp;

        if ( abs(x - last) < EPS ) {

            res_x = x;

            return true;

        }

        last = x;

    }

    return false;

}

vectsolveeqs(vect co) {

    vect root;

   

    if (co.size() <= 2) {

        if (co.size() == 2 && co[1] != 0) {

            root.push_back(-co[0] / co[1]);

        }

        return root;

    }

    vect dco = differentiate(co); // differentail coefficient

    vect droot = solveeqs(dco);

    droot.push_back(-INF);

    droot.push_back(INF);

    sort(droot.d, droot.d+droot.n);

    //sort(droot.begin(), droot.end());

   

    for (int i=1;i<droot.size(); ++i) {

        double initx = (droot[i-1] + droot)* 0.5, res;

        if (newton(co, dco, initx, res)) {

            root.push_back(res);

        }

    }

    return remove_dup(root);

}

//sigma(ax^n) <= 0

vector<seg> solvele(vect co) {

    vect root = solveeqs(co);

    vector<seg> res;

    bool neg = substitute(co, -INF) < 0 ;

    double last = -INF;

    for (int i=0;i<root.size();++i) {

        if (neg) {

            res.push_back(seg(last, root));

        }

        last = root;

        neg ^=1;

    }

    if (neg) res.push_back(seg(last, INF));

    return res;

}

vectnegpoly(vect co) {

    for (int i = 0; i < co.size(); ++i) {

        //co *= -1;

        co.d *= -1;

    }

    return co;

}

vector<seg> solvege(vect co) {

    return solvele(negpoly(co));

}

vectmultipoly(vect co1, vect co2) {

    vect res;

    for (int i = 0; i<co1.size()+co2.size()-1;++i) res.push_back(0);

    for (int i = 0; i < co1.size(); ++i) {

        for (int j = 0; j < co2.size(); ++j) {

            //res[i+j] += co1 * co2[j];

            res.d[i+j] += co1 * co2[j];

        }

    }

    return res;

}

vectsumpoly(vect co1, vect co2) {

    vect res;

    for (int i = 0; i < co1.size() || i < co2.size(); ++i) {

        double ta, tb;

        if (i < co1.size()) ta = co1; else ta = 0;

        if (i < co2.size()) tb = co2; else tb = 0;

        res.push_back(ta + tb);

    }

    return res;

}

// ax +b

vect lco(double a, double b) {

    vect res;

    res.push_back(b);

    res.push_back(a);

    return res;

}

segintersect(seg a, seg b) {

    seg res;

    if (a.l > b.l) res.l = a.l; else res.l = b.l;

    if (a.r < b.r) res.r = a.r; else res.r = b.r;

    return res;

}

void fix_time_limit(seg &ret) {

    if (ret.l < 0) ret.l = 0;

    if (ret.r > e) ret.r = e;

}

void calc_valid_seg( int a, int b, int c) {

    double x1 = x[b] - x[a], x2 = x[c] - x[a];

    double y1 = y[b] - y[a], y2 = y[c] - y[a];

    double vx1 = vx[b] - vx[a], vx2 = vx[c] - vx[a];

    double vy1 = vy[b] - vy[a], vy2 = vy[c] - vy[a];

    vect tmp;   

    vect co = sumpoly(multipoly(lco(vy1, y1), lco(vx2, x2)),negpoly(multipoly(lco(vy2, y2), lco(vx1, x1))));

    //  int start2 =clock();

    vector<seg> seg1 = solvege(co);

    tmp.push_back(-4*r*r);

    vect co2 = sumpoly(multipoly(co, co),multipoly(sumpoly(multipoly(lco(vx1,x1), lco(vx1,x1)), multipoly(lco(vx1,x1),lco(vx1,x1))), tmp));

    vector<seg> seg2 = solvele(co2);

    for (int i = 0; i < seg1.size(); ++i) {

        fix_time_limit(seg1);

        if (seg1.empty()) continue;

        for (int j = 0; j < seg2.size(); ++j){

            fix_time_limit(seg2[j]);

            if (seg2[j].empty()) continue;

            seg tmpseg = intersect(seg1,seg2[j]);

            if (!tmpseg.empty())add_seg(tmpseg);

        }

    }

}

int cmp(const point & l, const point & r) {

    if (l.pos != r.pos) return l.pos < r.pos;

    return l.mrk < r.mrk;

}

void sort_seg() {

    sort(points, points + pn, cmp);

}

void calc_max() {

    int i, tot = 0, maxx = 0;

    for (i = 0; i < pn; ++i) {

        if (points.mrk) ++tot; else --tot;

        if (tot > maxx) maxx = tot;

    }

    if (maxx + 2 > best) best = maxx + 2;

}

int main() {

    //freopen("D:\\Contest\\4+2\\2011\\ninjapro\\input1.txt","r", stdin);

    int t, i, j, k;

    scanf("%d", &t);

    while (t--) {

        scanf("%d%d%d", &n, &r, &e);

        for (i=0;i<n;++i) {

            scanf("%lf%lf%lf%lf", &x, &y,&vx, &vy);

        }

        best = 2;

        // enumerate the fixed line

        for (i=0;i<n;++i) {

            for (j=0;j<n;++j) {

                pn = 0;

                if (i == j) continue;

                for (k = 0; k < n; ++k) {

                    if (i == k || j == k) continue;

                    calc_valid_seg(i, j, k);

                }

                sort_seg();

                calc_max();

                //printf("%d %d\n", i,j);

            }

        }

        printf("%d\n", best);

    }

    return 0;

}