/*
** =====================================================================
** Demonstrate rounding behavior near the underflow limit.
**
** On AMD64, the output for the 32-bit version is
**
**  	u i s    e         x             y              x*y
**  	0 0 1 0x1.0p-01 0x1.800000p+0 0x1.000000p-127 0x1.800000p-127
**  	0 0 1 0x1.0p-02 0x1.400000p+0 0x1.800000p-127 0x1.e00000p-127
**  	0 0 1 0x1.0p-03 0x1.200000p+0 0x1.c00000p-127 0x1.f80000p-127
**  	0 0 1 0x1.0p-04 0x1.100000p+0 0x1.e00000p-127 0x1.fe0000p-127
**  	0 0 1 0x1.0p-05 0x1.080000p+0 0x1.f00000p-127 0x1.ff8000p-127
**  	0 0 1 0x1.0p-06 0x1.040000p+0 0x1.f80000p-127 0x1.ffe000p-127
**  	0 0 1 0x1.0p-07 0x1.020000p+0 0x1.fc0000p-127 0x1.fff800p-127
**  	0 0 1 0x1.0p-08 0x1.010000p+0 0x1.fe0000p-127 0x1.fffe00p-127
**  	0 0 1 0x1.0p-09 0x1.008000p+0 0x1.ff0000p-127 0x1.ffff80p-127
**  	0 0 1 0x1.0p-10 0x1.004000p+0 0x1.ff8000p-127 0x1.ffffe0p-127
**  	0 0 1 0x1.0p-11 0x1.002000p+0 0x1.ffc000p-127 0x1.fffff8p-127
**  	1 1 0 0x1.0p-12 0x1.001000p+0 0x1.ffe000p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-13 0x1.000800p+0 0x1.fff000p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-14 0x1.000400p+0 0x1.fff800p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-15 0x1.000200p+0 0x1.fffc00p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-16 0x1.000100p+0 0x1.fffe00p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-17 0x1.000080p+0 0x1.ffff00p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-18 0x1.000040p+0 0x1.ffff80p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-19 0x1.000020p+0 0x1.ffffc0p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-20 0x1.000010p+0 0x1.ffffe0p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-21 0x1.000008p+0 0x1.fffff0p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-22 0x1.000004p+0 0x1.fffff8p-127 0x1.000000p-126
**  	0 1 0 0x1.0p-23 0x1.000002p+0 0x1.fffffcp-127 0x1.000000p-126
**
** On SPARC64, the output differs in the u value in lines 13--23:
**
**  	u i s    e         x             y              x*y
**  	0 0 1 0x1.0p-01 0x1.800000p+0 0x1.000000p-127 0x1.800000p-127
**  	0 0 1 0x1.0p-02 0x1.400000p+0 0x1.800000p-127 0x1.e00000p-127
**  	0 0 1 0x1.0p-03 0x1.200000p+0 0x1.c00000p-127 0x1.f80000p-127
**  	0 0 1 0x1.0p-04 0x1.100000p+0 0x1.e00000p-127 0x1.fe0000p-127
**  	0 0 1 0x1.0p-05 0x1.080000p+0 0x1.f00000p-127 0x1.ff8000p-127
**  	0 0 1 0x1.0p-06 0x1.040000p+0 0x1.f80000p-127 0x1.ffe000p-127
**  	0 0 1 0x1.0p-07 0x1.020000p+0 0x1.fc0000p-127 0x1.fff800p-127
**  	0 0 1 0x1.0p-08 0x1.010000p+0 0x1.fe0000p-127 0x1.fffe00p-127
**  	0 0 1 0x1.0p-09 0x1.008000p+0 0x1.ff0000p-127 0x1.ffff80p-127
**  	0 0 1 0x1.0p-10 0x1.004000p+0 0x1.ff8000p-127 0x1.ffffe0p-127
**  	0 0 1 0x1.0p-11 0x1.002000p+0 0x1.ffc000p-127 0x1.fffff8p-127
**  	1 1 0 0x1.0p-12 0x1.001000p+0 0x1.ffe000p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-13 0x1.000800p+0 0x1.fff000p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-14 0x1.000400p+0 0x1.fff800p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-15 0x1.000200p+0 0x1.fffc00p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-16 0x1.000100p+0 0x1.fffe00p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-17 0x1.000080p+0 0x1.ffff00p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-18 0x1.000040p+0 0x1.ffff80p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-19 0x1.000020p+0 0x1.ffffc0p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-20 0x1.000010p+0 0x1.ffffe0p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-21 0x1.000008p+0 0x1.fffff0p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-22 0x1.000004p+0 0x1.fffff8p-127 0x1.000000p-126
**  	1 1 0 0x1.0p-23 0x1.000002p+0 0x1.fffffcp-127 0x1.000000p-126
**
** Thus, AMD64 detects underflow AFTER rounding, whereas SPARC64
** detects underflow BEFORE rounding.
**
**   [03-Aug-2021]
** =====================================================================
*/

#include <fenv.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static const float MINNORMAL = +0x1.0p-126F;

int
issubnormalf(float x)
{
    return (((-MINNORMAL < x) && (x < MINNORMAL)) && (x != 0.0F));
}

int
main(void)
{
    int i, s, u;
    volatile float e, w, x, y, z;

    (void)printf("MINNORMAL = %.1a\n\n", (double)MINNORMAL);

    (void)printf("u i s    e         x             y              x*y\n");

    for (e = 0.5F; (w = 1.0F + e, w != 1.0F); e *= 0.5F)
    {
	x = 1.0F + e;
	y = (1.0F - e) * MINNORMAL;
	(void)feclearexcept(FE_UNDERFLOW | FE_INEXACT);
	z = x * y;
	u = (fetestexcept(FE_UNDERFLOW) != 0);
	i = (fetestexcept(FE_INEXACT)   != 0);
	s = issubnormalf(z);
	(void)printf("%d %d %d %.1a %.6a %.6a %.6a\n", u, i, s, e, x, y, z);
    }

    return (EXIT_SUCCESS);
}