poj 2115 C Looooops 线性同余方程

发表于 2012-07-18 分类于算法，算法题阅读次数： Valine：
本文字数： 1.2k 阅读时长 ≈ 1 分钟

这个题目就是解线性同余方程，(a + nc) % 2的k次 = b % 2的k次。既然以前是学信安的，对数论本来就不排斥，最近还好好看了下算法导论。这个方程转换为nc = (b-a) % 2的k次。根据数论的知识， ax = b%n，需要保证gcd(a,n)|b，意思b是gcd(a,n)的倍数，这个一下子也很难解释清楚啊，不满足这个条件，就是没解了。还有，如果有解的话，解的个数就是d = gcd(a,n)。而且其中一个解是x0 = x’(b/ d)，其中x’是用扩展欧几里德算法求出来的，满足关系式ax’+ny’=d。
但是这个题不仅仅用到数论的这些知识，因为必须求满足条件的最小解，而如果有解的话是d个，而且满足解x = x0 + i(b/d)，（1<=i<=d)。既然要求最小的解，那么对解mod(n/d)即可了，因为它们之间的差都是n/d的倍数。

代码如下：

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

//扩展欧几里德算法
//d = a * x + b * y，d是a和b的最大公约数
long long egcd(long long a, long long b, long long x, long long y)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    else
    {
        long long nRet = egcd(b, a % b, x, y);
        long long t = x;
        x = y;
        y = t - (a / b) * y;
        return nRet;
    }
}

int main()
{
    long long nA, nB, nC, nK;

    while (scanf("%I64d%I64d%I64d%I64d", &nA, &nB, &nC, &nK),
            nA || nB || nC || nK)
    {
        long long x, y;
        long long n = pow((double)2, (double)nK) + 1e-8;
        long long d = egcd(n, nC, x, y);
        long long b = (nB - nA + n) % n;
        if (b % d)//如果d | b失败
        {
            printf("FOREVER\n");
        }
        else
        {
            //printf("y:%I64d, b:%I64d, d:%I64d n:%I64d\n", y, b, d, n);
            y = (y + n) % n;
            long long ans = (y * (b / d)) % (n / d);
            printf("%I64d\n", ans);
        }
    }

    return 0;
}

poj 1401 Factorial

发表于 2012-07-16 分类于算法，算法题阅读次数： Valine：
本文字数： 713 阅读时长 ≈ 1 分钟

这个题目是求N！后面有多少个0，注意N可能最大到10的9次。哈哈，直接枚举1-N有多少个2和5的因子，然后取小的值肯定会超时的。但是，我还是试了下，果断超时了。
那就只有想数学结论了，果断想到1-N中能被2整除的数字有N / 2。哈哈，再往后思考下，发现1-N中能被4整除的数字有N / 4个，再往后就是N / 8，一直到N 除以2的某个次方为0为止，那么把所有的值加起来就是2的因子的个数了。求5的因子的个数也是这样的方法了。
很明显，5的因子的个数一定会小于等于2的因子的个数。那么直接求5的因子的个数就行了。由于，N / 5的时候用到了向下取整，所以不能用等比数列求和公式，怎么把答案弄成一个公式，还不知道了。
PS：其实我这种思路的灵感来自于筛选素数的方法了。

代码如下：

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

int GetAns(int nN)
{
    int nAns = 0;
    while (nN)
    {
        nAns += nN / 5;
        nN /= 5;
    }
    return nAns;
}

int main()
{
    int nT;

    scanf("%d", &nT);
    while (nT--)
    {
        int nN;
        scanf("%d", &nN);
        printf("%d\n", GetAns(nN));
    }

    return 0;
}

poj 2407 Relatives 欧拉函数

发表于 2012-07-15 分类于算法，算法题阅读次数： Valine：
本文字数： 1.2k 阅读时长 ≈ 1 分钟

这个题一看就知道是求欧拉函数。欧拉函数描述的正式题意。欧拉函数的理解可以按照算法导论上面的说法，对0-N-1进行筛选素数。那么公式n∏(1-1/p)，其中p是n的素数因子，就可以得到直观的理解了。但是计算的时候，会将这个式子变形下，得到另外一个形式。
如图所示：

但是这个题，需要考虑下，有可能n是个大素数，直接进行因子分解的话会超时的。怎么办了，只能在分解的时候判断n是不是已经成为素数了，如果是素数，答案再乘以n-1就行了。为了加快判断，我用5mb的空间搞了个素数表，大于5000000的数字只能循环判断了。

代码如下，注意求欧拉函数的代码部分：

#include <stdio.h>
#include <math.h>
#define MAX (5000000)
bool bPrime[MAX];//false表示素数

