yx - C++博客

poj 1730 Perfect Pth Powers

通过这道题确实体会到A掉数学题确实还是需要经验了，不能猜对哪个地方会丧失精度的话，会一直wa的。其实，这道题我只想出了一半。
题意是 a的p次方 = n，其中n是32位整数，a和p都是整数，求满足条件的最大p。好吧，虽然我是在学数论，但是看到这题，我还是想起了
取对数。那么可以得到，p = ln(n) / ln(a)。既然要求最大的p，那么a最小即可了。那么直接从2开始枚举a不就可以了么。
可是直接枚举a的话肯定会超时的，因为a的范围太大了，比如n的是个大素数，a的范围就是2-n了，一定超时了。然后，我又想出另外一
种方法，对n分解因子，p就是所有因子的指数的最大公约数。呵呵，第二种方法更加会无情的超时，由于int范围很大，实现搞个素数表也不
可能。还是感觉时间不多了，就不多想了，然后搜了下，发现一句话，意识是枚举p。顿时觉得开朗起来，因为p最多是32。由前面可以得到
ln(a) = ln(n) / p。那么只要从32到1枚举p，保证a是整数即可。
后面发现这样精度难于控制，各种原因反正过不了题，看网上的代码，改成计算指数的形式了。因为 a = n的(1/p)次，这个可以用pow函
数算出来，如果a是整数，那么再计算pow(a,p)就会是n了。最难控制的是精度了，还有说n是负数的情况。不知道为什么直接处理负数答案
一直不对，只好把负数变为正数，同时判断p不能是偶数。

代码如下：

#include <stdio.h>
#include <math.h>

int main()
{
    double fN;//用double就不会溢出了,负数就可以直接转换为正数了

    while (scanf("%lf", &fN), fN)
    {
        bool bFlag = false;
        double fP = 31.0;
        if (fN < 0){fP = 32.0; fN = -fN; bFlag = true;};

        while (fP > 0)
        {
            //必须加上一个精度,防止往下误差
            double fA = pow(fN, 1.0 / fP) + 1e-8;
            //fA必须转换为int,因为一点点误差,pow之后就会放大很多
            double fTemp = pow((int)fA, fP);

            //必须对负数特殊判断,不可能出现偶数的p
            if (fabs(fN - fTemp) < 1e-8 && (!bFlag || ((int)fP) % 2))
            {
                printf("%.f\n", fP);
                break;
            }
            fP -= 1.0;
        }
    }

    return 0;
}

posted @ 2012-07-26 20:59 yx 阅读(918) | 评论 (0) | 编辑收藏

hdu 3627 Giant For

这个题是对可排序数据的实时增加删除查找，那天做比赛的时候一点都不会，想来想去觉得平衡树可以做，但是写平衡树是件很难的事情。
后面知道线段数可以做，虽然数据的范围很大，但是可以在全部读入数据后排序再离散化，然后进行线段树的操作，具体的代码没有写。
今天队友在网上发现一种用map和set可以水掉这题的方法。原来，这个方法最主要的使用了map和set里面的upper_bound操作，以前
居然忘记了这个东西了。既然这样，map和set也可以查前驱和后继了，但是注意low_bound查到的是小于等于的键。这个代码，注意是用
了一个map< int, set<int> > 集合把坐标都存起来了，进行添加删除和查找后继的操作。由于查找需要查找的元素是既比x大又比y大的元
素，就比较麻烦，需要循环x往后查找，但是这样就无情的超时了。然后，有一个优化，记录y的数目，那么当出现很大的y的时候，就不需要
查找了，然后才过了这个题。但是，数据变成很大的y对应的x很小的话，那么绝对过不了这个题了，只能用线段树做了。
现在觉得用map和set查找前驱和后继确实能水掉一些题啊。

代码如下：

#include <map>
#include <set>
#include <stdio.h>
using namespace std;

map< int, set<int> > ms;//存储x,y
map< int, set<int> >::iterator it;
map<int, int> my;//存储y的数目
set<int>::iterator msit;
int main()
{
    int nN;
    int nCase = 1;
    char szCmd[10];
    int nX, nY;
    int nTemp;

    while(scanf("%d", &nN), nN)
    {
        if (nCase > 1)
        {
            printf("\n");
        }

        printf("Case %d:\n", nCase++);
        ms.clear();
        my.clear();
        while (nN--)
        {
            scanf("%s", szCmd);
            scanf("%d%d", &nX, &nY);
            if (szCmd[0] == 'a')
            {
                if (my.find(nY) == my.end())
                {
                    my[nY] = 1;
                }
                else
                {
                    my[nY]++;
                }

                if (ms.find(nX) == ms.end())
                {
                    ms[nX].insert(nY);
                }
                else
                {
                    msit = ms[nX].find(nY);
                    if (msit == ms[nX].end())//会出现重复的数据
                    {
                        ms[nX].insert(nY);
                    }
                }
            }
            else if (szCmd[0] == 'r')
            {
                ms[nX].erase(nY);
                if(ms[nX].size() == 0)
                {
                    ms.erase(nX);
                }
                my[nY]--;
                if (my[nY] == 0)
                {
                    my.erase(nY);
                }
            }
            else if (szCmd[0] == 'f')
            {
                if (my.upper_bound(nY) == my.end())
                {
                    printf("-1\n");
                    continue;
                }
                while (true)
                {
                    it = ms.upper_bound(nX);
                    if (it == ms.end())//比nX大的不存在
                    {
                        printf("-1\n");
                        break;
                    }
                    nTemp = it->first;
                    msit = ms[nTemp].upper_bound(nY);
                    if (msit == ms[nTemp].end())//比nY大的不存在
                    {
                        nX = nTemp;
                        continue;//那么增加x,继续往后查
                    }
                    else
                    {
                        printf("%d %d\n", nTemp, *msit);
                        break;
                    }
                }
            }
        }
    }

    return 0;
}

