pnorm(z, log.p = TRUE) #R

Z値が大きすぎるとP値を計算するときに桁落ちする。
これを扱うのにRmpfrというパッケージが必要だ（と思い込んでいた）。
pnormだけでいけた。

pnorm(10000, lower.tail=F, log.p=T)
[1] -500007.8

log10 P-valueにする場合、以下の通り

log_pvalue   = log(2) + pnorm(10000, lower.tail=F, log.p=T) 
log10_pvalue = log_pvalue/log(10) 
mintissa     = 10^(log10_pvalue %% 1)
exponent     = log10_pvalue %/% 1
sprintf("%2fe%i", mintissa, exponent)
[1] "4.6e-217151"

参考
https://stackoverflow.com/questions/26983612/using-rmpfr-to-round-with-precision-in-r