#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <queue>
#include <map>
#include <set>
#include <utility>
#include <algorithm>
#include <complex>
#include <stdint.h>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <climits>
#include <cassert>

//#define NU_DEBUG

//#define STRICT_MPZ

#ifdef STRICT_MPZ
#include <gmpxx.h>
#endif

using namespace std;

#define REP(i,n) for(int i = 0; i < (int)(n); i++)
#define FOR(i,c) for(__typeof((c).begin()) i = (c).begin(); i != (c).end(); ++i)
#define ALLOF(c) (c).begin(), (c).end()

typedef double decimal;
const decimal EPS = 1e-8;

#define MAXN 1000  // ???

//typedef complex<int> P;
struct P {
  int x, y;
};

template <typename T>
struct F {
  int64_t i;
  T n;
  T d;
};
typedef F<int32_t> F1;
typedef F<int64_t> F2;

#ifdef STRICT_MPZ
struct Coeff {
  mpz_t n;
  mpz_t d;
};

// (a+sqrt(b))/c
struct ABC {
  mpz_t a, b, c;
};
typedef ABC Solution;
typedef ABC Area;
#else
typedef decimal Coeff;
typedef decimal Solution;
typedef decimal Area;
#endif

struct Quad {
  Coeff S, dS, ddS;
};

inline bool operator<(const F1& a, const F1& b){
  return a.i != b.i ? a.i < b.i : a.n * b.d < b.n * a.d;
}

inline F1 operator+(const F1& a, int x){
  F1 f;
  f.i = a.i + x;
  f.n = a.n;
  f.d = a.d;
  return f;
}

inline F1 operator*(const F1& a, int x){
  int64_t n = (int64_t)a.n * x;
  F1 f;
  f.i = a.i * x + n / a.d;
  f.n = n % a.d;
  if(f.n < 0){
    f.i--;
    f.n += a.d;
  }
  f.d = a.d;
  return f;
}

inline F2 operator+(const F2& a, const int64_t& x){
  F2 f;
  f.i = a.i + x;
  f.n = a.n;
  f.d = a.d;
  return f;
}

inline F2 operator+(const F1& a, const F1& b){
  F2 f;
  f.i = a.i + b.i;
  f.n = (int64_t)a.n * b.d + (int64_t)b.n * a.d;
  f.d = (int64_t)a.d * b.d;
  return f;
}

inline F2 operator-(const F1& a, const F1& b){
  F2 f;
  f.i = a.i - b.i;
  f.n = (int64_t)a.n * b.d - (int64_t)b.n * a.d;
  f.d = (int64_t)a.d * b.d;
  return f;
}

inline F2 operator-(const F2& a){
  F2 f;
  f.i = -a.i;
  f.n = -a.n;
  f.d = a.d;
  return f;
}

inline ostream& operator<<(ostream& os, const F1& a){
  return os << a.i << " + " << a.n << " / " << a.d;
}

inline ostream& operator<<(ostream& os, const F2& a){
  return os << a.i << " + " << a.n << " / " << a.d;
}

#ifdef STRICT_MPZ
#error mpz coeff reconstruct is not implemented!
#else
inline Coeff reconstruct(const vector<F2>& v){
  Coeff c0 = 0.0;
  Coeff c1 = 0.0;
  REP(i, v.size()){
    c0 += v[i].i;
    c1 += (Coeff)v[i].n / v[i].d;
  }
  return c0 + c1;
}
#endif

#ifdef STRICT_MPZ
#error mpz area reconstruct is not implemented!
#else
inline Area reconstruct(int64_t x){
  Area a = x;
  return a;
}
#endif

#ifdef STRICT_MPZ
#error mpz to_decimal is not implemented!
#else
inline decimal to_decimal(const Area& a){
  decimal ret = a;
  return ret;
}
#endif

int gcd(int x, int y){
  return y ? gcd(y, x%y) : x;
}

int N;

int scale;
P ps[MAXN+4];
F1 slope[MAXN+4];
int64_t area[MAXN+4][MAXN+4];
int64_t AREA;
vector<int> fp;
vector<int> same[MAXN+4];

