fork download
copy 
  1. #!/usr/bin/env python3
  2. from collections import namedtuple
  3. from math import ceil, floor
  4.  
  5. Circle = namedtuple('Circle', 'x, y, r')
  6.  
  7.  
  8. def count_lattice_points(r):
  9.     """Count the number of lattice points inside a circle of radius *r*.
  10.  
  11.     https://e...content-available-to-author-only...a.org/wiki/Gauss_circle_problem
  12.     """
  13.     assert isinstance(r, int)
  14.     return 1 + 4 * r + 4 * sum(int((r**2 - i**2)**.5) for i in range(1, r + 1))
  15.  
  16.  
  17. def count_lattice_points_intersection(a, b):
  18.     """The number of lattice points within intersection of circle *a* and *b*.
  19.  
  20.     x, y, r are all integers
  21.     """
  22.     # d = (a.x - b.x)**2 + (a.y - b.y)**2  # distance squared
  23.     # if d > (a.r + b.r)**2:     # no intersection
  24.     #     return 0
  25.     # elif d == (a.r + b.r)**2:  # touch in one point
  26.     #     # (y - a.y)/ a.r == (b.y - a.y) / (a.r + b.r)
  27.     #     # (x - a.x)/ a.r == (b.x - a.x) / (a.r + b.r)
  28.     #     return (a.r * (b.y - a.y) % (a.r + b.r) == 0 and
  29.     #             a.r * (b.x - a.x) % (a.r + b.r) == 0)
  30.     # elif d <= (a.r - b.r)**2:  # one circle within another
  31.     #     return count_lattice_points(min(a.r, b.r))
  32.  
  33.     # # circumferences intersect in two points
  34.     # assert (a.r - b.r)**2 < d < (a.r + b.r)**2
  35.     if a.r > b.r:
  36.         a, b = b, a  # make *a* be the smaller circle
  37.     assert a.r <= b.r
  38.  
  39.     # scan from bottom to top
  40.     count = 0
  41.     for y in range(a.y - a.r, a.y + a.r + 1):
  42.         if b.y - b.r <= y <= b.y + b.r:  # intersection possible
  43.             # find intersection boundaries for given y
  44.             ax1, ax2 = x_coordinates_intesect(a, y)
  45.             bx1, bx2 = x_coordinates_intesect(b, y)
  46.             if min(ax2, bx2) >= max(ax1, bx1):  # intersect
  47.                 count += floor(min(ax2, bx2)) - ceil(max(ax1, bx1)) + 1
  48.     return count
  49.  
  50.  
  51. def count_lattice_points_intersection_bruteforce(a, b):
  52.     return sum(((a.x - x)**2 + (a.y - y)**2 <= a.r**2 and
  53.                 (b.x - x)**2 + (b.y - y)**2 <= b.r**2)
  54.                for y in range(a.y - a.r, a.y + a.r + 1)
  55.                for x in range(a.x - a.r, a.x + a.r + 1))
  56.  
  57.  
  58. def x_coordinates_intesect(c, y):
  59.     """Get x-coordinates of intersection of circle *c* and horizontal line *y*."""
  60.     # solve quadratic equation
  61.     # (x - c.x)**2 + (y - c.y)**2 == c.r**2
  62.     assert c.r**2 >= (y - c.y)**2
  63.     D = (c.r**2 - (y - c.y)**2)**.5
  64.     return c.x - D, c.x + D
  65.  
  66. print(count_lattice_points_intersection(Circle(0, 0, 5), Circle(2, 2, 3)))
  67.  
  68. for a, b in [
  69.         (Circle(0, 0, 5), Circle(2, 2, 3)),
  70.         (Circle(0, 0, 4), Circle(2, 2, 3)),
  71.         (Circle(0, 0, 6), Circle(2, 2, 3)),
  72. ]:
  73.     got = count_lattice_points_intersection(a, b)
  74.     expected = count_lattice_points_intersection_bruteforce(a, b)
  75.     assert got == expected, (a, b, got, expected)
  76.  
Success #stdin #stdout 0.02s 28384KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
26
https://ideone.com/WQU2tV
language:
Python 3 (python  3.9.5)
created:
6 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!