2019 ICPC Asia-East Continent Final C. Dirichlet k-th root
定义 ∗ * ∗ 为狄利克雷卷积， ( f ∗ g ) ( n ) = ∑ d ∣ n f ( d ) g ( n d ) (f*g)(n)=\sum_{d|n}f(d)g(\frac nd) (f∗g)(n)=∑d∣nf(d)g(dn)，给出 n n n， f k f^k fk， k k k，求 f f f
n ≤ 1 0 5 , m o d = 998244353 , k < m o d n\le 10^5,\bmod =998244353, k< \bmod n≤105,mod=998244353,k<mod，保证 g ( 1 ) = 1 g(1)=1 g(1)=1
题解：
结论：对于
i
<
2
m
o
d
i<2^{mod}
i<2mod，
f
m
o
d
(
i
)
=
[
i
=
=
1
]
∗
f
(
1
)
f^{mod}(i)=[i==1]*f(1)
fmod(i)=[i==1]∗f(1)
所以 f = f k ∗ i n v ( k ) f=f^{k*inv(k)} f=fk∗inv(k)， i n v ( k ) inv(k) inv(k) 是 k k k 在 m o d 998244353 \bmod~998244353 mod 998244353 下的逆元。
证明：
对于 1 < i < 2 m o d 1<i<2^{mod} 1<i<2mod， f m o d ( i ) = ∑ i 1 ∗ i 2 ∗ . . . ∗ i m o d = i f ( i 1 ) f ( i 2 ) . . . f ( i m o d ) f^{mod}(i)=\sum\limits_{i_1*i_2*...*i_{mod}=i}f(i_1)f(i_2)...f(i_{mod}) fmod(i)=i1∗i2∗...∗imod=i∑f(i1)f(i2)...f(imod)，那么记 i x ≠ 1 i_x\neq1 ix=1 的 x x x 的个数为 p p p，选出这种组合的方案数是 C m o d p C_{mod}^p Cmodp，而对于 0 < p < m o d 0<p<mod 0<p<mod， C m o d p ≡ 0 ( m o d m o d ) C_{mod}^p\equiv 0\pmod {mod} Cmodp≡0(modmod)。
对于 i = 1 i=1 i=1， f m o d ( 1 ) = f ( 1 ) m o d ≡ f ( 1 ) f^{mod}(1)=f(1)^{mod}\equiv f(1) fmod(1)=f(1)mod≡f(1)
复杂度 O ( n ln n log m o d ) O(n\ln n\log \bmod) O(nlnnlogmod)
Code：
#include<bits/stdc++.h>
#define maxn 100005
using namespace std;
const int mod = 998244353;
int n,k,f[maxn],g[maxn];
int Pow(int a,int b){int s=1;for(;b;b>>=1,a=1ll*a*a%mod) b&1&&(s=1ll*s*a%mod);return s;}
void Mul(int *A,int *B,int *c){
static int a[maxn],b[maxn];
memcpy(a,A,(n+1)<<2),memcpy(b,B,(n+1)<<2),memset(c,0,(n+1)<<2);
for(int i=1;i<=n;i++) for(int j=1;i*j<=n;j++) c[i*j]=(c[i*j]+1ll*a[i]*b[j])%mod;
}
void Pow(int *a,int *s,int b){
s[1]=1;
for(;b;b>>=1,Mul(a,a,a)) if(b&1) Mul(a,s,s);
}
int main()
{
scanf("%d%d",&n,&k);
for(int i=1;i<=n;i++) scanf("%d",&g[i]);
Pow(g,f,Pow(k,mod-2));
for(int i=1;i<=n;i++) printf("%d%c",f[i],i==n?10:32);
}
然而此题还有更优的解法。
2018 ACM-ICPC World Finals H. Single Cut of Failure
一个 W ∗ H W*H W∗H 的矩形，给出一些端点在边上（不在交点且各不相同且不与矩形边重合）的线段，求最少用几条同样条件的线段可以切割所有给出的线段。输出方案。
每条线都与主对角线和副对角线中的某一条相交，所以答案最多为2，输出方案时偏离0.5即可。
判定是否能只用1条：将矩形看成排成一排的四条线段，给出的线段看做一些区间，相当于找到一个区间包含
n
n
n个属于不同线段的端点，two-pointer扫描即可。
因为不能让区间在矩形的同一边，所以右端点要尽量靠右，同时末端要加上一个点。
Code：
#include<bits/stdc++.h>
#define maxn 1000005
using namespace std;
int n,W,H;
bool vis[maxn];
struct node{
int x,id;
bool operator < (const node &p)const{return x<p.x;}
}a[maxn*2];
int turn(int x,int y){
return y==0?x:x==W?W+y:y==H?2*W-x+H:2*W+2*H-y;
}
void write(double x){
if(x<W) printf("%.1f %d",x,0);
else if(x<W+H) printf("%d %.1f",W,x-W);
else if(x<2*W+H) printf("%.1f %d",2*W+H-x,H);
else printf("%d %.1f",0,2*W+2*H-x);
}
void solve(){
for(int s=1,t=1;t<=2*n;s++){
while(t<=2*n&&!vis[a[t].id]) vis[a[t].id]=1,t++;
if(t-s==n) puts("1"),write(a[s].x-0.5),putchar(' '),write(a[t].x-0.5),exit(0);
vis[a[s].id]=0;
}
}
int main()
{
scanf("%d%d%d",&n,&W,&H);
for(int i=1,x,y;i<=n;i++){
scanf("%d%d",&x,&y),a[i*2-1]=(node){turn(x,y),i};
scanf("%d%d",&x,&y),a[i*2]=(node){turn(x,y),i};
}
a[2*n+1].x=2*W+2*H;
sort(a+1,a+1+2*n);
solve();
puts("2");
printf("%.1f %d %.1f %d\n",0.5,0,W-0.5,H);
printf("%d %.1f %d %.1f\n",W,0.5,0,H-0.5);
}
2018 ACM-ICPC World Finals I. Triangles
给出一个酱紫的图，数三角形个数。
先数
△
\triangle
△型，反过来同理。
数以每个点为右下角的
△
\triangle
△ 个数，记录每个点往左
L
L
L，往左上
U
U
U，往右上
D
D
D，能延伸多长。
对每一行从左往右扫，在
x
x
x处加1，扫到
x
+
D
x+D
x+D时将
x
x
x处的1减去，求
[
x
−
min
(
U
,
L
)
,
x
)
[x-\min(U,L),x)
[x−min(U,L),x)中为1的个数。
一条右上的斜线相当于“覆盖了”它往右的一部分点，可以对它们做贡献。
Code：
#include<bits/stdc++.h>
#define maxn 6005
#define maxm 12005
#define LL long long
using namespace std;
int n,m,U[maxn][maxm],D[maxn][maxm];
char a[maxn*2][maxm*2];
vector<int>G[maxm];
int arr[maxm];
void upd(int i,int v){for(;i;i-=i&-i) arr[i]+=v;}
int qsum(int i){int s=0;for(int lim=(m+1)/2;i<=lim;i+=i&-i) s+=arr[i];return s;}
int main()
{
freopen("1.in","r",stdin);
scanf("%d%d\n",&n,&m);
for(int i=1;i<=2*n-1;i++) fgets(a[i]+1,2*m+1,stdin);
LL ans=0;
for(int i=1;i<=2*n-1;i+=2){
int Mx=0;
for(int j=i/2&1?3:1,L=0,k=1;j<=2*m-1;j+=4,k++){
U[i][j] = a[i-1][j-1]=='\\' ? U[i-2][j-2]+1 : 0;
D[i][j] = a[i-1][j+1]=='/' ? D[i-2][j+2]+1 : 0;
L = a[i][j-1]=='-' ? L+1 : 0;
ans += qsum(k-min(U[i][j],L));
upd(k,1),G[k+D[i][j]].push_back(k),Mx=max(Mx,k+D[i][j]);
while(!G[k].empty()) upd(G[k].back(),-1),G[k].pop_back();
}
for(int k=m/2+1;k<=Mx;k++) while(!G[k].empty()) upd(G[k].back(),-1),G[k].pop_back();
}
for(int i=2*n-1;i>=1;i-=2){
int Mx=0;
for(int j=i/2&1?3:1,L=0,k=1;j<=2*m-1;j+=4,k++){
U[i][j] = a[i+1][j-1]=='/' ? U[i+2][j-2]+1 : 0;
D[i][j] = a[i+1][j+1]=='\\' ? D[i+2][j+2]+1 : 0;
L = a[i][j-1]=='-' ? L+1 : 0;
ans += qsum(k-min(U[i][j],L));
upd(k,1),G[k+D[i][j]].push_back(k),Mx=max(Mx,k+D[i][j]);
while(!G[k].empty()) upd(G[k].back(),-1),G[k].pop_back();
}
for(int k=m/2+1;k<=Mx;k++) while(!G[k].empty()) upd(G[k].back(),-1),G[k].pop_back();
}
printf("%lld\n",ans);
}
2018 ACM-ICPC World Finals C. Conquer the World
树上每个位置有一些军队，给出每个位置需要的军队数量，有边权，求移动军队使得满足条件的最小代价。
费用流最小权匹配，树上转化为模拟费用流的配对过程，用堆维护，详细题解见 cz_xuyixuan‘s blog
正确性：同一棵子树里面的对如果和小于0那么必然有一个会被选，所以提供反悔可能之后就可以提前选。
Code：
#include<bits/stdc++.h>
#define maxn 250005
#define maxp maxn*20
#define LL long long
using namespace std;
char cb[1<<20],*cs,*ct;
#define getc() (cs==ct&&(ct=(cs=cb)+fread(cb,1,1<<20,stdin),cs==ct)?0:*cs++)
void read(int &a){
char c;while(!isdigit(c=getc()));
for(a=c-'0';isdigit(c=getc());a=a*10+c-'0');
}
const LL inf = maxn*1e6;
int n,X[maxn],Y[maxn];
LL ans,dep[maxn];
int fir[maxn],nxt[maxn<<1],to[maxn<<1],w[maxn<<1],tot;
void line(int x,int y,int z){nxt[++tot]=fir[x],fir[x]=tot,to[tot]=y,w[tot]=z;}
struct node{
LL v; int s;
}a[maxp];
int rt0[maxn],rt1[maxn],lc[maxp],rc[maxp],sz,dis[maxp];
int newnode(node p){return a[++sz]=p,sz;}
void merge(int &i,int x,int y){
if(!x||!y) return void(i=x+y);
if(a[y].v<a[x].v) swap(x,y);
i=x,merge(rc[i],rc[x],y);
if(dis[lc[i]]<dis[rc[i]]) swap(lc[i],rc[i]);
dis[i]=dis[rc[i]]+1;
}
void Push(int &i,int x){merge(i,i,x);}
void Pop(int &i){merge(i,lc[i],rc[i]);}
void solve(int u,int v,LL d){
for(;rt0[u]&&rt1[v];){
LL x=a[rt0[u]].v,y=a[rt1[v]].v;
if((x+y-2*d)>=0) return;
int s=min(a[rt0[u]].s,a[rt1[v]].s);
ans+=(x+y-2*d)*s;
if(!(a[rt0[u]].s-=s)) Pop(rt0[u]);
if(!(a[rt1[v]].s-=s)) Pop(rt1[v]);
Push(rt0[u],newnode((node){2*d-y,s}));
Push(rt1[v],newnode((node){2*d-x,s}));
}
}
void dfs(int u,int ff){
if(X[u]>0) Push(rt0[u],newnode((node){dep[u],X[u]}));
if(Y[u]>0) Push(rt1[u],newnode((node){dep[u]-inf,Y[u]})),ans+=Y[u]*inf;
for(int i=fir[u],v;i;i=nxt[i]) if((v=to[i])!=ff){
dep[v]=dep[u]+w[i],dfs(v,u),solve(u,v,dep[u]),solve(v,u,dep[u]);
merge(rt0[u],rt0[u],rt0[v]),merge(rt1[u],rt1[u],rt1[v]);
}
}
int main()
{
freopen("1.in","r",stdin);
read(n);
for(int i=1,x,y,z;i<n;i++) read(x),read(y),read(z),line(x,y,z),line(y,x,z);
for(int i=1,t;i<=n;i++) read(X[i]),read(Y[i]),t=min(X[i],Y[i]),X[i]-=t,Y[i]-=t;
dfs(1,0);
printf("%lld\n",ans);
}
2018 ACM-ICPC World Finals D. Gem island
n n n个人，每个人初始有1个宝石，每过一天，在已有的宝石中等概率会有某个宝石变成两个，问 d d d 天后拥有宝石数前 r r r 大的人拥有的宝石数的和的期望。 n , d , r ≤ 500 n,d,r\le500 n,d,r≤500
去掉最开始的那一颗，假设每个人最后多获得了
a
i
a_i
ai 颗宝石，
∑
a
i
=
d
\sum a_i=d
∑ai=d，考虑这种情况发生的概率：
因为每天每个宝石分裂的概率是相同的，所以总情况数为
n
∗
(
n
+
1
)
∗
.
.
.
∗
(
n
+
d
−
1
)
n*(n+1)*...*(n+d-1)
n∗(n+1)∗...∗(n+d−1)
最后宝石的分布列是
{
a
i
}
\{a_i\}
{ai}，考虑
d
d
d 天中每天哪个人的宝石分裂，可重排列方案数是
d
!
∏
a
i
!
\frac {d!}{\prod a_i!}
∏ai!d!，然后每个人的宝石分裂的情况数是
1
∗
2
∗
.
.
.
∗
a
i
1*2*...*a_i
1∗2∗...∗ai，总情况数为
d
!
d!
d!
所以分布列为
{
a
i
}
\{a_i\}
{ai} 的概率为
d
!
(
n
+
d
−
1
)
!
(
n
−
1
)
!
\frac {d!(n+d-1)!}{(n-1)!}
(n−1)!d!(n+d−1)!，与具体的
a
i
a_i
ai 分布无关。
因为 d d d 是定值，所以最终每种分布列的概率都是相同的，我们只需要算出每种分布列的前 r r r 大之和累加之后除以总的分布方案数就可以了。
设 c n t i , j cnt_{i,j} cnti,j 表示给 i i i 个人分配 j j j 个宝石的方案数， s u m i , j sum_{i,j} sumi,j 表示给 i i i 个人分配 j j j 个宝石的所有方案中前 r r r 大之和的和。
因为要求前
r
r
r 大之和，所以我们DP的方式为“给一些数整体+1"，具体的，枚举给
k
k
k 个数+1，表示这
k
k
k个数比其它的
i
−
k
i-k
i−k 个数大：
c
n
t
i
,
j
=
∑
k
=
0
min
(
i
,
j
)
(
i
k
)
∗
c
n
t
k
,
j
−
k
s
u
m
i
,
j
=
∑
k
=
0
min
(
i
,
j
)
(
i
k
)
∗
(
s
u
m
k
,
j
−
k
+
min
(
k
,
r
)
∗
c
n
t
k
,
j
−
k
)
cnt_{i,j}=\sum_{k=0}^{\min(i,j)}\binom ik*cnt_{k,j-k}\\ sum_{i,j}=\sum_{k=0}^{\min(i,j)}\binom ik*(sum_{k,j-k}+\min(k,r)*cnt_{k,j-k})
cnti,j=k=0∑min(i,j)(ki)∗cntk,j−ksumi,j=k=0∑min(i,j)(ki)∗(sumk,j−k+min(k,r)∗cntk,j−k)
Code：
#include<bits/stdc++.h>
#define maxn 505
using namespace std;
int n,d,r;
double cnt[maxn][maxn],sum[maxn][maxn],c[maxn][maxn];
int main()
{
scanf("%d%d%d",&n,&d,&r);
for(int i=0,lim=max(n,d);i<=lim;i++)
for(int j=(c[i][0]=1);j<=i;j++)
c[i][j]=c[i-1][j-1]+c[i-1][j];
cnt[0][0]=1;
for(int i=1;i<=n;i++)
for(int j=0;j<=d;j++)
for(int k=0;k<=i&&k<=j;k++){
cnt[i][j]+=c[i][k]*cnt[k][j-k];
sum[i][j]+=c[i][k]*(sum[k][j-k]+min(k,r)*cnt[k][j-k]);
}
printf("%.8f\n",sum[n][d]/cnt[n][d]+r);
}
2018 ACM-ICPC World Finals E. Getting a Jump on Crime
一个
d
x
∗
d
y
d_x*d_y
dx∗dy 的网格图，每个格子的是底面为
w
∗
w
w*w
w∗w，高度为
h
x
,
y
h_{x,y}
hx,y 的长方体，给出速度大小
v
v
v，可以朝任何方向。起点在
l
x
,
l
y
l_x,l_y
lx,ly 的方格的中心，问到达每个格子最少跳几次，跳的抛物线不能穿过任何一个格子，重力加速度
g
=
9.80665
m/s
2
g=9.80665 \ \text{m/{s}}^2
g=9.80665 m/s2。
d
x
,
d
y
≤
20
,
d_x,d_y\le20,
dx,dy≤20,
思路比较简单，就是BFS，然后检验两个点能否到达。
因为横纵坐标之差和速度都是定的，所以抛物线可以解出来：
v
x
t
=
Δ
x
v
y
t
−
1
2
g
t
2
=
Δ
y
v
x
2
+
v
y
2
=
v
2
v_xt=\Delta_x\\ v_yt-\frac 12gt^2=\Delta_y\\ v_x^2+v_y^2=v^2
vxt=Δxvyt−21gt2=Δyvx2+vy2=v2
将
v
x
,
v
y
v_x,v_y
vx,vy用
t
t
t表示代入3式：
1
4
g
2
t
4
+
(
g
Δ
y
−
v
2
)
t
2
+
(
Δ
x
2
+
Δ
y
2
)
=
0
\frac 14g^2t^4+(g\Delta_y-v^2)t^2+(\Delta_x^2+\Delta_y^2)=0
41g2t4+(gΔy−v2)t2+(Δx2+Δy2)=0
t
t
t越大，
v
x
v_x
vx越小，
v
y
v_y
vy越大，抛物线曲线越高，所以解出
t
t
t 的最大值。
然后枚举横坐标，找到竖着的边界及其两边对应格子，算出穿过这个边界时的 y y y 即可，横着的边界同理。
Code：
#include<bits/stdc++.h>
#define maxn 25
using namespace std;
const double G = 9.80665;
int n,m,sx,sy,dis[maxn][maxn];
double W,V,h[maxn][maxn];
double sqr(double x){return x*x;}
bool check(int x,int y,int tx,int ty){
double d = W*sqrt(sqr(tx-x)+sqr(ty-y)), A = sqr(G)/4, B = (G*(h[tx][ty]-h[x][y])-V*V), C = sqr(d)+sqr(h[tx][ty]-h[x][y]);
if(B*B-4*A*C<0) return 0;
double t = sqrt((-B+sqrt(B*B-4*A*C))/(2*A));
double xv = (tx-x)/t, yv = (ty-y)/t, vx = d/t, vy = sqrt(V*V-vx*vx);
for(int i=min(tx,x);i<max(tx,x);i++){// -|-
t = (i-x+0.5)/xv;
double fy = h[x][y]+vy*t-G*t*t/2; int p = ceil(y-0.5+t*yv);
if(max(h[i][p],h[i+1][p])>=fy) return 0;
if(y-0.5+t*yv==p&&max(h[i][p+1],h[i+1][p+1])>=fy) return 0;
}
for(int i=min(ty,y);i<max(ty,y);i++){
t = (i-y+0.5)/yv;
double fy = h[x][y]+vy*t-G*t*t/2; int p = ceil(x-0.5+t*xv);
if(max(h[p][i],h[p][i+1])>=fy) return 0;
if(x-0.5+t*xv==p&&max(h[p+1][i],h[p+1][i+1])>=fy) return 0;
}
return 1;
}
int main()
{
scanf("%d%d%lf%lf%d%d",&n,&m,&W,&V,&sx,&sy);
for(int i=1;i<=m;i++) for(int j=1;j<=n;j++) scanf("%lf",&h[j][i]);
memset(dis,-1,sizeof dis);
static int q[maxn*maxn],l,r; q[l=r=1]=(sy-1)*n+sx,dis[sx][sy]=0;
while(l<=r){
int y=(q[l]-1)/n+1,x=(q[l]-1)%n+1; l++;
for(int i=1;i<=n;i++) for(int j=1;j<=m;j++) if(dis[i][j]==-1&&check(x,y,i,j))
q[++r]=(j-1)*n+i,dis[i][j]=dis[x][y]+1;
}
for(int i=1;i<=m;i++) for(int j=1;j<=n;j++)
if(~dis[j][i]) printf("%d%c",dis[j][i],j==n?10:32);
else printf("X%c",j==n?10:32);
}
2019 ICPC Asia-East Continent Final B. Black and White
n
∗
m
n*m
n∗m的网格，已黑白相间染色，左下角为白色，求从
(
0
,
0
)
(0,0)
(0,0)往右/上走到
(
n
,
m
)
(n,m)
(n,m)的轮廓线的左侧白格子减黑格子=k的方案数。
n
,
m
≤
1
0
5
n,m\le10^5
n,m≤105
考虑统计当前轮廓线上方白格子减去黑格子的和。
每次可以往右走或往上走，和不变/+1/-1。
为了解决+1/-1，我们一次走两步（需要
2
∣
(
n
+
m
)
2|(n+m)
2∣(n+m)），如果方向相同，那么对和以及分布都是没有影响的；如果方向不同，假设全都是先右再上，走了
p
p
p 步，那么最后的和等于
⌊
p
+
1
2
⌋
\lfloor \frac {p+1}2\rfloor
⌊2p+1⌋（奇数行白色多1，偶数行黑白相等），如果将其中的
q
q
q 步替换为先上再右，相当于白格子少了，和为
⌊
p
+
1
2
⌋
−
q
\lfloor \frac {p+1}2\rfloor-q
⌊2p+1⌋−q
然后枚举
p
p
p 组合数计算方案数即可。
如果 2 ∤ ( n + m ) 2\not|~~(n+m) 2∣ (n+m)，那么枚举最后一步往右还是往上，划归为 2 ∣ ( n + m ) 2|(n+m) 2∣(n+m) 的情况。
Code：
#include<bits/stdc++.h>
#define maxn 100005
using namespace std;
const int mod = 998244353;
int fac[maxn],inv[maxn];
int f(int n,int m,int k){
int ret=0,x=(n+m)/2;
for(int p=n&1,lim=min(n,m);p<=lim;p+=2){
int q=(p+1)/2-k;
if(0<=q&&q<=p) ret=(ret+1ll*fac[x]*inv[(n-p)>>1]%mod*inv[(m-p)>>1]%mod*inv[q]%mod*inv[p-q])%mod;
}
return ret;
}
int main()
{
fac[0]=fac[1]=inv[0]=inv[1]=1;
for(int i=2;i<maxn;i++) fac[i]=1ll*fac[i-1]*i%mod,inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
for(int i=2;i<maxn;i++) inv[i]=1ll*inv[i]*inv[i-1]%mod;
int T,n,m,k;
for(scanf("%d",&T);T--;){
scanf("%d%d%d",&n,&m,&k);
if((n+m)&1) printf("%d\n",(f(n,m-1,k)+f(n-1,m,k+(m&1)))%mod);
else printf("%d\n",f(n,m,k));
}
}
2018 ACM-ICPC World Finals G. Panda Preserve
简单多边形，求最小半径 r r r，使得以多边形顶点为圆心的半径为 r r r 的所有圆的并能够覆盖这个多边形。
对于多边形内部的某个点，离它最近的多边形顶点到它的距离记为 d p d_p dp，就是要求 d p d_p dp 的最大值。
对每个顶点用中垂线做半平面交求出最近范围，那么区域内离它最远的多边形内的点要么是半平面交的顶点，要么是多边形与半平面交的交点。
这个半平面交就是 voronoi \texttt{voronoi} voronoi 图，点数和边数都是 O ( n ) O(n) O(n) 级别的，所以可以枚举边与多边形的每条边相交，总复杂度 O ( n 2 log n + n 2 ) O(n^2\log n +n^2) O(n2logn+n2)
Code：
#include<bits/stdc++.h>
#define maxn 2015
using namespace std;
const double eps = 1e-10;
int n;
int sgn(double x){return x>eps?1:x<-eps?-1:0;}
struct Pt{
double x,y;
Pt(double x=0,double y=0):x(x),y(y){}
Pt operator + (const Pt &p)const{return Pt(x+p.x,y+p.y);}
Pt operator - (const Pt &p)const{return Pt(x-p.x,y-p.y);}
Pt operator * (const double &t)const{return Pt(x*t,y*t);}
double operator * (const Pt &p)const{return x*p.y-y*p.x;}
double len(){return sqrt(x*x+y*y);}
}a[maxn],p[maxn];
struct Line{
Pt p,v; double ang;
Line(Pt p=0,Pt v=0,bool flg=1):p(p),v(v){if(flg) ang=atan2(v.y,v.x);}
bool operator < (const Line &L)const{return sgn(ang-L.ang)?ang<L.ang:v*(L.p-p)<0;}
}L[maxn],q[maxn];
int cnt,l,r;
bool Onleft(Pt p,Line L){return L.v*(p-L.p)>0;}
Pt Inter(Line a,Line b){return a.p+a.v*((b.v*(a.p-b.p))/(a.v*b.v));}
void HalfPlane(){
sort(L+1,L+1+cnt),q[l=r=1]=L[1];
for(int i=2;i<=cnt;i++) if(sgn(L[i].ang-L[i-1].ang)){
while(l<r&&!Onleft(p[r-1],L[i])) r--;
while(l<r&&!Onleft(p[l],L[i])) l++;
q[++r]=L[i],p[r-1]=Inter(q[r-1],q[r]);
}
while(l<r-1&&!Onleft(p[r-1],q[l])) r--;
p[r]=Inter(q[r],q[l]);
}
bool InPolygon(Pt p){
bool crs=0;
for(int i=1;i<=n;i++){
if(!sgn((a[i]-p).len()+(a[i+1]-p).len()-(a[i+1]-a[i]).len())) return 1;
double u=a[i].y,v=a[i+1].y;
if(u==v||p.y<min(u,v)||p.y>=max(u,v)) continue;
crs^=(a[i].x+(a[i+1].x-a[i].x)*(p.y-u)/(v-u))>p.x;
}
return crs;
}
bool Intersect(Line a,Line b){
return ((b.p-a.p)*a.v) * ((b.p+b.v-a.p)*a.v) < 0
&& ((a.p-b.p)*b.v) * ((a.p+a.v-b.p)*b.v) < 0;
}
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;i++) scanf("%lf%lf",&a[i].x,&a[i].y); a[n+1]=a[1];
double ans=0;
for(int i=1;i<=n;i++){
L[cnt=1]=Line(Pt(-1e4,-1e4),Pt(1,0));
L[++cnt]=Line(Pt(1e4,-1e4),Pt(0,1));
L[++cnt]=Line(Pt(1e4,1e4),Pt(-1,0));
L[++cnt]=Line(Pt(-1e4,1e4),Pt(0,-1));
for(int j=1;j<=n;j++) if(i^j){
Pt o=a[i]+a[j],v=a[j]-a[i];
L[++cnt]=Line(Pt(o.x/2,o.y/2),Pt(-v.y,v.x));
}
HalfPlane(),p[r+1]=p[l];
for(int j=l;j<=r;j++){
if(InPolygon(p[j])) ans=max(ans,(p[j]-a[i]).len());
for(int k=1;k<=n;k++){
Line x=Line(p[j],p[j+1]-p[j],0),y=Line(a[k],a[k+1]-a[k],0);
if(Intersect(x,y)) ans=max(ans,(Inter(x,y)-a[i]).len());
}
}
}
printf("%.9f\n",ans);
}
353
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
