How to Efficiently Apply 3D Transformations to a Large Number of Points in Python?

I am working on a 3D graphics project in Python, and I have a list of thousands of 3D points that represent a complex object. I need to perform multiple transformations (translation, rotation, and scaling) on these points efficiently.

I have tried using loops to apply the transformations, but the performance is not satisfactory with a large number of points. Here's what I've done so far (A snippet from my recent project)

import math

def transform_points(points, dx, dy, dz, rx, ry, rz, sx, sy, sz):
    for point in points:
        # Translate
        point[0] += dx
        point[1] += dy
        point[2] += dz

        # Rotate around X-axis
        y = point[1] * math.cos(rx) - point[2] * math.sin(rx)
        z = point[1] * math.sin(rx) + point[2] * math.cos(rx)
        point[1] = y
        point[2] = z

        # Rotate around Y-axis
        x = point[0] * math.cos(ry) + point[2] * math.sin(ry)
        z = -point[0] * math.sin(ry) + point[2] * math.cos(ry)
        point[0] = x
        point[2] = z

        # Rotate around Z-axis
        x = point[0] * math.cos(rz) - point[1] * math.sin(rz)
        y = point[0] * math.sin(rz) + point[1] * math.cos(rz)
        point[0] = x
        point[1] = y

        # Scale
        point[0] *= sx
        point[1] *= sy
        point[2] *= sz

Is there a more efficient way to perform these transformations, possibly using a library or built-in functionality? The speed is crucial, as these transformations will be part of real-time rendering.

One of the main problem here is you are computing the cos and the sin twice in each iteration of the loop for each value.

One thing you could do is store them in a var, which would avoid to compute it every time, knowing that cos and sin are pretty expensive.

For example :

def transform_points(points, dx, dy, dz, rx, ry, rz, sx, sy, sz):
    cos_rx = math.cos(rx)
    sin_rx = math.sin(rx)

    cos_ry = math.cos(ry)
    sin_ry = math.sin(ry)

    cos_rz = math.cos(rz)
    sin_rz = math.sin(rz)

    for point in points:
        # Translate
        point[0] += dx
        point[1] += dy
        point[2] += dz

        # Rotate around X-axis
        y = point[1] * cos_rx - point[2] * sin_rx
        z = point[1] * sin_rx + point[2] * cos_rx
        point[1] = y
        point[2] = z

        # Rotate around Y-axis
        x = point[0] * cos_ry + point[2] * sin_ry
        z = -point[0] * sin_y + point[2] * cos_ry
        point[0] = x
        point[2] = z

        # Rotate around Z-axis
        x = point[0] * cos_rz - point[1] * sin_rz
        y = point[0] * sin_rz + point[1] * cos_rz
        point[0] = x
        point[1] = y

        # Scale
        point[0] *= sx
        point[1] *= sy
        point[2] *= sz

In this way, you would have to compute once the sin and cos for rx, ry and rz, so only 6 computations, instead of twice for each iteration of the loop for each value, so points.length * 6 * 2 ,which chan be big if the number of points increases.