The authors describe a method of random generation of the integrals A1,2(t,t+h) = integral-t+h/t integral-s/t dw1(r)dw2(s) - integral-t+h/t integral-s/t dw2(r)dw1(s), together with the increments w1(t+h) - w1(t) and w2(t+h) - w2(t) of a two-dimensional Brownian path (w1(t), w2(t)). The method chosen is based on Marsaglia's ''rectangle-wedge-tail'' method, generalised to higher dimensions. The motivation is the need for a numerical scheme for simulation of strong solutions of general multidimensional stochastic differential equations with an order of convergence O(h), where h is the stepsize. Previously, no method has obtained an order of convergence better than O(square-root h) in the general case.