posted @ 2012-07-26 12:22 yx 阅读(877) | 评论 (0) | 编辑收藏

poj 2551 Ones and poj 2262 Goldbach's Conjecture

第一个题用到了同余的性质，这是数论里面最基本的性质，但是做题时候不一定能够自己发现。题意是n*m = 11111...，给出n，
用一个m乘以n得到的答案全是1组成的数字，问1最小的个数是多少。可以转换为n*m=(k*10+1)，那么可以得到(k*10+1)%n==0。
当然最开始的k是1，那么我们不断的增长k = （10*k+1）。看增长多少次，就是有多少个1了。因为要避免溢出，所以需要不断%n。
因为同余的性质，所以可以保证%n之后答案不变。
第二个用到素数筛选法。素数筛选法的原理是筛去素数的倍数，由于是从小循环到大的，所以当前的值没被筛掉的话，则一定是素数，
这个判断导致复杂度不是n的平方。

poj 2551 代码：

#include <stdio.h>

int main()
{
    int nN;

    while (scanf("%d", &nN) == 1)
    {
        int nCnt = 1;
        int nTemp = 1;
        while (1)
        {
            if (nTemp % nN == 0)break;
            else nTemp = (nTemp * 10 + 1) % nN;
            ++nCnt;
        }
        printf("%d\n", nCnt);
    }

    return 0;
}

poj 2262 代码：

#include <stdio.h>
#include <string.h>
#include <math.h>

#define MAX (1000000 + 10)
bool bPrime[MAX];
void InitPrime()
{
    memset(bPrime, true, sizeof(bPrime));
    bPrime[0] = bPrime[1] = false;
    for (int i = 2; i <= MAX; ++i)
    {
        if (bPrime[i])
        for (int j = 2 * i; j <= MAX; j += i)
        {
            bPrime[j] = false;
        }
    }
}

int main()
{
    int nN;

    InitPrime();
    while (scanf("%d", &nN), nN)
    {
        int i;
        for (i = 2; i < nN; ++i)
        {
            if (i % 2 && (nN - i) % 2 && bPrime[i] && bPrime[nN - i])
            {
                printf("%d = %d + %d\n", nN, i, nN - i);
                break;
            }
        }
        if (i == nN)
        {
            printf("Goldbach's conjecture is wrong.\n");
        }
    }

    return 0;
}

posted @ 2012-07-25 22:35 yx 阅读(812) | 评论 (0) | 编辑收藏

poj 2187 Beauty Contest

这个题我是按照discussion里面的说法，先求凸包，然后枚举过的。因为开始先把求凸包算法里面的用到了数组名搞混了，无故wa了好
多次。后面求凸包换了种算法过了。结果发现bug是搞混了数组名，然后把前面wa掉的代码下载下来，改好之后也都过了。
这个题主要是凸包算法需要处理有重复点，有多点共线之类的情况。那个按极角排序后，再求凸包的算法，对共点共线处理的不是很好，
不过那个算法也过了这个题。有个直接按坐标排序后，再求上凸包和下凸包的算法，可以处理共点共线的情况。这个算法比较优美啊，既不
需要找y坐标最小的点，也不需要按极角排序，直接按坐标排序下，然后求凸包即可。
这个算法的一点解释：http://www.algorithmist.com/index.php/Monotone_Chain_Convex_Hull
另外，演算法笔记：http://www.csie.ntnu.edu.tw/~u91029/ConvexHull.html#a3上也有提到这个算法，我也是从这上面看到的。
这个算法可以假设是Graham排序基准点在无限远处，于是夹角大小的比较可以直接按水平坐标比较。

代码如下：

