Problem: Magic Circle

Problem Statement

~ How to draw a magic circle ~

1. Prepare a sufficiently large and white floor.
2. Draw a circle with radii 1, 2, ..., R centered on the point at the floor coordinates (0,0).
3. Paint blue between the circle with radius 1 and the circle with radius 2, between the circle with radius 3 and the circle with radius 4, ... Do not paint outside the circle with radius R.
4. Draw a polygon with N vertices on the floor.
5. Repaint the white area inside the polygon with blue and the blue area with white.



It is said that the power of the magic circle becomes stronger as the magic circle contains more blue.
Therefore, I would like you to find the area of ​​the blue area included in the magic circle.

The following figure is an example of a magic circle that can be drawn by this procedure.

<image>

Constraints

* 3 ≤ N ≤ 100
* 1 ≤ R ≤ 100
* Polygons are included inside or on the circumference of a circle with radius R.
* Polygons do not have self-intersections.
* The vertices of the polygon are given in a counterclockwise order.

Input

Input follows the following format. All given numbers are integers.


N R
x_1 y_1
...
x_N y_N

(x_i, y_i) is the i-th vertex of the polygon.

Output

Output the area of ​​the blue area on one line.
The output value must have an absolute or relative error of less than 10 ^ {-8} with the true value.

Examples

Input

3 1
0 0
1 0
0 1


Output

0.500000000


Input

3 2
0 0
2 0
0 2


Output

8.995574288


Input

3 2
1 -1
1 1
-2 0


Output

11.123246567


Input

4 3
1 1
-1 1
-1 -1
1 -1


Output

11.707963268

Solution (in CPP):

Solution:
#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-10)
#define equals(a,b) (fabs((a)-(b))<EPS)

using namespace std;

static const int COUNTER_CLOCKWISE = 1;
static const int CLOCKWISE = -1;
static const int ONLINE_BACK = 2;
static const int ONLINE_FRONT = -2;
static const int ON_SEGMENT = 0;


class Point{
public:
  double x,y;

  Point(double x = 0,double y = 0): x(x),y(y){}

  Point operator + (Point p){return Point(x+p.x,y+p.y);}
  Point operator - (Point p){return Point(x-p.x,y-p.y);}
  Point operator * (double a){return Point(a*x,a*y);}
  Point operator / (double a){return Point(x/a,y/a);}
  Point operator * (Point a){ return Point(x * a.x - y * a.y, x * a.y + y * a.x); }

  bool operator < (const Point& p) const{ return !equals(x,p.x)?x<p.x:y<p.y; }

  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }

};

struct Segment{
  Point p1,p2;
  Segment(Point p1 = Point(),Point p2 = Point()):p1(p1),p2(p2){}
  bool operator == (const Segment& p)const { return p.p1 == p1 && p.p2 == p2; }
};

typedef Point Vector;
typedef Segment Line;
typedef vector<Point> Polygon;

ostream& operator << (ostream& os,const Point& a){ os << "(" << a.x << "," << a.y << ")"; }

ostream& operator << (ostream& os,const Segment& a){ os << "( " << a.p1 << " , " << a.p2 << " )"; }

double dot(Point a,Point b){ return a.x*b.x + a.y*b.y; }

double cross(Point a,Point b){ return a.x*b.y - a.y*b.x; }

double norm(Point a){ return a.x*a.x+a.y*a.y; }

double abs(Point a){ return sqrt(norm(a)); }

//rad は角度をラジアンで持たせること
Point rotate(Point a,double rad){ return Point(cos(rad)*a.x - sin(rad)*a.y,sin(rad)*a.x + cos(rad)*a.y); }

// 度をラジアンに変換
double toRad(double agl){ return agl*M_PI/180.0; }

// a => prev, b => cur, c=> next
// prev から cur へ行って next へ行く際の角度を求める
double getArg(Point a,Point b,Point c){
  double arg1 = atan2(b.y-a.y,b.x-a.x);
  double arg2 = atan2(c.y-b.y,c.x-b.x);
  double arg = fabs( arg1 - arg2 );
  while( arg > M_PI ) arg -= 2.0 * M_PI;
  return fabs(arg);
}

int ccw(Point p0,Point p1,Point p2){
  Point a = p1-p0;
  Point b = p2-p0;
  if(cross(a,b) > EPS)return COUNTER_CLOCKWISE;
  if(cross(a,b) < -EPS)return CLOCKWISE;
  if(dot(a,b) < -EPS)return ONLINE_BACK;
  if(norm(a) < norm(b))return ONLINE_FRONT;
  return ON_SEGMENT;
}

Point projection(Line l,Point p) {
  double t = dot(p-l.p1, l.p1-l.p2) / norm(l.p1-l.p2);
  return l.p1 + (l.p1-l.p2)*t;
}

double distanceLP(Line l, Point p) {
  return abs(p - projection(l, p));
}

double cross3p(Point p,Point q,Point r){ return (r.x-q.x) * (p.y -q.y) - (r.y - q.y) * (p.x - q.x); }
  
bool collinear(Point p,Point q,Point r) { return fabs(cross3p(p,q,r)) < EPS; }
  
bool ccwtest(Point p,Point q,Point r){ return cross3p(p,q,r) > 0; }
 