// global retval
// scale
// ps: normalized points
// fp: feature points (id)
// same[N]: same x-coodinate list
void read_and_init(){
  // line, ese-normalize
  int X, Y;
  cin >> X >> Y;
  int g = gcd(X, Y);
  assert(g != 0);
  X /= g;
  Y /= g;

  // S-scale
  scale = 2 * (X * X + Y * Y);

  // calc normalized coordinates
  vector<P> tmp(N);
  multimap<int, int> ids;
  int miny = INT_MAX;
  REP(i, N){
    int x, y;
    cin >> x >> y;
    // rotate
    P& p = tmp[i];
    p.x = Y * x - X * y;
    p.y = X * x + Y * y;
    ids.insert(make_pair(p.x, i));
    if(p.y < miny) miny = p.y;
  }
  // force p.y > 0
  REP(i, N){
    tmp[i].y -= miny - 1;
  }

  // shift index to make ps[0] leftmost
  pair<int, int> i0 = *ids.begin();
  int offset = i0.second;
  REP(i, N){
    int x = i - offset;
    if(x < 0) x += N;
    ps[x] = tmp[i];
  }

  // feature points
  fp.clear();
  REP(i, N){
    same[i].clear();
  }
  int prevx = INT_MIN;
  int previd = -1;
  FOR(it, ids){
    int x = it->first;
    int id = it->second - offset;
    if(id < 0) id += N;
    if(x != prevx){
      fp.push_back(id);
      prevx = x;
      previd = id;
    }else{
      same[previd].push_back(id);
    }
  }
  assert((int)fp.size() > 1);
  
  // calc slope and area
  REP(i, N){
    int j = i + 1;
    if(j == N) j = 0;
    int dx = ps[i].x - ps[j].x;
    int dy = ps[i].y - ps[j].y;
    if(dx != 0){
      int g = gcd(dx, dy);
      if(dx / g < 0) g = -g;
      slope[i].i = 0;
      slope[i].n = dy / g;
      slope[i].d = dx / g;
    }

    area[i][i] = (int64_t)0;
    area[i][j] = (int64_t)dx * (ps[i].y + ps[j].y);
  }
  // area i->j
  for(int di = 2; di < N; di++){
    REP(i, N){
      int j = i + di - 1;
      if(j >= N) j -= N;
      int k = j + 1;
      if(k == N) k = 0;
      area[i][k] = area[i][j] + area[j][k];
    }
  }
  AREA = area[0][1] + area[1][0];

#if 0
  cerr << "scale = " << scale << endl;
  REP(i, N){
    cerr << "P" << i << ": (" << ps[i].x << ", " << ps[i].y << ")" << endl;
    cerr << "slope = " << slope[i] << endl;
    cerr << "area =";
    REP(j, N){
      cerr << " " << area[i][j];
    }
    cerr << endl;
  }
  REP(i, fp.size()){
    cerr << "fp" << i << ": " << fp[i] << endl;
    cerr << "same =";
    REP(j, same[fp[i]].size()){
      cerr << " " << same[fp[i]][j];
    }
    cerr << endl;
  }
#endif
}

inline Area substitute(Quad& q, int x){
  return q.ddS * x * x + 2 * q.dS * x + q.S;
}

inline Area substitute(Quad& q, Solution x){
  return q.ddS * x * x + 2 * q.dS * x + q.S;
}

Area solve_one_equation(Quad& q0, Quad& q1, int lo, int hi){
  Coeff a = q1.ddS - q0.ddS;
  Coeff b = q1.dS - q0.dS;
  Coeff c = q1.S - q0.S;
#if 0
  cerr << "coeff: " << a << " " << b << " " << c << endl;
#endif
  Solution s;
  if(a == 0){
    s = -c / (b*2);
  }else{
    s = (-b - sqrt(b*b - a*c)) / a;
    if(s < lo || s > hi){
      s = (-b - sqrt(b*b - a*c)) / a;
    }
  }
#if 0
  cerr << "s = " << s << endl;
#endif
  return substitute(q0, s);
}

