{
  checkoest(L, a)=

  m = floor(log(L) / log(sqrtq));
  u = (sqrtq^(m+1) - L) / (L*sqrtq - sqrtq^m);
  t0 = solve(x=Pi/(m+1), Pi/m, cos((m+1)*x/2) + u * cos((m-1)*x/2));
  a = cos(t0)*sqrtq;
  ((L-1)*a + q - L) / (q + 1 - 2*a) <= g - eps
}

{
  oesterle(q,g)=

  eps = 10^(-10);
  sqrtq = sqrt(q);
  b = q + g * floor(2*sqrtq);
  a = -1;
  while (b - a > 1,
    L = floor((a+b)/2);
    if (checkoest(L),
      a = L
    ,
      b = L
    );
  );
  if (checkoest(b), b+1, a+1)
}