We modified J. G. Williams' formulation of the precession and the nutation by using the 3-1-3-1 rotation [1] so as to express them in an arbitrary inertial frame of reference. It gives the precession-nutation matrix as the product of four rotational matrices as NP = R₁(-ε)·R₃·(-ψ)·R₁·(φ)·R₃·(γ) and the precession one similarly as P = R₁(-ε)·R₃·(-ψ)·R₁·(φ)·R₃·(γ) . Here φ and &gamma are the angles to specify the location of the ecliptic pole of date in the given inertial frame, ψ and ψ are the true and mean ecliptic angles of precession, respectively, and ψ and ψ are the true and mean obliquities of the ecliptic, respectively. As a result, the pole coordinates of the true and mean equators are explicitly given in terms of the newly introduced precession angles. Although the expression of nutation matrix is unchanged, we recommend the usage of the above form of NP instead of preparing P and N separately, because of faster evaluation. The formulation is robust in the sense it avoids a singularity caused by finite pole offsets near the epoch. Facing the singularity is inevitable in the current IAU formulation. By using a recent theory of the forced nutation of the non-rigid Earth, SF2001 [2], we converted the true pole offsets referred to the ICRF, observed by VLBI for 1979-2000, and complied by USNO, to the offsets in the above three angles of precession, ψ,φ and γ , while we fixed ε as the combination of the linear part provided in SF2001 and the quadratic and higher terms derived by Williams (1994). From the converted offsets, we determined the best-fit polynomial expressions of the three precession angles in the ICRF by a weighted least square method where we kept the quadratic and cubic terms as the same as in Williams (1994). These constitute a new set of fundamental expressions of the precessional quantities. The combination of the new precession formula and the periodic part of SF2001 serves a good approximation of the precession-nutation matrix in the ICRF.