/***********************************************************************
Report sizeof() values, and the machine epsilon and smallest
floating-point number, to investigate floating-point support.

[21-Nov-2001]
***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <float.h>

float fstore(float x)
{
    return (x);
}

double dstore(double x)
{
    return (x);
}

#if defined(HAVE_LONGDOUBLE)
long double ldstore(long double x)
{
    return (x);
}

long double fppow2(int n)
{
    long double x, power;
    x = (n < 0) ? (1.0L/2.0L) : 2.0L;
    n = (n < 0) ? -n : n;
    power = 1.0L;
    while (n-- > 0)
	power *= x;
    return (power);
}
#else
double fppow2(int n)
{
    double x, power;
    x = (n < 0) ? (1.0L/2.0L) : 2.0L;
    n = (n < 0) ? -n : n;
    power = 1.0L;
    while (n-- > 0)
	power *= x;
    return (power);
}
#endif

void test_float(void)
{
    float x;

    (void)printf("float:\n");

    (void)printf("\tsizeof(float) =            %2u\n", (unsigned int)sizeof(float));
    (void)printf("\tFLT_MANT_DIG =             %2d\n", FLT_MANT_DIG);

    x = 1.0F;
    while (fstore(1.0F + x/2.0F) != 1.0F)
	x /= 2.0F;
    (void)printf("\tmachine epsilon =           %13.5e  ", x);

    if (x == (float)fppow2(-23))
	(void)printf("[IEEE 754 32-bit macheps]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");

    x = 1.0F;
    while (fstore(x / 2.0F) != 0.0F)
	x /= 2.0F;
    (void)printf("\tsmallest positive number =  %13.5e  ", x);

    if (x == (float)fppow2(-149))
	(void)printf("[IEEE 754 smallest 32-bit subnormal]\n");
    else if (x == (float)fppow2(-126))
	(void)printf("[IEEE 754 smallest 32-bit normal]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");
    (void)printf("\n");
}

void test_double(void)
{
    double x;

    (void)printf("double:\n");

    (void)printf("\tsizeof(double) =           %2u\n", (unsigned int)sizeof(double));
    (void)printf("\tDBL_MANT_DIG =             %2d\n", DBL_MANT_DIG);

    x = 1.0;
    while (dstore(1.0 + x/2.0) != 1.0)
	x /= 2.0;
    (void)printf("\tmachine epsilon =           %13.5le  ", x);

    if (x == (double)fppow2(-52))
	(void)printf("[IEEE 754 64-bit macheps]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");

    x = 1.0;
    while (dstore(x / 2.0) != 0.0)
	x /= 2.0;
    (void)printf("\tsmallest positive number =  %13.5le  ", x);

    if (x == (double)fppow2(-1074))
	(void)printf("[IEEE 754 smallest 64-bit subnormal]\n");
    else if (x == (double)fppow2(-1022))
	(void)printf("[IEEE 754 smallest 64-bit normal]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");
    (void)printf("\n");
}

void test_long_double(void)
{
#if defined(HAVE_LONGDOUBLE)
    long double x;

    (void)printf("long double:%s\n",
		 (sizeof(double) == sizeof(long double)) ?
		 "  [**WARNING**: long double is really double]" : "");

    (void)printf("\tsizeof(long double) =      %2u\n", (unsigned int)sizeof(long double));
    (void)printf("\tLDBL_MANT_DIG =            %2d\n", LDBL_MANT_DIG);

    x = 1.0L;
    while (ldstore(1.0L + x/2.0L) != 1.0L)
	x /= 2.0L;
    (void)printf("\tmachine epsilon =           %13.5Le  ", x);

    if (x == (long double)fppow2(-52))
	(void)printf("[IEEE 754 64-bit macheps]\n");
    else if (x == (long double)fppow2(-63))
	(void)printf("[IEEE 754 80-bit macheps]\n");
    else if (x == (long double)fppow2(-112))
	(void)printf("[IEEE 754 128-bit macheps]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");

    x = 1.0L;
    while (ldstore(x / 2.0L) != 0.0L)
	x /= 2.0L;
    (void)printf("\tsmallest positive number =  %13.5Le  ", x);

#if LDBL_MANT_DIG == 53
    if (x == (long double)fppow2(-1074))
	(void)printf("[IEEE 754 smallest 64-bit subnormal]\n");
    else if (x == (long double)fppow2(-1022))
	(void)printf("[IEEE 754 smallest 64-bit normal]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");
#elif LDBL_MANT_DIG == 64
    if (x == (long double)fppow2(-16446))
	(void)printf("[IEEE 754 smallest 80-bit subnormal]\n");
    else if (x == (long double)fppow2(-16382))
	(void)printf("[IEEE 754 smallest 80-bit normal]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");
#elif LDBL_MANT_DIG == 113
    if (x == (long double)fppow2(-16494))
	(void)printf("[IEEE 754 smallest 128-bit subnormal]\n");
    else if (x == (long double)fppow2(-16382))
	(void)printf("[IEEE 754 smallest 128-bit normal]\n");
    else
	(void)printf("[not IEEE 754 conformant]\n");
#else
    (void)printf("[not IEEE 754 conformant]\n");
#endif

#endif
    (void)printf("\n");
}


int
main(int argc, char* argv[])
{
    test_float();
    test_double();
    test_long_double();

    return (EXIT_SUCCESS);
}