You have a very large square wall and a circular dartboard placed on the wall. You have been challenged to throw darts into the board blindfolded. Darts thrown at the wall are represented as an array of points on a 2D plane. 
Return the maximum number of points that are within or lie on any circular dartboard of radius r.
 
Example 1:

Input: points = [[-2,0],[2,0],[0,2],[0,-2]], r = 2
Output: 4
Explanation: Circle dartboard with center in (0,0) and radius = 2 contain all points.

Example 2:

Input: points = [[-3,0],[3,0],[2,6],[5,4],[0,9],[7,8]], r = 5
Output: 5
Explanation: Circle dartboard with center in (0,4) and radius = 5 contain all points except the point (7,8).

Example 3:
Input: points = [[-2,0],[2,0],[0,2],[0,-2]], r = 1
Output: 1

Example 4:
Input: points = [[1,2],[3,5],[1,-1],[2,3],[4,1],[1,3]], r = 2
Output: 4

 
Constraints:

1 <= points.length <= 100
points[i].length == 2
-10^4 <= points[i][0], points[i][1] <= 10^4
1 <= r <= 5000

Difficulty: interview

Starter code:

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        

Solutions:
class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        ans = 1
        for x, y in points: 
            angles = []
            for x1, y1 in points: 
                if (x1 != x or y1 != y) and (d:=sqrt((x1-x)**2 + (y1-y)**2)) <= 2*r: 
                    angle = atan2(y1-y, x1-x)
                    delta = acos(d/(2*r))
                    angles.append((angle-delta, +1)) #entry
                    angles.append((angle+delta, -1)) #exit
            angles.sort(key=lambda x: (x[0], -x[1]))
            val = 1
            for _, entry in angles: 
                ans = max(ans, val := val+entry)
        return ans 



        
#https://www.geeksforgeeks.org/angular-sweep-maximum-points-can-enclosed-circle-given-radius/    
#        class Solution {
#    public int numPoints(int[][] points, int r) {
#        int count = 1;
#        for(int i = 0; i < points.length; i++) {
#            Map<Double, Integer> angles = new HashMap<>();
#            for(int j = 0; j < points.length; j++) {
#                if (i != j) {
#                    if (points[i][0] != points[j][0] || points[i][1] != points[j][1]) {
#                        int d = (points[i][0] - points[j][0]) * (points[i][0] - points[j][0]) + (points[i][1] - points[j][1]) * (points[i][1] - points[j][1]);
#                        double angle = Math.atan2(points[j][1] - points[i][1], points[j][0] - points[i][0]);
#                        double delta = acos(d / (2 * r));
#                        double entry = angle - delta;
#                        double exit = angle + delta;
#                        map.put(entry, map.getOrDefault(entry, 0) + 1);
#                        map.put(exit, map.getOrDefault(exit, 0) - 1);
#                    }
#                }
#            }
#            Map<String, Integer> result = angles.entrySet().stream()
#                .sorted(Map.Entry.comparingByValue(Comparator.reverseOrder()))
#                .collect(Collectors.toMap(Map.Entry::getKey, Map.Entry::getValue,
#                        (oldValue, newValue) -> oldValue, LinkedHashMap::new));
#            
#        }
#        return count;
#    }
#}


class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        res = 1
        for x, y in points:
            order = []
            for x1, y1 in points:
                if x != x1 or y != y1:
                    d = math.sqrt((x1 - x)**2 + (y1 - y)**2)
                    if d <= 2 * r:
                        delta = math.acos(d / (2 * r))
                        angle = math.atan2(y1 - y, x1 - x)
                        order.append((angle - delta, 1))
                        order.append((angle + delta, -1))
            order.sort(key=lambda x: (x[0], -x[1]))
            val = 1
            for _, entry in order: 
                res = max(res, val := val+entry)
            # count = 1
            # for _, entry in order:
            #     count += entry
            #      res = max(res, count)
        return res

from math import sqrt

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:

        n = len(points)
        if n <= 2:
            return n

        best = 1

        rr = r ** 2
        for i in range(n):
            for j in range(i + 1, n):
                [x1, y1] = points[i]
                [x2, y2] = points[j]

                xm = (x1 + x2) / 2
                ym = (y1 + y2) / 2

                q = sqrt((x2 - x1)**2 + (y2 - y1)**2)
                qq = (q / 2)**2
                rq = rr - qq
                if rq < 0:
                    continue

                xc = xm + sqrt(rq) * (y1 - y2) / q
                yc = ym + sqrt(rq) * (x2 - x1) / q

                curr = 2
                for k in range(n):
                    if k == i or k == j:
                        continue

                    [x, y] = points[k]
                    if sqrt((x - xc)**2 + (y - yc)**2) <= r:
                        curr += 1

                best = max(best, curr)

        return best





class Solution:
    def numPoints(self, A: List[List[int]], r: int) -> int:
            # def numPoints(self, A, r):
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
            d = ((x1 - x2) ** 2 + (y1 - y2) ** 2)/4.0
            if d > r * r: continue
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d) **0.5 / (d * 4) ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d) ** 0.5/(d * 4) **0.5
            #     x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            #     y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
            # res = max(res, sum((x-x0)**2 + (y - y0)**2 <= r * r + 0.0001 for x, y in A))
        return res
        # res = 1
        # for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
        #     d = ((x1 - x2)**2 + (y1 - y2)**2) / 4.0
        #     if d > r * r: continue
        #     x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
        #     y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
        #     res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
        # return res


