We report on accurate, long-term numerical simulations of the orbits of the major planets in our solar system. The equations of motion are directly integrated by a Störmer multistep scheme, which is optimized to reduce round-off errors. The physical models are successively refined to include corrections due to general relativity and the finite size of the lunar orbit. In one case, the Earth-Moon system is resolved as two separate bodies, and the results are compared with those based on analytically averaging the lunar orbit. Through this comparison, a better analytical model is obtained. The computed orbits are in good agreement with those of previous studies for the past 5 Myr but not for earlier times. The inner planets exhibit chaotic behavior with a Lyapunov time of exponential separation of nearby orbits equal to about 4 Myr. Modeling uncertainties and chaos in the inner solar system restrict the accuracy of the computations beyond the past 50 Myr. We do not observe marked chaos in the motion of the Jovian planets in our 90 Myr integration, and we infer that the Lyapunov time for those planets is at least 30 Myr.