[RFC] Quad-float support, round 4

[Steve Kargl]

On Tue, Sep 14, 2010 at 07:16:53AM -0700, Tim Prince wrote:
>  On 9/13/2010 6:51 PM, Steve Kargl wrote:
> >On Mon, Sep 13, 2010 at 08:17:56PM +0200, FX wrote:
> >>(Yeah, Intel compiler doesn't like real(kind=10)...)
> >>
> >>Again, comparing __float128 to double:
> >>SQRT is slower by 160, SIN is slower by only 6, COS by 10,
> >>ASINH by merely 10 and ERFC by 6 again. I'm actually amazed
> >>at how it does that!
> >>
> >First, I agree that worrying about performance at this point
> >should not be a priority.  Perhaps, someone like Tim Prince
> >might step forward to lend a hand in optimizing the code.
> >
>
> I haven't figured out yet where you did this in the source code, nor 
> what ground rules you are applying.
> If you have decided to use an iterative improvement for sqrt(), the 
> native x87 sqrt in 64-bit precision mode followed by a single Heron's 
> step should be sufficient.  I don't think people would normally change 
> precision mode when running under cygwin, but there's no guarantee.
> 

It is done in libquad/src/math/sqrtq.c.  FX uses 1 or 2 Heron
iterations depending on the availability of sqrt() and sqrtl(),
and depending on the magnitude of the argument he may use
frexpq() and scalbnq() for argument reduction.  sqrtq() 
isn't the interesting case because FX reported the 
__float128 to double comparisons as follows the ratios:

        gfc   ifc
sqrtq   200   160   (ie, __float128 is XXX times slower than double)
sinq    100     6
cosq    150    10
asinhq  120    10
erfcq   150     6

on the same hardware.  Note, gfc = gfortran, and ifc = ifort.

The question then becomes what are the algorithms that Intel uses
to compute the last four functions.  I purposed one possible
algorithm for sinq() where an argument (after argument reduction
to a suitable range) is split into 3 double pieces.  Here, I am
assuming that at most a double with 53-bits of precision is
available on all targets.  One then uses trig identities to
replace sinq(x) by products on sin() and cos() of the 3 doubles.

Another possible algorithm would be to use a lookup tables
with a 100 or so entries and use a Taylor's series about
the table entries to determine the desired value.  One then
has

sinq(x) = sin(x_i) + dx * cos(x_i) + O(dx**2)

with dx = x - x_i, 0 <= x_i < pi / 4, and i = 0, ..., N.  Here,
one needs to chose the x_i and N to give the desired accuracy
where sin(x_i) and cos(x_i) are represented by 2 parts (the high
53-bits and low 53-bit).

-- 
Steve