Leon's shooting range problem
Area-based algorithm, \$O(k^2+n\log(k))\$
Our general strategy is to instead select points (outside the sticker rectangles) from the square (-1,-1) to (1,1), then repeat if we don't get a point within the circle. On average, we will have to select \$\frac{4}{\pi}\$ points from the square for each point in the circle.
The runtime is not depend on rectangle coverage, but it does depend on measures of rectangle intersection (denoted L in the complexity analysis in the code). For the listed complexity, I use \$L=O(k)\$, which is probably not tight. If \$L=O(\log(k))\$ can be proven, then the complexity is \$O((k+n)\log(k))\$ (see get_heights for an algorithmic alternative).
Define &quot;open area&quot; to be the area that points could be generated in. Recall that we are using a square instead of a circle, so the total area is piece-wise linear in x.

Precompute the total open height function \$h(x)\$ (this is piece-wise constant, so store the change at each transition: either the start or end of a new rectangle). Similarly, pre-compute the cumulative areas \$\int_{-1}^{x} h(x)\$, which is piece-wise linear, so once again store the value at each transition. This can be done in \$O(k \log(k))\$ time

Precompute the intervals on which y coordinates can be distributed for a given x (again, only storing this interval on each transition). This might be able to be done in \$O(k \log(k))\$, but this algorithm only succeeds at \$O(k^2)\$ time

