klee.cpp

 #include <bits/stdc++.h>
#include <cassert>
using namespace std;



struct Box {
    vector<double> lo, hi;
};

struct Bounds {
    vector<double> lo, hi;
};

static int D;

double volume(const Bounds& g) {
    double v = 1.0;
    for (int i = 0; i < D; i++) {
        double w = g.hi[i] - g.lo[i];
        if (w <= 0.0) return 0.0;
        v *= w;
    }
    return v;
}



void simplify(vector<Box>& B, Bounds& g) {

for(int a=0;a<D;a++){

  

  vector<Box> non_slab;
  vector<pair<double,int>> events;

  for(int b=0;b<B.size();b++){

    bool slab=true;

    for(int d=0;d<D;d++){

      if(d==a)continue;

      if(B[b].lo[d] > g.lo[d] || B[b].hi[d] < g.hi[d]){

       slab=false;

       break;
      }

    }

    if(slab){

      double lft = max(B[b].lo[a] , g.lo[a]);
      double rgt = min(B[b].hi[a] , g.hi[a]);

      if(lft >= rgt)continue;

      events.emplace_back(lft,1);
      events.emplace_back(rgt,-1);

    }
    else non_slab.push_back(B[b]);


  }

  sort(events.begin(),events.end(),[](auto& a,auto& b){
      return a.first != b.first ? a.first < b.first : a.second > b.second;
      });

  vector<pair<double,double>> rem;

  int active=0;
  double beg=0;

  for(auto& e : events){
    if(active == 0 && e.second == 1)beg=e.first;
    active+=e.second;

    if(active==0)rem.emplace_back(beg,e.first);
  }

  vector<double> pref(rem.size()+1,0.0);

  for(int i=0;i<rem.size();i++){
    pref[i+1]=pref[i]+(rem[i].second-rem[i].first);
  }

  auto shift=[&](double x){
    if(rem.empty())return x;

    int low=0,high=rem.size()-1,idx=-1;

    while(low<=high){
      int mid=low+(high-low)/2;

      if(rem[mid].first<x){
        idx=mid;
        low=mid+1;
      }
      else high=mid-1;
    }

    if(idx==-1)return x;

    double netshft=pref[idx];

    if(x>rem[idx].first)netshft+=min(x,rem[idx].second)-rem[idx].first;

    return x-netshft;
  };

  for(auto& box: non_slab){
    box.lo[a] = shift(box.lo[a]);
    box.hi[a] = shift(box.hi[a]);
  }
  

  g.lo[a]=shift(g.lo[a]);
  g.hi[a]=shift(g.hi[a]);

  B=move(non_slab);
  }

}

double cut(vector<Box>& B, Bounds& g, int axis) {

  vector<pair<double,double>> coords;

  double netwt=0.0;

  for(int i=0;i<B.size();i++){
    double currw=0.0;
    
     double lft=max(B[i].lo[axis],g.lo[axis]);
     double rgt=min(B[i].hi[axis],g.hi[axis]);

   
     
     if(lft>=rgt)continue;

    for(int j=0;j<D;j++){
      if(j!=axis){

        if(B[i].lo[j]>g.lo[j] || B[i].hi[j]<g.hi[j]){
     
        currw=pow(2.0,((double)((axis+1)+(j+1))/D));

        if(lft>g.lo[axis] && lft<g.hi[axis]){
        
          coords.emplace_back(lft,currw);
          netwt+=currw;
        }

        if(rgt>g.lo[axis] && rgt<g.hi[axis]){
        
          coords.emplace_back(rgt,currw);
          netwt+=currw;
        }
        }
      }

    }

  }

  if(coords.empty()){
    return (g.lo[axis]+g.hi[axis])/2.0;
  }

  sort(coords.begin(),coords.end());

  double pref=0.0;

  for(auto& [crd,cwt]:coords){
    pref+=cwt;
    if(pref>=netwt/2.0)return crd;
  }

  return coords.back().first;
}

double measure(vector<Box> B, Bounds g, int axis,int depth) {

  assert(depth<500);

  //now we will write the recursion logic using the helper functions 
  simplify(B,g);

    if(B.empty()){
    return volume(g);
  }

  if(B.size()==1){
    //This is the base case if there is one single box we simply get the volume by clipping against the boundaries of the domain
    
    vector<double> hi_clip,lo_clip;

    hi_clip.resize(D),lo_clip.resize(D);

    for(int i=0;i<D;i++){

      lo_clip[i] = max(B[0].lo[i] , g.lo[i]);

      hi_clip[i] = min(B[0].hi[i] , g.hi[i]);
    }

    double clipped_box=1.0;

    for(int i=0;i<D;i++){

      double curr=hi_clip[i]-lo_clip[i];

      if(curr<=0){
        //If the box is out of bounds of domain in any dimeonsion, it is completely outside

        clipped_box=0.0;
        break;
      }

      clipped_box*=curr;
    }

    return volume(g)-clipped_box;
  }
    
  double midpt=cut(B,g,axis);

  vector<Box> lftB,rgtB;

  Bounds lftG=g,rgtG=g;

  lftG.hi[axis]=midpt;
  rgtG.lo[axis]=midpt;

  for(auto& bx:B){
    if(bx.lo[axis]<midpt){
      Box tmp=bx;
      tmp.hi[axis]=min(tmp.hi[axis],midpt);
      lftB.push_back(tmp);
    }
    if(bx.hi[axis]>midpt){
      Box tmp=bx;
      tmp.lo[axis]=max(tmp.lo[axis],midpt);
      rgtB.push_back(tmp);
    }
  }

  

  double lft_void = measure(lftB,lftG,(axis+1)%D,depth+1);
  double rgt_void = measure(rgtB,rgtG,(axis+1)%D,depth+1);

  return lft_void+rgt_void;

}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    cout << setprecision(10);

    int n;

    //Taking the dimensions and number of boxes as input 
    cin>>D>>n;

    vector<Box> boxes(n);
    //Declaring a AoS to keep track of the boxes 

    for (int b = 0; b < n; b++) {
        boxes[b].lo.resize(D);
        boxes[b].hi.resize(D);
        for (int i = 0; i < D; i++) cin >> boxes[b].lo[i] >> boxes[b].hi[i];
        //Per box we need to take the bounds of each box as the input 
    }

    Bounds domain;
    //This signifies the domain of our space 

    domain.lo.resize(D);
    domain.hi.resize(D);

    for(int i=0;i<D;i++)cin>>domain.lo[i]>>domain.hi[i];
   
    double empty = measure(boxes, domain, 0,0);
    double uni = volume(domain) - empty;
    cout << uni << "\n";
    return 0;
}