Quaternions are popular because they are usually the fastest technique but are notoriously non-intuitive. This is the best guide I have found explaining them in a practical way: https://danceswithcode.net/engineeringnotes/quaternions/quat...
I used that guide to write quaternion functions for one of my projects. It is three times faster than the (also gimbal safe) trigonometric method I used previously. The quaterion functions can be found at the top of: https://github.com/kevinloch/bsrender/blob/main/src/process-.... The rotation quaternion is constructed at the bottom of https://github.com/kevinloch/bsrender/blob/main/src/init-sta.... The vector quaternion is constructed by simply copying the Euclidian x,y,z into i,j,k while setting r=0, and the output x,y,z is copied from i,j,k after rotation.
Beware my program uses a non-standard coordinate system where +y is to the left, being derived from equatorial coordinates used in astronomy where the x->y angle (Right Ascension) increases to the East. Removing the y-axis inversion/reversion steps in the quaternion functions should make them work with conventional coordinates.
The speed advantage of quaternions depends on the use case. If the points being rotated are natively Euclidian and you are performing the same rotation on a billion points then it's a clear winner as the expensive trig functions are only done once while the actual rotation steps are all + and *.