def METHOD_NAME(f,gx,gy,gz,function_type, eps = 1.e-6):
  """
  Exercise linear, quadratic and tricubic interpolation. Thorough.
  Behaviors at corners and boundares are not exercised.
  """
  for n in [25,50,100]:
    O = set_map(f=f, n=n)
    n_trials = 100
    tl,  tq,  tt  = 0,0,0
    gxl, gxq, gxt = 0,0,0
    gyl, gyq, gyt = 0,0,0
    gzl, gzq, gzt = 0,0,0
    for trial in range(n_trials):
      sx = random.randrange(1,30)/10.*random.choice([-1,0,1])
      sy = random.randrange(1,30)/10.*random.choice([-1,0,1])
      sz = random.randrange(1,30)/10.*random.choice([-1,0,1])
      x,y,z = [5+sx, 5+sy, 5+sz]
      e = get_exact(f=f, gx=gx,gy=gy,gz=gz, x=x,y=y,z=z)
      if(abs(e.t)<eps or abs(e.gx)<eps or abs(e.gy)<eps or abs(e.gz)<eps):
        continue
      site_frac = [x/O.box,y/O.box,z/O.box]
      i = get_interpolated(site_frac=site_frac, O=O)
      tl += re(i.l[0], e.t)
      tq += re(i.q[0], e.t)
      tt += re(i.t[0], e.t)
      #
      gxl += re(i.l[1], e.gx)
      gyl += re(i.l[2], e.gy)
      gzl += re(i.l[3], e.gz)
      #
      gxq += re(i.q[1], e.gx)
      gyq += re(i.q[2], e.gy)
      gzq += re(i.q[3], e.gz)
      #
      gxt += re(i.t[1], e.gx)
      gyt += re(i.t[2], e.gy)
      gzt += re(i.t[3], e.gz)
    r = []
    for it in [tl,gxl,gyl,gzl, tq,gxq,gyq,gzq, tt,gxt,gyt,gzt]:
      r.append(it/n_trials)
    tl,gxl,gyl,gzl, tq,gxq,gyq,gzq, tt,gxt,gyt,gzt = r
    #
    s = "(%10.6f %10.6f %10.6f %10.6f)"
    fmt = " ".join(["%5.3f", s, s, s])
    print(fmt%(O.step[0], tl,gxl,gyl,gzl, tq,gxq,gyq,gzq, tt,gxt,gyt,gzt))
    #
    cntr=0
    if(function_type=="L"):
      cntr+=1
      for it in [tl,gxl,gyl,gzl, tq,gxq,gyq,gzq, tt,gxt,gyt,gzt]:
        assert approx_equal(it, 0.0, 1.e-4)
    if(function_type=="Q"):
      cntr+=1
      for it in [tq,gxq,gyq,gzq, tt,gxt,gyt,gzt]:
        assert approx_equal(it, 0.0, 1.e-4)
    if(function_type=="C"):
      cntr+=1
      for it in [tt,gxt,gyt,gzt]:
        assert approx_equal(it, 0.0, 1.e-4)
    assert cntr==1