Go: package main

import (
	"fmt"
	"math"
	"math/cmplx"
)

func main() {
	for {
		var n int
		fmt.Scan(&n)

		if n == 0 {
			return
		}

		arr := make([]complex128, n)

		for i := 0; i < n; i++ {
			var x, y float64

			fmt.Scanf("%f%f", &x, &y)

			arr[i] = complex(x, y)
		}

		max := 1
		for i := 0; i < n; i++ {
			for j := i + 1; j < n; j++ {
				sub := arr[j] - arr[i]

				if cmplx.Abs(sub) >= 2 {
					continue
				}
				sum := arr[j] + arr[i]

				v := 1
				{
					halfd := cmplx.Abs(sub / 2)

					center := sum/2 + sub*cmplx.Rect(math.Sqrt(1-halfd*halfd)/cmplx.Abs(sub), float64(v)*math.Pi/2)

					counter := 0
					for k := 0; k < n; k++ {
						if i == k || j == k || cmplx.Abs(arr[k]-center) <= 1 {
							counter++
						}
					}
					if max < counter {
						max = counter
					}
				}

				v = -1
				{
					halfd := cmplx.Abs(sub / 2)

					center := sum/2 + sub*cmplx.Rect(math.Sqrt(1-halfd*halfd)/cmplx.Abs(sub), float64(v)*math.Pi/2)

					counter := 0
					for k := 0; k < n; k++ {
						if i == k || j == k || cmplx.Abs(arr[k]-center) <= 1 {
							counter++
						}
					}
					if max < counter {
						max = counter
					}
				}

			}
		}

		fmt.Println(max)
	}
}

C++: #include <iostream>
#include <cmath>
#include <vector>
#include <utility>
#include <stdio.h>

#if 0 
  #define SPAM(a) printf a
#else
  #define SPAM(a) (void)0
#endif

using namespace std;
const double EPS= 1e-5;
const double MAXDIST = 2.0;
const double MAXRAD  = 1.0;


typedef pair<double,double> Point;

//TODO check all function good for double
//DONE check printf good for double
vector<Point > pointList;

using namespace std;

int inRange(Point A, Point B, double num){
    double dx, dy, d;
    dx = A.first - B.first;
    dy = A.second - B.second;
    d = hypot(dx,dy);
    SPAM(("                    d = %lf\n", d));
    if (abs(num-d) < 0.00001 ){
            SPAM(("REALLY CLOSE\n"));
    }
    if ( EPS > d -num ){
        return 1;
    }
    return 0;
}

int get_Centre(Point pA,
               Point pB,
               Point &pC, 
               Point &pD )
{
   // SPAM(("in get Centre\n"));
  double x0 = pA.first;
  double y0 = pA.second;
  double x1 = pB.first;
  double y1 = pB.second;
  
  
  double a, dx, dy, d, h, rx, ry;
  double x2, y2;

  /* dx and dy are the vertical and horizontal distances between
   * the circle centers.
   */
  dx = x1 - x0;
  dy = y1 - y0;
  //SPAM(("dx = %lf, dy = %lf\n", dx, dyi));
  
  /* Determine the straight-line distance between the centers. */
  //d = sqrt((dy*dy) + (dx*dx));
  d = hypot(dx,dy); // Suggested by Keith Briggs
  //SPAM("         d = %lf\n", d);

  /* Check for solvability. */
  if (EPS > MAXDIST -d)
  {
    /* no solution. circles do not intersect. */
    return 0;
  }
    //SPAM("check solvability\n");

  /* 'point 2' is the point where the line through the circle
   * intersection points crosses the line between the circle
   * centers.  
   */

  /* Determine the distance from point 0 to point 2. */
  a = ( d / 2.0 ) ;
  //SPAM("a = %lf\n", a);
    
  /* Determine the coordinates of point 2. */
  x2 = (x0 + x1)/2.0 ;
  y2 = (y0 + y1)/2.0 ;

  /* Determine the distance from point 2 to either of the
   * intersection points.
   */
    //SPAM("before sqrt\n");
  //h = sqrt(1 - a*a);
  //SPAM("a*a = %lf\n", a*a);
  
  h = sqrtl(MAXRAD - a*a);
  //SPAM("h = %lf\n", h);
  //SPAM("after sqrt\n");
  /* Now determine the offsets of the intersection points from
   * point 2.
   */
  //SPAM("get Offset\n");
  //SPAM("          rx = -dy * (h / d) = %lf * (%lf /  %lf)\n", -dy, h, d);
  //SPAM("          ry = dx * (h / d) = %lf * (%lf /  %lf)\n", dx, h, d);
  rx = -dy * (h/d);
  ry = dx * (h/d);
  //SPAM("          rx = %lf, ry = %lf\n", rx, ry);
  //SPAM("get int points\n");
  /* Determine the absolute intersection points. */
  pC.first = x2 + rx;
  pC.second = y2 + ry;
  pD.first = x2 - rx;
  pD.second = y2 - ry;
  //SPAM("in getCentre pC = (%lf, %lf)\n", pC.first, pC.second);

  return 1;
}