class Solution:
    def numPoints(self, A: List[List[int]], r: int) -> int:
            # def numPoints(self, A, r):
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
            d = ((x1 - x2) ** 2 + (y1 - y2) ** 2)/4.0
            if d > r * r: continue
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d) **0.5 / (d * 4) ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d) **0.5 / (d * 4) ** 0.5
            #     x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            #     y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
            # res = max(res, sum((x-x0)**2 + (y - y0)**2 <= r * r + 0.0001 for x, y in A))
        return res
        # res = 1
        # for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
        #     d = ((x1 - x2)**2 + (y1 - y2)**2) / 4.0
        #     if d > r * r: continue
        #     x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
        #     y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
        #     res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
        # return res


class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        max_pts = 500
        dis = [[0 for _ in range(max_pts)] for _ in range(max_pts)]

        def getPointsInside(i, r, n):
            # This vector stores alpha and beta and flag
            # is marked true for alpha and false for beta
            angles = []

            for j in range(n):
                if i != j and dis[i][j] <= 2 * r:
                    # acos returns the arc cosine of the complex
                    # used for cosine inverse
                    B = math.acos(dis[i][j] / (2 * r))

                    # arg returns the phase angle of the complex
                    x1, y1 = points[i]
                    x2, y2 = points[j]
                    A = math.atan2(y1 - y2, x1 - x2)
                    alpha = A - B
                    beta = A + B
                    angles.append((alpha, False))
                    angles.append((beta, True))

            # angles vector is sorted and traversed
            angles.sort()
            # count maintains the number of points inside
            # the circle at certain value of theta
            # res maintains the maximum of all count
            cnt, res = 1, 1
            for angle in angles:
                # entry angle
                if angle[1] == 0: cnt += 1
                # exit angle
                else: cnt -= 1

                res = max(cnt, res)

            return res

        # Returns count of maximum points that can lie
        # in a circle of radius r.
        # dis array stores the distance between every
        # pair of points
        n = len(points)
        for i in range(n - 1):
            for j in range(i + 1, n):
                # abs gives the magnitude of the complex
                # number and hence the distance between
                # i and j
                x1, y1 = points[i]
                x2, y2 = points[j]
                dis[i][j] = dis[j][i] = sqrt((x1 - x2)**2 + (y1 - y2)**2)

        # This loop picks a point p
        ans = 0
        # maximum number of points for point arr[i]
        for i in range(n):
            ans = max(ans, getPointsInside(i, r, n))

        return ans

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        ans = 1
        for (x1, y1), (x2, y2) in itertools.combinations(points, 2):
            q = math.sqrt((x1-x2)**2+(y1-y2)**2)
            if q > 2*r + 0.00001:
                continue
            # only need to calculate a circle on one side
            # consider 1-D case
            x0 = (x1+x2)/2 + math.sqrt(r**2-(q/2)**2)*(y1-y2)/q
            y0 = (y1+y2)/2 + math.sqrt(r**2-(q/2)**2)*(x2-x1)/q
            #print(x0, x1)
            ans = max(ans, sum(1 for x, y in points if (x-x0)**2+(y-y0)**2 <= r**2 + 0.00001))
        return ans

import math
from typing import List


