I tried to implement my proposal in a C-code. That is a mixture of analytic and numeric integration. It does $10^6$ rectangles with half-percent relative precision in about 16 seconds, which is a bit better than the corresponding Iosif's 30 minutes. You can play with parameters to trade speed for precision and vice versa too. The code should be self-explanatory but feel free to ask questions if something is unclear.

*Edit:* This is the best and the fastest version. $n$ is gone now and the guaranteed relative precision is $1/N^2$ (the constant $1$ is correct, so if you want $10^{-3}$ accuracy (to compare with Mathematica time), just set $N=34$ and get $10^6$ pairs in under 10 seconds. The time is essentially proportional to $N$ For $10^{-5}$ accuracy $N=340$ and 83 seconds suffice. I'll explain the algorithm a bit later; now it makes sense :-)
  

    #include <stdlib.h>
    #include <stdio.h>
    #include <math.h>
    #include <time.h>
    
    const double pi=3.141592653589, ppi=pi/57.6, pi2=pi/2, dl=sqrt(0.6)/2;
    
    double gghh(double a, double b, double c, double d, double x, double t)
    {
    double y=fabs(t), g=y*a, h=2*d;
    if(y<=a-b) g=(a*a+y*y+b*b/3)*0.5; 
    else if(y<a+b) {double z=a+b-y; g+=z*z*z/(12*b);}
    
    y=fabs(t-x); 
    if(y>c-d) h-=(y-c+d);
    return g*h;
    }
    
    double F(double a, double b, double c, double d, double x, double aa, double bb)
    {
    double t2=(aa+bb)*0.5, bbaa=bb-aa, dt=dl*bbaa;
    return bbaa*(gghh(a,b,c,d,x,t2-dt)+gghh(a,b,c,d,x,t2+dt)+1.6*gghh(a,b,c,d,x,t2))/(a*c*d);
    }
    
    double D(double a1,double b1, double c1, double d1, double a2,double b2, double c2, double d2, int N)
    {
    double s=0.0;
    double X1=b1-a1, Y1=d1-c1, X2=b2-a2, Y2=d2-c2, S1=(a2+b2-a1-b1), S2=(c2+d2-c1-d1);
    
    double t0=pi2/N, cs=cos(t0), ss=sin(t0), dcs=2*cs*cs-1, dss=2*cs*ss; 
    double SS=fabs(S1)+fabs(S2)+fabs(X1)+fabs(X2)+fabs(Y1)+fabs(Y2);
    SS*=0.00000001;
    for(int k=0; k<N;++k)
    { 
    double csnew=cs*dcs-ss*dss;
    ss=ss*dcs+cs*dss; cs=csnew;
    double U[4]={fabs(X1*cs)+SS, fabs(Y1*ss)+SS, fabs(X2*cs)+SS, fabs(Y2*ss)+SS};
    double x=-fabs(S1*cs+S2*ss);
    
    
    for(int kk=0;kk<3;++kk)
    {
    int kkk=3-kk;
    for(int j=0;j<kkk;++j)
    if(U[j]>U[j+1]) {double u=U[j]; U[j]=U[j+1]; U[j+1]=u;}
    }
    
    
    double U0=U[0], U1=U[1], U2=U[2], U3=U[3];
    
    double V[4]={-U3-U2,-U3+U2,U3-U2,U3+U2}, 
    VV[4]={x-U1-U0,x-U1+U0,x+U1-U0,x+U1+U0};
    
    double W[8]; 
    int i=0, ii=0, kstart=-1, kfinish=-1;
    while(ii<4)
    {
    ++kfinish; 
    if(V[i]<VV[ii]) {W[kfinish]=V[i]; ++i;}
    else {W[kfinish]=VV[ii]; if(ii==0) kstart=kfinish; ++ii;} 
    }
    
    for(int kk=kstart;kk<kfinish;++kk)
    s+=F(U3,U2,U1,U0,x,W[kk],W[kk+1]);
    }
    return ppi*s/N;
    }
    
    
    
    double unitrand()
    {
    return (rand()+0.0)/RAND_MAX;
    }
    
    
    int main()
    {
    time_t now=time(0);
    srand(now);	
    
    int N=1000;
    
    double m=100,M=0;
    
    for(int k=0; k<1000000;++k)
    {
    if(k%10000==0) {printf("%d %.12f %.12f\n",k/10000,m,M);}
    double 
    a1=unitrand(),b1=a1+unitrand(),
    a2=unitrand(),b2=a2+unitrand(),
    c1=unitrand(),d1=c1+unitrand(),
    c2=unitrand(),d2=c2+unitrand();
    
    double r=D(a1,b1,c1,d1,a2,b2,c2,d2,N);
    
    if(k%1000==0)
    {
    r/=D(a1,b1,c1,d1,a2,b2,c2,d2,600);
    if(r<m) m=r;
    if(r>M) M=r;
    }
    }
    printf("\n%.12f",D(1,2,3,5,4,6,7,8,4000));
    printf("\n%.12f",D(1,2,3,5,4,6,7,8,N)/D(1,2,3,5,4,6,7,8,2000)-1);
    printf("\n%.12f",D(0,2,0,2,0,2,0,2,N)/D(0,2,0,2,0,2,0,2,2000)-1);
    printf("\n%.12f",D(0,3,0,0.0001,0,3,0,0.0001,N)-1);
    printf("\n%.12f",D(0,2,0,0,0,0,0,2,N)/D(0,2,0,0,0,0,0,2,2000)-1);
    printf("\n%.12f",D(0,0,0,0,3,3,4,4,N)/5-1);
    return 0;
    }