- C20250002's blog
discuss搬运
-
C20250002 LV 5 @ 2026-1-12 17:56:53
由于discuss必须login才能看
故在此搬运:
原地址：Here
格点凸多边形 (lattice convex polygon)
(提什么意见都可以)
能不能用代码simply实现不重要，重要的是想法
work similar to this: James Paterson's blog
输入N, N<262144
输出格点凸N边形面积的最小值
不要求最佳，按你的输出与的真正的最小值 (记为F(N)) 的相对误差记分
(相对误差为x, 给-log x分)
References:
oeis:
paper:
这里不截取了，因为网上真的很难搜到，要专门去问人，搜到了也不是好结果
目前最好结论是由一位日本数学家得到的
He showed:
存在常数极限，大概为0.180
我之前的construction能够做到
大概是0.183吧
其实也是迄今最优的
做法是取一个原点为圆心的圆盘，去掉其中横纵坐标不互素的点，剩余的向量按斜率首尾相接得到这个多边形，形状大概是个圆形
(你猜我为什么放数论🧐)
我没想到什么better idea，就调整吧，规定方向优先级
但是这个调整有很多可能性：
-
一是我们可以生成一个大一点的，有一个去点策略（比如每次调整到局部最优去掉一个点，还可以使这个点位于没怎么被调整过的区域；或者一开始就均匀去点；或者两者兼具，比如每次调整到局部最优去掉一些点）
-
二是我们可以调整adjust函数，我写的是只考虑相邻5个点的，改进可以考虑相邻6个点然后再对中间两个找最优等等，performances are expected to be better
-
三是我们可以增加容忍策略，比如连续调整面积不变若干次后查询步长增加，直到调整使面积变化（因为我之前自己做过很多的奇奇妙妙的数论优化，很多都有这种improvement，容忍步长大概在log N到log log N间比较适配）
Update at 2026/1/12 17:32
实现了“做法是取一个原点为圆心的圆盘，去掉其中横纵坐标不互素的点”，而且应该已经很优了，详情可见make_circle及相关代码，to improve speed，不得不有些看起来像石山的东西
Update at 2026/1/16 17:39
实现了“剩余的向量按斜率”，想了一些优化出来(最好去看我代码sort_circle部分)，已经尽可能去掉浮点运算和利用图形结构减少比较次数了
在校图书馆机上sort 262144个顶点都能做到36ms 🤣
(空格和缩进很乱)
#include<math.h>
#include<iostream>
#include<time.h>
#include<iomanip>
#define sup 262144
#define version "v1.2.0"
#define debug 1
#define overload_num 3
#define tolerence_ring 1
using namespace std;
#define sp2 0.6079271018540267// 6 / pi^2
#define p2s 1.644934066848226// pi^2 / 6
#define ps 0.5235987755982988// pi / 6
#define ps_sqrt 0.723601254558267659363// sqrt(pi / 6)
#define half_sqrt 0.7071067811865475// sqrt(0.5)
#define x_len 371//1 + floor(sqrt(sup) * ps_sqrt)
#define y_len 262//1 + floor(x_len / sqrt(2))
long real_x_len = 0;
long real_y_len = 0;
long polygon[sup][2] = {{0}};//unused
long store_arrow[y_len] = {0};//start from 1
long sorted_vector[sup][2] = {{0}};//start from 0
bool circle[x_len][y_len] = {{0}};//both start from 1
int vert = 0;
long stable_amount = 0;
long adjust_amount = 0;
long circle_num = 0;
long circle_num_store = 0;
long long gcd(long long a, long long b) {//faster gcd
long long i = __builtin_ctz(a);
long long j = __builtin_ctz(b);
long long k = min(i, j);
long long d;
b >>= j;
while (a) {
a >>= i;
d = b - a;
i = __builtin_ctz(d);
if (a < b) b = a;
if (d < 0) a = -d;
else a = d;
}
return b << k;
}
bool is_coprime(long a, long b){//just to make it faster
if(min(a, b) == 1){return 1;}
if((a & 1) == 0 && (b & 1) == 0){return 0;}
if((a % 3) == 0 && (b % 3) == 0){return 0;}
if((a % 5) == 0 && (b % 5) == 0){return 0;}
return (1 == gcd(a, b));
}
int near(long double q, bool is_up){
if(is_up){return floor(q) + 1;}
else{return ceil(q) - 1;}
}
long double diff(long double q, bool is_up){
return abs(q - near(q, is_up));
}
short triangle_area(long a[2], long c[2], long e[2]){
// | 1 a b |
// | 1 c d |
// | 1 e f |
long long s = c[0]*e[1]-c[1]*e[0]-a[0]*e[1]+a[1]*e[0]+a[0]*c[1]-a[1]*c[0];
if(s > 0){return 1;}
if(s < 0){return -1;}
else{return 0;}
}
long double triangle_area_n(long a[2], long c[2], long e[2]){
return c[0]*e[1]-c[1]*e[0]-a[0]*e[1]+a[1]*e[0]+a[0]*c[1]-a[1]*c[0];
}
int modnum(int i){//addon
return (i + vert) % vert;
}
bool greater_than(long long a, long long b, long long c, long long d){
//a/b > c/d?
return ((a * d) > (b * c));
//integers are great
//i love integers uwu
}
void make_circle(){
long double rad = sqrt(vert) * ps_sqrt;
long double rad_sq = vert * ps;
long y = floor(rad / sqrt(2));
real_y_len = y;
real_x_len = rad;
long current_x = y;
for(int i = y; i > 0; i--){
while(current_x * current_x + i * i <= rad_sq){current_x ++;}
current_x--;
store_arrow[i] = current_x;
circle[current_x][i] = 1;
circle[i][i] = 1;
}
bool have_met = 0;
for(int j = real_y_len - 1; j > 0; j--){//here real_y_len-1 is to prevent that the first row have only one 1
for(int i = 1; i <= real_x_len; i++){
if(circle[i][j] && have_met){
circle[i][j] = is_coprime(i, j);
have_met = 0;
continue;
}
if(circle[i][j]){
have_met = 1;
circle[i][j] = is_coprime(i, j);
}
if(have_met && !circle[i][j] && is_coprime(i, j)){circle[i][j] = 1;}
}
}
circle[real_y_len][real_y_len] = (real_y_len <= 1);//forget to erase it lol
for(int i = 1; i < real_x_len; i++){
for(int j = 1; j < real_y_len; j++){
circle_num += (short)circle[i][j];
}
}
circle_num_store = circle_num;
}
void count_circle(){
long long cnt = 0;
cnt = circle_num_store << 3;
cout << "generated vertex count: " << cnt << endl;
}
void sort_circle(){
for(int i = 1; i < real_y_len; i++){//for safety lol
while(!circle[store_arrow[i]][i]){store_arrow[i]--;}
}
//it is a must to change float computation to int one
//or it will be too slow
//the flood-like algorithm, shifting from right to left
while(circle_num >= 0){
int min_y = 1;
for(int i = 1; i < real_y_len; i++){//where are the arrows? sort the flattest one out
if(!greater_than(i, store_arrow[i], min_y, store_arrow[min_y])){min_y = i;}
}
sorted_vector[circle_num_store - circle_num][0] = store_arrow[min_y];
sorted_vector[circle_num_store - circle_num][1] = min_y;
store_arrow[min_y] -= 1;
if((store_arrow[min_y] == min_y) && (min_y > 1)){
store_arrow[min_y] = 0;
circle_num--;
continue;
}
while(!circle[store_arrow[min_y]][min_y]){store_arrow[min_y]--;}//refresh state
circle_num--;
continue;
}
}
//now the major promotion must be drawn here and in the func auto_adjust one
void make_convex(/*int vert*/){
int quad = ceil(vert / 4.0 + overload_num);//overload protect
long len = quad * (quad + 1) / 2;
for(int i = 0; i < quad; i++){
polygon[i][0] = i * (i + 1) / 2;
polygon[i][1] = quad - i;
polygon[quad + i][0] = len + i;
polygon[quad + i][1] = i * (i + 1) / 2;
polygon[2 * quad + i][0] = len + quad - i * (i + 1) / 2;
polygon[2 * quad + i][1] = len + i;
polygon[3 * quad + i][0] = quad - i;
polygon[3 * quad + i][1] = len + quad - i * (i + 1) / 2;
}
for(int i = vert; i < 4 * quad; i++){//cancel overload
//something can be propelled here forwards the montgomery method
//by picking random points to eliminate then do the adjustments
polygon[i][0] = 0;
polygon[i][1] = 0;
}
}
void adjust(int vertnum){
//the trivial adjustment
//but it seemes always to work well awa
int v1 = modnum(vertnum - 2);
int v2 = modnum(vertnum - 1);
int v3 = modnum(vertnum);
int v4 = modnum(vertnum + 1);
int v5 = modnum(vertnum + 2);
long double k;
long double b;
long st[2];
//a slightly better soloution but it seems to be nonsense
//find the shorter edge of v2v4
long dist = min(abs(polygon[v2][0] - polygon[v4][0]), abs(polygon[v2][1] - polygon[v4][1]));
bool b_p = (dist != abs(polygon[v2][0] - polygon[v4][0]));
if(dist != 0){
//little confusing
//unkind for future debug
k = ((polygon[v2][!b_p] - polygon[v4][!b_p]) * 1.0) / ((polygon[v2][b_p] - polygon[v4][b_p]) * 1.0);
normal:
b = polygon[v2][!b_p] - polygon[v2][b_p] * k;
//[!b_p] = k * [b_p] + b
//is v3 upon v2v4?
long double len = (polygon[v3][!b_p] - k * polygon[v3][b_p] - b);
bool facing = (len > 0);
short test_1 = triangle_area(polygon[v1], polygon[v2], polygon[v3]);
short test_2 = triangle_area(polygon[v5], polygon[v4], polygon[v3]);
len = abs(len);
long left = min(polygon[v2][b_p], polygon[v4][b_p]);
long double area = triangle_area_n(polygon[v2], polygon[v3], polygon[v4]);
//long right = max(polygon[v2][b_p], polygon[v4][b_p]);
for(int i = 0; i < dist; i++){
st[b_p] = left + i;
st[!b_p] = near(k * st[b_p] + b, facing);
if((diff(k * st[b_p] + b, facing) < len)){
if((triangle_area(polygon[v1], polygon[v2], st) == test_1) && (triangle_area(polygon[v5], polygon[v4], st) == test_2)){
polygon[v3][0] = st[0];
polygon[v3][1] = st[1];
}
}
}
if(triangle_area_n(polygon[v2], polygon[v3], polygon[v4]) == area){//when auto stop
stable_amount++;
}
else{stable_amount = 0;}
}
if(dist == 0){
dist = max(abs(polygon[v2][0] - polygon[v4][0]), abs(polygon[v2][1] - polygon[v4][1]));
b_p = !b_p;
k = 0;
goto normal;
}
}
void ring_adjust(){
for(int i = 0; i < vert; i++){
adjust(i);
}
long stx = polygon[0][0], sty = polygon[0][1];
for(int i = 0; i < vert; i++){
stx = min(stx, polygon[i][0]);
sty = min(sty, polygon[i][1]);
}//do the shift
for(int i = 0; i < vert; i++){
polygon[i][0] -= stx;
polygon[i][1] -= sty;
}
}
void auto_adjust(){
while(stable_amount < tolerence_ring * vert){//no idea but to type sth like this uwu
adjust_amount++;
ring_adjust();
}
}
void output_polygon(){
for(int i = 0; i < vert; i++){
cout << polygon[i][0] << "," << polygon[i][1] << " ";
}
}
unsigned long long polygon_area() {
unsigned long long area = 0;
for (int i = 0; i < vert; i++) {
int j = (i + 1) % vert;
area += polygon[i][0] * polygon[j][1] - polygon[i][1] * polygon[j][0];
}
return area;
}
int main(){
if(!debug){
cout << "ver: " << version << endl;
cout << "input the vertex number of the polygon: " << endl;
cin >> vert;
make_convex();
// output_polygon();
long long tic = clock();
auto_adjust();
long long toc = clock();
// cout << endl;
output_polygon();
cout << endl;
cout << "adjust amount: " << adjust_amount << endl;
unsigned long long area = polygon_area();
cout << "polygon area: " << fixed << setprecision(1) << area / 2.0 << endl;
cout << "ratio: " << fixed << setprecision(8) << area / (2 * pow(vert, 3)) << endl;
cout << "time used: " << (toc - tic) / 1000 << " ms" << endl;
return 0;
}
else{
cout << "ver: " << version << endl;
cout << "testing sort_circle" << endl;
cout << "input the vertex number of the polygon: " << endl;
cin >> vert;
make_circle();
/*
for(int j = real_y_len; j > 0; j--){
for(int i = 1; i <= real_x_len; i++){
cout << circle[i][j] << " ";
}
cout << endl;
}
*/
count_circle();
long long tic = clock();
sort_circle();
long long toc = clock();
cout << "sort time: " << (toc - tic) / 1000 << " ms" << endl;
for(int i = 0; i <= circle_num_store; i++){
cout << "(" << sorted_vector[i][0] << ", " << sorted_vector[i][1] << "); ";
}
cout << endl;
return 0;
}
}
额，之所以有一段生成凸多边形的代码被我注释了，是因为它的performance太差了，经常调着调着就局部最优了，使我红温
马上就会有新的了！~
并不是很有时间把那段make_convex改成我说的那样了，for这可能很annoying
(真的没时间啦)
Anyway, 速度和结果都可能会存在大问题，这个确实不是很好写（吧？）
另：我不会模拟退火，别喷，而且感觉退火的代价也挺大的
欢迎有兴趣的前来看看! (或者comment some优化below也是极好的)
Link:
看第二块代码也可以