[tel]

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.

[amock]

It seems like you'd be IO bound when scanning the array so using larger than machine precision integers shouldn't slow things down. Particularly since you know the number of bits you'll need beforehand.