FZU1120 A Pilot in Danger! [计算几何+数论]

毕设需要,在蛮早之前就写了这题,没时间记录一下,在这里学到的。

这个是大家耳熟能详的一个方法,用射线法去判断点是否在一个简单多边形内(边不自交,不论凹凸),如果是凸的话,可以有有向面积法,对于凹的还有累积度数法、编码法什么的,具体见计算机图形或是算法相关书籍,不赘述了。

后面加了个数论题,关乎于一个结论,在这里学到的,问说有多少个正整数不能表示成px+py的形式(x>=0,y>=0,p,q是素数且p!=q)。

1 点在多边形内的判断(射线法)

对于相交,考虑下面几个图:

Polygon

射线经过多边形的情况就这四种,第一种与点相交,第二种与边相交,与点相交可能还在内部,与边相交可能和边重合。对于第1、3、4个,可以用一种巧妙的技巧来处理下,可以使得判断变得简单,就是枚举边的时候,如果边的起始点在射线上的话,就不考虑该边,因此对于有公共顶点的两条边,就不会重复计算两次了,对于有点在射线上的话,如果两边的点在射线异端的话,就统计该店,否则不统计,累计后,判断统计的次数是奇数次,就说明是在多边形内部,否则就不在。对于第四中情况,上边的处理,等同于把中间的那些边都砍掉,退化成第一种的情况了。所以,说白了,就第一种和第二种两种情况,大致算法如下:

1. 枚举边

1.1. if 如果射线和线段直接相交,统计

1.2. elseif 如果射线和有向线段的起始点相交,就continue

1.3. elseif 如果射线和有向线段的终点相交的话,判断两端的点,如果在异侧,就统计

2. 判断统计数的奇偶

题目在这里,代码如下:

/*
判断点是否在多边形内
*/
#include <iostream>
#include <string>
#include <algorithm>
#include <vector>
#include <set>
#include <stdio.h>
#include <math.h>
using namespace std;
 
const int MAX = 100;
const double INF = 1e9;
const double EPS = 1e-6;
 
//是点也是向量
class CPoint
{
public:
    double x, y;
public:
    CPoint() {}
    CPoint(double _x, double _y) 
        : x(_x), y(_y) {}
    CPoint operator - (CPoint& t)  { return CPoint(t.x - x, t.y - y); }
    CPoint operator + (CPoint& t)  { return CPoint(t.x + x, t.y + y); }
 
    double ChaJi(CPoint t)  { return x * t.y - t.x * y; }
    double NeiJi(CPoint t)  { return x * t.x + y * t.y; }
};
 
class CLine
{
public:
    CPoint x, y;
public:
    CLine() {}
    CLine(CPoint _x, CPoint _y) 
        : x(_x), y(_y) {}
    double GetMinX()  { return min(x.x, y.x); }
    double GetMaxX()  { return max(x.x, y.x); }
    double GetMinY()  { return min(x.y, y.y); }
    double GetMaxY()  { return max(x.y, y.y); }
};
 
class CPolygon
{
public:
    CPoint point[MAX];
    CLine line[MAX];
    int line_num;
public:
    void SetLineByPoint()
    {
        for(int i = 0; i < line_num; i++)
        {
            line[i] = CLine(point[i], point[(i + 1) % line_num]);
        }
    }
    bool JudgePointIn(CPoint x)
    {
        int cnt = 0;
        CLine dq = CLine(x, CPoint(INF, x.y));
        
        for(int i = 0; i < line_num; i++)
        {
            if(PointOnLine(x, line[i]))  return true;
 
            if(LineXiangJiao(line[i], dq))  cnt++;
            else
            {
                if(PointOnLine(line[i].x, dq))  continue;
                else if(PointOnLine(line[i].y, dq)
                    && !OnTheSameSide(line[i].x, line[(i + 1) % line_num].y, dq))  cnt++;
                //else if(PointOnLine(line[i].y, dq)
                //    && PointOnLine(line[(i + 1) % line_num].y, dq)
                //    && !OnTheSameSide(line[i].x, line[(i + 2) % line_num].y, dq))  cnt++;
            }
        }
        
        if(cnt & 1)  return true;
        else  return false;
    }
    //精度控制
    bool IsZero(double x)  { return fabs(x) < EPS ? true : false; }
    bool BigThanZero(double x)  { return x >= EPS ? true : false; }
    bool SmaThanZero(double x)  { return x <= -EPS ? true : false; }
    //判断点是否在线段上
    bool PointOnLine(CPoint x, CLine L)
    {
        CPoint da = L.x - x;
        CPoint db = L.y - x;
 
        if(IsZero(da.ChaJi(db)) && SmaThanZero(da.NeiJi(db)))  return true;
        else  return false;
    }
    //两条线段相交,不包括端点
    bool LineXiangJiao(CLine x, CLine y)
    {
        //包围盒测试
        if(x.GetMaxX() < y.GetMinX() ||
            x.GetMaxY() < y.GetMinY() || 
            x.GetMinX() > y.GetMaxX() ||
            x.GetMinY() > y.GetMinY())  return false;
 
        if(!OnTheSameSide(x.x, x.y, y) &&
            !OnTheSameSide(y.x, y.y, x))  return true;
        else  return false;
    }
    //判断x,y是否在L的同侧
    bool OnTheSameSide(CPoint x, CPoint y, CLine L)
    {
        CPoint da, db, dc;
        
        da = L.y - L.x;
        db = x - L.x;
        dc = y - L.x;
 
        if(SmaThanZero(da.ChaJi(db) * da.ChaJi(dc)))  return false;
        else  return true;
    }
};
 
