The best and most concise answer is here
The "gold standard" is the SPICE toolkit. You can read more about using spice in this excellent answer.
The JPL Development Ephemerides are lists of the Chebyshev coefficients you seek, along with some additional information. They can be downloaded as text files or binary. If you are familliar with Python you can have a look inside the Skyfield package (Github) to see how it's done there.
Here's a quick way to get DE421: https://pypi.org/project/de421/
Here they are in ascii format: ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii
Further reading:
according to ftp://ssd.jpl.nasa.gov/pub/eph/planets/ascii/ascii_format.txt:
INTERNAL FORMAT OF THE EPHEMERIS FILES
--------------------------------------
On the first record of an export binary file or in the "GROUP 1050 of the ascii
"header", there are 3 sets of 15 integers each. (Older versions have only the first 13
integers in each set)
The 15 triplets give information about the location, order and time-coverage of
the chebychev polynomials corresponding to the following 15 items:
Mercury
Venus
Earth-Moon barycenter
Mars
Jupiter
Saturn
Uranus
Neptune
Pluto
Moon (geocentric)
Sun
Earth Nutations in longitude and obliquity (IAU 1980 model)
Lunar mantle libration
Lunar mantle angular velocity
TT-TDB (at geocenter)
Word (1,i) is the starting location in each data record of the chebychev
coefficients belonging to the ith item. Word (2,i) is the number of chebychev
coefficients per component of the ith item, and Word (3,i) is the number of
complete sets of coefficients in each data record for the ith item.
Items not stored on the ascii files have 0 coefficients [Word (3,i)].
Data Records ("GROUP 1070")
These records contain the actual ephemeris data in the form of chebychev
polynomials.
The first two double precision words in each data record contain
Julian date of earliest data in record.
Julian date of latest data in record.
The remaining data are chebychev position coefficients for each component of
each body on the tape. The chebychev coefficients for the planets represent
the solar system barycentric positions of the centers of the planetary systems.
There are three Cartesian components (x, y, z), for each of the items #1-11;
there are two components for the 12th item, nutations : d(psi) and d(epsilon);
there are three components for the 13th item, librations : phi, theta, psi;
there are three components for the 14th item, mantle omega_x,omega_y,omega_z;
there is one component for the 17th item, TT-TDB.
Planetary positions are stored in units of kilometers (TDB-compatible).
The nutations and librations are stored in units of radians.
The mantle angular velocities are stored in radians/day.
TT-TDB is stored in seconds.