Except that the sum and sum of squares could both be too large for machine precision. I'd suggest: Scan over the array accumulating two quantities. P <- P

Except that the sum and sum of squares could both be too large for machine precision. I'd suggest: Scan over the array accumulating two quantities. P <- P*arr[i]/i; % prod(arr)/N! Q <- Q + i - arr[i]; % (N/2)(N+1) - sum(arr); If A is the missing value and B is the extra value, after doing the whole array P = B/A and Q = A-B. Therefore, A = Q/(1-P) and B = AP. You could conceivably compute P in log terms as well to make it a bit more stable if needed.