#include <stdio.h>
#include <algorithm>
#include <math.h>
using namespace std;
struct Point
{
    int x, y;
    bool operator<(const Point& p) const
    {
        return x < p.x || x == p.x && y < p.y;
    }
};
Point pts[50100];
Point pcs[50100];
int nN;
int nM;
inline int SDis(const Point& a, const Point& b)
{
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

void Convex()
{
    sort(pts, pts + nN);

    nM = 0;
    for (int i = 0; i < nN; ++i)
    {
        while(nM >= 2 && Cross(pcs[nM - 2], pcs[nM - 1], pts[i]) <= 0)
        {
            nM--;
        }
        pcs[nM++] = pts[i];
    }
    for (int i= nN - 2, t = nM + 1; i >= 0; --i)
    {
        while (nM >= t && Cross(pcs[nM - 2], pcs[nM - 1], pts[i]) <= 0)
        {
            nM--;
        }
        pcs[nM++] = pts[i];
    }
    nM--;//起点会被重复包含
}

int main()
{
    while (scanf("%d", &nN) == 1)
    {
        for (int i = 0; i < nN; ++i)
        {
            scanf("%d%d", &pts[i].x, &pts[i].y);
        }
        Convex();
        int nMax = -1;
        for (int i = 0; i < nM; ++i)
        {
            for (int j = i + 1; j < nM; ++j)
            {
                nMax = max(nMax, SDis(pcs[i], pcs[j]));
            }
        }
        printf("%d\n", nMax);
    }

    return 0;
}

也可以用旋转卡壳算法来求最远点对，此题的完整代码如下：

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <algorithm>
using namespace std;

struct Point
{
    int x, y;
    bool operator < (const Point& p)const
    {
        return x < p.x || x == p.x && y < p.y;
    }
};

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

//输入点集合,输出凸包
void Convex(vector<Point>& in, vector<Point>& out)
{
    int nN = in.size();
    int nM = 0;

    sort(in.begin(), in.end());
    out.resize(nN);

    for (int i = 0; i < nN; ++i)
    {
        while (nM >= 2 && Cross(out[nM - 2], out[nM - 1], in[i]) <= 0)
        {
            nM--;
        }
        out[nM++] = in[i];
    }

    for (int i = nN - 2, t = nM + 1; i >= 0; --i)
    {
        while (nM >= t && Cross(out[nM - 2], out[nM - 1], in[i]) <= 0)
        {
            nM--;
        }
        out[nM++] = in[i];
    }
    out.resize(nM);
    out.pop_back();//起始点重复
}

int SDis(Point a,Point b)
{
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

int RC(vector<Point>& vp)
{
    int nP = 1;
    int nN = vp.size();
    vp.push_back(vp[0]);
    int nAns = 0;
    for (int i = 0; i < nN; ++i)
    {
        while (Cross(vp[i], vp[i + 1], vp[nP + 1]) > Cross(vp[i], vp[i + 1], vp[nP]))
        {
            nP = (nP + 1) % nN;
        }
        nAns = max(nAns, max(SDis(vp[i], vp[nP]), SDis(vp[i + 1], vp[nP + 1])));
    }
    vp.pop_back();
    return nAns;
}

int main()
{
    int nN;
    vector<Point> in, out;
    Point p;
    while (scanf("%d", &nN) == 1)
    {
        in.clear(), out.clear();
        while (nN--)
        {
            scanf("%d%d", &p.x, &p.y);
            in.push_back(p);
        }
        Convex(in, out);

        printf("%d\n", RC(out));
    }

    return 0;
}

关于旋转卡壳的算法描述，网上有很多资料，比如，http://www.cppblog.com/staryjy/archive/2010/09/25/101412.html
尤其关于这个求最远点对的。

posted @ 2012-07-23 22:18 yx 阅读(949) | 评论 (0) | 编辑收藏

poj 3525 Most Distant Point from the Sea

这个题的题意是给定一个凸多边形表示的海岛，求海岛离大海最远的距离。可以转化为一个凸多边形内部最多能够放入一个多大的圆。
显然可以对圆的半径进行二分，但是怎么确定圆心了。确定是否存在圆心，可以把原来的凸多边形往内部移动r（圆的半径）的距离之后，
再对新的多边形求半平面交，如果半平面交存在（是个点即可），那么当前大小的圆能够放入。
求半平面交的算法可以用上一篇中的N*N复杂度的基本算法。本题还涉及到一个知识，就是如何把一条直线往逆时针方向或者顺时针方向
移动R的距离。其实，可以根据单位圆那种思路计算。因为相当于以原来直线上的一点为圆心，以r为半径做圆，而且与原来的直线成90的夹
角，那么后来点的坐标是（(x0 + cos(PI / 2 +θ ))，（y0 + sin(PI / 2 + θ))），转化一下就是(x0 - sinθ，y0 + cosθ)。那么直接可以
求出dx = / (vp[i].y - vp[(i + 1) % nN].y) * fR / fDis，dy = (vp[(i + 1) % nN].x - vp[i].x) * fR / fDis，fDis是线段的长度。

代码如下：

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <algorithm>
#include <vector>
using namespace std;

const double fPre = 1e-8;

struct Point
{
    double x,y;
    Point(){}
    Point(const Point& p){x = p.x, y = p.y;}
    Point(double fX, double fY):x(fX), y(fY){}
    Point& operator+(const Point& p)
    {
        x += p.x;
        y += p.y;
        return *this;
    }
    Point& operator+=(const Point& p)
    {
        return *this = *this + p;
    }

    Point& operator-(const Point& p)
    {
        x -= p.x;
        y -= p.y;
        return *this;
    }
    Point& operator*(double fD)
    {
        x *= fD;
        y *= fD;
        return *this;
    }
};
typedef vector<Point> Polygon;
int DblCmp(double fD)
{
    return fabs(fD) < fPre ? 0 : (fD > 0 ? 1 : -1);
}

double Cross(Point a, Point b)
{
    return a.x * b.y - a.y * b.x;
}

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

Point Intersection(Point a1, Point a2, Point b1, Point b2)
{
    Point a = a2 - a1;
    Point b = b2 - b1;
    Point s = b1 - a1;
    return a1 + a * (Cross(b, s) / Cross(b, a));
}

Polygon Cut(Polygon& pg, Point a, Point b)
{
    Polygon pgRet;
    int nN = pg.size();

    for (int i = 0; i < nN; ++i)
    {
        double fC = Cross(a, b, pg[i]);
        double fD = Cross(a, b, pg[(i + 1) % nN]);

        if (DblCmp(fC) >= 0)
        {
            pgRet.push_back(pg[i]);
        }
        if (DblCmp(fC * fD) < 0)
        {
            pgRet.push_back(Intersection(a, b, pg[i], pg[(i + 1) % nN]));
        }
    }
    //printf("pgRet number:%d\n", pgRet.size());
    return pgRet;
}

double Dis(Point a, Point b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
//返回半平面的顶点个数
int HalfPlane(Polygon& vp, double fR)
{
    Polygon pg;
    pg.push_back(Point(-1e9, -1e9));
    pg.push_back(Point(1e9, -1e9));
    pg.push_back(Point(1e9, 1e9));
    pg.push_back(Point(-1e9, 1e9));
    int nN = vp.size();
    for (int i = 0; i < nN; ++i)
    {
        double fDis = Dis(vp[i], vp[(i + 1) % nN]);
        double dx = (vp[i].y - vp[(i + 1) % nN].y) * fR / fDis;
        double dy = (vp[(i + 1) % nN].x - vp[i].x) * fR / fDis;
        Point a = vp[i], b = vp[(i + 1) % nN], c(dx, dy);
        a += c;
        b += c;
        //printf("%f %f %f %f\n", a.x, a.y, b.x, b.y);
        pg = Cut(pg, a, b);
        if (pg.size() == 0)
        {
            return 0;
        }
    }
    return pg.size();
}

int main()
{
    int nN;
    vector<Point> vp;

    while (scanf("%d", &nN), nN)
    {
        vp.clear();
        Point p;
        for (int i = 0; i < nN; ++i)
        {
            scanf("%lf%lf", &p.x, &p.y);
            vp.push_back(p);
        }
        double fMin = 0.0, fMax = 10000.0;
        while (DblCmp(fMin - fMax))
        {
            double fMid = (fMin + fMax) / 2;
            int nRet = HalfPlane(vp, fMid);
            //printf("fMid:%f, nRet:%d\n", fMid, nRet);
            if (nRet == 0)
            {
                fMax = fMid;
            }
            else
            {
                fMin = fMid;
            }
        }
        printf("%.6f\n", fMax);
    }

    return 0;
}

posted @ 2012-07-23 16:45 yx 阅读(847) | 评论 (0) | 编辑收藏

poj 3130 How I Mathematician Wonder What You Are!

半平面交的一个题，也是求多边形的核心。求出这个好像也可以用于解决一些线性规划问题。我用的是N*N的基本算法，每加入一条直线，
就对原来求出的半平面交进行处理，产生新的核心。
代码参照台湾的一个网站演算法笔记上的内容和代码。表示这个网站巨不错，求凸包的算法也参照了这个网站上的内容和代码。
半平面交的地址：http://www.csie.ntnu.edu.tw/~u91029/Half-planeIntersection.html#a4

代码思路主要是：先读入所有的多边形顶点，放入一个vector（vp）里面，然后对多边形的每条边求一个半平面。刚开始的时候，用一个
vector（Polygon）保存代表上下左右四个无限远角的四个点，表示原始的半平面。然后，用读入的多边形的每条边去切割原来的半平面。
切割的过程是，如果原来（Polygon）中的点在当前直线的指定一侧，那么原来的点还是有效的。如果原来的点和它相邻的下一个点与当前
直线相交，那么还需要把交点加入Polygon集合。
还有求交点的方法比较奇葩，类似于黑书上面的那种根据面积等分的方法。

代码如下：

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <vector>
#include <algorithm>
using namespace std;

double fPre = 1e-8;
struct Point
{
    double x;
    double y;
    Point(){}
    Point(double fX, double fY)
    {
        x = fX, y = fY;
    }
};
typedef vector<Point> Polygon;
typedef pair<Point, Point> Line;
Point operator+(const Point& a, const Point& b)
{
    Point t;
    t.x = a.x + b.x;
    t.y = a.y + b.y;
    return t;
}

Point operator-(const Point& a, const Point& b)
{
    Point t;
    t.x = a.x - b.x;
    t.y = a.y - b.y;
    return t;
}

Point operator*(Point a, double fD)
{
    Point t;
    t.x = a.x * fD;
    t.y = a.y * fD;
    return t;
}

int DblCmp(double fD)
{
    return fabs(fD) < fPre ? 0 : (fD > 0 ? 1 : -1);
}

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}
//3点叉积
double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}
//向量叉积
double Cross(Point a, Point b)
{
    return a.x * b.y - a.y * b.x;
}

