/*
** Floating point relationals test program
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

typedef struct fpCCR {
    int _exception_bits;	
    int _trap_enable_bits;	
    int _sticky_bits;
    int _rounding_mode;
    unsigned int _unused :16;
    unsigned int _op2	:3;
    unsigned int _dst	:3;
    unsigned int _op1	:3;
    unsigned int _last_operation : 7;
    union {
	float sf;
	double df;
	long double xf;
    } _operand1;
    union {
	float sf;
	double df;
	long double xf;
    } _operand2;
    union {
	int si;
	float sf;
	double df;
	long double xf;
    } _result;
} fpCCR_t;

#define INEXACT 1
#define UNDER	2
#define OVER	4
#define	DIV0	8
#define INVOP	16

#define ROUND_TO_NEAREST	0 // round result to nearest representable value
#define ROUND_TO_ZERO		1 // round result towards zero
#define ROUND_TO_PLUS		2 // round result towards plus infinity
#define ROUND_TO_MINUS		3 // round result towards minus infinity

extern fpCCR_t _fpCCR; 

typedef union {
    long double ld;
    struct {
	unsigned int s:1;
	unsigned int e:15;
	unsigned int z:16;
	int u;
	int l;
    } sei2;
    struct {
	unsigned int s:1;
	unsigned int e:15;
	unsigned int z:16;
	long long ll;
    } ll;
    struct {
	int se;
	int u;
	int l;
    } i3;
} ldoub;

typedef union {
    long long ll;
    struct {
	int u;
	int l;
    } i2;
} llong;

typedef union {
    double d;
    struct {
	unsigned long long s:1;
	unsigned long long e:11;
	unsigned long long f:52;
    } sef;
    struct {
	unsigned int s:1;
	unsigned int e:11;
	unsigned int u:20;
	int l;
    } sei;
    struct {
	int u;
	int l;
    } i2;
    long long ll;
} doub;

typedef union {
    float f;
    struct {
	unsigned long s:1;
	unsigned long e:8;
	unsigned long f:23;
    } sef;
    long l;
} flt;

extern void ran_init(unsigned init);
extern unsigned ran_nxt(void);

static double dftype(int nxt);
static float sftype(int nxt);
static long double xftype(int nxt);

static void __attribute__ ((interrupt)) fptrap(void);

static int trapped;

#define MAX_TYPES 16

int main(void) {
    char inps[132];
    double da,db;
    float sa,sb;
    long double xa,xb;
    int r,s;
    int i,j,k;
    int *dav,*dbv,*sav,*sbv,*xav,*xbv;
    char *op;
    
    // Enable the 5206e instruction cache:
    __asm__ __volatile__("movec %0,%/cacr" : : "d" (0x81000103));
    
    printf("Relational operator test\n");
    
    dav = (int *)&da;
    dbv = (int *)&db;
    sav = (int *)&sa;
    sbv = (int *)&sb;
    xav = (int *)&xa;
    xbv = (int *)&xb;
    
    ran_init(0);
    
//    _fpCCR._trap_enable_bits = INEXACT | UNDER | OVER | DIV0 | INVOP;

    // Set vector for TRAP #14
//    *(int *)(0x300000b8) = (int)fptrap;

    for (j = 0; j < MAX_TYPES; j += 1) {
	for (k = 0; k < MAX_TYPES; k += 1) {
	    da = dftype(j);
	    db = dftype(k);
	    for (i = 0; i < 6; i += 1) {
		_fpCCR._sticky_bits = 0;
		trapped = 0;
		switch (i) {
		case 0:
		    r = da == db;
		    op = "==";
		    break;
		case 1:
		    r = da != db;
		    op = "!=";
		    break;
		case 2:
		    r = da > db;
		    op = "> ";
		    break;
		case 3:
		    r = da >= db;
		    op = ">=";
		    break;
		case 4:
		    r = da < db;
		    op = "< ";
		    break;
		default:
		case 5:
		    r = da <= db;
		    op = "<=";
		    break;
		}
		s = _fpCCR._sticky_bits;
		printf("%08x%08x %s %08x%08x = %d   %c%c%c%c%c %c"
		    "   (%e %s %e = %d)\n",
		    dav[0],dav[1],op,dbv[0],dbv[1],r,(s&INVOP)?'I':'i',
		    (s&DIV0)?'D':'d',(s&OVER)?'O':'o',
		    (s&UNDER)?'U':'u',(s&INEXACT)?'X':'x',(trapped)?'T':'t',
		    da,op,db,r);
	    }
	    gets(inps);
	}
    }
    
    for (j = 0; j < MAX_TYPES; j += 1) {
	for (k = 0; k < MAX_TYPES; k += 1) {
	    sa = sftype(j);
	    sb = sftype(k);
	    for (i = 0; i < 6; i += 1) {
		_fpCCR._sticky_bits = 0;
		trapped = 0;
		switch (i) {
		case 0:
		    r = sa == sb;
		    op = "==";
		    break;
		case 1:
		    r = sa != sb;
		    op = "!=";
		    break;
		case 2:
		    r = sa > sb;
		    op = "> ";
		    break;
		case 3:
		    r = sa >= sb;
		    op = ">=";
		    break;
		case 4:
		    r = sa < sb;
		    op = "< ";
		    break;
		case 5:
		    r = sa <= sb;
		    op = "<=";
		    break;
		}
		s = _fpCCR._sticky_bits;
		printf("%08x %s %08x = %d   %c%c%c%c%c %c"
		    "   (%e %s %e = %d)\n",
		    sav[0],op,sbv[0],r,(s&INVOP)?'I':'i',
		    (s&DIV0)?'D':'d',(s&OVER)?'O':'o',
		    (s&UNDER)?'U':'u',(s&INEXACT)?'X':'x',(trapped)?'T':'t',
		    sa,op,sb,r);
	    }
	    gets(inps);
	}
    }

    for (j = 0; j < MAX_TYPES; j += 1) {
	for (k = 0; k < MAX_TYPES; k += 1) {
	    xa = xftype(j);
	    xb = xftype(k);
	    for (i = 0; i < 6; i += 1) {
		_fpCCR._sticky_bits = 0;
		trapped = 0;
		switch (i) {
		case 0:
		    r = xa == xb;
		    op = "==";
		    break;
		case 1:
		    r = xa != xb;
		    op = "!=";
		    break;
		case 2:
		    r = xa > xb;
		    op = "> ";
		    break;
		case 3:
		    r = xa >= xb;
		    op = ">=";
		    break;
		case 4:
		    r = xa < xb;
		    op = "< ";
		    break;
		default:
		case 5:
		    r = xa <= xb;
		    op = "<=";
		    break;
		}
		s = _fpCCR._sticky_bits;
		printf("%08x%08x%08x %s %08x%08x%08x = %2d   %c%c%c%c%c %c\n",
		    xav[0],xav[1],xav[2],op,xbv[0],xbv[1],xbv[2],r,
		    (s&INVOP)?'I':'i',(s&DIV0)?'D':'d',(s&OVER)?'O':'o',
		    (s&UNDER)?'U':'u',(s&INEXACT)?'X':'x',(trapped)?'T':'t');
	    }
	    gets(inps);
	}
    }
    
    printf("Done!\n");
    return(0);
}

/*
** Generate various types of doubles
*/
static double dftype(int nxt) {
    double q;
    int exp;
    unsigned u,l;
    int sgn = nxt & 1;
    
    switch (nxt >> 1) {
    case 0:			// Normal finites
	exp = ((2046 * (ran_nxt() >> 12)) >> 20) + 1;
	u = (ran_nxt() & 0x000fffff) | (exp << 20);
	l = ran_nxt();
	break;
	
    case 1:			// Signed zero
	u = 0;
	l = 0;
	break;
	
    case 2:			// Denorm
	do {
	    exp = (52 * (ran_nxt() >> 6)) >> 26;
	    u = 0;
	    if (exp >= 20) {
		l = ran_nxt() >> (exp - 20);
	    } else {
		l = ran_nxt();
		u = ran_nxt() >> (exp + 12);
	    }
	} while (!l && !u);
	break;
	
    case 3:			// Power of 2
	exp = ((2046 * (ran_nxt() >> 12)) >> 20) + 1;
	u = exp << 20;
	l = 0;
	break;
	
    case 4:			// Signed max
	u = 0x7fefffff;
	l = 0xffffffff;
	break;
	
    case 5:			// Signed infinity
	u = 0x7ff00000;
	l = 0;
	break;
	
    case 6:			// Quiet NaN
	u = (ran_nxt() >> 1) | 0x7ff80000;
	l = ran_nxt();
	break;
	
    case 7:			// Signaling NaN
	u = (ran_nxt() >> 13) | 0x7ff00000;
	l = ran_nxt();
	break;
	
    default:
	u = l = 0;
	break;
    }
    
    u |= sgn << 31;
    
    *(int *)&q = u;
    *((int *)&q + 1) = l;
    return(q);
}