int main()
{
    int n, c = 0;
    while(scanf("%d", &n), n)
    {
        CPolygon m_Polygon;
        m_Polygon.line_num = n;
        for(int i = 0; i < n; i++)
        {
            double x, y;
            scanf("%lf%lf", &x, &y);
            m_Polygon.point[i] = CPoint(x, y);
        }
        m_Polygon.SetLineByPoint();
 
        //CPoint o;
        //scanf("%lf%lf", &o.x, &o.y);
 
        int p, q;
        scanf("%d%d", &p, &q);
 
        printf("Pilot %d\n", ++c);
 
        if(m_Polygon.JudgePointIn(CPoint(0, 0))) 
            printf("The pilot is in danger!\nThe secret number is %d.\n", (p - 1) * (q - 1) / 2);
        else  printf("The pilot is safe.\n");
        
        printf("\n");
    }
}

2 有多少个正整数不能表示为px+py的形式

趁这次机会,把以前学习的一些数论知识重新看了下,觉得比之前深入理解了些。

2.1 最大公约数

求最大公约数GCD(a, b),现在都能背着写出来了,但让我去推导,因为没有了解到本质,所以也不会对复杂度分析,就是常说的“只知其然,不知其所以然”,了解到算法的本质,才可以做到较好的灵活应用,为己所用,现在能记得以前学的东西里面,很大很大的一部分都是来自于推导也说明了这点。

GCD(a, b) = GCD(b, a % b) ,这里假设a > b,因为如果a < b的话,一次GCD之后,就变成了a > b,然后就一直这样下去。把该等式证明后,相信代码就好写了,就是把这两个偏序的数,一致的不断的变小(引自Crazyb0y的话),把规模缩小到可解,然后回朔就是。

我证明这个的方法比较笨,显得臃肿,就是分别证明GCD(a, b) | GCD(b, a % b)和GCD(b, a % b) | GCD(a, b)。这个过程就不详述了,说明其中一个证明过程,记g = GCD(a, b), a > b,然后证g | a % b,因为a = b * x + r, 0 <= r < b,r就是a % b,因为g | a, g | b,所以a = g * d1, b = g * d2,d1, d2 >= 1,a % b = r = a – b * x = g * d1 – g * d2 * x = g * (d1 – d2 * x),所以g | r = (a % b)。

这个求最大公约数的算法就叫做辗转相除法(Euclidean Algorithm),然后这里简单分析它的复杂度(来自Crazyb0y的集训队论文),考虑(a, b) –> (b, a % b)这样一个步骤,它把a变成了a % b,然后下次用同样的手法处理了b,如果b <= a / 2,则a % b < a / 2,如果b > a / 2,则a % b <= a – b < a / 2,所以它把一个数给砍半了,所以复杂度就为log(max(a, b)),就是logN级别的了。

2.2 ax+by=g

将辗转相除法扩展一下就是扩展欧几里得算法(Extended Euclidean Algorithm),下面说明下,如何求出该类方程的一组解。

ax + by = g, g = gcd(a, b),对该方程做下变形,目标是把系数做成辗转相除法中的方法,让它变小,到可解的时候,解出解后回朔,下面用到的除法是正除以。

(a / b * b + a % b)x + by = g  =>  b(a / b * x + y) + (a % b)x +  = g

设x’ = a / b * x + y, y’ = x然后递归求出右边等式bx’ + (a % b)y’ = g的解,根据x’,y’求出x,y即可,当a % b = 0的时候,(x’,y’)=(1,0)便是一组解。