//求直线交点的一种简便方法
//平行四边形面积的比例等于高的比例
Point Intersection(Point a1, Point a2, Point b1, Point b2)
{
    Point a = a2 - a1;
    Point b = b2 - b1;
    Point s = b1 - a1;

    return a1 + a * (Cross(b, s) / Cross(b, a));
}

Polygon HalfPlane(Polygon& pg, Point a, Point b)
{
    Polygon pgTmp;
    int nN = pg.size();
    for (int i = 0; i < nN; ++i)
    {
        double fC = Cross(a, b, pg[i]);
        double fD = Cross(a, b, pg[(i + 1) % nN]);
        if (DblCmp(fC) >= 0)
        {
            pgTmp.push_back(pg[i]);
        }
        if (fC * fD < 0)
        {
            pgTmp.push_back(Intersection(a, b, pg[i], pg[(i + 1) % nN]));
        }
    }
    return pgTmp;
}

int main()
{
    int nN;
    Point p;
    vector<Point> vp;
    Polygon pg;

    while (scanf("%d", &nN), nN)
    {
        vp.clear();
        for (int i = 0; i < nN; ++i)
        {
            scanf("%lf%lf", &p.x, &p.y);
            vp.push_back(p);
        }
        pg.clear();
        pg.push_back(Point(-1e9, 1e9));
        pg.push_back(Point(-1e9, -1e9));
        pg.push_back(Point(1e9, -1e9));
        pg.push_back(Point(1e9, 1e9));
        for (int i = 0; i < nN; ++i)
        {
            pg = HalfPlane(pg, vp[i], vp[(i + 1) % nN]);
            if (pg.size() == 0)
            {
                printf("0\n");
                break;
            }
        }
        if (pg.size())
        {
            printf("1\n");
        }
    }

    return 0;
}

