The supernova impostor SN 2010da located in the nearby galaxy NGC 300, later identified as a likely supergiant B[e] high-mass X-ray binary, was simultaneously observed by NuSTAR and XMM-Newton between 2016 December 16 and 20, over a total time span of 310 ks. We report the discovery of a strong periodic modulation in the X-ray flux with a pulse period of 31.6 s and a very rapid spin-up, and confirm therefore that the compact object is a neutron star. We find that the spin period is changing from 31.71 s to 31.54 s over that period, with a spin-up rate of -5.56 x 10-7 s s-1, likely the largest ever observed from an accreting neutron star. The spectrum is described by a power-law and a disk black-body model, leading to a 0.3-30 keV unabsorbed luminosity of 4.7 x 10^39 erg s-1. Applying our best-fit model successfully to the spectra of an XMM-Newton observation from 2010, suggests that the lower fluxes of NGC 300 ULX1 reported from observations around that time are caused by a large amount of absorption, while the intrinsic luminosity was similar as seen in 2016. A more constant luminosity level is also consistent with the long-term pulse period evolution approaching an equilibrium value asymptotically. We conclude that the source is another candidate for the new class of ultraluminous X-ray pulsars.