What is the format of the data from the JPL's HORIZONS system? I'm looking for accurate positions and velocities of the planets in the solar system over several decades. I want to simulate their trajectories using Newtonian laws of motion, and compare the precession of Mercury's perihelion with observed data. I used the JPL's HORIZONS system but I can't manage to find these positions and velocities in the output data. 

Could you explain how the data is formatted? (The HORIZONS web-interface Tutorial didn't help.)
WebGeocalc (http://wgc.jpl.nasa.gov:8080/webgeocalc/) is actually a 
fairly useful tool, and the next best thing to downloading the CSPICE 
libraries, so let me provide a slightly more detailed answer to your 
question. 

When you first visit WebGeocalc, you will see something like this: 

 

Because you are looking for perihelions, which involve distance, 
scroll down and choose "distance finder": 

 

Since you want to find the perihelion of Mercury, fill out the top 
half of the form like this: 

 

For this example, I'll find all perihelions from 1900 to 2201. A 
perihelion is a minimum of distance, and since you're looking for all 
perihelions, choose "is local minimum" on the lower half of the form, 
so it looks as follows: 



Scroll down once again (you can leave the remaining form values as is) 
and click the "Calculate" button: 



When you get your results, click "Save All Intervals" to save the 
times of Mercury's perihelions. Note that each result is an interval, 
but the interval is 0 seconds long, so it's actually a point in time: 

 

Now click on "Calculation Menu" to get back to the main menu: 



We'll now find Mercury's position at all of these perihelions so you 
can see how the perihelion position changes. 

To find Mercury's position at these times, first click on "State 
Vector", which will take you to a form. Fill out the top half of the 
form as follows: 



Fill out the lower half of the form as follows, dragging the results 
window (for the times of Mercury's perihelions) into the input times 
window: 



Click "Calculate" to receive your results: 



See also http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/Tutorials/pdf/individual_docs/36_webgeocalc.pdf

NOTE: if you're going to do extensive calculations, you'll probably 
want to download the C SPICE library and kernels: 
http://naif.jpl.nasa.gov/naif/tutorials.html

EDIT: I emailed Jon D. Giorgini (jdg@tycho.jpl.nasa.gov) re the open-sourceness of HORIZONS, and got this answer (which I'm not sure is helpful, but...):



Horizons output is copyrighted in the sense one can't take a
table output by Horizons and re-publish it as one's own.
JPL is a contractor organization and so the output is not prepared
by an officer or employee of the US government, and it is not
transferred to the US government, so copyright remains relevant.

The data values themselves, say the RA/DEC of Pluto at
some instant, are not legally copyrighted, though it would 
be unethical to publish it without citation.

Legal copyright protection is for fixed representations
analogous to a book; the book has copyright protection, but 
an individual word extracted from the book (i.e., into a 
different representation) is not itself protected.  Up to 
a point

Note that Horizons reads SPICE files for planets, natural
satellites, and spacecraft, then derives results it outputs 
from scratch from that basis. But for comets and asteroids, 
it creates all information on demand from scratch.

Further, we (JPL Solar System Dynamics) put the planet/satellite/
spacecraft information into the SPICE format to begin with. So 
a SPICE file might be like the blank pages of a book; not 
copyrighted until the original work is written into it, but not
itself copyrighted (beyond the source code, of course).

But if you need a more trustworthy legal opinion for a specific 
situation, I can put you in touch with legal staff. If this is an 
actual issue and not just a point of curiosity, that's the way to
since perhaps other rules would apply to your actual situation.