  The Cassini spacecraft's Grand Finale orbits provided a unique opportunity to
probe Saturn's gravity field and interior structure. Doppler measurements
yielded unexpectedly large values for the gravity harmonics J_6, J_8, and J_10
that cannot be matched with planetary interior models that assume uniform
rotation. Instead we present a suite of models that assume the planet's
interior rotates on cylinders, which allows us to match all the observed even
gravity harmonics. For every interior model, the gravity field is calculated
self-consistently with high precision using the Concentric Maclaurin Spheroid
(CMS) method. We present an acceleration technique for this method, which
drastically reduces the computational cost, allows us to efficiently optimize
model parameters, map out allowed parameter regions with Monte Carlo sampling,
and increases the precision of the calculated J_2n gravity harmonics to match
the error bars of the observations, which would be difficult without
acceleration. Based on our models, Saturn is predicted to have a dense central
core of 15-18 Earth masses and an additional 1.5-5 Earth masses of heavy
elements in the envelope. Finally, we vary the rotation period in the planet's
deep interior and determine the resulting oblateness, which we compare with the
value from radio occultation measurements by the Voyager spacecraft. We predict
a rotation period of 10:33:34 h +- 55s, which is in agreement with recent
estimates derived from ring seismology.
