/*
** =====================================================================
** 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>

#define PRINTF	(void)printf

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

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

extern int m68k_fpcr(void);
extern int m68k_fpsr(void);

void
printfenv(void)
{
    fenv_t fenv;	

    (void)fegetenv(&fenv);	
    PRINTF("fenv.__control_register     = 0x%08x\n", fenv.__control_register);
    PRINTF("fenv.__status_register      = 0x%08x\n", fenv.__status_register);
    PRINTF("fenv.__instruction_address  = 0x%08x\n", fenv.__instruction_address);
    (void)fesetenv(&fenv);
}

int
fetestexcept(int excepts)
{
    int fpcr, fpsr, result;

    fpcr = m68k_fpcr();
    fpsr = m68k_fpsr();
    printfenv();

    PRINTF("DEBUG: sizeof(int) = %d\n", (int)sizeof(int));
    PRINTF("DEBUG: fpcr = 0x%08x\n", fpcr);
    PRINTF("DEBUG: fpsr = 0x%08x\n", fpsr);

    result = 0;

    if ((excepts & FE_INEXACT) && (fpsr & 0x0080))
	result |= FE_INEXACT;

    if ((excepts & FE_DIVBYZERO) && (fpsr & 0x0100))
	result |= FE_DIVBYZERO;

    if ((excepts & FE_UNDERFLOW) && (fpsr & 0x0200))
	result |= FE_UNDERFLOW;

    if ((excepts & FE_OVERFLOW) && (fpsr & 0x0400))
	result |= FE_OVERFLOW;

    if ((excepts & FE_INVALID) && (fpsr & 0x0800))
	result |= FE_INVALID;
}

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

    PRINTF("MINNORMAL = %.1a\n\n", (double)MINNORMAL);
    printfenv();	

    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);
	PRINTF("%d %d %d %.1a %.6a %.6a %.6a\n", u, i, s, e, x, y, z);
    }

    printfenv();	

    return (EXIT_SUCCESS);
}