posted @ 2012-07-23 10:41 yx 阅读(1025) | 评论 (0) | 编辑收藏

poj 1584 A Round Peg in a Ground Hole

这个题需要多个计算几何算法。第一个是判断一系列点是否能够构成凸多边形，第二个是判断一个点是否在一个简单多边形内部，
第三个是求一个点到一条线段（或者说直线）的距离，第四个是判断一个圆是否则一个凸多边形内部。
其实，我是要判断一个圆是否则一个凸多边形内部而用到算法二和三。其实，有不需要判断圆心是否则多边形内部的算法。
算法一的思想，求所有边的偏转方向，必须都是逆时针或者顺时针偏转。算法二则是我前面发的那篇改进弧长法判断点和多边形的关系，
算法三尤其简单，直线上面取2点，用叉积求出这三点构成的三角形面积的2倍，再除以底边。算法四则是先判断圆心在多边形内部，然后
判断圆心到所有边的距离要大于圆的半径。
贴出代码，纯粹为了以后作为模版使用等，防止遗忘，方便查找，其实现在也能手敲出来了。

代码如下：

#include <stdio.h>

#include <string.h>
#include <math.h>
#include <algorithm>
#include <vector>
using namespace std;

const double fPre = 1e-8;
int DblCmp(double fD)
{
    if (fabs(fD) < fPre)
    {
        return 0;
    }
    else
    {
        return fD > 0 ? 1 : -1;
    }
}

struct Point
{
    double x, y;
    bool operator == (const Point& p)
    {
        return DblCmp(x - p.x) == 0 && DblCmp(y - p.y) == 0;
    }
};

Point operator-(const Point& a, const Point& b)
{
    Point p;
    p.x = a.x - b.x;
    p.y = a.y - b.y;
    return p;
}

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

bool IsConvexPolygon(vector<Point>& vp)
{
    int nN = vp.size();
    int nDirection = 0;
    bool bLine = true;//避免所有点共线
    for (int i = 0; i < nN; ++i)
    {
        int nTemp = DblCmp(Cross(vp[i], vp[(i + 1) % nN], vp[(i + 2) % nN]));
        if (nTemp)
        {
            bLine = false;
        }
        //这次的方向和上次的方向必须是相同的或者是3点和3点以上共线的情况
        if (nDirection * nTemp < 0)
        {
            return false;
        }
        nDirection = nTemp;
    }
    return bLine == false;
}

int GetQuadrant(Point p)
{
    return p.x >= 0 ? (p.y >= 0 ? 0 : 3) : (p.y >= 0 ? 1 : 2);
}

