Equations for Milankovitch Factors?
There are many models that estimate orbital parameters at long timescales. All of them rely in rather complicated formulas that account for the perturbations generated by other planets, the Moon, and many other factors. Some of the models even consider relativistic effects as described in this answer.

An orbital solution that I often see cited and used in papers dealing with climatic forcing is LA93, computed by Laskar et al. (1993) and described in the paper Orbital, precessional, and insolation quantities for the Earth from -20 MYR to +10 MYR. However, more accurate solutions are available by now, like LA2004 and more notably LA2010. I'll describe below how to download output data and codes for each of these solutions.

For LA93, the codes (in fortran) and outputs of that model are available at the VizieR catalogue, where you can query the model output, or just go into this FTP folder and download the files for the code and long tables with the model outputs every one thousand years.

If by "the three Milankovitch factors" you mean Obliquity, Eccentricity, and   Precession index (Eccentricity $\times$ sin(longitude of perihelion)). That data is ready to download and use in these two files:

CLIVAR0N.ASC (-20Myr to 0Myr)
CLIVAR0P.ASC (0Myr to +10Myr)

The byte-per-byte description of these files as described in the README is

--------------------------------------------------------------------------------
   Bytes Format  Units   Label     Explanations
--------------------------------------------------------------------------------
   4- 11  F8.0   1000yr  t         Time from J2000 in 1000 years
  13- 31  D19.12 ---     e         Eccentricity
  33- 51  D19.12 rad     eps       obliquity (radians)
  53- 71  D19.12 ---     CP        e * sin(longitude of perihelion
                                           from moving equinox)
--------------------------------------------------------------------------------


With the code you can also compute common climatic proxies such as isolation at 65°N or any other latitude.

For LA2004 (-50 Myr to +20 Myr) the same research team implemented a very useful online calculator that allows to compute insolation quantities derived from the orbital and precessional quantities for any latitude of interest. It outputs the results in tabular form. The calculator interface looks like this so you get an idea of the available options.



The LA2010 solution is described in the paper La2010: a new orbital solution for the long-term motion of the Earth, and goes 250 Myr back in time. You can find similar data files and codes for this solution at the Astronomical Solutions for Earth Paleoclimates page. To my understanding, LA2010 represent the state-of-the-art solution for orbital parameters and it does include relativistic corrections as well.