rem: bilinear_interp(old_lp, new_lp);
add: bi(old_lp, new_lp);
context: int multi(laplace& lp){  int min_size = 5;  int iter = 0;  double res;		// residual  int I = lp.I();  int J = lp.J();  int i, j;  int mult;  int depth;  pair<double, double> x(-4,4);  pair<double, double> y(-4,4);    if(I != J)    return 0;  if(I < min_size)    return 0;  // find depth of recursion  mult = I;  depth = 0;  while(mult != min_size && mult != 0){    mult >>= 1;    depth++;  }  if(mult == 0)    return 0;  residual_norm norm(lp, EPS);  int size = mult;  boundary top, right, bottom, left;  top = boundary(boundary::DIRICHLET, size, 0.0);  right = boundary(top);  bottom= boundary(top);  left= boundary(boundary::DIRICHLET, size, 1.0);  laplace base(size, size, x, y, top, right, bottom, left);  gauss_seidel gs(20000);  iter = gs.solve(base, norm);  cout << "# iter " << min_size << ": " << iter << endl;  laplace old_lp(base), new_lp(base);  size *= 2;  int iter_i =0;  while(size < I){    top = boundary(boundary::DIRICHLET, size, 0.0);    right = boundary(top);    bottom= boundary(top);    left = boundary(boundary::DIRICHLET, size, 1.0);    old_lp = new_lp;    new_lp = laplace(size, size, x, y, top, right, bottom, left);    bilinear_interp(old_lp, new_lp);    iter_i =gs.solve(new_lp, norm);     iter += iter_i;    cout << "# iter " << size << ": " << iter_i << endl;    size *= 2;  }  if(min_size < I){    bilinear_interp(new_lp, lp);    iter_i= gs.solve(lp, norm);    cout << "# iter " << size << ": " << iter_i << endl;    iter += iter_i;  }  else    lp = base;  return iter;}