bool IsPtInPolygon(vector<Point>& vp, Point p)
{
    int nN = vp.size();
    int nA1, nA2, nSum = 0;
    int i;

    nA1 = GetQuadrant(vp[0] - p);
    for (i = 0; i < nN; ++i)
    {
        int j = (i + 1) % nN;
        if (vp[i] == p)
        {
            break;
        }
        int nC = DblCmp(Cross(p, vp[i], vp[j]));
        int nT1 = DblCmp((vp[i].x - p.x) * (vp[j].x - p.x));
        int nT2 = DblCmp((vp[i].y - p.y) * (vp[j].y - p.y));
        if (!nC && nT1 <= 0 && nT2 <= 0)
        {
            break;
        }
        nA2 = GetQuadrant(vp[j] - p);
        switch ((nA2 - nA1 + 4) % 4)
        {
            case 1:
                nSum++;
                break;
            case 2:
                if (nC > 0)
                {
                    nSum += 2;
                }
                else
                {
                    nSum -= 2;
                }
                break;
            case 3:
                nSum--;
                break;
        }
        nA1 = nA2;
    }

    if (i < nN || nSum)
    {
        return true;
    }
    return false;
}

double PtDis(Point a, Point b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (b.y - a.y) * (b.y - a.y));
}
//点p到直线ab的距离
//h = (2 * Spab) / |ab|
double GetDis(Point a, Point b, Point p)
{
    return fabs(Cross(a, b, p)) / PtDis(a, b);
}

bool IsCircleInPolygon(vector<Point>& vp, Point p, double fR)
{
    if (!IsPtInPolygon(vp, p))
    {
        return false;
    }

    int nN = vp.size();
    for (int i = 0; i < nN; ++i)
    {
        if (GetDis(vp[i], vp[(i + 1) % nN], p) < fR)
        {
            return false;
        }
    }
    return true;
}

int main()
{
    int nN;
    double fR, fPx, fPy;
    vector<Point> vp;
    Point p;

    while (scanf("%d%lf%lf%lf", &nN, &fR, &fPx, &fPy), nN >= 3)
    {
        vp.clear();
        for (int i = 0; i < nN; ++i)
        {
            scanf("%lf%lf", &p.x, &p.y);
            vp.push_back(p);
        }

        if (IsConvexPolygon(vp))
        {
            p.x = fPx;
            p.y = fPy;
            if (IsCircleInPolygon(vp, p, fR))
            {
                printf("PEG WILL FIT\n");
            }
            else
            {
                printf("PEG WILL NOT FIT\n");
            }
        }
        else
        {
            printf("HOLE IS ILL-FORMED\n");
        }
    }

    return 0;
}

posted @ 2012-07-20 21:41 yx 阅读(1031) | 评论 (0) | 编辑收藏

poj 1265 Area

此题用到了几个知识，一个是求多边形面积的公式。然后是，根据顶点都在整点上求多边形边界上的顶点数目的方法。最后一个是pick
定理。根据前面2个信息和pick定理算出在多边形内部的整点的个数。
求多边形面积的方法还是叉积代表有向面积的原理，把原点看做另外的一个点去分割原来的多边形为N个三角形，然后把它们的有向面
积加起来。
判断边界上点的个数是根据Gcd（dx,dy)代表当前边上整数点的个数的结论。这个结论的证明其实也比较简单，假设dx = a，dy = b。
初始点是x0，y0，假设d = Gcd（a，b）。那么边上的点可以被表示为（x0 + k * (a / d)，y0 + k * （b / d))。为了使点是整数点，
k必须是整数，而且0<= k <=d，所以最多有d个这个的点。
求多边形内部点的个数用的是pick定理。面积 = 内部点 + 边界点 / 2 - 1。

代码如下：

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
using namespace std;
#define MAX (100 + 10)

struct Point
{
    double x, y;
};
Point pts[MAX];

int nN;
const int IN = 1;
const int EAGE = 2;
const int OUT = 3;
const double fPre = 1e-8;

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

double GetArea()
{
    double fArea = 0.0;
    Point ori = {0.0, 0.0};

    for (int i = 0; i < nN; ++i)
    {
        fArea += Cross(ori, pts[i], pts[(i + 1) % nN]);
    }
    return fabs(fArea) * 0.5;
}

int gcd(int nX, int nY)
{
    if (nX < 0)
    {
        nX = -nX;
    }
    if (nY < 0)
    {
        nY = -nY;
    }
    if (nX < nY)
    {
        swap(nX, nY);
    }
    while (nY)
    {
        int nT = nY;
        nY = nX % nY;
        nX = nT;
    }
    return nX;
}

int main()
{
    int nT;
    int nI, nE;
    double fArea;

    scanf("%d", &nT);
    int dx ,dy;

    for (int i = 1; i <= nT; ++i)
    {
        scanf("%d", &nN);
        nI = nE = 0;
        pts[0].x = pts[0].y = 0;
        for (int j = 1; j <= nN; ++j)
        {
            scanf("%d%d", &dx, &dy);
            pts[j].x = pts[j - 1].x + dx;
            pts[j].y = pts[j - 1].y + dy;
            nE += gcd(dx, dy);
        }
        fArea = GetArea();
        nI = (fArea + 1) - nE / 2;
        printf("Scenario #%d:\n%d %d %.1f\n\n", i, nI, nE, fArea);
    }

    return 0;
}