/*
** Generate various types of single floats
*/
static float sftype(int nxt) {
    float q;
    int exp;
    unsigned u;
    int sgn = nxt & 1;
    
    switch (nxt >> 1) {
    case 0:			// Normal finites
	exp = ((254 * (ran_nxt() >> 12)) >> 20) + 1;
	u = (ran_nxt() & 0x000fffff) | (exp << 23);
	break;
	
    case 1:			// Signed zero
	u = 0;
	break;
	
    case 2:			// Denorm
	do {
	    exp = (23 * (ran_nxt() >> 6)) >> 26;
	    u = ran_nxt() >> (exp + 9);
	} while (!u);
	break;
	
    case 3:			// Power of 2
	exp = ((254 * (ran_nxt() >> 12)) >> 20) + 1;
	u = exp << 23;
	break;
	
    case 4:			// Signed max
	u = 0x7f7fffff;
	break;
	
    case 5:			// Signed infinity
	u = 0x7f800000;
	break;
	
    case 6:			// Quiet NaN
	u = (ran_nxt() >> 1) | 0x7fc00000;
	break;
	
    case 7:			// Signaling NaN
	u = (ran_nxt() >> 10) | 0x7f800000;
	break;
	
    default:
	u = 0;
	break;
    }
    
    u |= sgn << 31;
    
    *(int *)&q = u;
    return(q);
}

