通過這道題確實體會到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;
}
這個題是對可排序數據的實時增加刪除查找，那天做比賽的時候一點都不會，想來想去覺得平衡樹可以做，但是寫平衡樹是件很難的事情。
后面知道線段數可以做，雖然數據的范圍很大，但是可以在全部讀入數據后排序再離散化，然后進行線段樹的操作，具體的代碼沒有寫。
今天隊友在網上發現一種用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;
}
第一個題用到了同余的性質，這是數論里面最基本的性質，但是做題時候不一定能夠自己發現。題意是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;
}
這個題我是按照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 <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://m.shnenglu.com/staryjy/archive/2010/09/25/101412.html 尤其關于這個求最遠點對的。
#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;
}
這個題的題意是給定一個凸多邊形表示的海島，求海島離大海最遠的距離。可以轉化為一個凸多邊形內部最多能夠放入一個多大的圓。
顯然可以對圓的半徑進行二分，但是怎么確定圓心了。確定是否存在圓心，可以把原來的凸多邊形往內部移動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;
}
半平面交的一個題，也是求多邊形的核心。求出這個好像也可以用于解決一些線性規劃問題。我用的是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;
}
這個題需要多個計算幾何算法。第一個是判斷一系列點是否能夠構成凸多邊形，第二個是判斷一個點是否在一個簡單多邊形內部，
第三個是求一個點到一條線段（或者說直線）的距離，第四個是判斷一個圓是否則一個凸多邊形內部。
其實，我是要判斷一個圓是否則一個凸多邊形內部而用到算法二和三。其實，有不需要判斷圓心是否則多邊形內部的算法。
算法一的思想，求所有邊的偏轉方向，必須都是逆時針或者順時針偏轉。算法二則是我前面發的那篇改進弧長法判斷點和多邊形的關系，
算法三尤其簡單，直線上面取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;
}
此題用到了幾個知識，一個是求多邊形面積的公式。然后是，根據頂點都在整點上求多邊形邊界上的頂點數目的方法。最后一個是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;
}
轉角法判斷點和多邊形的關系大家都知道，原理比較簡單，在多邊形內掃過的轉角一定是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;
}
這個題用到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;
}