bool onSegment(Point p,Point q,Point r){
  return collinear(p,q,r) && equals(sqrt(pow(p.x-r.x,2)+pow(p.y-r.y,2)) + sqrt(pow(r.x-q.x,2) + pow(r.y-q.y,2) ),sqrt(pow(p.x-q.x,2)+pow(p.y-q.y,2)) ) ;
}

#define pow2(a) ((a)*(a))
typedef pair<double,double> dd;
const double DINF = 1e20;
bool LT(double a ,double b){ return !equals(a,b) && a < b; }
bool LTE(double a,double b){ return  equals(a,b) || a < b; }


dd calc(double x1,double y1,double vx1,double vy1,
        double x2,double y2,double vx2,double vy2,double r){
  double VX = (vx1-vx2), X = (x1-x2), VY = (vy1-vy2), Y = (y1-y2);
  double a = pow2(VX) + pow2(VY), b = 2*(X*VX+Y*VY), c = pow2(X) + pow2(Y) - pow2(r);
  dd ret = dd(DINF,DINF);
  double D = b*b - 4 * a * c;

  if( LT(D,0.0) ) return ret;

  if( equals(a,0.0) ) {
    if( equals(b,0.0) ) return ret;
    if( LT(-c/b,0.0)  ) return ret;
    ret.first = - c / b;
    return ret;
  }

  if( equals(D,0.0) ) D = 0;
  ret.first  = ( -b - sqrt( D ) ) / ( 2 * a );
  ret.second = ( -b + sqrt( D ) ) / ( 2 * a );
  if( !equals(ret.first,ret.second) && ret.first > ret.second ) swap(ret.first,ret.second);
  return ret;
}


const Point ZERO = Point(0,0);
//線分ABと中心がcp,半径rの円との共通部分の面積を求める
inline double calculator_TypeA(Point A,Point B,Point cp,double r){
  A = A - cp, B = B - cp; 
  if( A == ZERO || B == ZERO ) return 0;
  double cross_value = cross(A,B);
  if( equals(cross_value,0.0) ) return 0;
  double sig = LT(cross_value,0.0) ? -1 : 1;
  Segment AB = Segment(A,B);
  double nearest_distance = distanceLP(AB,ZERO);
  double distance_OA = abs(A);
  double distance_OB = abs(B);

  if( LTE(0.0,r-distance_OA) && LTE(0.0,r-distance_OB) && LTE(0.0,r-nearest_distance) ) {
    return sig * fabs( cross_value / 2.0 );
  } else if( LTE(0.0,distance_OA-r) && LTE(0.0,distance_OB-r) && LTE(0.0,nearest_distance-r) ) {
    return sig * ( r * r * (M_PI-getArg(A,ZERO,B)) ) / 2.0;
  } else if( LTE(0.0,distance_OA-r) && LTE(0.0,distance_OB-r) && LT(0.0,r-nearest_distance) ) {

    Point proj_p = projection(AB,ZERO);
    if( onSegment(AB.p1,AB.p2,proj_p) ) {
      Vector e = ( A - B ) / abs( A - B );
      dd tmp = calc(A.x,A.y,e.x,e.y,0,0,0,0,r);
      Point r_p1 = A + e * tmp.first;
      Point r_p2 = A + e * tmp.second; 
      double ret = r * r * (M_PI-getArg(B,ZERO,A)) / 2.0;
      double subtract = r * r * (M_PI-getArg(r_p1,ZERO,r_p2)) / 2.0 - fabs(cross(r_p1,r_p2))/2.0 ;
      return sig * ( ret - subtract );
    } else {
      return sig * ( r * r * (M_PI-getArg(B,ZERO,A)) ) / 2.0; 
    }
  } else {

    if( LT(distance_OB-r,0.0) ) swap(A,B);
    Vector e = ( A - B ) / abs( A - B );
    dd tmp = calc(A.x,A.y,e.x,e.y,0,0,0,0,r);
    Point r_p1 = A + e * tmp.first;
    Point r_p2 = A + e * tmp.second; 
    if( onSegment(A,B,r_p2) ) r_p1 = r_p2;
    double ret = fabs(cross(r_p1,A)) * 0.5;
    ret += r * r * (M_PI-getArg(r_p1,ZERO,B)) * 0.5;
    return sig * ret;
  }
  assert(false);
}

double calc(const Polygon &poly,Point cp,double r){
  double sum = 0;
  rep(i,(int)poly.size()){  
    sum += calculator_TypeA(poly[i],poly[(i+1)%(int)poly.size()],cp,r);
  }
  return fabs(sum);
}


int N,R;
Polygon ps;

double compute(int R){ return calc(ps,ZERO,R) - ((R-1>0)?calc(ps,ZERO,R-1):0); }

int main(){
  cin >> N >> R;
  ps.resize(N);
  rep(i,N) cin >> ps[i].x >> ps[i].y;
  double area = 0;
  for(int r=2;r<=R;r+=2) area += ( r * r * M_PI ) - ( (r-1) * (r-1) * M_PI );
  for(int r=1;r<=R;r++){
    if( r & 1 ) area += compute(r);
    else        area -= compute(r);
  }
  printf("%.10lf\n",area);
  return 0;
}