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

/* Load/store fp control word, stolen from ieeefp.h. */
#define __fldcw(addr)  __asm __volatile("fldcw %0" : : "m" (*(addr)))
#define __fnstcw(addr) __asm __volatile("fnstcw %0" : "=m" (*(addr)))

/* Floating point precision control. */
#define FP_PC_MASK 0x300
#define FP_PC_SHIFT 8
enum fp_pc {
	FP_PC_REAL4 = 0,
	FP_PC_REAL8 = 2,
	FP_PC_REAL10 = 3
};

static const char * const FP_PC_STR[] = {
	[FP_PC_REAL4] = "REAL4",
	[FP_PC_REAL8] = "REAL8",
	[FP_PC_REAL10] = "REAL10"
};

static void
fp_pc_check(enum fp_pc pc)
{

	switch (pc) {
	case FP_PC_REAL4:
	case FP_PC_REAL8:
	case FP_PC_REAL10:
		break;
	default:
		abort();
	}
}

static const char *
fp_pc_str(enum fp_pc pc)
{

	fp_pc_check(pc);
	return (FP_PC_STR[pc]);
}

static enum fp_pc
fp_get_precision(void)
{
	unsigned cw;
	enum fp_pc pc;

	__fnstcw(&cw);
	pc = (enum fp_pc)((cw & FP_PC_MASK) >> FP_PC_SHIFT);
	return (pc);
}

static void
fp_set_precision(enum fp_pc pc)
{
	unsigned cw;

	fp_pc_check(pc);

	__fnstcw(&cw);
	cw &= ~FP_PC_MASK;
	cw |= (unsigned)pc << FP_PC_SHIFT;
	__fldcw(&cw);
}

/*
 * Floating point operations.
 *
 * Access via volatile pointers to force load from and store to memory.
 */

static void
fsubl(volatile long double *r, volatile long double *x,
    volatile long double *y)
{

	*r = *x - *y;
}

static void
fdivl(volatile long double *r, volatile long double *x,
    volatile long double *y)
{

	*r = *x / *y;
}

static void
demo(void)
{
	long double max, epsilon, two, r;
	enum fp_pc pc;

	pc = fp_get_precision();
	printf("LDBL_MAX demo, mode %s (%#x)\n", fp_pc_str(pc), (unsigned)pc);

	max = LDBL_MAX;
	epsilon = LDBL_EPSILON;
	printf("max                   = %La\n", max);
	printf("epsilon               = %La\n", epsilon);

	fsubl(&r, &max, &epsilon);
	printf("max - epsilon         = %La\n", r);

	fsubl(&r, &max, &max);
	printf("max - max             = %La\n", r);

	two = 2;
	fdivl(&r, &max, &two);
	printf("max / 2               = %La\n", r);

	r = nextafterl(max, INFINITY);
	printf("nextafterl(max, +inf) = %La\n", r);

	r = nextafterl(max, -INFINITY);
	printf("nextafterl(max, -inf) = %La\n", r);
}

int
main(void)
{

	demo();

	printf("\n");
	fp_set_precision(FP_PC_REAL10);
	demo();

	return (0);
}