On Tue, Sep 3, 2013 at 2:46 PM, Glenn Fowler <[email protected]> wrote:
> On Tue, 3 Sep 2013 02:50:31 +0200 Tina Harriott wrote:
>> I'm currently looking into a complex mathematical simulation (fluid
>> dynamics) and found references to NaN (Not A Number) with payloads.
>> The Nan values are used to represent 'no data' in the datasets and the
>> payloads are used for context information. The code has almost 2
>> million lines of code and removing the feature would be a major
>> effort.
>
>> As we are in the process of using ksh93 as "interface" to the
>> simulation using libshell we need access to Nan values and their
>> payloads.
>
>> Can anyone give an estimate how hard it would be to add support for
>> Nan payloads to ksh93?
>
> I'm aware of NaN payload but have not seen C code that
>         * has MIN/MAX macros for the min/max payload value for a give FP type

See code attached code example... the maximum is the mask-1 ... the
minimum is either 0 or 1... I don't remember anymore...

>         * sets a NaN value with a specific payload
>         * extracts a payload from a NaN
> are there standard api's for this?

Yes... the function |nan()| - see
http://pubs.opengroup.org/onlinepubs/009695399/functions/nan.html
-- snip --
       double nan(const char *tagp);
       float nanf(const char *tagp);
       long double nanl(const char *tagp);
-- snip --

The trouble is... the exact string input is
compiler/libc/platforms-specific... and beyond highend HPC platforms
most platforms do not implement |nan()| with payload support.

