  We perform Monte Carlo calculation of correlation functions in 4d N=4 super
Yang-Mills theory on R*S^3 in the planar limit. In order to circumvent the
well-known problem of lattice SUSY, we adopt the idea of a novel large-N
reduction, which reduces the calculation to that of corresponding correlation
functions in the plane-wave matrix model or the BMN matrix model. This model is
a 1d gauge theory with 16 supersymmetries, which can be simulated in a manner
similar to the recent studies of the D0-brane system. We study two-point and
three-point functions of chiral primary operators at various coupling constant,
and find that they agree with the free theory results up to overall constant
factors. The ratio of the overall factors for two-point and three-point
functions agrees with the prediction of the AdS/CFT correspondence.