Generate points. For each point: (average \$\frac{4}{\pi} n\$ times for a total runtime of \$O(n\log(k))\$
a. Choose a random value \$r\$ for the area (\$O(1)\$)
b. Find \$x\$ such that \$\int_{-1}^{x} h(t) dt = r\$. This is accomplished through binary search on the pre-computed heights \$h\$ and cumulative areas \$\int h\$, so this is done in \$O(\log(k))\$ time
c. Choose a random value \$y\$ for the given x such that (x,y) is in an open area. Through similar binary search as 3b, this is done in \$O(\log(k))\$ time
d. Reject the point if it is outside the circle (constant factor 4/pi)


The total runtime is \$O(k\log(k)) + O(k^2) + O(n\log(k)) = O(k^2+n\log(k))\$
In the following example images, I used the same n=10000, so more rectangle coverage meant a higher density of points.



Full python code below, including matplotlib display:
import math
import random
# for binary search
import bisect
# matplotlib for plotting the output points
import matplotlib.pyplot as plt
import matplotlib

def find_le(ls, x):
    # find_le is O(log(|ls|))
    # Find rightmost value less than or equal to x in ls
    i = bisect.bisect_right(ls, x)
    return i-1, ls[i-1]

class Interval:
    def __init__(self, lo, hi):
        # each range is tuple (lo, hi)
        # ranges shall be sorted in ascending order
        # each range does not overlap with other ranges
        self.ranges = [(lo, hi)]

    def remove(self, lo, hi):
        new_ranges = []
        for self_range in self.ranges:
            if lo &lt;= self_range[0] &lt;= hi &lt; self_range[1]:
                new_ranges.append((hi, self_range[1]))
            elif self_range[0] &lt; lo &lt;= self_range[1] &lt;= hi:
                new_ranges.append((self_range[0], lo))
            elif self_range[0] &lt; lo &lt; hi &lt; self_range[1]:
                new_ranges.append((self_range[0], lo))
                new_ranges.append((hi, self_range[1]))
            elif lo &lt;= self_range[0] &lt; self_range[1] &lt;= hi:
                pass
            elif hi &lt;= self_range[0] or self_range[1] &lt;= lo:
                # no overlap
                new_ranges.append(self_range)
        self.ranges = new_ranges

    def cache_lengths(self):
        # If L is the number of subintervals of this Interval
        # then cache_lengths is O(L) ⊂ O(k)

        # in cum_lengths, (0.6, 0.7) means that the total length to the
        # left of x=0.7 in this interval is 0.6
        self.cum_lengths = [(0, -1)]
        last_length = 0
        for self_range in self.ranges:
            length = self_range[1] - self_range[0]
            self.cum_lengths.append((last_length, self_range[0]))
            last_length += length
        self.total_length = last_length

    def random_point(self):
        # If L is the number of subintervals of this Interval
        # then random_point is O(log(L)) ⊂ O(log(k))
        r = random.random()*self.total_length
        # This find_le is O(log(L))
        i, (length_left, x_left) = find_le(self.cum_lengths, (r,))
        return (r - length_left) + x_left

    def __repr__(self):
        return &quot;I&quot; + str(self.ranges)

def interval_complement(rects):
    # O(L) where L is the length of rects list
    interval = Interval(-1, 1)
    for rect in rects:
        interval.remove(rect.bottom, rect.top)
    return interval

class Rect:
    def __init__(self, left, right, bottom, top):
        self.left = min(left, right)
        self.right = max(left, right)
        self.bottom = min(bottom, top)
        self.top = max(bottom, top)
        self.width = self.right - self.left
        self.height = self.top - self.bottom

    def __repr__(self):
        return f&quot;Rect({self.left}, {self.right}, {self.bottom}, {self.top})&quot;

    def overlap(self, others):
        # returns vertical (y) overlap between self and other rects
        total_overlap = 0
        # https://pyinterval.readthedocs.io/en/latest/guide.html#operations is
        # looking so tempting right now, but unions/merges would
        # require Principle Inclusion-Exclusion which can be k^2 worst case
        starts = sorted(others, key=lambda rect: rect.bottom, reverse=True)
        ends = sorted(others, key=lambda rect: rect.top, reverse=True)
        active_rects = set()
        overlap_start = self.bottom
        while ends:
            if starts and starts[-1].bottom &lt; ends[-1].top:
                active_rect = starts.pop()
                if not active_rects and active_rect.bottom &gt; self.bottom:
                    # This is the start of an overlap interval
                    overlap_start = active_rect.bottom
                active_rects.add(active_rect)
            else:
                active_rect = ends.pop()
                active_rects.remove(active_rect)

                if not active_rects:
                    # This is the end of an overlap interval
                    # Rare cases of negative = started and ended below start of self
                    total_overlap += max(0, min(active_rect.top, self.top) - overlap_start)
        return total_overlap

    def as_matplotlib_rect(self):
        return matplotlib.patches.Rectangle((self.left, self.bottom), self.width, self.height)

def get_heights(rects):
    # get_heights is O(k log(k)) total
    # Sorting: O(k log(k))
    # Use these lists as a queue
    starts = sorted(rects, key=lambda rect: rect.left, reverse=True)
    ends = sorted(rects, key=lambda rect: rect.right, reverse=True)
    # list of tuples (open height starting at x, x, open intervals starting at x)
    heights = [(2, -1, Interval(-1,1))]
    active_rects = set()
    while ends:
        # inner loop performed once for k left edges and k right edges,
        # so this is O(k)
        if starts and starts[-1].left &lt; ends[-1].right:
            # next x-pos is the start of a new rect
            rect = starts.pop()
            height_diff = rect.height - rect.overlap(active_rects)
            height = heights[-1][0] - height_diff
            active_rects.add(rect)
            # interval_complement is O(L) so this is O(k^2), not sure if it is O(k log(k))
            heights.append((height, rect.left, interval_complement(active_rects)))
        else:
            # next x-pos is the end of an old rect
            rect = ends.pop()
            active_rects.remove(rect)
            # rect.overlap is O(L) as well
            height_diff = rect.height - rect.overlap(active_rects)
            # Two choices to get O(k log(k)):
               # 1. prove that |active_rects| is O(log(k)) under the given approaches
               # 2. avoid using interval_complement on each transition using some increments approach, and avoid doing rect.overlap on each transition
            height = heights[-1][0] + height_diff
            heights.append((height, rect.right, interval_complement(active_rects)))
    heights.append((2, 1, Interval(-1,1)))
    # cum_areas is the cumulative areas from -1 to x (of not rects)
    # Not subtracting curve because we will be sampling from square
    cum_areas = [(0, -1, Interval(-1, 1))]
    last_height = 2
    for height, x, interval in heights[1:]:
        last_area, last_x, _ = cum_areas[-1]
        new_area = (x-last_x)*last_height
        cum_areas.append((last_area + new_area, x, interval))
        last_height = height
    return (heights, cum_areas)

def random_pick_points(rects, num_points):
    # random_pick_points is O((n+k) log(k))
    # Step 1. Precompute heights: O(k log (k))
    heights, cum_areas = get_heights(rects)
    total_area = cum_areas[-1][0]
    points = []
    # Step 2. Precompute vertical intervals: O(k)
    for _, _, y_interval in cum_areas:
        y_interval.cache_lengths()
    # Step 3. Sample points outside rects and within circle
    # Total time complexity of this loop is O(n log(k))
    while len(points) &lt; num_points:
        # This loop is reached average (4/pi)*n = O(n) times
        # Each occurence has complexity O(log(k))
        # Step 3a. Choose a random area value
        area = random.random()*total_area
        # Step 3b. Choose the x value such that area is the total area
        # outside rects and within 2x2 square to the left of x
        # There are 2k+1 total cum_areas so this find_le is O(log(k))
        i, (area_left, x_left, y_interval) = find_le(cum_areas, (area,))
        height, _, _ = heights[i]
        x = x_left + (area - area_left) / height
        # Step 3c. Choose a y value uniformly in the open spaces of that x value
        # Interval.random_point() is O(log(k))
        y = y_interval.random_point()
        # Step 3d. Rejection sample until points are inside the circle
        # On average, takes constant 4/pi random points to land within the circle
        if x*x + y*y &lt; 1:
            points.append((x, y))
    # Done!
    return points

def display_points(rects, points):
    fig, ax = plt.subplots(1)
    ax.plot(*zip(*points), 'ko', markersize=1, alpha=0.5)
    ax.axis([-1, 1, -1, 1])
    ax.add_collection(
        matplotlib.collections.PatchCollection(
            [rect.as_matplotlib_rect() for rect in rects],
            facecolor = 'r', edgecolor='None', alpha=0.5
        )
    )
    plt.show()

def display_random_points(rects, num_points):
    points = random_pick_points(rects, num_points)
    display_points(rects, points)

def random_rects(k):
    return [Rect(*[random.random()*2-1 for _ in range(4)]) for i in range(k)]

k = 20
# rects = random_rects(k)
# print(&quot;Rects:&quot;, rects)
n = 10000
# print(&quot;Points:&quot;, random_pick_points(rects, n))
display_random_points(random_rects(k), n)

Try it online I guess, if you can read a list of points.