/*
** Generate various types of long doubles
*/
static long double xftype(int nxt) {
    ldoub q;
    llong x;
    long long b;
    int p;
    int sgn = nxt & 1;
    
    q.sei2.z = 0;
    switch (nxt >> 1) {
    case 0:			// Normalized finites
	x.i2.u = 0;
	x.i2.l = ran_nxt();
	x.ll *= 32767;//1022;
	q.sei2.e = x.i2.u;// + 0x3c00;
	q.sei2.u = ran_nxt() | 0x80000000;
	q.sei2.l = ran_nxt();
	break;
	
    case 1:			// Signed zero
	q.sei2.e = 0;
	q.sei2.u = 0;
	q.sei2.l = 0;
	break;
	
    case 2:			// Unnormalized number
	x.i2.u = 0;
	x.i2.l = ran_nxt();
	x.ll *= 32767;
	q.sei2.e = x.i2.u;
	p = (64 * (ran_nxt() >> 6)) >> 26;
	q.sei2.l = ran_nxt();
	q.sei2.u = ran_nxt();
	b = 1;
	b <<= (p);
	q.ll.ll &= (b - 1);
	q.ll.ll |= b;
	break;
	
    case 3:			// Power of 2
	x.i2.u = 0;
	x.i2.l = ran_nxt();
	x.ll *= 32767;
	q.sei2.e = x.i2.u;
	q.sei2.l = 0;
	q.sei2.u = 0x80000000;
	break;
	
    case 4:			// Signed max
	q.sei2.e = 0x7ffe;
	q.sei2.u = 0xffffffff;
	q.sei2.l = 0xffffffff;
	break;
	
    case 5:			// Signed infinity
	q.sei2.e = 0x7fff;
	q.sei2.u = ran_nxt() & 0x80000000;
	q.sei2.l = 0;
	break;
	
    case 6:			// Quiet NaN
	q.sei2.e = 0x7fff;
	q.sei2.u = ran_nxt() | 0x40000000;
	q.sei2.l = ran_nxt();
	break;
	
    case 7:			// Signaling NaN
	q.sei2.e = 0x7fff;
	q.sei2.u = ran_nxt() & 0xbfffffff;
	q.sei2.l = ran_nxt();
	break;
	
    default:
	break;
    }
    q.sei2.s = sgn;
    
    return(q.ld);
}

static void __attribute__ ((interrupt)) fptrap(void) {
    trapped = 1;
    return;
}