I'm not sure if this may be helpful, but for the particular case, using the attached script the result is correct.
The idea is to write the correlation as

exp (log (corr))

and evaluate log (corr) using logarithm properties.

(file #47601)

This has a very long history in numerical mathematics.
Simple application of the formula for cor(x,y) using
finite representation of x and y will result in cases
where cor(x,y) evaluates to a result outside [-1,1].

Therefore, the cor(x,y) should always return a value
on [-1,1], but it should also indicate a loss of
precision in any case where the most precisely evaluated
result is outside the bound.

A large literature on this is available.

corr(x,y) - 1 = 2.2204e-16

numerical precision limit is hit in these lines:

    c = cov (x);
    s = sqrt (diag (c));
    retval = c ./ (s * s');

results could maybe be constrained to [-1,1]

    retval = min(max(c ./ (s * s'), -1), 1);

However the result could still be 1-eps. The user of corr needs to anticipate that the result is subject to numerical errors.