Magnetic fields play an important role in star formation in local as well as high-redshift galaxies, and recent studies of dynamo amplification in the first dark matter halos suggest that they were also likely present during the formation of the first stars in the Universe at redshifts 15 and above. We systematically study for the first time how magnetic fields can impact the primordial initial mass function. We carry out 200 3D magnetohydrodynamic simulations using the adaptive mesh refinement code FLASH with different initial turbulent magnetic field strengths as inspired from the theory of small-scale dynamo, producing a suite of more than 1100 first stars in total. We detect a strong statistical signature of suppressed fragmentation in the presence of strong magnetic fields, leading to a dramatic reduction in the number of primordial stars with masses low enough that they might be expected to survive to the present day. However, individual simulations are highly chaotic: around one-third of the runs result in the formation of a single massive star regardless of magnetic field strength, whereas there is a high diversity of clustered systems in others. While the origin of primordial magnetic fields is still not clear, we further show that they are crucial to derive a realistic initial mass function for the first stars. In fact, the role of primordial magnetic fields in setting the Population III initial mass function seems more important than that of their counterpart for the present-day initial mass function.