abstract class Eds
{
  public final static int EULER     = 0;
  public final static int MILSHTEIN = 1;
  public final static int RUNKUT4   = 2;
  public final static int HEUN      = 3;

  public int method = EULER;

  float oneStep(float x, float h, float t, float dw)
  {
    switch (method)
      {
        case EULER:     return euler(x,h,t,dw);
        case MILSHTEIN: return milshtein(x,h,t,dw);
        case RUNKUT4:   return runkut4(x,h,t,dw);
        case HEUN:      return heun(x,h,t,dw);
        default:        return euler(x,h,t,dw);
      }
  }

  private float euler(float x, float h, float dw, float t)
  {
    return x + b(x, t)*h + sig(x, t)*dw;
  }

  private float milshtein(float x, float h, float dw, float t)
  {
    float     g, dgg;
    g = sig(x, t);
    dgg = (float) 0.5 * dsigdx(x, t) * g;
    return x + (b(x, t) - dgg) * h + g * dw + dgg * dw * dw;
  }

  private float runkut4(float x, float h, float t, float dw)
  {
    float     f0, f1, f2, f3, g0, g1, g2, g3, x1, x2, x3, t1, t2;
    t1 = t + ((float) 0.5) * h;
    t2 = t + h;
    f0 = a(x, t);
    g0 = sig(x, t);
    x1 = x + ((float) 0.5) * (f0 * h + g0 * dw);
    f1 = a(x1, t1);
    g1 = sig(x1, t1);
    x2 = x + ((float) 0.5) * (f1 * h + g1 * dw);
    f2 = a(x2, t1);
    g2 = sig(x2, t1);
    x3 = x + f2 * h + g2 * dw;
    f3 = a(x3, t2);
    g3 = sig(x3, t2);
    return x + ((f0 + 2 * (f1 + f2) + f3) * h + 
		(g0 + 2 * (g1 + g2) + g3) * dw) / 6;
  }
  
  private float heun(float x, float h, float t, float dw)
  {
    float     xx;
    xx = x + a(x, t) * h + sig(x, t) * dw;
    return x + ((float) 0.5) * (a(x, t) + a(xx, t + h)) * h
      + ((float) 0.5) * (sig(x, t) + sig(xx, t + h)) * dw;
  }


  abstract float b(float x, float t); //Ito drift 
  abstract float a(float x, float t); //Stratonovitch drift
  abstract float sig(float x, float t); //Diffusion
  abstract float dsigdx(float x, float t); //Derivative of diff.

}