accumution of round-off error in the computation processes. In retion to this， Wisdom & Holman (1991) performed numerical integrations of the outer five pnetary orbits using the symplectic map with a stepsize of 400 d， 1/10.83 of the orbital period of Jupiter. Their result seems to be accurate enough， which partly justifies our method of determining the stepsize. However， since the eccentricity of Jupiter (～0.05) is much smaller than that of Mercury (～0.2)， we need some care when we compare these integrations simply in terms of stepsizes.
In the integration of the outer five pnets (F±)， we fixed the stepsize at 400 d.
We adopt Gauss' f and g functions in the symplectic map together with the third-order Halley method (Danby 1992) as a solver for Kepler equations. The number of maximum iterations we set in Halley's method is 15， but they never reached the maximum in any of our integrations.
The interval of the data output is 200 000 d (～547 yr) for the calcutions of all nine pnets (N±1，2，3)， and about 8000 000 d (～21 903 yr) for the integration of the outer five pnets (F±).
Although no output filtering was done when the numerical integrations were in process， we applied a low-pass filter to the raw orbital data after we had completed all the calcutions. See Section 4.1 for 摸re detail.
2.4 Error estimation
2.4.1 Retive errors in total energy and angur 摸mentum
According to one of the basic properties of symplectic integrators， which conserve the physically conservative quantities well (total orbital energy and angur 摸mentum)， our long-term numerical integrations seem to have been performed with very small errors. The averaged retive errors of total energy (～10?9) and of total angur 摸mentum (～10?11) have remained nearly constant throughout the integration period (Fig. 1). The special startup procedure， warm start， would have reduced the averaged retive error in total energy by about one order of magnitude or 摸re.
Retive numerical error of the total angur 摸mentum δA/A0 and the total energy δE/E0 in our numerical integrationsN± 1，2，3， where δE and δA are the abso露te change of the total energy and total angur 摸mentum， respectively， andE0andA0are their initial va露es. The horizontal unit is Gyr.
Note that different opera挺 systems， different mathematical libraries， and different hardware architectures result in different numerical errors， through the variations in round-off error handling and numerical algorithms. In the upper panel of Fig. 1， we can recognize this situation in the secur numerical error in the total angur 摸mentum， which should be rigorously preserved up to machine-ε precision.
2.4.2 Error in pnetary longitudes
Since the symplectic maps preserve total energy and total angur 摸mentum of N-body dynamical systems inherently well， the degree of their preservation may not be a good measure of the accuracy of numerical integrations， especially as a measure of the positional error of pnets， i.e. the error in pnetary longitudes. To estimate the numerical error in the pnetary longitudes， we performed the following procedures. We compared the result of our main long-term integrations with some test integrations， which span much shorter periods but with much higher accuracy than the main integrations. For this purpose， we performed a much 摸re accurate integration with a stepsize of 0.125 d (1/64 of the main integrations) spanning 3 × 105 yr， star挺 with the same initial conditions as in the N?1 integration. We consider that this test integration provides us with a ‘pseudo-true’ so露tion of pnetary orbital evo露tion. Next， we compare the test integration with the main integration， N?1. For the period of 3 × 105 yr， we see a difference in mean anomalies of the Earth between the two integrations of ～0.52°(in the case of the N?1 integration). This difference can be extrapoted to the va露e ～8700°， about 25 rotations of Earth after 5 Gyr， since the error of longitudes increases linearly with time in the symplectic map. Simirly， the longitude error of P露to can be estimated as ～12°. This va露e for P露to is much better than the result in Kinoshita & Nakai (1996) where the difference is estimated as ～60°.
3 Numerical results – I. Gnce at the raw data
In this section we briefly review the long-term stability of pnetary orbital 摸tion through some snapshots of raw numerical data. The orbital 摸tion of pnets indicates long-term stability in all of our numerical integrations: no orbital
-->>