We present analysis of the light curves (LCs) of 77 hydrogen-poor superluminous supernovae (SLSNe I) discovered during the Zwicky Transient Facility Phase I operation. We find that the majority (67%) of the sample can be fit equally well by both magnetar and ejecta-circumstellar medium (CSM) interaction plus $^56$Ni decay models. This implies that LCs alone cannot unambiguously constrain the physical power sources for an SLSN I. However, 23% of the sample show inverted V-shape, steep-declining LCs or features of long rise and fast post-peak decay, which are better described by the CSM+Ni model. The remaining 10% of the sample favors the magnetar model. Moreover, our analysis shows that the LC undulations are quite common, with a fraction of 18%-44% in our gold sample. Among those strongly undulating events, about 62% of them are found to be CSM-favored, implying that the undulations tend to occur in the CSM-favored events. Undulations show a wide range in energy and duration, with median values (and 1σ errors) being as $1.7 % -0.7 % ^+1.5 % ,E_mathrmrad,mathrmtotal$ and $28.8-9.1^+14.4$ days, respectively. Our analysis of the undulation timescales suggests that intrinsic temporal variations of the central engine can explain half of the undulating events, while CSM interaction (CSI) can account for the majority of the sample. Finally, all of the well-observed He-rich SLSNe Ib either have strongly undulating LCs or the LCs are much better fit by the CSM+Ni model. These observations imply that their progenitor stars have not had enough time to lose all of the He-envelopes before supernova explosions, and H-poor CSM are likely to present in these events.