gdal - Get data from a raster using two coordinate arrays with Python

I'm trying to find a function in Python, that given two arrays with geographic coordinates, returns a list or array of raster values in those positions. that is to say:

import numpy as np
"""geog coordinates on x axis """
xc = [878571.041, 878583.54109659, 878596.04119319]
""" geographical coordinates on Y axis """
Yc = [187189.14, 187176.64015077, 187164.14030153]

"""list to array"""
X = np.array (Xc)
Y = np.array (Yc)
"""raster values at those coordinates"""
z = function (X, Y) 
or
z = function (Xc, Yc)

The output I expect is an array with raster values in the x, y positions. Is there a function in GDAL, Python or any other library that allows me to do this? Or must I necessarily use a for cycle to obtain these values.