The other issue is that often compiler/libc/platforms-combinations to
not preserve the payload when a value is casted between _different_
floating-point types. That problem isn't new and has haunted HPC
programmers since eternity... but the solution is to simply test with
C99 |isnan()| and then add some copy-payload code. Since |isnan()| is
basically a bit-pattern test the overhead is insignificant (on SX-8X
that's around 8 instructions more compared to a plain cast).

If we implement NaN's with payloads we have to make sure that _every_
cast between |Sfdouble| (basically a |long double|) and another
floating-point type is guarded with such copy-payload code... except
when we pass the value to a arithmetricfunction like |sin()|/|cos()|
etc. which (intentionally) discard the payload if they pass the Nan
through.

Attached (as "nan_payload_test1_20130904_001.c.txt") is some prototype
code which shows now to handle NaN payloads...

* Notes:
- The code explicitly requires ISO C99 semantics
- The code relies on |__int128_t| to handle |long double|. Originally
I thought about handling the code using byte-by-byte masking... but it
turned out that this creates a lot of mess and is endian-specific. The
final code may have to resort to that to be portable but for now I
use |__int128_t| to keep the code readable

----

Bye,
Roland

-- 
  __ .  . __
 (o.\ \/ /.o) [email protected]
  \__\/\/__/  MPEG specialist, C&&JAVA&&Sun&&Unix programmer
  /O /==\ O\  TEL +49 641 3992797
 (;O/ \/ \O;)
/***********************************************************************
*                                                                      *
*               This software is part of the ast package               *
*            Copyright (c) 2013 AT&T Intellectual Property             *
*                      and is licensed under the                       *
*                  Common Public License, Version 1.0                  *
*                    by AT&T Intellectual Property                     *
*                                                                      *
*                A copy of the License is available at                 *
*            http://www.opensource.org/licenses/cpl1.0.txt             *
*         (with md5 checksum 059e8cd6165cb4c31e351f2b69388fd9)         *
*                                                                      *
*              Information and Software Systems Research               *
*                            AT&T Research                             *
*                           Florham Park NJ                            *
*                                                                      *
*              Roland Mainz <[email protected]>                 *
*                                                                      *
***********************************************************************/

#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>

/*
 * Notes about portablity:
 * - The ISO C99 standards supports a portable way to write the
 *   payload of a quiet NaN in ISO C.
 * - The |nan()|/|nanf()|/|nanl()| functions (from <math.h>, see
 *   section 7.12.11.2 in the ISO C99 spec) accept strings as
 *   arguments.
 * - |strtof()|/|strtod()|/|strtold()| functions (from <stdlib.h>,
 *   see section 7.20.1.3) accept strings in the form
 *   "NAN(character sequence)".
 * - |fscanf()|+|sscanf()| follow |strtod()|
 * - the character sequence in "NAN(character sequence)" is
 *   interpreted in an (libc+compiler) implementation-specific
 *   way (=not portable).
 */

void print_bits_32(uint32_t x)
{
        int i;
        for (i=0 ; i < 32 ; i++)
        {
                /* fixme: endian! */
                printf("%d", (x & (1<<31))?1:0);
                x = x << 1;
        }
}
void print_bits_64(uint64_t x)
{
        int i;
        for (i=0 ; i < 64 ; i++)
        {
                /* fixme: endian! */
                printf("%d", (x & (1LL<<63))?1:0);
                x = x << 1;
        }
}
void print_bits_128(__int128_t x)
{
        int i;
        for (i=0 ; i < 128 ; i++)
        {
                /* fixme: endian! */
                printf("%d", (x & (((__int128_t)1)<<127))?1:0);
                x = x << 1;
        }
}

void set_nan_payload_binary32(float *val, unsigned int payload)
{
        if (!isnanf(*val))
                return;

        volatile uint32_t *bitval_ptr = ((volatile uint32_t *)val);

        const uint32_t exp_bits = log2f(FLT_MAX_EXP)+1;
        const uint32_t frac_bits = 32-1-exp_bits;
        const uint32_t frac_mask = powf(2, frac_bits-1)-1;

        *bitval_ptr &= ~frac_mask;
        *bitval_ptr |= (payload & frac_mask);
}

unsigned int get_nan_payload_binary32(float x)
{
        unsigned int payload;

        if(!isnanf(x))
                return 0;

        const uint32_t exp_bits = log2(FLT_MAX_EXP)+1;
        const uint32_t frac_bits = 32-1-exp_bits;
        const uint32_t frac_mask = pow(2, frac_bits-1)-1;

        volatile uint32_t *bitval_ptr = ((volatile uint32_t *)&x);

        payload = *bitval_ptr & frac_mask;

        return payload;
}

void set_nan_payload_binary64(double *val, unsigned long payload)
{
        if (!isnan(*val))
                return;

        volatile uint64_t *bitval_ptr = ((volatile uint64_t *)val);

        const uint64_t exp_bits = log2(DBL_MAX_EXP)+1;
        const uint64_t frac_bits = 64-1-exp_bits;
        const uint64_t frac_mask = pow(2, frac_bits-1)-1;

        *bitval_ptr &= ~frac_mask;
        *bitval_ptr |= (payload & frac_mask);
}

unsigned long get_nan_payload_binary64(double x)
{
        unsigned long payload;

        if(!isnan(x))
                return 0UL;

        const uint64_t exp_bits = log2(DBL_MAX_EXP)+1;
        const uint64_t frac_bits = 64-1-exp_bits;
        const uint64_t frac_mask = pow(2, frac_bits-1)-1;

        volatile uint64_t *bitval_ptr = ((volatile uint64_t *)&x);

        payload = *bitval_ptr & frac_mask;

        return payload;
}

const int ldbl_bits=80/*128*/;

void set_nan_payload_binary80(long double *val, unsigned long payload)
{
        if (!isnanl(*val))
                return;

        volatile __int128_t *bitval_ptr = ((volatile __int128_t *)val);

        /* what about Intel i387 80bit |long double|? ? */

        const __int128_t exp_bits = log2l(LDBL_MAX_EXP)+1;
        const __int128_t frac_bits = ldbl_bits-1-exp_bits;
        const __int128_t frac_mask = ((__int128_t)powl(2, frac_bits-1))-1;

        *bitval_ptr &= ~frac_mask;
        *bitval_ptr |= (payload & frac_mask);
}

unsigned long get_nan_payload_binary80(long double x)
{
        unsigned long payload;

        if(!isnanl(x))
                return 0;

        const __int128_t exp_bits = log2l(LDBL_MAX_EXP)+1;
        const __int128_t frac_bits = ldbl_bits-1-exp_bits;
        const __int128_t frac_mask = ((__int128_t)powl(2, frac_bits-1))-1;

        volatile __int128_t *bitval_ptr = ((volatile __int128_t *)&x);

        payload = *bitval_ptr & frac_mask;

        return payload;
}

float cast_binary64_to_binary32(double d)
{
        float f;

        if (isnan(d))
        {
                unsigned long payload;

                payload = get_nan_payload_binary64(d);
                f = (float)d;

                set_nan_payload_binary32(&f, payload);
        }
        else
        {
                f = (float)d;
        }

        return f;
}


float cast_binary80_to_binary32(long double ld)
{
        float f;

        if (isnanl(ld))
        {
                unsigned long payload;

                payload = get_nan_payload_binary80(ld);
                f = (float)ld;

                set_nan_payload_binary32(&f, payload);
        }
        else
        {
                f = (float)ld;
        }

        return f;
}


double cast_binary80_to_binary64(long double ld)
{
        double d;

        if (isnanl(ld))
        {
                unsigned long payload;

                payload = get_nan_payload_binary80(ld);
                d = (double)ld;

                set_nan_payload_binary64(&d, payload);
        }
        else
        {
                d = (double)ld;
        }

        return d;
}

int main(int ac, char *av[])
{
        volatile float val=NAN;

        long i;

        uint32_t *xp=(uint32_t *)&val;

        printf("orig:\t" ); print_bits_32(*xp); printf("\tval=%f\n", val);

        const uint32_t exp_bits = log2(FLT_MAX_EXP)+1;
        const uint32_t frac_bits = 32-1-exp_bits;
        const uint32_t frac_mask = pow(2, frac_bits-1)-1;

        printf("exp_bits=%d, frac_bits=%d, frac_mask=%lx\n", exp_bits, 
frac_bits, (unsigned long)frac_mask);

        printf("frcmsk:\t" ); print_bits_32(frac_mask); printf("\n");

        *xp |= 0x01;
        printf("0x01:\t" ); print_bits_32(*xp); printf("\tval=%f\n", val);

        for (i=0 ; i <= frac_mask ; i++)
        {
                volatile uint32_t *bitval_ptr = ((volatile uint32_t *)&val);

                val = NAN;
                *bitval_ptr &= ~frac_mask;
                *bitval_ptr |= (i & frac_mask);
                if (!isnan(val))
                        printf("%f not a nan for %ld\n", val, (long)i);
                if (isinf(val))
                        printf("%f is a inf for %ld\n", val, (long)i);
        }

        {
                float fval;
                fval = -NAN;
                set_nan_payload_binary32(&fval, 667);
                printf("flt: val with -NAN667=%f\n", fval);
                printf("flt: payload of -NAN667=%d\n", 
(int)get_nan_payload_binary32(fval));
        }

        {
                double dval;
                dval = -NAN;
                set_nan_payload_binary64(&dval, 667);
                printf("dbl: val with -NAN667=%f\n", dval);
                printf("dbl: payload of -NAN667=%lu\n", (unsigned 
long)get_nan_payload_binary64(dval));
                printf("dbl: payload of sys_cast float(-NAN667)=%lu\n", 
(unsigned long)get_nan_payload_binary32((float)dval));
                printf("dbl: payload of own_cast float(-NAN667)=%lu\n", 
(unsigned long)get_nan_payload_binary32(cast_binary64_to_binary32(dval)));
        }

        {
                long double ldval;
                ldval = -NAN;
                set_nan_payload_binary80(&ldval, 667);
                printf("ldbl: val with -NAN667=%Lf\n", ldval);
                printf("ldbl: payload of -NAN667=%lu\n", (unsigned 
long)get_nan_payload_binary80(ldval));
                printf("ldbl: payload of sys_cast double(-NAN667)=%lu\n", 
(unsigned long)get_nan_payload_binary64((double)ldval));
                printf("ldbl: payload of own_cast double(-NAN667)=%lu\n", 
(unsigned long)get_nan_payload_binary64(cast_binary80_to_binary64(ldval)));
        }

        return EXIT_SUCCESS;
}

_______________________________________________
ast-developers mailing list
[email protected]
http://lists.research.att.com/mailman/listinfo/ast-developers

Reply via email to