#include <iostream>
#include <iomanip>
#include "nr.h"
using namespace std;

// Driver for routine bcuint

int main(void)
{
        const DP xx_d[4]={0.0,2.0,2.0,0.0};
        const DP yy_d[4]={0.0,0.0,2.0,2.0};
        int i;
        DP ansy,ansy1,ansy2,ey,ey1,ey2;
        DP x1,x1l,x1u,x1x2,x2,x2l,x2u,xxyy;
        Vec_DP y(4),y1(4),y12(4),y2(4);
        Vec_DP xx(xx_d,4),yy(yy_d,4);

        x1l=xx[0];
        x1u=xx[1];
        x2l=yy[0];
        x2u=yy[3];
        for (i=0;i<4;i++) {
          xxyy=xx[i]*yy[i];
          y[i]=xxyy*xxyy;
          y1[i]=2.0*yy[i]*xxyy;
          y2[i]=2.0*xx[i]*xxyy;
          y12[i]=4.0*xxyy;
        }
        cout << endl << setw(6) << "x1" << setw(9) << "x2";
        cout << setw(8) << "y" << setw(12) << "expect";
        cout << setw(7) << "y1" << setw(11) << "expect";
        cout << setw(7) << "y2" << setw(11) << "expect" << endl << endl;
        cout << fixed << setprecision(4);
        for (i=0;i<10;i++) {
          x2=(x1=0.2*(i+1));
          NR::bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,ansy1,ansy2);
          x1x2=x1*x2;
          ey=x1x2*x1x2;
          ey1=2.0*x2*x1x2;
          ey2=2.0*x1*x1x2;
          cout << setw(8) << x1 << setw(9) << x2 << setw(9) << ansy;
          cout << setw(9) << ey << setw(9) << ansy1 << setw(9) << ey1;
          cout << setw(9) << ansy2 << setw(9) << ey2 << endl;
        }
        return 0;
}