Короче, вот этот код до 10^6 точек вроде бы считает, если больше, то возникает переполнение: // POLYGON.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define POINTS 1000000
struct point {
double x;
double y;
};
point polygon[POINTS];
int _tmain(int argc, _TCHAR* argv[])
{
double pi = 3.1415926535897932384626433832795;
double radius;
for (int i = 0; i < POINTS; i++)
{
radius = 2.0 + (double)rand()/RAND_MAX;
polygon[i].x = radius*cos(2*pi*i/POINTS);
polygon[i].y = radius*sin(2*pi*i/POINTS);
}
point intra = {0.5,0.5};
point extra = {3.5,4.5};
double delta;
double angle = 0.0;
for (int i = 0; i < POINTS-1; i++)
{
delta = acos(((polygon[i].x-intra.x)*(polygon[i+1].x-intra.x)+(polygon[i].y-intra.y)*(polygon[i+1].y-intra.y))/(sqrt((polygon[i].x-intra.x)*(polygon[i].x-intra.x)+(polygon[i].y-intra.y)*(polygon[i].y-intra.y))*sqrt((polygon[i+1].x-intra.x)*(polygon[i+1].x-intra.x)+(polygon[i+1].y-intra.y)*(polygon[i+1].y-intra.y))));
if (0 < (polygon[i].x-intra.x)*(polygon[i+1].y-intra.y)-(polygon[i+1].x-intra.x)*(polygon[i].y-intra.y))
{
angle += delta;
}
else
{
angle -= delta;
}
}
delta = acos(((polygon[POINTS-1].x-intra.x)*(polygon[0].x-intra.x)+(polygon[POINTS-1].y-intra.y)*(polygon[0].y-intra.y))/(sqrt((polygon[POINTS-1].x-intra.x)*(polygon[POINTS-1].x-intra.x)+(polygon[POINTS-1].y-intra.y)*(polygon[POINTS-1].y-intra.y))*sqrt((polygon[0].x-intra.x)*(polygon[0].x-intra.x)+(polygon[0].y-intra.y)*(polygon[0].y-intra.y))));
if (0 < (polygon[POINTS-1].x-intra.x)*(polygon[0].y-intra.y)-(polygon[0].x-intra.x)*(polygon[POINTS-1].y-intra.y))
{
angle += delta;
}
else
{
angle -= delta;
}
printf("area = %f\n",angle);
if (pi < angle) printf("intra point found!\n");
else printf("extra point found!\n");
angle = 0.0;
for (int i = 0; i < POINTS-1; i++)
{
delta = acos(((polygon[i].x-extra.x)*(polygon[i+1].x-extra.x)+(polygon[i].y-extra.y)*(polygon[i+1].y-extra.y))/(sqrt((polygon[i].x-extra.x)*(polygon[i].x-extra.x)+(polygon[i].y-extra.y)*(polygon[i].y-extra.y))*sqrt((polygon[i+1].x-extra.x)*(polygon[i+1].x-extra.x)+(polygon[i+1].y-extra.y)*(polygon[i+1].y-extra.y))));
if (0 < (polygon[i].x-extra.x)*(polygon[i+1].y-extra.y)-(polygon[i+1].x-extra.x)*(polygon[i].y-extra.y))
{
angle += delta;
}
else
{
angle -= delta;
}
}
delta = acos(((polygon[POINTS-1].x-extra.x)*(polygon[0].x-extra.x)+(polygon[POINTS-1].y-extra.y)*(polygon[0].y-extra.y))/(sqrt((polygon[POINTS-1].x-extra.x)*(polygon[POINTS-1].x-extra.x)+(polygon[POINTS-1].y-extra.y)*(polygon[POINTS-1].y-extra.y))*sqrt((polygon[0].x-extra.x)*(polygon[0].x-extra.x)+(polygon[0].y-extra.y)*(polygon[0].y-extra.y))));
if (0 < (polygon[POINTS-1].x-extra.x)*(polygon[0].y-extra.y)-(polygon[0].x-extra.x)*(polygon[POINTS-1].y-extra.y))
{
angle += delta;
}
else
{
angle -= delta;
}
printf("area = %f\n",angle);
if (pi < angle) printf("intra point found!\n");
else printf("extra point found!\n");
return 0;
}