// 左側の連結成分(面積は単調増加)と右側の連結成分(面積は単調減少)に分けて考える
// 0<=x<=dXにおいてmax(f_i(x))をプロットすると、極小値を持つ単峰になる
// 左側グラフと右側グラフに交点がない場合、極小値は自明
// 交点が存在する場合、左側×右側の交点のy座標の最大値が求める値となる
//
// 計算すべき交点の数が単純計算で連結成分の数の半分の二乗になるので、
// オーダー表記的には変わらないが、高速になることが期待される
//
// 区間内に交点があるかどうかの判定はx=0とx=dXのときの値を見てやればよい
// 交点がある場合、二次方程式の解の計算はそれなりに遅くてO(sqrt)かかるとすると、
// この関数の計算量はO(交点の数 * sqrt)ぐらいだと見積もることができる
//
// 区間を一つ右に移動したときに面積の計算式が変わる連結成分は高々定数個
// 残りの連結成分同士の交点は定義域全体で高々二つずつであることを考えると、
// ならし計算量はO(N^2 * sqrt)ぐらいになっていると嬉しいなぁ
// O(N^3 * sqrt)になるような入力つくれません＞＜
//
// 今までの最小値をもらってきて枝刈りしてるので、
// 最小値をとるような区間を優先的に持ってくると速くなる
// O(N)ぐらいで最小値をとりそうな順に並べることができるだろうか?
Area solve_equations(vector<Quad>& left, vector<Quad>& right, int dX, Area& mini){
  Area maxi = reconstruct((int64_t)0);
  int nl = (int)left.size();
  int nr = (int)right.size();
  vector<Area> lS0(nl), lS1(nl);
  vector<Area> rS0(nr), rS1(nr);
  REP(i, nl){
    Quad& q = left[i];
    lS0[i] = substitute(q, 0);
    lS1[i] = substitute(q, dX);
    if(maxi < lS0[i]) maxi = lS0[i];
  }
  REP(i, nr){
    Quad& q = right[i];
    rS0[i] = substitute(q, 0);
    rS1[i] = substitute(q, dX);
    if(maxi < rS1[i]) maxi = rS1[i];
  }
  if(!(maxi < mini)) return maxi;
  REP(i, nl){
    if(lS1[i] < maxi) continue;
    REP(j, nr){
      if(rS0[j] < maxi) continue;
      if(!(lS0[i] < rS0[j] && rS1[j] < lS1[i])) continue;
#if 0
      cerr << i << " vs " << j << endl;
      cerr << left[i].S << " " << left[i].dS << " " << left[i].ddS << endl;
      cerr << right[j].S << " " << right[j].dS << " " << right[j].ddS << endl;
#endif
      Area a = solve_one_equation(left[i], right[j], 0, dX);
#if 0
      cerr << "a = " << a << endl;
#endif
      if(maxi < a){
        maxi = a;
        if(!(maxi < mini)) return maxi;
      }
    }
  }
  return maxi;
}

