import numpy as np
# Gordon-Webb-Wolpert pair, each a union of 7 right-isosceles half-squares (legs = 1)
D1 = np.array([(0,0),(0,1),(2,3),(2,2),(3,2),(2,1),(1,1),(1,0)],float)
D2 = np.array([(1,0),(0,1),(0,2),(2,2),(2,3),(3,2),(2,1),(1,1)],float)
def area(P):
    x,y=P[:,0],P[:,1]; return 0.5*(np.dot(x,np.roll(y,-1))-np.dot(y,np.roll(x,-1)))
def inside(P,pts):
    # even-odd ray casting
    x,y=pts[:,0],pts[:,1]; c=np.zeros(len(pts),bool); n=len(P)
    for i in range(n):
        x1,y1=P[i]; x2,y2=P[(i+1)%n]
        m=((y1>y)!=(y2>y))
        with np.errstate(divide='ignore',invalid='ignore'):
            xi=x1+(y-y1)*(x2-x1)/(y2-y1)
        c^= m & (x<xi)
    return c
def on_boundary(P,pts,tol=1e-9):
    b=np.zeros(len(pts),bool); n=len(P)
    for i in range(n):
        a=P[i]; d=P[(i+1)%n]-a; L2=d@d
        t=np.clip(((pts-a)@d)/L2,0,1); q=a+t[:,None]*d
        b|= np.linalg.norm(pts-q,axis=1)<tol
    return b