[Gnuplot-info] Configurable FP precision (Re: accuracy in calculation)

[Nikita Z. <coo...@ma...>]

On Wed, 10 Apr 2019 09:56:07 -0700
sfeam <eam...@gm...> wrote:

> On Wednesday, 10 April 2019 18:11:19 Patrick Dupre wrote:
> > Hello,
> >  
> > I guess that it may be a good illustration of what I get.
> > running under gnuplot
> >  
> > I={0,1}
> > set sample 1001
> > sqrt_pi = sqrt (pi)
> > sqrt_ln2 = sqrt (log(2))
> > z_ (x, wL, wG) = (x + I * abs (wL)) / wG
> > w_ (z, wG) = faddeeva (z) / wG / sqrt_pi
> > w_prim (z, wG) = (-2 * z * w_ (z, wG) + 2 * I / pi / wG) / wG
> > dd_Fadd (x, wL, wG) = (wGG = abs (wG) / sqrt_ln2, z = z_ (x, wL,
> > wGG), -2 * (w_ (z, wGG) / wGG + z * w_prim (z, wGG)) / wGG) plot
> > [-1000:1000] real (dd_Fadd (x, 270, 0.206)) plot [-1000:1000] real
> > (dd_Fadd (x, 270, 0.106))  
> 
> Since I am not familiar with these functions, I do not know what to
> look for or what to expect.
>   
> You say "the difference should never be zero even if it is close to
> 1e-10" That indeed is what seems to result from plotting
> 
>     abs( real (dd_Fadd (x, 270, 0.206)) -  real (dd_Fadd (x, 270,
> 0.106)))
> 
> Have I misunderstood the description?
> 
> 	Ethan
> 
> 

This side is possible to improve by adding build configuration option
to change FP type. Of course some options may be meaningless or
unavailable on some platforms (gcc targets), so i describe it with
x86_64 in mind.

I did own experiment while writing audio spectrum analyzer. For now it
is some playground for my experiments, and among them i tried
configurable FP type in hope it will enhance fftw analyzis precision,
since this library has all implementations - for float, double, long
double, quad (using gcc libquadmath).

Now about numeric FP type configuration.
config.h:
#define FFT_PRECISION <N>
where <N> - integer in range 1..3.

common.h:
#if FFT_PRECISION == 1
    typedef float real_t;
    #define real_rsize (sizeof(real_t))
    #define M_PIf M_PI
    #define REAL_EPSILON FLT_EPSILON
#elif FFT_PRECISION == 2
    typedef double real_t;
    #define real_rsize (sizeof(sample_t))
    #define fabsf      fabs
    #define fmodf      fmod
......
    #define crealf     creal
    #define cimagf     cimag
    #define nextafterf nextafter

    #define M_PIf M_PI
    #define REAL_EPSILON DBL_EPSILON
......
#elif FFT_PRECISION == 3
    typedef long double real_t;
.....
    #define cimagf     cimagl
    #define nextafterf nextafterl

    #define FLT_EPSILON LDBL_EPSILON

    #ifdef _GNU_SOURCE
        #define M_PIf M_PIl
    #else
        #define M_PIf M_PI
    #endif
#elif FFT_PRECISION == 4
    typedef float128 real_t;
    /* or typedef __float128 real_t */
.....
    #define cimagf     cimagq
    #define nextafterf nextafterq

    #ifdef _GNU_SOURCE
        #define M_PIf M_PIq
    #else
        #define M_PIf M_PI
    #endif
#end /* FFT_PRECISION */

This way proposes to use all math functions and other type-dependant
names with mandatory prefix. Prefix 'f' (sinf, absf) is selected
because name 'float' is more common than 'double' (d).

One implementation drawback here is that all float-specific names
(sinf, absf) become affected by configuration even when they should not
(if you eventually happen to use float type). It is ok, if no math used
with other than real_t data, but it might be worth to use own
'r' (real) prefix: absr, sinr, M_PIr...

I'm sure, this would be give huge precision boost, remembering that each
single bit adds 2x precision. __float128 of course is probably going to
be slower than 'long double', yet may be supported not everywhere, but
it looks excellent compromise between fixed-width language types and
multiple precision frameworks.

Once this configurable FP type system is established, it is possible in
future to easily add wider types, as they appear (if not yet exist).
For case, if one such new type will be called 'real' (using 'r'
prefix), would be better, if separate prefix is preferable, use
something more specific, e.g. gpfloat ('gpf' prefix).