Constraining cosmology using weak gravitational lensing consists of comparing a measured feature vector of dimension N b with its simulated counterpart. An accurate estimate of the N b × N b feature covariance matrix C is essential to obtain accurate parameter confidence intervals. When C is measured from a set of simulations, an important question is how large this set should be. To answer this question, we construct different ensembles of Nr realizations of the shear field, using a common randomization procedure that recycles the outputs from a smaller number Ns ≤ Nr of independent ray-tracing N -body simulations. We study parameter confidence intervals as a function of (Ns, Nr) in the range 1 ≤ Ns ≤ 200 and 1 ≤ Nr 10 5 . Previous work [1] has shown that Gaussian noise in the feature vectors (from which the covariance is estimated) lead, at quadratic order, to an O(1/Nr) degradation of the parameter confidence intervals. Using a variety of lensing features measured in our simulations, including shear-shear power spectra and peak counts, we show that cubic and quartic covariance fluctuations lead to additional O(1/N 2 r ) error degradation that is not negligible when Nr is only a factor of few larger than N b . We study the large Nr limit, and find that a single, 240Mpc/h sized 512 3 -particle N -body simulation (Ns = 1) can be repeatedly recycled to produce as many as Nr = few × 10 4 shear maps whose power spectra and high-significance peak counts can be treated as statistically independent. As a result, a small number of simulations (Ns = 1 or 2) is sufficient to forecast parameter confidence intervals at percent accuracy.PACS numbers: 95.36.+x, 95.30.Sf, 98.62.Sb