[J3] Questions about sqrt -- mostly for vendors

[Steven G. Kargl]

On Tue, Jan 29, 2019 at 02:01:50PM -0800, Van Snyder via J3 wrote:
> Many applications need the reciprocal square root.  For example,
> celestial mechanics, electrodynamics, coulomb calculations, molecular
> dynamics....
> 
> I have been told that many processors offer a reciprocal square root
> instruction, because it's faster to compute it than to compute the
> square root.
> 
> Rolf Strebel's PhD thesis points out that there are ten steps in the
> Goldschmidt algorithm, but there are bubbles in the pipeline.  He shows
> how to compute a vector of 1/sqrt(x) in such a way as to fill the
> pipeline.
> 
> Do Fortran processors exploit instructions to compute 1/sqrt(x) when
> sqrt(x) appears in the denominator of a fraction?
> 
> Does elemental square root arrange interleaved calculations to fill
> bubbles in the pipeline?
> 
> Would the standard benefit from an RSQRT intrinsic function, or can most
> current optimizers use the RSQRT instruction directly, instead of
> multiplying by x to get sqrt(x), and then dividing by sqrt(x)?

There is a old bug report for gfortran/gcc from 2007 that addresses
this very issue.  Some of the analysis there might be informative

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=31723

I suspect other compiler vendors has considered this issue.

For at least x86 class hardware, a rsqrt instruction is only available
for single precision.

gfortran will normally use sqrt(x) and a division for y = 1 / sqrt(x).
If, however, -ffast-math or -funsafe-math-optimization is specified
on the command line, gfortran will use rsqrt(x) and a Newton-Raphson
step to ensure precision.  This applies to REAL x.  For DOUBLE PRECISION
x, there is no rsqrt so gfortran uses sqrt and does the division.

If you're looking to add new intrinsic subprograms to Fortran, I
think it would be more profitable to add SINPI, COSPI, TANPI,
and SINCOS.  Each has an obvious benefit to Fortran programmers.

-- 
Steve