class Solution:
  def numPoints(self, points: List[List[int]], r: int) -> int:
    max_pts = 500
    dis = [[0 for _ in range(max_pts)] for _ in range(max_pts)]

    def getPointsInside(i, r, n):
      # This vector stores alpha and beta and flag
      # is marked true for alpha and false for beta
      angles = []

      for j in range(n):
        if i != j and dis[i][j] <= 2 * r:
          # acos returns the arc cosine of the complex
          # used for cosine inverse
          B = math.acos(dis[i][j] / (2 * r))

          # arg returns the phase angle of the complex
          x1, y1 = points[i]
          x2, y2 = points[j]
          A = math.atan2(y1 - y2, x1 - x2)
          alpha = A - B
          beta = A + B
          angles.append((alpha, False))
          angles.append((beta, True))

      # angles vector is sorted and traversed
      angles.sort()
      # count maintains the number of points inside
      # the circle at certain value of theta
      # res maintains the maximum of all count
      cnt, res = 1, 1
      for angle in angles:
        # entry angle
        if angle[1] == 0:
          cnt += 1
        # exit angle
        else:
          cnt -= 1

        res = max(cnt, res)

      return res

    # Returns count of maximum points that can lie
    # in a circle of radius r.
    # dis array stores the distance between every
    # pair of points
    n = len(points)
    for i in range(n - 1):
      for j in range(i + 1, n):
        # abs gives the magnitude of the complex
        # number and hence the distance between
        # i and j
        x1, y1 = points[i]
        x2, y2 = points[j]
        dis[i][j] = dis[j][i] = math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)

    # This loop picks a point p
    ans = 0
    # maximum number of points for point arr[i]
    for i in range(n):
      ans = max(ans, getPointsInside(i, r, n))

    return ans

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(points, 2):
            d = (x2 - x1) ** 2 + (y2 - y1) ** 2
            if d > 4 * r ** 2 or d == 0:
                continue
            
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d / 4) ** 0.5 / d ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d / 4) ** 0.5 / d ** 0.5
            res = max(res, sum((x - x0) ** 2 + (y - y0) ** 2 <= r * r + 0.00001 for x, y in points))
            
            x0 = (x1 + x2) / 2.0 - (y2 - y1) * (r * r - d / 4) ** 0.5
            y0 = (y1 + y2) / 2.0 + (x2 - x1) * (r * r - d / 4) ** 0.5
            res = max(res, sum((x - x0) ** 2 + (y - y0) ** 2 <= r * r + 0.00001 for x, y in points))
        return res

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        
        def count(x0, y0):
            nonlocal ans
            
            cnt = 0
            for x, y in points:
                if (x - x0) * (x - x0) + (y - y0) * (y - y0) <= r * r + 1e-6:
                    cnt += 1
            
            ans = max(ans, cnt)
            
        
        n = len(points)
        ans = 1
        for i in range(n):
            for j in range(i + 1, n):
                x1, y1 = points[i]
                x2, y2 = points[j]
                x_mid, y_mid = (x1 + x2) / 2, (y1 + y2) / 2
                sq = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)
                d = r * r - sq / 4
                if d < 0:
                    continue
                d = sqrt(d)
                sq = sqrt(sq)
                sin = abs(x1 - x2) / sq
                cos = abs(y1 - y2) / sq
                if (x2 - x1) * (y1 - y2) < 0:
                    cos = -cos
                x0, y0 = x_mid - d * cos, y_mid - d * sin
                count(x0, y0)
                x0, y0 = x_mid + d * cos, y_mid + d * sin
                count(x0, y0)
                
        return ans

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:        
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(points, 2):
            d = (x2 - x1) ** 2 + (y2 - y1) ** 2
            if d > 4 * r ** 2 or d == 0:
                continue
            
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d / 4) ** 0.5 / d ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d / 4) ** 0.5 / d ** 0.5
            res = max(res, sum((x - x0) ** 2 + (y - y0) ** 2 <= r * r + 0.00001 for x, y in points))
            
            x0 = (x1 + x2) / 2.0 - (y2 - y1) * (r * r - d / 4) ** 0.5
            y0 = (y1 + y2) / 2.0 + (x2 - x1) * (r * r - d / 4) ** 0.5
            res = max(res, sum((x - x0) ** 2 + (y - y0) ** 2 <= r * r + 0.00001 for x, y in points))
        return res        

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        res = 1
        for x, y in points:
            order = []
            for x1, y1 in points:
                if x != x1 or y != y1:
                    d = math.sqrt((x1 - x)**2 + (y1 - y)**2)
                    if d <= 2 * r:
                        delta = math.acos(d / (2 * r))
                        angle = math.atan2(y1 - y, x1 - x)
                        order.append((angle - delta, 1))
                        order.append((angle + delta, -1))
            order.sort(key=lambda x: (x[0], -x[1]))
            count = 1
            for _, entry in order:
                count += entry
                res = max(res, count)
        return res

