Celestial Programming : Precession IAU 1976

$$ \begin{align*} t &=\frac{JD-JD0}{36525} \\ \\ T &= \frac{JD0 - 2451545.0}{36525} \\ \\ \zeta &= (2306.2181 + 1.39656T - 0.000139T^2)t + (0.30188 - 0.000344T)t^2 + 0.017998t^3 \\ z &= (2306.2181 + 1.39656T - 0.000139T^2)t + (1.09468 + 0.000066T)t^2 + 0.018203t^3 \\ \theta &= (2004.3109 -0.85330T - 0.000217T^2)t + (-0.42665 - 0.000217T)t^2 - 0.041833t^3 \\ \\ \tan(\alpha - z) &= \frac{\cos(\delta_0) \sin(\alpha_0 + z)}{\cos(\theta) \cos(\delta_0) \cos(\alpha_0 + z) - \sin(\theta) \sin(\delta_0)} \\ \\ \sin(\delta) &= \sin(\theta) \cos(\delta_0) \cos(\alpha_0 + z) + \cos(\theta) \sin(\delta_0) \\ \end{align*} $$ JD - JD to Precess to.
JD0 - Starting epoch (2451545 for J2000)
\(\alpha_0\) - Starting RA
\(\delta_0\) - Starting Dec
\(\alpha\) - Ending RA
\(\delta\) - Ending RA

Alternate matrix method

Using rotation matrixes $$ \vec r=R_3(-z)R_2(-\theta)R_3(-\zeta)\vec r_0 $$ or explaning out the rotations into a single matrix:
$$ \bf P = \begin{bmatrix} \cos z \cos \theta \cos \zeta - \sin z \sin \zeta & -\cos z \cos \theta \sin \zeta - \sin z \cos \zeta & -\cos z \sin \theta \\ \sin z \cos \theta \cos \zeta + \cos z \sin \zeta & -\sin z \cos \theta \sin \zeta + \cos z \cos \zeta & -\sin z \sin \theta \\ \sin \theta \cos \zeta & -\sin \theta \sin \zeta & \cos \theta \end{bmatrix} $$

Simplified, if starting from J2000

$$ \begin{align*} \zeta &=2306.2181"t + 0.30188"t^2 + 0.017998"t^3 \\ z &=2306.2181"t + 1.09468"t^2 + 0.018203"t^3 \\ \theta &=2004.3109"t - 0.42665"t^2 - 0.041833"t^3 \\ \end{align*} $$