We consider quantum dynamics of the order parameter in the discrete pairing model (Richardson model) in thermodynamic equilibrium. The integrable Richardson Hamiltonian is represented as a direct sum of Hamiltonians acting in different Hilbert spaces of single-particle and paired/empty states. This allows us to factorize the full thermodynamic partition function into a combination of simple terms associated with real spins on singly occupied states and the partition function of the quantum XY model for Anderson pseudospins associated with the paired/empty states. Using coherent-state path integral, we calculate the effects of superconducting phase fluctuations exactly. The contribution of superconducting amplitude fluctuations to the partition function in the broken-symmetry phase is shown to follow from the Bogoliubov-de Gennes equations in imaginary time. These equations in turn allow several interesting mappings, e.g., they are shown to be in a one-to-one correspondence with the one-dimensional Schrödinger equation in supersymmetric quantum mechanics. However, the most practically useful approach to calculate functional determinants is found to be via an analytical continuation of the quantum order parameter to real time, Δ (τ→it), such that the problem maps onto that of a driven two-level system. The contribution of a particular dynamic order parameter, Δ (τ), to the partition function is shown to correspond to the sum of the Berry phase and dynamic phase accumulated by the pseudospin. We also examine a family of exact solutions for two-level-system dynamics on a class of elliptic functions and suggest a compact expression to estimate the functional determinants on such trajectories. The possibility of having quantum soliton solutions coexisting with classical BCS mean field is discussed.