title: Getting high precision values from qnorm in the tail
body: <h2>The problem</h2>
  <p>I am looking for high precision values for the normal distribution in the tail <code>(1e-10 and 1 - 1e-10)</code>, as the R package that I am using sets any number which is out of this range to these values and then calls the <code>qnorm</code> and <code>qt</code> function.</p>
  <p>What I have noticed is that the <code>qnorm</code> implementation in R is not symmetric when looking at the tails. This is quite surprising to me, as it is well known that this distribution is symmetric, and I have seen implementations in other languages that are symmetric. I have checked the <code>qt</code> function and it is also not symmetric in the tails.</p>
  <p>Here are the results from the qnorm function:</p>
  <pre><code>x       qnorm(x)                qnorm(1-x)              qnorm(1-x) + qnorm(x)
  1e-2    -2.3263478740408408     2.3263478740408408      0.0 (i.e &lt; machine epsilon)
  1e-3    -3.0902323061678132     3.0902323061678132      0.0 (i.e &lt; machine epsilon)
  1e-4    -3.71901648545568       3.7190164854557084      2.8421709430404007e-14
  1e-5    -4.2648907939228256     4.2648907939238399      1.014299755297543e-12
  1e-10   -6.3613409024040557     6.3613408896974208      -1.2706634855419452e-08
  </code></pre>
  <p>It is quite clear that at a value of <code>x</code> close to 0 or 1, this function breaks down. Yes, in &quot;normal&quot; use this isn't a problem, but I am looking at fringe cases and multiplying small probabilities by very large values, in which case the error <code>(1e-08)</code> becomes a large value.</p>
  <p>Note: I have tried this with <code>1-x</code> and with entering the actual number <code>0.00001</code> and <code>0.99999</code> and the accuracy issue is still there.</p>
  <h2>The questions</h2>
  <p>Firstly, is this a known problem with the <code>qnorm</code> and <code>qt</code> implementations? I could not find anything in the documentation, the algorithm is supposed to be accurate 16 digits for p values from <code>10^-314</code>as described in the <a href="https://www.jstor.org/stable/2347330?seq=1#page_scan_tab_contents" rel="noreferrer">Algorithm AS 241</a> paper.</p>
  <p>Quote from R doc:</p>
  <blockquote>
  <p>Wichura, M. J. (1988) Algorithm AS 241: The percentage points of the normal distribution. Applied Statistics, 37, 477–484.</p>
  <p>which provides precise results up to about 16 digits.</p>
  </blockquote>
  <p>If the R code implements the 7 digit version, why does it claim 16 digits? Or is it &quot;accurate&quot; but the original algorithm is not symmetric and wrong?</p>
  <p>If R does implement both versions of <a href="https://www.jstor.org/stable/2347330?seq=1#page_scan_tab_contents" rel="noreferrer">Algorithm AS 241</a> can I turn the 16 digit version on?</p>
  <p>Or, is there a more accurate version of <code>qnorm</code> in R?
  Or, another solution to my problem where I need high precision in the tails for quantile functions.</p>
  <h3>R version</h3>
  <pre><code>&gt;version 
  platform       x86_64-w64-mingw32          
  arch           x86_64                      
  os             mingw32                     
  system         x86_64, mingw32             
  status                                     
  major          3                           
  minor          3.2                         
  year           2016                        
  month          10                          
  day            31                          
  svn rev        71607                       
  language       R                           
  version.string R version 3.3.2 (2016-10-31)
  nickname       Sincere Pumpkin Patch   
  </code></pre>