  We devise an optimal method to measure the temporal power spectrum of the
lensing and intrinsic fluctuations of multiply-imaged strongly lensed quasar
light curves, along with the associated time delays. The method is based on a
Monte-Carlo Markov Chain (MCMC) sampling of a putative gaussian likelihood, and
accurately recovers the input properties of simulated light curves, as well as
the "Time Delay Challenge". We apply this method to constrain the dimensionless
cosmological (non-linear) matter power spectrum on milliparsec scales
(comparable to the size of the solar system), to $\Delta_{\rm NL}^2< 4 \times
10^7$ at $k_{\rm NL} \sim 10^3 {\rm pc}^{-1}$. Using a semi-analytic nonlinear
clustering model which is calibrated to simulations, the corresponding
constraint on the primordial (linear) scalar power spectrum is ${\cal P}_{\cal
R} < 3 \times 10^{-9}$ at $k_{\rm L}\sim$ 3 pc$^{-1}$. This is the strongest
constraint on primordial power spectrum at these scales, and is within an order
of magnitude from the standard $\Lambda$CDM prediction. We also report
measurements of temporal spectra for intrinsic variabilities of quasar light
curves, which can be used to constrain the size of the emitting region in
accretion disks. Future cadenced optical imaging surveys, such as LSST, should
increase the number of observed strongly lensed quasars by 3 orders of
magnitude and significantly improve these measurements, even though
improvements in modelling quasar accretion and stellar microlensing are
necessary.
