Using 1200 CPUs of the National Supercomputer TH-A1 and a parallel integral algorithm based on the 3500th-order Taylor expansion and the 4180-digit multiple precision data, we have done a reliable simulation of chaotic solution of Lorenz equation in a rather long interval 0<=t<=10000 LTU (Lorenz time unit). Such a kind of mathematically reliable chaotic simulation has never been reported. It provides us a numerical benchmark for mathematically reliable long-term prediction of chaos. Besides, it also proposes a safe method for mathematically reliable simulations of chaos in a finite but long enough interval. In addition, our very fine simulations suggest that such a kind of mathematically reliable long-term prediction of chaotic solution might have no physical meanings, because the inherent physical micro-level uncertainty due to thermal fluctuation might quickly transfer into macroscopic uncertainty so that trajectories for a long enough time would be essentially uncertain in physics.