// もう日本語でいいや
// 左側(x座標の小さい側)から右側へ直線を動かしていく
// 各特徴点(x座標)と次の特徴点との間における最大値を求める
// 直線を特徴点からΔ右にずれたところに置けば、直線上に点が来ることはない
// Δずらしたことにするには、直線上の点を左側にあると思って計算すればよい
// 直線が特徴点にあるときの面積を求め、直線の傾きから増分を計算する
//
// 面積の計算は線をたどりながら積分計算
// ここでx座標は整数だがy座標が分数になるので、
// 結局全体では多倍長分数型が必要になるのではないか疑惑
// というか二次方程式を解くところでa+b*sqrt(c)型が要る
//
// 積分計算は、前計算で求めた頂点間の符号つき面積と、
// 直線から一番近い二頂点への辺の符号つき面積の和としてO(1)で求める
// 一回のiterationでO(直線と多角形の交点の数)になる
// 全体でならせばO(N^2)でナイーブな方法と変わらないが、定数倍速いはず
// 前計算で求める面積は64bit型に収まるので、1M個あっても高々8MB
Area solve(){
  vector<bool> isleft(N, false);
  Area mini = reconstruct(AREA);
  // Δ右にずらすと本来くっつかないとこがくっつくので、
  // そこのとこは別に計算する
  Area exact_maxi = reconstruct((int64_t)0);
  REP(_, fp.size()-1){
    int lfp = fp[_];
    int rfp = fp[_+1];
    isleft[lfp] = true;
    REP(i, same[lfp].size()){
      isleft[same[lfp][i]] = true;
    }
    int X0 = ps[lfp].x;
    int dX = ps[rfp].x - ps[lfp].x;

    // 直線と交差する辺IDを列挙
    // 0番の点を最左に置いたので、横切り方は右向き、左向き、右向き、…の順
    // 右から左に移動した点に関係するところでしか変化がないが、
    // どうせ交点を再計算するので全部調べる
    vector<F1> cp;
    map<F1, int> m[2];
    vector<int> ids;
    REP(i, N){
      int j = i + 1;
      if(j == N) j = 0;
      int l, r;
      if(isleft[i] && !isleft[j]){
        ids.push_back(i);
        l = i, r = j;
      }else if(!isleft[i] && isleft[j]){
        ids.push_back(i);
        l = j, r = i;
      }else{
        continue;
      }
      int ldx = X0 - ps[l].x;
      int rdx = ps[r].x - X0;
      int64_t num = (int64_t)ps[l].y * rdx + (int64_t)ps[r].y * ldx;
      int denom = ldx + rdx;
      assert(num > 0 && denom > 0 && num / denom < INT_MAX);
      F1 f;
      f.i = num / denom;
      f.n = num % denom;
      f.d = ldx + rdx;
      cp.push_back(f);
      int id = (int)ids.size()-1;
      m[id&1][f] = id;
    }

    // 多角形をたどっていって直線にぶつかって切られたとき、
    // 次にどの辺に遷移すればいいか(adj)を求める
    // 直線との交点のy座標が小さいものから順に、右向き、左向き、…が並んでいるので
    // 二つずつペアにして組み合わせる
    //
    // ついでに直線の切れ端の長さと右に動かしたときの増分を求める
    // 普通に長さはy座標の差、増分は傾きの差
    // 左と右で符号を逆にしておく
    int n = (int)ids.size();
    vector<int> adj(n);
    vector<F2> deltaS(n);
    vector<F2> deltadS(n);
    for(__typeof(m[0].begin()) it = m[0].begin(), jt = m[1].begin();
        it != m[0].end(); ++it, ++jt){
      F1 f0 = it->first;
      int lid = it->second;
      F1 f1 = jt->first;
      int rid = jt->second;
      F2 ds = f1 - f0;
      F2 dds = slope[ids[rid]] - slope[ids[lid]];

      adj[lid] = rid;
      deltaS[lid] = ds;
      deltadS[lid] = dds;
      adj[rid] = lid;
      deltaS[rid] = -ds;
      deltadS[rid] = -dds;
    }

    // 各連結成分の二次関数表現を求める
    // visitedでない辺からスタートして順にたどっていくだけ
    // ここから多倍長有理数型が必要になる
    // とりあえずvector<F2>からdecimalに変換する関数でお茶を濁す
    vector<Quad> objs[2];
    vector<bool> visited(n, false);
    REP(i, n){
      if(visited[i]) continue;
      vector<F2> Ss, dSs, ddSs;
      // 直線に沿って移動した先がp
      // そこから辺上を移動して直線にぶつかった点がq
      int q = i;
      do{
        int p = adj[q];
        int pv = ids[p] + 1;
        if(pv == N) pv = 0;
        q = p + 1;
        if(q == n) q = 0;
        int qv = ids[q];
#if 0
        cerr << "S0: " << cp[p] << " : " << ps[pv].y << " : " << X0 << " : " << ps[pv].x << endl;
        cerr << "S1: " << cp[q] << " : " << ps[qv].y << " : " << ps[qv].x << " : " << X0 << endl;
        cerr << "S2: " << area[pv][qv] << endl;
#endif
        F2 S = (cp[p] + ps[pv].y) * (X0 - ps[pv].x)
             + (cp[q] + ps[qv].y) * (ps[qv].x - X0)
             + area[pv][qv];
        F2 dS = deltaS[q];
        F2 ddS = deltadS[q];
        Ss.push_back(S);
        dSs.push_back(dS);
        ddSs.push_back(ddS);
        visited[q] = true;
      }while(q != i);
      Quad obj;
      obj.S = reconstruct(Ss);  // abs = int64, int64 + int64/int32 + int64/int32
      obj.dS = reconstruct(dSs);  // abs = int32, int32/int32 - int32/int32
      obj.ddS = reconstruct(ddSs);  // abs = int32, int32/int32 - int32/int32
      objs[i&1].push_back(obj);
#if 0
      cerr << "obj-" << ids[i] << ": " << obj.S << " " << obj.dS << " " << obj.ddS << endl;
#endif
    }

    // x = 0の場合を特別に計算
    REP(i, objs[1].size()){
      Area a = substitute(objs[1][i], 0);
      if(a > exact_maxi) exact_maxi = a;
    }
    if(exact_maxi < mini){
      mini = exact_maxi;
    }
    Area res = solve_equations(objs[0], objs[1], dX, mini);
    if(res < mini){
      mini = res;
    }
    exact_maxi = reconstruct((int64_t)0);
    REP(i, objs[0].size()){
      Area a = substitute(objs[0][i], dX);
      if(a > exact_maxi) exact_maxi = a;
    }
  }
  return mini;
}

int main(){
  while(cin >> N){
    if(N == 0) break;
    read_and_init();
    Area res = solve();
    cout << fixed << showpoint << setprecision(10)
         << to_decimal(res) / scale << endl;
  }

  return 0;
}