class Solution:
    def numPoints(self, A, r):
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
            d = ((x1 - x2) ** 2 + (y1 - y2) ** 2) / 4.0
            if d > r * r:
                continue
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d) ** 0.5 / (d * 4) ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d) ** 0.5 / (d * 4) ** 0.5
            res = max(res, sum((x - x0) ** 2 + (y - y0) ** 2 <= r * r + 0.00001 for x, y in A))

            x0 = (x1 + x2) / 2.0 - (y2 - y1) * (r * r - d) ** 0.5
            y0 = (y1 + y2) / 2.0 + (x2 - x1) * (r * r - d) ** 0.5
            res = max(res, sum((x - x0) ** 2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
        return res

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        
        n=len(points)
        ans=1
        for i in range(n-1):
            x1,y1=points[i][0],points[i][1]
            for j in range(i+1,n):
                x2,y2=points[j][0],points[j][1]
                d=(x1-x2)**2+(y1-y2)**2
                
                if d>4*r*r:
                    continue
                
                cx,cy=(x1+x2)/2,(y1+y2)/2
                h=math.sqrt(r*r-d/4)
                
                
                
                if y1==y2:
                    dx,dy=0,h
                
                else:
                    m=-(x1-x2)/(y1-y2)
                    dx,dy=h/math.sqrt(1+m*m),h*m/(math.sqrt(1+m*m))
                    
                for sign in [1,-1]:
                    count=0
                    px,py=cx+sign*dx,cy+sign*dy
                    
                    for k,(x,y) in enumerate(points):
                        if (x-px)**2+(y-py)**2<=r*r+1e-7:
                            count+=1
                    ans=max(ans,count)
                    
        return ans 
                    
                    
                    
                    
                    
                


class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        max_pts = 500
        dis = [[0 for _ in range(max_pts)] for _ in range(max_pts)]

        def getPointsInside(i, r, n):
            # This vector stores alpha and beta and flag
            # is marked true for alpha and false for beta
            angles = []

            for j in range(n):
                if i != j and dis[i][j] <= 2 * r:
                    # acos returns the arc cosine of the complex
                    # used for cosine inverse
                    B = math.acos(dis[i][j] / (2 * r))

                    # arg returns the phase angle of the complex
                    x1, y1 = points[i]
                    x2, y2 = points[j]
                    A = math.atan2(y1 - y2, x1 - x2)
                    alpha = A - B
                    beta = A + B
                    angles.append((alpha, False))
                    angles.append((beta, True))

            # angles vector is sorted and traversed
            angles.sort()
            # count maintains the number of points inside
            # the circle at certain value of theta
            # res maintains the maximum of all count
            cnt, res = 1, 1
            for angle in angles:
                # entry angle
                if angle[1] == 0: cnt += 1
                # exit angle
                else: cnt -= 1

                res = max(cnt, res)

            return res

        # Returns count of maximum points that can lie
        # in a circle of radius r.
        # dis array stores the distance between every
        # pair of points
        n = len(points)
        for i in range(n - 1):
            for j in range(i + 1, n):
                # abs gives the magnitude of the complex
                # number and hence the distance between
                # i and j
                x1, y1 = points[i]
                x2, y2 = points[j]
                dis[i][j] = dis[j][i] = sqrt((x1 - x2)**2 + (y1 - y2)**2)

        # This loop picks a point p
        ans = 0
        # maximum number of points for point arr[i]
        for i in range(n):
            ans = max(ans, getPointsInside(i, r, n))

        return ans        

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        best = 1
        for i in range(len(points)):
            p1 = complex(*points[i])
            for j in range(i+1, len(points)):
                p2 = complex(*points[j])
                d = abs(p1 - p2)
                if d > 2 * r:
                    continue
                h = math.sqrt(r * r - d * d / 4)
                c = (p1 + p2) / 2 + (p1 - p2) * h / d * 1j
                
                count = 0
                for x, y in points:
                    if (x - c.real) ** 2 + (y - c.imag) ** 2 <= r ** 2 + 1e-6:
                        count += 1
                best = max(best, count)
        return best

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        eps = 1e-9
        
        ret = 1
        n = len(points)
        if n == 1:
            return ret
        
        points.sort()
        
        def isect(x1, y1, x2, y2):
            dx2 = (x2 - x1) ** 2
            dy2 = (y2 - y1) ** 2
            if dx2 + dy2 > 4 * r * r:
                return []
            cx = (x1 + x2) / 2
            cy = (y1 + y2) / 2
            if dx2 + dy2 == 4 * r * r:
                return [(cx, cy)]
            ss2 = (dx2 + dy2) / 4
            ol = math.sqrt(r ** 2 - ss2)
            diffx = ol * math.sqrt(dy2 / (dx2 + dy2))
            diffy = ol * math.sqrt(dx2 / (dx2 + dy2))
            return [(cx - diffx, cy + diffy), (cx + diffx, cy - diffy)]
        
        for i in range(n-1):
            a, b = points[i]
            for j in range(i+1, n):
                c, d = points[j]
                l = isect(a, b, c, d)
                for x, y in l:
                    cur = 0
                    lb = bisect.bisect_left(points, [x-r-eps,y-r-eps])
                    ub = bisect.bisect_right(points, [x+r+eps,y+r+eps])
                    for k in range(lb, ub):
                        dx = points[k][0] - x
                        dy = points[k][1] - y
                        if dx ** 2 + dy ** 2 <= (r + eps) ** 2:
                            cur += 1
                    ret = max(ret, cur)
        return ret


class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        \"\"\"
        :type points: List[List[int]]
        :type r: int
        :rtype: int
        \"\"\"
        res = 1
        for i in range(len(points)):
            for j in range(i + 1, len(points)):
                x1 = points[i][0]
                y1 = points[i][1]
                x2 = points[j][0]
                y2 = points[j][1]
                d = ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) ** 0.5
                if d > r * 2:
                   continue
                _x = (y2 - y1) * ((r * r - d * d/ 4.0) ** 0.5) / d
                _y = (x2 - x1) * ((r * r - d * d/ 4.0) ** 0.5) / d
                res = max(res, self.count(points, (x1 + x2) / 2.0 + _x, (y1 + y2) / 2.0 - _y, r))
                res = max(res, self.count(points, (x1 + x2) / 2.0 - _x, (y1 + y2) / 2.0 + _y, r))
        return res
        
    def count(self, points, x, y, r):
        count = 0
        for a, b in points:
            d = ((a - x)*(a - x) + (b - y)*(b - y))**0.5
            if d <= r + 0.0001 :
                count += 1
        return count

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        @functools.lru_cache(None)
        def find(cx, cy, r):
            ans = 0
            for x, y in points:
                if (x-cx) ** 2 + (y-cy) ** 2 <= r ** 2 + 1e-6:
                    ans += 1
            return ans
        
        ans = 1
        for i in range(len(points)):
            x1, y1 = points[i]
            for j in range(i+1, len(points)):
                x2, y2 = points[j]
                deno = (y1-y2) ** 2 + (x1-x2) ** 2
                if deno > 4 * r * r:
                    continue
                k = math.sqrt(r * r / deno - 0.25)
                ans = max(ans, max(find((x1+x2)/2 + t*(y1-y2), (y1+y2)/2 - t * (x1-x2), r) for t in (k, -k)))
        return ans

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        def find(cx, cy, r):
            ans = 0
            for x, y in points:
                if (x-cx) ** 2 + (y-cy) ** 2 <= r ** 2 + 1e-6:
                    ans += 1
            return ans
        
        ans = 1
        for i in range(len(points)):
            x1, y1 = points[i]
            for j in range(i+1, len(points)):
                x2, y2 = points[j]
                deno = (y1-y2) ** 2 + (x1-x2) ** 2
                if deno > 4 * r * r:
                    continue
                k = math.sqrt(r * r / deno - 0.25)
                ans = max(ans, max(find((x1+x2)/2 + t*(y1-y2), (y1+y2)/2 - t * (x1-x2), r) for t in (k, -k)))
        return ans

import math as m
class Solution(object):
    def dist(self, p1, p2):
        dx = p1[0] - p2[0]
        dy = p1[1] - p2[1]
        return m.sqrt(dx * dx + dy * dy)
    
    def intersect(self, p1, p2, r):
        res = []
        d = self.dist(p1, p2) - 2 * r
        if (d > 0):
            res = []
        elif (d == 0):
            res = [[(p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2]] 
        else:
            mid_x = (0.0 + p1[0] + p2[0]) / 2
            mid_y = (0.0 + p1[1] + p2[1]) / 2
            # ep ^ 2 + fun ^ 2 == r ^ 2 - dist(mid, p1) ^ 2;
            # ep * (p1[0] - p2[0]) + fun * (p1[1] - p2[1]) == 0 -> fun = -(p1[0] - p2[0]) / (p1[1] - p2[1])  * eps;
            if (p1[1] != p2[1]):
                ratio = - (0.0 + p1[0] - p2[0]) / (0.0 + p1[1] - p2[1])
                eps = m.sqrt((r ** 2 - self.dist([mid_x, mid_y], p1) ** 2) / (ratio ** 2 + 1))
                fun =  eps * ratio
            else:
                eps = 0
                fun = m.sqrt(r ** 2 - self.dist([mid_x, mid_y], p1) ** 2)
            # res update    
            res  =[[mid_x + eps, mid_y + fun], [mid_x - eps, mid_y - fun]]
        return res    
            
    def numPoints(self, points, r):
        l = len(points)
        result = 1
        for i in range(l):
            for j in range(i + 1, l):
                c = self.intersect(points[i], points[j], r)
                for p in c:
                    au = 0
                    for k in range(l):
                        if (self.dist(p, points[k]) <= r): au += 1
                    result = max(result, au)
        return result
            


import math
class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        ##We have that the circle that will contain the most points can be chosen so that at least one point lies on the circumference. Indeed, suppose the best circle contains no points on circumference. Then, we can do a contraction of the circle at center until it hits the first point (still be maximum number of points). Then we can do a dilation from that first point till we have a circle of original radius, and we see that it any point that was originally in the smaller circle will be in the larger dilated circle.
        ##So we need to essentially ask the question if a circle of radius r passes through each given point, what is maximum number of points the circle can contain then. Then over all points we take the maximum.
        self.total = len(points)
        maxPoints = 0
        for i in range(self.total):
            maxPoints = max(maxPoints,self.numPointsInCircle(points,r,i))
        return maxPoints
    def numPointsInCircle(self,points,r,chosenIndex):
        maxPoints = 0
        enterList = []
        exitList = []
        numberPoints = 0
        for i in range(self.total):
            answer = self.computeRange(points,r,chosenIndex,i)
            if answer != None:
                if len(answer) == 0:
                    numberPoints += 1
                else:
                    enterAngle,exitAngle = answer
                    enterList.append(enterAngle)
                    exitList.append(exitAngle)
                    if enterAngle > exitAngle:##includes case where the enterAngle = pi and exitAngle = -pi, because at exact angle of -pi we have the point included.
                        numberPoints += 1
        exitList.sort()
        enterList.sort()
        maxPoints = numberPoints
        exitCounter = 0
        enterCounter = 0
        pointsToCheck = len(exitList)
        while max(exitCounter,enterCounter) < pointsToCheck:
            currentExit = exitList[exitCounter]
            currentEnter = enterList[enterCounter]
            if currentEnter <= currentExit:
                numberPoints += 1
                maxPoints = max(maxPoints,numberPoints)
                enterCounter += 1
            else:
                numberPoints -= 1
                exitCounter += 1
        while exitCounter < pointsToCheck:
            numberPoints -= 1
            exitCounter += 1
        while enterCounter < pointsToCheck:
            numberPoints += 1
            maxPoints = max(maxPoints,numberPoints)
            enterCounter += 1
        return maxPoints
        
        
        
    def computeRange(self,points,r,chosenIndex,otherIndex):##
        chosenPair = points[chosenIndex]
        otherPair = points[otherIndex]
        if chosenPair == otherPair:
            return []##this will be used to indicate that this will not be considered in computeEnterExit
        else:
            distance = self.distanceBetweenPoints(chosenPair,otherPair)
            if distance > 2*r:
                return None
            else:
                angleOne = self.computeAngle(chosenPair,otherPair)
                angleTwo = math.acos(distance/(2*r))##remember which objects and modules you are using.
                exit = angleOne+angleTwo
                exit = exit-2*math.pi if exit >= math.pi else exit##assures exit is in range of [-math.pi,math.pi)
                enter = angleOne-angleTwo
                enter = enter + 2*math.pi if enter < -math.pi else enter ##assures enter is in range of (-math.pi,math.pi]
                return [enter,exit]
    def computeAngle(self,pairOne,pairTwo):
        x1,y1 = pairOne
        x2,y2 = pairTwo
        return math.atan2((y2-y1),(x2-x1))
    def numberPointsStart(self,enterList,exitList):
        pass
    def distanceBetweenPoints(self,pairOne,pairTwo):
        x1,y1 = pairOne
        x2,y2 = pairTwo
        return math.sqrt((x1-x2)**2+(y1-y2)**2)

import math as m
class Solution(object):
    def dist(self, p1, p2):
        dx = p1[0] - p2[0]
        dy = p1[1] - p2[1]
        return m.sqrt(dx * dx + dy * dy)
    
    def intersect(self, p1, p2, r):
        res = []
        d = self.dist(p1, p2) - 2 * r
        if (d > 0):
            res = []
        elif (d == 0):
            res = [[(p1[0] + p2[0]) / 2, (p1[1] + p2[1]) / 2]] 
        else:
            mid_x = (0.0 + p1[0] + p2[0]) / 2
            mid_y = (0.0 + p1[1] + p2[1]) / 2
            # ep ^ 2 + fun ^ 2 == r ^ 2 - dist(mid, p1) ^ 2;
            # ep * (p1[0] - p2[0]) + fun * (p1[1] - p2[1]) == 0 -> fun = -(p1[0] - p2[0]) / (p1[1] - p2[1])  * eps;
            if (p1[1] != p2[1]):
                ratio = - (0.0 + p1[0] - p2[0]) / (0.0 + p1[1] - p2[1])
                eps = m.sqrt((r ** 2 - self.dist([mid_x, mid_y], p1) ** 2) / (ratio ** 2 + 1))
                fun =  eps * ratio
            else:
                eps = 0
                fun = m.sqrt(r ** 2 - self.dist([mid_x, mid_y], p1) ** 2)
            # res update    
            res  =[[mid_x + eps, mid_y + fun], [mid_x - eps, mid_y - fun]]
        return res    
            
    def numPoints(self, points, r):
        l = len(points)
        result = 1
        for i in range(l):
            for j in range(i + 1, l):
                c = self.intersect(points[i], points[j], r)
                for p in c:
                    au = 0
                    for k in range(l):
                        if (self.dist(p, points[k]) <= r): au += 1
                    result = max(result, au)
        return result

import math
from decimal import Decimal as D
class Solution:
#     def numPoints(self, points: List[List[int]], r: int) -> int:
#         def find_intersect(point1,point2,r):
#             x1 = point1[0]
#             y1 = point1[1]
#             x2 = point2[0]
#             y2 = point2[1]
#             d = math.sqrt((x1-x2)**2 + (y1-y2)**2)/2
#             ans = set([])
#             if d <= r:
#                 x = (x1+x2)/2
#                 y = (y1+y2)/2
#                 de = math.sqrt(r**2-d**2)
#                 #ans.add((x + de*(y1-y)/d, y + de*(x-x1)/d))
#                 #ans.add((x - de*(y1-y)/d, y - de*(x-x1)/d))
#                 tmp = sqrt(4*(r**2-d**2)/((x1-x2)**2 + (y1-y2)**2))
#                 print(tmp)
#                 ans.add((x + tmp*(y1-y), y + tmp*(x-x1)))
#                 ans.add((x - tmp*(y1-y), y - tmp*(x-x1)))
                
#             return ans
        
#         record = {}
#         for i in range(len(points)-1):
#             for j in range(i+1,len(points)):
#                 point1 = points[i]
#                 point2 = points[j]
#                 ans = find_intersect(point1,point2,r)
#                 #print(ans,point1,point2)
#                 for a in ans:
#                     if a not in record:
#                         record[a] = set([i,j])
#                     else:
#                         record[a].add(i)
#                         record[a].add(j)
       
        
#         return max([len(record[k]) for k in record])
    def numPoints(self, points: List[List[int]], r: int) -> int:        
        def circles_from_p1p2r(p1, p2, r):
            (x1, y1), (x2, y2) = p1, p2

            # delta x, delta y between points
            dx, dy = x2 - x1, y2 - y1

            # dist between points
            q = math.sqrt(dx**2 + dy**2)

            # if two points are too far away, there is no such circle
            if q > 2.0*r:
                return []

            # find the halfway point
            x3, y3 = (x1+x2)/2, (y1+y2)/2

            # distance along the mirror line
            d = math.sqrt(r**2-(q/2)**2)

            # One circle
            c1 = [x3 - d*dy/q, y3 + d*dx/q]

            # The other circle
            c2 = [x3 + d*dy/q, y3 - d*dx/q]
            return [c1, c2]

        # now the main function
        res = 0
        n = len(points)

        for p in range(n):
            for q in range(p+1, n):

                # Find the two candidate circle for each pair of points
                TwoCirs = circles_from_p1p2r(points[p], points[q], r)

                # count how many dots are inside the circle
                for center in TwoCirs: # center = TwoCirs[1]
                    cnt = 0

                    for dots in points: # dots =points[0]
                        if (dots[0]-center[0])**2 + (dots[1]-center[1])**2 <= r**2+10**-6:
                            cnt += 1
                    res = max(res, cnt)

        return res if res else 1

class Solution:
    def getC(self, p1, p2, r):
        [x1, y1], [x2, y2] = p1, p2
        dist = sqrt(((x1-x2)**2)+((y1-y2)**2))
        if(dist>2.0*r):
            return []
        x3, y3 = (x1+x2)/2, (y1+y2)/2
        magnitude = sqrt((r**2)-((dist/2)**2))
        c1 = [x3 - magnitude*(y1-y2)/dist, y3 + magnitude*(x1-x2)/dist]
        c2 = [x3 + magnitude*(y1-y2)/dist, y3 - magnitude*(x1-x2)/dist]
        return [c1,c2]
        
    def numPoints(self, points: List[List[int]], r: int) -> int:
        ans=0
        for i in range(len(points)):
            for j in range(i+1, len(points)):
                C = self.getC(points[i], points[j], r)
                for c in C:
                    cnt=0
                    for pt in points:
                        if(((pt[0]-c[0])**2)+((pt[1]-c[1])**2)<=(r**2)+1e-6):
                            cnt+=1
                    ans = max(ans, cnt)
        return ans if ans else 1


class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        
        def circle_centers(p1, p2):
            [x1,y1], [x2,y2] = p1, p2
            half_dist_2 = ((x2-x1) ** 2 + (y2-y1) ** 2) / 4
            
            offset2 = r**2 - half_dist_2
            midx, midy = (x1+x2)/2, (y1+y2)/2
            
            if abs(offset2) < 1e-4:
                return [[midx, midy]]
            
            if offset2 < 0:
                return []
            
            offset = offset2 ** 0.5
            dist = ((x2-x1) ** 2 + (y2-y1) ** 2) ** .5
            
            
            dx, dy = offset * (y1-y2) / dist, offset * (x2-x1) / dist
            return [[midx+dx,midy+dy], [midx-dx,midy-dy]]
        
        
        def count_members_around(cx, cy):
            print(cx, cy)
            return len([0 for px, py in points if (px-cx)**2 + (py-cy)**2 < r**2 + 1e-4])
        
        result = [count_members_around(cx, cy)\\
                    for i in range(len(points))\\
                    for j in range(i+1, len(points))\\
                    for cx, cy in circle_centers(points[i], points[j])]
        
        if len(result) == 0:
            return 1
        print(result)
        return max(result)



                    

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        def circles_from_p1p2r(p1, p2, r):
            (x1, y1), (x2, y2) = p1, p2


            # delta x, delta y between points
            dx, dy = x2 - x1, y2 - y1


            # dist between points
            q = math.sqrt(dx**2 + dy**2)


            # if two points are too far away, there is no such circle
            if q > 2.0*r:
                return []


            # find the halfway point
            x3, y3 = (x1+x2)/2, (y1+y2)/2


            # distance along the mirror line
            d = math.sqrt(r**2-(q/2)**2)


            # One circle
            c1 = [x3 - d*dy/q, y3 + d*dx/q]


            # The other circle
            c2 = [x3 + d*dy/q, y3 - d*dx/q]
            return [c1, c2]


        # now the main function
        res = 0
        n = len(points)


        for p in range(n):
            for q in range(p+1, n):


                # Find the two candidate circle for each pair of points
                TwoCirs = circles_from_p1p2r(points[p], points[q], r)


                # count how many dots are inside the circle
                for center in TwoCirs: # center = TwoCirs[1]
                    cnt = 0


                    for dots in points: # dots =points[0]
                        if (dots[0]-center[0])**2 + (dots[1]-center[1])**2 <= r**2+10**-6:
                            cnt += 1
                    res = max(res, cnt)

        return res if res else 1

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        def valid(x1,y1,x2,y2,r):
            L=math.sqrt((y2-y1)**2+(x2-x1)**2)
            if 2*r<L:
                return False
            else:
                return True
        def getcenter(x1,y1,x2,y2,r):
            L=math.sqrt((y2-y1)**2+(x2-x1)**2)
            if r**2-(L/2)**2<0:
                return None
            d=math.sqrt(r**2-(L/2)**2)
            a,b=d*(y1-y2)/L, d*(x2-x1)/L
            midx=(x1+x2)*0.5
            midy=(y1+y2)*0.5
            return midx+a,midy+b,midx-a,midy-b
        def check(x,y,r):
            count=0
            for u in points:
                if (u[0]-x)**2+(u[1]-y)**2<=r**2+0.001:
                    count+=1
            return count
        L=len(points)
        if L<=1:
            return L
        ans=1
        for i in range(L-1):
            for j in range(i+1,L):
                if valid(points[i][0],points[i][1],points[j][0],points[j][1],r):
                    cx1,cy1,cx2,cy2=getcenter(points[i][0],points[i][1],points[j][0],points[j][1],r)
                    print(cx1,cy1,cx2,cy2)
                    ans=max(ans,check(cx1,cy1,r))
                    ans=max(ans,check(cx2,cy2,r))
        return ans

class Solution:
    # O(n^3) time, O(1) space
    def numPoints(self, points: List[List[int]], r: int) -> int:
        n = len(points)
        result = 1
        
        def calc_dist(pt1, pt2):
            x1, y1 = pt1
            x2, y2 = pt2
            return math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
        
        def get_centers(i, j):
            d_xy = calc_dist(points[i], points[j])
            if d_xy > r * 2:
                return []
            delta_x = points[j][0] - points[i][0]
            delta_y = points[j][1] - points[i][1]
            x_mid = points[i][0] + delta_x / 2
            y_mid = points[i][1] + delta_y / 2
            if d_xy == r * 2:
                return [(x_mid, y_mid)]
            
            d = math.sqrt(r ** 2 - (d_xy / 2) ** 2)
            return [
                (x_mid + d * delta_y / d_xy, y_mid - d * delta_x / d_xy),
                (x_mid - d * delta_y / d_xy, y_mid + d * delta_x / d_xy)
            ]
        
        def count_points_inside(center, i, j):
            count = 2
            for idx in range(n):
                if idx == i or idx == j:
                    continue
                count += calc_dist(center, points[idx]) <= r
            return count
        
        for i in range(n):  # O(n)
            for j in range(i + 1, n):  # O(n)
                centers = get_centers(i, j)
                for center in centers:
                    result = max(result, count_points_inside(center, i, j))  # O(n)
        
        return result


class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:

        def within(pt, cen, r2):
            a1, a2 = pt
            p1, p2 = cen
            return (a1 - p1)**2 + (a2 - p2)**2 <= r2

        ans = 1

        for idx, a in enumerate(points):
            for b in points[:idx]:
                a1, a2 = a
                b1, b2 = b

                if (a1-b1)**2+(a2-b2)**2 <= 4*r**2:
                    pos_cir = 0
                    neg_cir = 0

                    mid = ((a1+b1)/2, (a2+b2)/2)
                    perp = (b2-a2, a1-b1)
                    perp_dis = (perp[0]**2 + perp[1]**2)**.5
                    p_dis = (r**2 - ((a1-b1)**2+(a2-b2)**2)/4)**.5

                    cen = (mid[0]+perp[0]*p_dis/perp_dis,
                           mid[1]+perp[1]*p_dis/perp_dis)

                    ocen = (mid[0]+perp[0]*p_dis/perp_dis,
                            mid[1]+perp[1]*p_dis/perp_dis)

                    for c in points:

                        if within(c, cen, r**2+0.00000001):
                            pos_cir += 1
                        if within(c, ocen, r**2+0.00000001):
                            neg_cir += 1

                    ans = max(ans, pos_cir, neg_cir)

        return ans

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        def getpoints(l,i,r,n,dis):
            d={}
            pi=3.1415626
            angles=[]
            for j in range(n):
                if i!=j and dis[i][j]<=2*r:
                    #print(dis[i][j]/(2*r))
                    B=math.acos(dis[i][j]/(2*r))
                    A=math.acos((l[j][0]-l[i][0])/dis[i][j])
                    if l[j][1]-l[i][1]<0:
                        A=2*pi-A
                    a=A-B
                    b=A+B
                    #print(i,j,(A/3.14)*180,(B*180/3.14))
                    angles.append((a,0))  #start
                    angles.append((b,1))  #end

            angles.sort()
            #print(l[i],angles)
            ct=1
            res=1
            for a in angles:
                if a[1]==0:
                    ct+=1
                else:
                    ct-=1
                res=max(res,ct)
            #print(i,res)
            return res
        n=len(points)
        dis=[[0 for j in range(n)] for i in range(n)]
        
        for i in range(n-1):
            for j in range(i+1,n):
                dis[i][j]=((points[i][0]-points[j][0])**2+(points[i][1]-points[j][1])**2)**(0.5)
                dis[j][i]=((points[i][0]-points[j][0])**2+(points[i][1]-points[j][1])**2)**(0.5)
        ans=0
        #print(dis)
        #print(dis[0][3],dis[3][0])
        for i in range(n):
            ans=max(ans,getpoints(points,i,r,n,dis))
        return ans

class Solution:
     def numPoints(self, A, r):
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
            d = ((x1 - x2)**2 + (y1 - y2)**2) / 4.0
            if d > r * r: continue
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
        return res

class Solution:
    def numPoints(self, points: List[List[int]], r: int) -> int:
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(points, 2):
            d = ((x1 - x2)**2 + (y1 - y2)**2) / 4.0
            if d > r * r: continue
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in points))
        return res
                