int main(){
    //for each iteration until 0
    int N;
    scanf("%d", &N);
    while(N !=0){
        pointList.clear();
        //for each set of points
        for(int i=0; i < N; i++){
            //take in points 
            double X,Y;
            scanf("%lf %lf", &X, &Y);
            pointList.push_back(make_pair(X,Y));
        }        
            
        //SPAM("put all points in\n");
        //setMax = 0
        int setMax = 0;
        //for each point
        for(int i =0; i < N-1; i++){
            Point pA;
            pA = pointList[i];
            SPAM(("pA = (%lf, %lf)\n", pA.first, pA.second));
            //for each second point (point1 -> fin)
            for (int j=i+1; j < N; j++){
                Point pB = pointList[j];
                
                //SPAM("finding centre\n");
                //find centre
                Point pC, pD;
                //SPAM("--------> pB = (%lf,%lf)\n", pB.first, pB.second);
                //continue if distance > 2
                if( get_Centre(pA, pB, pC, pD) == 0){
                    //SPAM("         not in range\n");
                    continue;
                }
                //SPAM("out getCentre pC = (%lf, %lf)\n", pC.first, pC.second);
                //SPAM("after get_Centre\n");
                //circleMax = 0;
                int circleCMax = 0;
                //SPAM("          going through points\n");
                //go through all points 
                for(auto p: pointList){
                    //see if centre-> point < 1, circleMax++
                    SPAM(("                    p = (%lf, %lf)\n", p.first, p.second));
                    if ( inRange(p, pC, MAXRAD)){
                        circleCMax++;
                        SPAM(("                        in range\n"));
                    }
                }
                SPAM(("        circleCMax = %d\n\n", circleCMax));
                setMax = max(setMax, circleCMax);
                int circleDMax = 0;
                for(auto p: pointList){
                    //see if centre-> point < 1, circleMax++
                    SPAM(("                    p = (%lf, %lf)\n", p.first, p.second));
                    if ( inRange(p, pD, MAXRAD)){
                        circleDMax++;
                        SPAM(("                        in range\n"));
                    }
                }
                SPAM(("        circleDMax = %d\n\n", circleDMax));
                setMax = max(setMax, circleDMax);   
            }
                
        }
        //print setMax
        if(setMax == 0){
            printf("1\n");
        } else {
            printf("%d\n", setMax);
        }
        
        //scan next N in
        scanf("%d", &N);
    }
}

C++: #include <iostream>
#include <algorithm>
#include <cmath>
#include <complex>

using namespace std;

typedef complex<double> Point;

const double EPS = 1e-8;

int n;
Point data[300];

int calc(Point a, Point b){
  int res = 0;
  double d = abs(b - a) / 2.0;
  Point ab = (b - a) / (2.0 * d);
  if(d < 1){
    Point o = ab * d;
    d = sqrt(1.0 - d * d);
    o += Point(d * ab.imag() + a.real(), -d * ab.real() + a.imag());
    for(int i=0;i<n;i++){
      if(abs(o - data[i]) < 1.0 + EPS) res++;
    }
  }
  return res;
}

main(){
  while(cin >> n && n){
    for(int i=0;i<n;i++){
      cin >> data[i].imag() >> data[i].real();
    }
    int ans = 1;
    for(int i=0;i<n;i++){
      for(int j=i+1;j<n;j++){
        ans = max(ans, calc(data[i], data[j]));
        ans = max(ans, calc(data[j], data[i]));
      }
    }
    cout << ans << endl;
  }
}