下面贴一下欧几里得算法和扩展欧几里得算法的代码:

#include <iostream>
#include <string>
#include <algorithm>
#include <vector>
#include <set>
#include <stdio.h>
#include <math.h>
using namespace std;
 
int Euclidean_Algorithm(int a, int b)    //辗转相除法
{
    if(b == 0)  return a;
    else  return Euclidean_Algorithm(b, a % b);        //使得a > b
}
 
void Extended_Euclidean_Algorithm(int a, int b, int& x, int& y)    //ax+by=gcd(a,b)=g
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return;
    }
 
    int xx, yy;
    Extended_Euclidean_Algorithm(b, a % b, xx, yy);
 
    x = yy;
    y = xx - a / b * yy;
}
 
int main()
{
    int a, b;
    while(scanf("%d%d", &a, &b))
    {
        int x, y;
        Extended_Euclidean_Algorithm(a, b, x, y);
        printf("%d %d %d\n", Euclidean_Algorithm(a, b),
            x, y);
    }
}

下面就是关于这道数论的证明,参考这儿

证明,如果GCD(p, q) = 1,p>=1,q>=1,则不能表示为px+qy,x>=0,y>=0的非负整数有(p-1)*(q-1)/2个,下面一步一步证明,下面的变量如无特别说明,都是整数,

(1) 不能表示为px+qy,x>=0,y>=0,(p,q)=1,p>=1,q>=1的最大正整数是pq-p-q;

先证明pq-p-q不能表示为所述形式:

  假设能表示,则pq-p-q=px+qy,

  则p(q-1-x)=q(y+1),

  则q | (q-1-x), p | (y+1)

  则q | (1+x),因为q(y+1)>0,所以q-1-x>0,所以1<=x+1<=q-1,所以产生矛盾,pq-p-q不能表示为上述形式。

再证明pq-p-q是最大的:

  先看这样的一个序列1p,2p,3p,…,(q-1)p,这q个数%p后的值是在{0,1,2,q-1}中,因为如果有两个是相等的,设ip=jp(mod q),0<=i<j<=q-1,则(j-i)p%q=0,则(j-i)%q=0,而0<j-i<=q-1,所以产生矛盾,说明上述是正确的。

  同理n-1p,n-2p,n-3p,…,n-(q-1)p,这个序列也是一样,所以设(n-up)%q=0,则n=up+vq,u∈[0,q-1],

  n=up+vq>pq-p-q,则vq>pq-p-q-up>=pq-p-q-pq+p,推出v>-1,即v>=0,说明对于n是能表示为上述形式的。

(2) 对于任意一个非负整数n,都可表示为n=px+qy,的形式,其中x∈[0, q – 1],y为整数即可;

这个在(1)中已经给出证明,按照那种方法就可以构造出来,这里可以如此调整出来,

知道n=px+qy=p(x-tq)+q(y+tp),可以通过调整x-tq,把它调整至[0,q-1]这个区间内。

(3) m=pq-p-q,对于n和m-n,n∈[0, m],n和m-n这两个数,有且仅有一个能表示为px+qy的形式,(p,q)=1,p>=1,q>=1,x>=0,y>=0;

对于n和m-n都想表示成那种形式,为使u1,u2,v1,v2属于[0,+∞)为此,

n=pu1+qv1, m-n=pu2+qv2,这里u1,u2∈[0, q-1],因为pu1+qv1<=pq-p-q,所以pu1<=pq-p-q-qv1<=pq-p,所以u1<=q-1,同理u2。

所以pq-p-q=pu1+pu2+qv1+qv2,则p(q-1-u1-u2)=q(v1+v2+1),

所以p | (v1+v2+1), q | (q-1-u1-u2),

所以q | (u1+u2+1),因为1<=u1+u2<=2q-1,所以u1+u2+1=q,所以v1+v2=-1,所以v1,v2不能同时>=0,也不能同时<0。只能一个符合要求,一个不符合要求。

(4) 不能表示为px+qy,x>=0,y>=0,(p,q)=1,p>=1,q>=1的整数有(p-1)*(q-1)/2个。

对于[0, pq-p-q]这个区间的pq-p-q+1个数(偶数),能进行和为pq-p-q的配对,共有(p-1)*(q-1)/2对,这些对里面,有且仅有一个是能表示为,另一个是不能表示的,而对于n>pq-p-q,前面已经证明是可以表示为的,所以不能表示为的,共有(p-1)*(q-1)/2个。

原文地址:https://www.cnblogs.com/litstrong/p/2022497.html