void InitPrime()
{
    bPrime[0] = bPrime[1] = true;
    int nMax = sqrt((double)MAX) + 1;
    for (int i = 2; i <= nMax; ++i)
    {
        if (!bPrime[i])
            for (int j = i * 2; j < MAX; j += i)
            {
                bPrime[j] = true;
            }
    }
}

bool IsPrime(int nN)
{
    if (nN < MAX)
    {
        return !bPrime[nN];
    }
    else
    {
        int nMax = sqrt((double)nN) + 1;
        for (int i = 2; i <= nMax; ++i)
        {
            if (nN % i == 0)
            {
                return false;
            }
        }
        return true;
    }
}

int main()
{
    int nN;

    InitPrime();
    while (scanf("%d", &nN), nN)
    {
        if (nN == 1)
        {
            printf("0\n");
            continue;
        }
        int nAns = 1;
        for (int i = 2; i <= nN; ++i)
        {
            if (IsPrime(nN))
            {
                nAns *= nN - 1;
                break;
            }
            if (nN % i == 0)
            {
                nAns *= i - 1;
                nN /= i;
                while (nN % i == 0)
                {
                    nAns *= i;
                    nN /= i;
                }
            }
        }
        printf("%d\n", nAns);
    }

    return 0;
}

poj 1730 Perfect Pth Powers

发表于 2012-07-13 分类于算法，算法题阅读次数： Valine：
本文字数： 1.1k 阅读时长 ≈ 1 分钟

通过这道题确实体会到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;
}

hdu 3627 Giant For

发表于 2012-07-11 分类于算法，算法题阅读次数： Valine：
本文字数： 2.1k 阅读时长 ≈ 2 分钟

这个题是对可排序数据的实时增加删除查找，那天做比赛的时候一点都不会，想来想去觉得平衡树可以做，但是写平衡树是件很难的事情。
后面知道线段数可以做，虽然数据的范围很大，但是可以在全部读入数据后排序再离散化，然后进行线段树的操作，具体的代码没有写。
今天队友在网上发现一种用map和set可以水掉这题的方法。原来，这个方法最主要的使用了map和set里面的upper_bound操作，以前居然忘记了这个东西了。既然这样，map和set也可以查前驱和后继了，但是注意low_bound查到的是小于等于的键。这个代码，注意是用了一个map< int, set > 集合把坐标都存起来了，进行添加删除和查找后继的操作。由于查找需要查找的元素是既比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;
}

poj 2551 Ones and poj 2262 Goldbach's Conjecture

发表于 2012-07-09 分类于算法，算法题阅读次数： Valine：
本文字数： 1.3k 阅读时长 ≈ 1 分钟

第一个题用到了同余的性质，这是数论里面最基本的性质，但是做题时候不一定能够自己发现。题意是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;
}

poj 2187 Beauty Contest

发表于 2012-07-07 分类于算法，算法题阅读次数： Valine：
本文字数： 3.6k 阅读时长 ≈ 3 分钟

这个题我是按照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尤其关于这个求最远点对的。

poj 3525 Most Distant Point from the Sea

发表于 2012-07-05 分类于算法，算法题阅读次数： Valine：
本文字数： 3.3k 阅读时长 ≈ 3 分钟

这个题的题意是给定一个凸多边形表示的海岛，求海岛离大海最远的距离。可以转化为一个凸多边形内部最多能够放入一个多大的圆。
显然可以对圆的半径进行二分，但是怎么确定圆心了。确定是否存在圆心，可以把原来的凸多边形往内部移动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;
}

poj 3130 How I Mathematician Wonder What You Are!

发表于 2012-07-03 分类于算法，算法题阅读次数： Valine：
本文字数： 2.7k 阅读时长 ≈ 2 分钟

半平面交的一个题，也是求多边形的核心。求出这个好像也可以用于解决一些线性规划问题。我用的是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;
}

poj 1584 A Round Peg in a Ground Hole

发表于 2012-06-30 分类于算法，算法题阅读次数： Valine：
本文字数： 3.4k 阅读时长 ≈ 3 分钟

这个题需要多个计算几何算法。第一个是判断一系列点是否能够构成凸多边形，第二个是判断一个点是否在一个简单多边形内部，第三个是求一个点到一条线段（或者说直线）的距离，第四个是判断一个圆是否则一个凸多边形内部。
其实，我是要判断一个圆是否则一个凸多边形内部而用到算法二和三。其实，有不需要判断圆心是否则多边形内部的算法。算法一的思想，求所有边的偏转方向，必须都是逆时针或者顺时针偏转。算法二则是我前面发的那篇改进弧长法判断点和多边形的关系，算法三尤其简单，直线上面取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;
}