import itertools

class Solution:
  def numPoints(self, points: List[List[int]], r: int) -> int:
    # http://www4.comp.polyu.edu.hk/~csbxiao/paper/2004/ISPAN04_cover.pdf
    # TC: O(N^3), possible O(N^2logN), SC: O(1)
    ans = 1
    for (x1, y1), (x2, y2) in itertools.combinations(points, 2):
      d = ((x1 - x2) ** 2 + (y1 - y2) ** 2) / 4.0
      if d <= r * r:
        x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d) ** 0.5 / (d * 4) ** 0.5
        y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d) ** 0.5 / (d * 4) ** 0.5
        ans = max(ans, sum((x - x0) ** 2 + (y - y0) ** 2 <= r * r + 0.00001 for x, y in points))
    return ans

class Solution:
\tdef numPoints(self, P: List[List[int]], r: int) -> int:
\t\tans = 1
\t\tfor (x1, y1), (x2, y2) in itertools.combinations(P, 2):
\t\t\td = ((x1 - x2)**2 + (y1 - y2)**2) / 4.0
\t\t\tif d > r * r:
\t\t\t\tcontinue
\t\t\tx0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d) ** 0.5 / (d * 4) ** 0.5
\t\t\ty0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d) ** 0.5 / (d * 4) ** 0.5
\t\t\tans = max(ans, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in P))
\t\treturn ans
\t\t\t


class Solution:
    def numPoints(self, A: List[List[int]], r: int) -> int:
            # def numPoints(self, A, r):
        res = 1
        for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
            d = ((x1 - x2) ** 2 + (y1 - y2) ** 2)/4.0
            if d > r * r: continue
            x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
            res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
            # res = max(res, sum((x-x0)**2 + (y - y0)**2 <= r * r + 0.0001 for x, y in A))
        return res
        # res = 1
        # for (x1, y1), (x2, y2) in itertools.combinations(A, 2):
        #     d = ((x1 - x2)**2 + (y1 - y2)**2) / 4.0
        #     if d > r * r: continue
        #     x0 = (x1 + x2) / 2.0 + (y2 - y1) * (r * r - d)**0.5 / (d * 4) ** 0.5
        #     y0 = (y1 + y2) / 2.0 - (x2 - x1) * (r * r - d)**0.5 / (d * 4) ** 0.5
        #     res = max(res, sum((x - x0)**2 + (y - y0)**2 <= r * r + 0.00001 for x, y in A))
        # return res