posted @ 2012-07-20 15:36 yx 阅读(1031) | 评论 (0) | 编辑收藏

zoj 1081 Points Within - 改进弧长法判断点和简单多边形的关系

转角法判断点和多边形的关系大家都知道，原理比较简单，在多边形内扫过的转角一定是360度，在边界上和外面则不一定。
实现起来也比较麻烦，浮点误差比较大，而且还要考虑些特殊情况。
在网上找到一种叫做改进弧长法的算法，原理和转角法类似，但是做了很多重要的改进。比如，计算转角改成了计算叉积，根据叉积决定
旋转方向，还要根据计算下一个点的象限决定偏转多少，每次偏转的都是90度的倍数。
该算法可以方便判断出点在多边形内，还是边界上，还是在多边形外面。

摘自别人对该算法的描述如下：
首先从该书中摘抄一段弧长法的介绍：“弧长法要求多边形是有向多边形，一般规定沿多边形的正向，边的左侧为多边形的内侧域。
以被测点为圆心作单位圆，将全部有向边向单位圆作径向投影，并计算其中单位圆上弧长的代数和。若代数和为0，则点在多边形外部；
若代数和为2π则点在多边形内部；若代数和为π，则点在多边形上。”
按书上的这个介绍，其实弧长法就是转角法。但它的改进方法比较厉害：将坐标原点平移到被测点P，这个新坐标系将平面划分为4个
象限，对每个多边形顶点P ，只考虑其所在的象限，然后按邻接顺序访问多边形的各个顶点P，分析P和P[i+1]，有下列三种情况：
（1）P[i+1]在P的下一象限。此时弧长和加π/2；
（2）P[i+1]在P的上一象限。此时弧长和减π/2；
（3）P[i+1]在Pi的相对象限。首先计算f=y[i+1]*x-x[i+1]*y（叉积），若f=0，则点在多边形上；若f<0，弧长和减π；若f>0，弧长和加π。
最后对算出的代数和和上述的情况一样判断即可。
实现的时候还有两点要注意，第一个是若P的某个坐标为0时，一律当正号处理；第二点是若被测点和多边形的顶点重合时要特殊处理。

以上就是书上讲解的内容，其实还存在一个问题。那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错。
如边(3,0)-(-3,0)，按以上的处理，象限分别是第一和第二，这样会使代数和加π/2，有可能导致最后结果是被测点在多边形外。而实际上
被测点是在多边形上（该边穿过该点）。
对于这点，我的处理办法是：每次算P和P[i+1]时，就计算叉积和点积，判断该点是否在该边上，是则判断结束，否则继续上述过程。
这样牺牲了时间，但保证了正确性。
具体实现的时候，由于只需知道当前点和上一点的象限位置，所以附加空间只需O(1)。实现的时候可以把上述的“π/2”改成1，“π”改成2，
这样便可以完全使用整数进行计算。不必考虑顶点的顺序，逆时针和顺时针都可以处理，只是最后的代数和符号不同而已。整个算法编写
起来非常容易。

代码如下：

#include <stdio.h>
#include <math.h>

#define MAX (100 + 10)
struct Point
{
    double x,y;
};

Point pts[MAX];
const int OUT = 0;
const int IN = 1;
const int EDGE = 2;
const double fPre = 1e-8;

int DblCmp(double fD)
{
    if (fabs(fD) < fPre)
    {
        return 0;
    }
    else
    {
        return fD > 0 ? 1 : -1;
    }
}

int GetQuadrant(Point p)
{
    return DblCmp(p.x) >= 0 ? (DblCmp(p.y) >= 0 ? 0 : 3) :
               (DblCmp(p.y) >= 0 ? 1 : 2);
}

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

int PtInPolygon(Point* pts, int nN, Point p)
{
    int i, j, k;
    for (j = 0; j < nN; ++j)
    {
        pts[j].x -= p.x;
        pts[j].y -= p.y;
    }
    int nA1, nA2;
    int nSum = 0;
    nA1 = GetQuadrant(pts[0]);
    for (i = 0; i < nN; ++i)
    {
        k = (i + 1) % nN;
        if (DblCmp(pts[k].x) == 0 && DblCmp(pts[k].y) == 0)
        {
            break;//与顶点重合
        }
        int nC = DblCmp(Det(pts[i].x, pts[i].y,
                            pts[k].x, pts[k].y));
        if (!nC && DblCmp(pts[i].x * pts[k].x) <= 0
                && DblCmp(pts[i].y * pts[k].y) <= 0)
        {
            break;//边上
        }
        nA2 = GetQuadrant(pts[k]);
        if ((nA1 + 1) % 4 == nA2)
        {
            nSum += 1;
        }
        else if ((nA1 + 2) % 4 == nA2)
        {
            if (nC > 0)
            {
                nSum += 2;
            }
            else
            {
                nSum -= 2;
            }
        }
        else if ((nA1 + 3) % 4 == nA2)
        {
            nSum -= 1;
        }
        nA1 = nA2;
    }

    for (j = 0; j < nN; ++j)
    {
        pts[j].x += p.x;
        pts[j].y += p.y;
    }

    if (i < nN)
    {
        return EDGE;
    }
    else if (nSum)//逆时针nSum == 4, 顺时针nSum == -4
    {
        return IN;
    }
    else
    {
        return OUT;
    }
}

