unexpected behaviour of rnorm(): summary

The original script was

f <- function(n){max(rnorm(n))}
plot(sapply(rep(5000,4000),f))   #[this takes my PC about 30 seconds]

Received wisdom is that "The Marsaglia-Multicarry and Kinderman-Ramage
options don't play nicely together in the extreme tails of the Normal
distribution" (PD) but I discovered last night that

f.uppertail <- function(n,upper=3){x <- rnorm(n);x[x>upper]}
plot(f.uppertail(1e7))

shows the effect much more directly...and viewed this way, the default
RNG would seem to be more seriously flawed.  The second script above
considers the upper (ie Z>3) tail of the Gaussian: and it has gaps.
This is borderline "extreme tails" because pnorm(3) is only about
0.998.

In contrast, Bug #1664 refers to the maximal value of a series of
rnorm() values which, as PBR points out, _is_ a rather severe test of
any RNG; the difference between f.uppertail() and #1664 would be that
#1664 simulates

dmaxnorm <- function(x,n){n* dnorm(x)*pnorm(x)^(n-1)}

which involves PDFs and CDFs; f.uppertail() is solely a rnorm() issue.
If nothing else, f.uppertail() is a simpler test for rnorm().