答案由$3$部分构成：
$1$.抛物线的极值。
$2$.询问区间的左端点在抛物线上的值。
$3$.询问区间的右端点在抛物线上的值。
对于$1$，就是某个矩形范围内最大值查询，使用KD-Tree可以在$O(sqrt{n})$的时间内完成询问。
对于$2$和$3$，对编号建线段树，每个节点维护该区间内最上方的轮廓线，那么每个节点的轮廓线可以由左右儿子的轮廓线通过归并得到。
归并时不断对两条抛物线求交点，在相邻关键点之间取较高的一段。
为了减小常数，可以在归并完成后将连续的所属同一条抛物线的轮廓线连起来。
查询时在线段树上找到$O(log n)$个节点，然后在那些节点的轮廓线上二分查找即可。
总时间复杂度$O(nlog n+m(log^2n+sqrt{n}))$。
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<vector>
#define PB push_back
using namespace std;
const int N=50010;
const double eps=1e-6;
namespace KD{
int root,cmp_d,ans,X1,X2,Y1,Y2;
struct node{int d[2],l,r,Max[2],Min[2],val,sum;}t[N];
inline bool cmp(const node&a,const node&b){return a.d[cmp_d]<b.d[cmp_d];}
inline void umax(int&a,int b){if(a<b)a=b;}
inline void umin(int&a,int b){if(a>b)a=b;}
inline void up(int x){
if(t[x].l){
umax(t[x].Max[0],t[t[x].l].Max[0]);
umin(t[x].Min[0],t[t[x].l].Min[0]);
umax(t[x].Max[1],t[t[x].l].Max[1]);
umin(t[x].Min[1],t[t[x].l].Min[1]);
umax(t[x].sum,t[t[x].l].sum);
}
if(t[x].r){
umax(t[x].Max[0],t[t[x].r].Max[0]);
umin(t[x].Min[0],t[t[x].r].Min[0]);
umax(t[x].Max[1],t[t[x].r].Max[1]);
umin(t[x].Min[1],t[t[x].r].Min[1]);
umax(t[x].sum,t[t[x].r].sum);
}
}
int build(int l,int r,int D){
int mid=(l+r)>>1;
cmp_d=D;
nth_element(t+l+1,t+mid+1,t+r+1,cmp);
t[mid].Max[0]=t[mid].Min[0]=t[mid].d[0];
t[mid].Max[1]=t[mid].Min[1]=t[mid].d[1];
t[mid].sum=t[mid].val;
if(l!=mid)t[mid].l=build(l,mid-1,!D);
if(r!=mid)t[mid].r=build(mid+1,r,!D);
return up(mid),mid;
}
void ask(int x){
if(t[x].Min[0]>X2||t[x].Max[0]<X1||t[x].Min[1]>Y2||t[x].Max[1]<Y1||t[x].sum<=ans)return;
if(t[x].Min[0]>=X1&&t[x].Max[0]<=X2&&t[x].Min[1]>=Y1&&t[x].Max[1]<=Y2){
ans=t[x].sum;
return;
}
if(t[x].d[0]>=X1&&t[x].d[0]<=X2&&t[x].d[1]>=Y1&&t[x].d[1]<=Y2)umax(ans,t[x].val);
if(t[x].l)ask(t[x].l);
if(t[x].r)ask(t[x].r);
}
inline int query(int A,int B,int C,int D){
ans=0;
X1=A,X2=B,Y1=C,Y2=D;
ask(root);
return ans;
}
inline void init(int i,int x,int y){t[i].d[0]=i,t[i].d[1]=x,t[i].val=y;}
}
inline void read(int&a){char c;while(!(((c=getchar())>='0')&&(c<='9')));a=c-'0';while(((c=getchar())>='0')&&(c<='9'))(a*=10)+=c-'0';}
int n,m,i,A,B,C,D,a[N][3];double ans;
struct E{
double a,b,c;long long X,Y;
E(){}
E(double p,double x,double y,int _X,int _Y){
a=-y/(p-x)/(p-x);
b=-2*a*x;
c=y+a*x*x;
X=1LL*_X*_X,Y=_Y;
}
inline double f(double x){return a*x*x+b*x+c;}
}e[N];
inline double pos(const E&a,const E&b,double l,double r){
double A=a.a-b.a,B=a.b-b.b,C=a.c-b.c;
if(a.X*b.Y==a.Y*b.X){
if(fabs(B)<eps)return r;
double x=-C/B;
if(x>l+eps&&x<r)return x;
return r;
}
double d=B*B-4*A*C;
if(d<-eps)return r;
d=sqrt(d);
double x=(-B-d)/(2*A),y=(-B+d)/(2*A),ret=r;
if(x>l+eps&&x<r)ret=x;
if(y>l+eps&&y<r)ret=min(ret,y);
return ret;
}
struct P{
double l,r;int p;
P(){}
P(double _l,double _r,int _p){l=_l,r=_r,p=_p;}
};
vector<P>v[131100],c;
vector<double>q;
inline void merge(vector<P>&f,const vector<P>&a,const vector<P>&b){
int na=a.size(),nb=b.size(),i=0,j=0,k,top=0;
c.clear();
q.clear();
q.PB(-N);
while(i<na&&j<nb)if(a[i].r<b[j].r||j==nb){
if(a[i].r>q[top]+eps)q.PB(a[i].r),top++;
i++;
}else{
if(b[j].r>q[top]+eps)q.PB(b[j].r),top++;
j++;
}
for(;i<na;i++)if(a[i].r>q[top]+eps)q.PB(a[i].r),top++;
for(;j<nb;j++)if(b[j].r>q[top]+eps)q.PB(b[j].r),top++;
for(i=j=k=0;k<top;k++){
double l=q[k],r=q[k+1];
while(i<na&&a[i].r+eps<r)i++;
while(j<nb&&b[j].r+eps<r)j++;
if(i==na||a[i].l>l+eps){
c.PB(P(l,r,b[j].p));
continue;
}
if(j==nb||b[j].l>l+eps){c.PB(P(l,r,a[i].p));continue;}
if(!a[i].p||!b[j].p){c.PB(P(l,r,a[i].p+b[j].p));continue;}
double m,mid;
while(l+eps<r){
m=pos(e[a[i].p],e[b[j].p],l,r);
mid=(l+m)/2;
if(e[a[i].p].f(mid)>e[b[j].p].f(mid)){
c.PB(P(l,m,a[i].p));
}else{
c.PB(P(l,m,b[j].p));
}
l=m;
}
}
for(i=0,k=c.size();i<k;i=j){
for(j=i;j<k&&c[i].p==c[j].p;j++);
f.PB(P(c[i].l,c[j-1].r,c[i].p));
}
}
void build(int x,int a,int b){
if(a==b){
v[x].PB(P(-N,::a[a][0],0));
v[x].PB(P(::a[a][0],2*::a[a][1]-::a[a][0],a));
v[x].PB(P(2*::a[a][1]-::a[a][0],N,0));
return;
}
int mid=(a+b)>>1;
build(x<<1,a,mid),build(x<<1|1,mid+1,b);
merge(v[x],v[x<<1],v[x<<1|1]);
}
inline void query(const vector<P>&v,int o){
int l=0,mid,r=v.size()-1;
while(l<=r){
mid=(l+r)>>1;
if(v[mid].r+eps<o)l=mid+1;
else if(v[mid].l-eps>o)r=mid-1;
else{
ans=max(ans,e[v[mid].p].f(o));
return;
}
}
}
void ask(int x,int a,int b){
if(A<=a&&b<=B){
query(v[x],C);
query(v[x],D);
return;
}
int mid=(a+b)>>1;
if(A<=mid)ask(x<<1,a,mid);
if(B>mid)ask(x<<1|1,mid+1,b);
}
int main(){
read(n);
for(i=1;i<=n;i++){
int p,x,y;
read(p),read(x),read(y);
KD::init(i,x,y);
a[i][0]=p,a[i][1]=x,a[i][2]=y;
e[i]=E(p,x,y,p-x,y);
}
KD::root=KD::build(1,n,0);
build(1,1,n);
read(m);
while(m--){
read(A),read(B),read(C),read(D);
ans=KD::query(A,B,C,D);
ask(1,1,n);
printf("%.6f
",ans);
}
return 0;
}