int main()
{
    int nN, nM;
    int nCase = 1;

    while (scanf("%d%d", &nN, &nM), nN)
    {
        if (nCase > 1)
        {
            printf("\n");
        }

        for (int i = 0; i < nN; ++i)
        {
            scanf("%lf%lf", &pts[i].x, &pts[i].y);
        }
        printf("Problem %d:\n", nCase++);
        for (int i = 0; i < nM; ++i)
        {
            Point p;
            scanf("%lf%lf", &p.x, &p.y);
            if (PtInPolygon(pts, nN, p))
            {
                printf("Within\n");
            }
            else
            {
                printf("Outside\n");
            }
        }
    }

    return 0;
}

posted @ 2012-07-20 10:30 yx 阅读(1505) | 评论 (0) | 编辑收藏

poj 3348 Cows

这个题用到2个计算几何算法。求解凸包和简单多边形面积。凸包算法详细解释见算法导论。求解多边形面积的思想是将多边形分解为三角
形，一般是假设按顺序取多边形上面连续的2点与原点组合成一个三角形，依次下去用叉积求有向面积之和，最后取绝对值即可。注意，这些
点必须是在多边形上逆时针或者顺时针给出的，而求出凸包刚好给了这样的一系列点。
凸包代码，其实先找出一个y坐标最小的点，再对剩下的所有点按极角排序。然后对排序后的点进行一个二维循环即可。二维循环的解释是
当加入新的点进入凸包集合时候，如果与以前加入的点形成的偏转方向不一致，那么前面那些点都需要弹出集合。

代码如下：

#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <math.h>
using namespace std;
#define MAX_N (10000 + 10)

struct Point
{
    double x, y;
    bool operator <(const Point& p) const
    {
        return y < p.y || y == p.y && x < p.x;
    }
};

Point pts[MAX_N];
int nN;
Point ans[MAX_N];
int nM;

double Det(double fX1, double fY1, double fX2, double fY2)
{
    return fX1 * fY2 - fX2 * fY1;
}

double Cross(Point a, Point b, Point c)
{
    return Det(b.x - a.x, b.y - a.y, c.x - a.x, c.y - a.y);
}

bool CrossCmp(const Point& a, const Point& b)
{
    double cs = Cross(pts[0], a, b);
    return cs > 0 || cs == 0 && a.x < b.x;
}

void Scan()
{
    nM = 0;
    sort(pts + 1, pts + nN, CrossCmp);//对所有点按极角排序,逆时针偏转

    //第0-2个点，其实不会进入第二重循环的
    //从第3个点开始，就依次检查与凸包中前面点形成的边的偏转方向
    //如果不是逆时针偏转，则弹出该点
    for (int i = 0; i < nN; ++i)
    {
        while (nM >= 2 && Cross(ans[nM - 2], ans[nM - 1], pts[i]) <= 0)
        {
            nM--;
        }
        ans[nM++] = pts[i];
    }
}

double GetArea()
{
    double fAns = 0.0;
    Point ori = {0.0, 0.0};
    for (int i = 0; i < nM; ++i)
    {
        fAns += Cross(ori, ans[i], ans[(i + 1) % nM]);
    }
    return fabs(fAns) * 0.5;
}

int main()
{
    while (scanf("%d", &nN) == 1)
    {
        for (int i = 0; i < nN; ++i)
        {
            scanf("%lf%lf", &pts[i].x, &pts[i].y);
            if (pts[i] < pts[0])
            {
                swap(pts[i], pts[0]);//pts[0]是y坐标最小的点
            }
        }

        Scan();//扫描出凸包
        double fArea = GetArea();
        printf("%d\n", (int)(fArea / 50));
    }

    return 0;
}

posted @ 2012-07-18 20:28 yx 阅读(1050) | 评论 (0) | 编辑收藏

poj 1730 Perfect Pth Powers

hdu 3627 Giant For

poj 2551 Ones and poj 2262 Goldbach's Conjecture

poj 2187 Beauty Contest

poj 3525 Most Distant Point from the Sea

poj 3130 How I Mathematician Wonder What You Are!

poj 1584 A Round Peg in a Ground Hole

poj 1265 Area

zoj 1081 Points Within - 改进弧长法判断点和简单多边形的关系

poj 3348 Cows

导航

统计

公告

常用链接

留言簿(3)

随笔分类

随笔档案

me

好友

同学

网友

搜索

最新评论

阅读排行榜

评论排行榜