本篇题解参考 @djwj223，但是本人太弱了，对着代码看了 INF 年才看懂。因此写一篇具体讲解实现方面的题解。在此先表示感谢。
思路
考虑最终的答案圆会是怎样的。第一种是半径达到了 \(R\) 的上界，不能继续扩充。显然这种情况可以把圆移动一下使某个点在圆上，以此进行计算。第二种是楼上所说的与两个圆相切并过一个点。此时可以先钦定与一个圆相切且过某个点。然后进行二分。
第一种情况
考虑枚举每一个点 \(A\)，我们钦定答案圆的半径为 \(R\) 且过这个点。相当于可能成为答案的圆（或者说半径）绕着这个点旋转（见下图）。
然后考虑每个其他的点和圆，我们可以求出答案圆旋转到哪个位置区间时可以与它们有交。可以利用余弦定理求出这个区间。
如图，利用余弦定理求出 \(\angle BAC\) 和 \(\angle BAD\)，再使用 atan2 计算 \(AB\) 的旋转角（\(x\) 正方向为 \(0\)，转回去就是 \(2\pi\)），然后即可求出 \(AC\) 和 \(AD\) 的旋转角，调整到 \([0,2\pi)\)，带上贡献（圆为 \(-\infty\)，点为 \(1\)）排序后前缀和即可。
注意像图中的 \(AC\) 和 \(AD\) 横跨 \(x\) 正半轴，调整后的旋转角加比减小，需要把计算的变量初始值加上贡献。这样可以保证在 \(AC\) 处把贡献减掉，在 \(AD\) 处加回来。
这部分的代码（pw 是平方，\(d\) 是预处理的圆心距离）：
struct cir{
int x,y;
ld r;
}p[N];
struct s{
ld pos;
int del;
bool operator<(s b){
return pos<b.pos;
}
}q[N<<1];
inline ld fun(ld a){ //调整
while(a>=2*pi){
a-=2*pi;
}
while(a<-eps){
a+=2*pi;
}
return max(a,(ld)0.0);
}
inline void solve_r(int x){
ld A,B,C,l,r,deg,f;
int res=1,del,cnt=0;
for(int i=1;i<=n+m;i++){ //前 n 个是圆，后面是点。
if(x==i){
continue;
}
if(d[i][x]>pw(maxr*2+p[i].r)){ //根本不交
continue;
}
A=sqrt(d[i][x]);
B=maxr;
C=maxr+p[i].r;
deg=acos((B*B+d[i][x]-C*C)/(2*A*B));
f=atan2(p[i].y-p[x].y,p[i].x-p[x].x);
l=fun(f-deg);
r=fun(f+deg);
if(fabs(l-r)<eps){ //减小精度影响
r=l+eps;
}
del=(i<=n?-2100:1); //贡献
if(l>r){ //横跨
res+=del;
}
q[++cnt]={l,del};
q[++cnt]={r,-del};
}
ans=max(ans,res);
sort(q+1,q+cnt+1);
for(int i=1;i<=cnt;i++){
res+=q[i].del;
ans=max(ans,res);
}
}
第二种情况
第二种考虑是先钦定与一个圆相切且过某个点。此时可能的答案圆如下图（mid 是二分时的半径）：
如何求出这个点 \(B\) 的坐标呢？考虑以 \(A\) 为圆心，\(mid\) 为半径的圆，和以 \(C\) 为圆心 \(mid+r\) 为半径的圆。两者的（两个）交点即为 \(B_1\) 和 \(B_2\)。求出 \(B\) 后，判断当前答案圆与第二个圆的位置关系即可继续二分。二分结束后可以直接计算并更新答案。
这个求交点的过程与第一种情况雷同，唯一不同是这里需要坐标而不是角度。但是使用上述的余弦定理方法速度非常慢，因此这里可以直接计算三角函数值并使用和角公式。
这部分的代码：
inline void get_intersect(cir a,cir b,int qq){ //计算交点
int dis=pw(a.x-b.x)+pw(a.y-b.y);
ld cosar=(dis+pw(a.r)-pw(b.r))/(2*sqrt(dis)),sinar=sqrt(max(pw(a.r)-pw(cosar),(ld)0.0)),cosb=(b.x-a.x)/sqrt(dis),sinb=(b.y-a.y)/sqrt(dis);
if(qq){
sinar=-sinar; //某一个交点。另外一个角度加上负值。
}
xx=a.x+cosar*cosb-sinar*sinb;
yy=a.y+sinar*cosb+sinb*cosar;
}
inline int get(ld px,ld py,ld r){ //计算答案
int ans=0;
for(int i=1;i<=n;i++){
if(pw(px-p[i].x)+pw(py-p[i].y)-pw(r+p[i].r)<-eps){
return 0;
}
}
for(int i=1;i<=m;i++){
ans+=pw(px-p[i+n].x)+pw(py-p[i+n].y)-pw(r)<eps;
}
return ans;
}
inline void solve_smaller(int x,int y,int z){
ld l,r,mid;
if(d[x][y]>pw(maxr*2+p[x].r+p[y].r)){ //至少要刚好能卡上
return;
}
for(int i=0;i<=1;i++){
l=(sqrt(d[x][z])-p[x].r)/2; //下界是仅有一个交点的时候。
r=maxr;
if(l>r+eps){
return;
}
for(int kk=1;kk<=100;kk++){
mid=(l+r)/2;
get_intersect({p[x].x,p[x].y,p[x].r+mid},{p[z].x,p[z].y,mid},i);
if(pw(xx-p[y].x)+pw(yy-p[y].y)>=pw(mid+p[y].r)){
l=mid; //没有交
}
else{
r=mid;
}
}
get_intersect({p[x].x,p[x].y,p[x].r+l},{p[z].x,p[z].y,l},i);
ans=max(ans,get(xx,yy,l));
}
}
这道题不怎么卡精度。这是好的。但是时间要卡，需要开 O2。这好吗？这不好。
Code
#include<bits/stdc++.h>
#define ld long double
using namespace std;
const ld eps=1e-7,pi=acos(-1);
const int N=2020;
int n,m,maxr,d[N][N],ans;
ld xx,yy;
struct cir{
int x,y;
ld r;
}p[N];
struct s{
ld pos;
int del;
bool operator<(s b){
return pos<b.pos;
}
}q[N<<1];
inline ld fun(ld a){
while(a>=2*pi){
a-=2*pi;
}
while(a<-eps){
a+=2*pi;
}
return max(a,(ld)0.0);
}
inline int pw(int x){
return x*x;
}
inline ld pw(ld x){
return x*x;
}
inline void solve_r(int x){
ld A,B,C,l,r,deg,f;
int res=1,del,cnt=0;
for(int i=1;i<=n+m;i++){
if(x==i){
continue;
}
if(d[i][x]>pw(maxr*2+p[i].r)){
continue;
}
A=sqrt(d[i][x]);
B=maxr;
C=maxr+p[i].r;
deg=acos((B*B+d[i][x]-C*C)/(2*A*B));
f=atan2(p[i].y-p[x].y,p[i].x-p[x].x);
l=fun(f-deg);
r=fun(f+deg);
if(fabs(l-r)<eps){
r=l+eps;
}
del=(i<=n?-2100:1);
if(l>r){
res+=del;
}
q[++cnt]={l,del};
q[++cnt]={r,-del};
}
ans=max(ans,res);
sort(q+1,q+cnt+1);
for(int i=1;i<=cnt;i++){
res+=q[i].del;
ans=max(ans,res);
}
}
inline void get_intersect(cir a,cir b,int qq){
int dis=pw(a.x-b.x)+pw(a.y-b.y);
ld cosar=(dis+pw(a.r)-pw(b.r))/(2*sqrt(dis)),sinar=sqrt(max(pw(a.r)-pw(cosar),(ld)0.0)),cosb=(b.x-a.x)/sqrt(dis),sinb=(b.y-a.y)/sqrt(dis);
if(qq){
sinar=-sinar;
}
xx=a.x+cosar*cosb-sinar*sinb;
yy=a.y+sinar*cosb+sinb*cosar;
}
inline int get(ld px,ld py,ld r){
int ans=0;
for(int i=1;i<=n;i++){
if(pw(px-p[i].x)+pw(py-p[i].y)-pw(r+p[i].r)<-eps){
return 0;
}
}
for(int i=1;i<=m;i++){
ans+=pw(px-p[i+n].x)+pw(py-p[i+n].y)-pw(r)<eps;
}
return ans;
}
inline void solve_smaller(int x,int y,int z){
ld l,r,mid;
if(d[x][y]>pw(maxr*2+p[x].r+p[y].r)){
return;
}
for(int i=0;i<=1;i++){
l=(sqrt(d[x][z])-p[x].r)/2;
r=maxr;
if(l>r+eps){
return;
}
for(int kk=1;kk<=100;kk++){
mid=(l+r)/2;
get_intersect({p[x].x,p[x].y,p[x].r+mid},{p[z].x,p[z].y,mid},i);
if(pw(xx-p[y].x)+pw(yy-p[y].y)>=pw(mid+p[y].r)){
l=mid;
}
else{
r=mid;
}
}
get_intersect({p[x].x,p[x].y,p[x].r+l},{p[z].x,p[z].y,l},i);
ans=max(ans,get(xx,yy,l));
}
}
int main(){
scanf("%d%d%d",&n,&m,&maxr);
for(int i=1;i<=n;i++){
scanf("%d%d%Lf",&p[i].x,&p[i].y,&p[i].r);
}
for(int i=1;i<=m;i++){
scanf("%d%d",&p[i+n].x,&p[i+n].y);
}
for(int i=1;i<=n+m;i++){
for(int j=i+1;j<=n+m;j++){
d[i][j]=d[j][i]=pw(p[i].x-p[j].x)+pw(p[i].y-p[j].y);
}
}
for(int i=n+1;i<=n+m;i++){
solve_r(i);
}
for(int i=1;i<=n;i++){
for(int j=i+1;j<=n;j++){
for(int k=1;k<=m;k++){
solve_smaller(i,j,k+n);
}
}
}
printf("%d",ans);
}