Aims. We aim to measure the mass of the supermassive black hole (SMBH) in the S0 galaxy NGC 4751 using CO J:3-2 emission from the 5 ''-scale nuclear rotating molecular disk. Methods. We imaged the kpc-scale molecular gas disk in NGC 4751 at 0.'' 22 (similar to 28 pc) spatial and 28 km s-1 spectral resolution in the CO J:3-2 emission line and neighboring continuum, with the Atacama Large Millimeter Array (ALMA). We used Hubble Space Telescope (HST) imaging and stellar and ionized gas kinematics at 100 pc to kpc-scales, derived from integral field spectroscopy, to determine the galaxy morphology and the circular velocity attributed to the stellar potential. We used the Markov chain Monte Carlo (MCMC) algorithm in the KINematic Molecular Simulation (KinMS) package to obtain the model parameters that best fit the observed molecular gas kinematics in the ALMA datacube. Results. Strong CO emission was detected over radii of similar to 0.'' 2 to 5 '', with isolated CO clumps detected out to 7 ''. The molecular disk kinematics is axisymmetric and rotation-dominated, with radial velocities between 400 km s(-1) and 660 km s(-1), a kinematic major axis position angle (PA) of 355 degrees, and an inclination close to 78 degrees. The intrinsic velocity dispersion is similar to 16 km s(-1), and there is no evidence for significant non-rotational kinematics. The kinematic center of the disk coincides with the compact nuclear 345 GHz source. The SMBH sphere of influence (SOI) is well resolved along all position angles. The (rotation) velocity curve due to the stellar potential ( Vradialmax similar to 430 V radial max similar to 430 V (max)(radia) km s-1) is determined by fitting the luminosity profile of NGC 4751 in an (H-band) image from the Wide Field Camera 3 (WFC3) aboard HST, and constraining the mass-to-light ratio (M/L) at this waveband using the molecular- and ionized-gas kinematics at radii greater than or similar to 4 '', outside the SMBH SOI. Several KinMS fits, all using a distance (D) of 26.3 Mpc, but with variations in other input quantities, resulted in SMBH masses of 3.22 - 4.33 x 10(9) M-circle dot and M/L values of 1.1-2.3 in the F160W band. In each fit, the statistical errors of these values are on the level of a few percent. Conclusions. Based on the results of the multiple KinMS fits, we argued for and adopted a value of 3.3 x 109 (D/26.3) M circle dot for the black hole mass, along with a (constant with radius) M/L of 2.28/(26.3/D)2 in the F160W band. We estimated the (one sigma) errors to be 20% in each of these. We find that the primary driver of the uncertainty (apart from distance) is the stellar potential in this dusty S0 galaxy. This CO-based mass is similar to 2.4 times higher than a previous stellar-dynamics based SMBH mass measurement using the same distance. We argue that this new value is more robust given the clear and well resolved Keplerian-rotation dominated signature in the molecular disk, as well as its robust values of inclination (78 degrees) and PA (355 degrees), further supported by the consistency among derived values across different datasets and methods.