/*
   This file is part of Fatgrind, a floating point accuracy tracking tool.

   Copyright (C) 2009 Stephane Marchesin

   This program is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public License as
   published by the Free Software Foundation; either version 2 of the
   License, or (at your option) any later version.

   This program is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307, USA.

   The GNU General Public License is contained in the file COPYING.
*/

// XXX Il y a un bug a cause de ce qui semble etre un probleme d'interaction GMP <-> valgrind quand on utilise les optimisations ASM.
// desactiver ces optimisations semble regler le probleme
// pour l'instant l'hypothese la plus probable est une corruption de l'etat de la FPU du programme hote par l'utilisation du MMX dans GMP:
// - execution du programme hote
// - context switch valgrind vers les callback fatgrind qui ne sauve pas la FPU
// - execution GMP qui utilise la FPU et donc change son etat
// - retour au programme hote
// - le programme hote n'apprecie pas le changement d'etat de sa FPU
// -> a creuser

#include "pub_tool_basics.h"
#include "pub_tool_tooliface.h"
#include "pub_tool_libcassert.h"
#include "pub_tool_libcbase.h"
#include "pub_tool_libcprint.h"
#include "pub_tool_mallocfree.h"
#include "pub_tool_machine.h"
#include "pub_tool_debuginfo.h"
#include "pub_tool_vki.h"
#include "gmp-4.2.4/gmp.h"

// 32 bits out of a 128 bit vector
#define Ity_V32 (Ity_V128+1)
// 64 bits out of a 128 bit vector
#define Ity_V64 (Ity_V128+2)

#define FILE_LEN              256

static int is_64bit=0;

//#define DEBUG

#define panic() panic_real(__FILE__,__func__, __LINE__)
static void panic_real(const char * file, const char* func, int line)
{
	VG_(printf)("\n");
	VG_(printf)("*************************************************\n");
	VG_(printf)("Panique a bord, droit dans le decor !\n");
	VG_(printf)("%s:%d  %s()\n",file,line,func);
	VG_(printf)("*************************************************\n");
	VG_(tool_panic)("");
}

static IRExpr** mkIRExprVec_11 ( IRExpr* arg1, IRExpr* arg2, IRExpr* arg3,
                         IRExpr* arg4, IRExpr* arg5, IRExpr* arg6,
                         IRExpr* arg7, IRExpr* arg8, IRExpr* arg9,
			 IRExpr* arg10, IRExpr* arg11) {
   IRExpr** vec = LibVEX_Alloc(12 * sizeof(IRExpr*));
   vec[0] = arg1;
   vec[1] = arg2;
   vec[2] = arg3;
   vec[3] = arg4;
   vec[4] = arg5;
   vec[5] = arg6;
   vec[6] = arg7;
   vec[7] = arg8;
   vec[8] = arg9;
   vec[9] = arg10;
   vec[10] = arg11;
   vec[11] = NULL;
   return vec;
}


typedef void (*GMP_1FUNC) (mpf_t, mpf_t);
typedef void (*GMP_2FUNC) (mpf_t, mpf_t, mpf_t);
typedef void (*GMP_3FUNC) (mpf_t, mpf_t, mpf_t, mpf_t);
typedef void (*GMP_4FUNC) (mpf_t, mpf_t, mpf_t, mpf_t, mpf_t);


static void mpf_recip(mpf_t r, mpf_t a)
{
	mpf_t one;

	mpf_init_set_ui(one, 1);
	mpf_div(r,one,a);

	mpf_clear(one);
}

static void mpf_rsqrt(mpf_t r, mpf_t a)
{
	mpf_t i,one;

	mpf_init(i);
	mpf_init_set_ui(one, 1);
	mpf_sqrt(i,a);
	mpf_div(r,one,i);

	mpf_clear(i);
	mpf_clear(one);
}

static void mpf_sqrt_round(mpf_t r, mpf_srcptr roundmode, mpf_srcptr  a)
{
	mpf_sqrt(r,a);
}

static void mpf_identity_round(mpf_t r, mpf_srcptr roundmode, mpf_srcptr  a)
{
	mpf_init_set(r,a);
}

static void mpf_add_round(mpf_t r, mpf_srcptr roundmode, mpf_srcptr  a, mpf_srcptr b)
{
	mpf_add(r,a,b);
}

static void mpf_sub_round(mpf_t r, mpf_srcptr roundmode, mpf_srcptr  a, mpf_srcptr b)
{
	mpf_sub(r,a,b);
}

static void mpf_mul_round(mpf_t r, mpf_srcptr roundmode, mpf_srcptr  a, mpf_srcptr b)
{
	mpf_mul(r,a,b);
}

static void mpf_div_round(mpf_t r, mpf_srcptr roundmode, mpf_srcptr  a, mpf_srcptr b)
{
	mpf_div(r,a,b);
}

static void print_float(float f)
{
	int i = (int)f;
	int i2 = (int)((f - (float)i)*100000.0);
	VG_(printf)("%d.%05d ",i,i2);
//	VG_(printf)("++%x -> %d.%05d++ ",*((unsigned*)(&f)),i,i2);
}

static void print_double(double f)
{
	int i = (int)f;
	int i2 = (int)((f - (float)i)*100000.0);
	VG_(printf)("%d.%05d ",i,i2);
//	VG_(printf)("++%llx -> %d.%05d++ ",*((unsigned long long*)(&f)),i,i2);
}

static struct 
{
	IROp op;
	char name[16];
	IRType t[4];
	IRType ret;
	void (*gmp_exec);
} fp_instr_list[] = 
{

{Iop_AddF64,		"AddF64",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_add_round},
{Iop_SubF64,		"SubF64",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_sub_round}, 
{Iop_MulF64,		"MulF64",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_mul_round}, 
{Iop_DivF64,		"DivF64",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_div_round},

{Iop_AddF64r32,		"AddF64r32",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_add_round},
{Iop_SubF64r32,		"SubF64r32",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_sub_round},
{Iop_MulF64r32,		"MulF64r32",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_mul_round},
{Iop_DivF64r32,		"DivF64r32",	{Ity_I32, 	Ity_F64, 	Ity_F64,	Ity_INVALID},	Ity_F64,	mpf_div_round},

//{Iop_NegF64,		"NegF64",	{Ity_F64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_INVALID,	mpf_neg},
//{Iop_AbsF64,		"AbsF64",	{Ity_F64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_INVALID,	mpf_abs},

{Iop_SqrtF64,		"SqrtF64",	{Ity_I32, 	Ity_F64, 	Ity_INVALID,	Ity_INVALID},	Ity_F64,	mpf_sqrt_round},
{Iop_SqrtF64r32,	"SqrtF64r32",	{Ity_I32, 	Ity_F64, 	Ity_INVALID,	Ity_INVALID},	Ity_F64,	mpf_sqrt_round},

{Iop_F64toI16,		"F64toI16",	{Ity_I32, 	Ity_F64, 	Ity_INVALID,	Ity_INVALID},	Ity_I16,	mpf_identity_round},
{Iop_F64toI32,		"F64toI32",	{Ity_I32, 	Ity_F64, 	Ity_INVALID,	Ity_INVALID},	Ity_I32,	mpf_identity_round},
{Iop_F64toI64,		"F64toI64",	{Ity_I32, 	Ity_F64, 	Ity_INVALID,	Ity_INVALID},	Ity_I64,	mpf_identity_round},

//{Iop_I16toF64,		"I16toF64",	{Ity_I16, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_F64,	mpf_identity_round},
//{Iop_I32toF64,		"I32toF64",	{Ity_I32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_F64,	mpf_identity_round},
{Iop_I64toF64,		"I64toF64",	{Ity_I32, 	Ity_I64, 	Ity_INVALID,	Ity_INVALID},	Ity_F64,	mpf_identity_round},

//{Iop_F32toF64,		"F32toF64",	{Ity_F32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_INVALID,	mpf_identity_round},
{Iop_F64toF32,		"F64toF32",	{Ity_I32, 	Ity_F64, 	Ity_INVALID,	Ity_INVALID},	Ity_F32,	mpf_identity_round},

{Iop_Add32F0x4,		"Add32F0x4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_add},
{Iop_Sub32F0x4,		"Sub32F0x4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_sub},
{Iop_Mul32F0x4,		"Mul32F0x4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_mul},
{Iop_Div32F0x4,		"Div32F0x4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_div},
{Iop_Recip32F0x4,	"Recip32F0x4",	{Ity_V32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_recip},
{Iop_Sqrt32F0x4,	"Sqrt32F0x4",	{Ity_V32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_sqrt},
{Iop_RSqrt32F0x4,	"RSqrt32F0x4",	{Ity_V32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_rsqrt},

{Iop_Add64F0x2,		"Add64F0x2",	{Ity_V64, 	Ity_V64, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_add},
{Iop_Sub64F0x2,		"Sub64F0x2",	{Ity_V64, 	Ity_V64, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_sub},
{Iop_Mul64F0x2,		"Mul64F0x2",	{Ity_V64, 	Ity_V64, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_mul},
{Iop_Div64F0x2,		"Div64F0x2",	{Ity_V64, 	Ity_V64, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_div},
{Iop_Recip64F0x2,	"Recip64F0x2",	{Ity_V64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_recip},
{Iop_Sqrt64F0x2,	"Sqrt64F0x2",	{Ity_V64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_sqrt},
{Iop_RSqrt64F0x2,	"RSqrt64F0x2",	{Ity_V64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_rsqrt},

{Iop_Add32Fx4,		"Add32Fx4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_add},
{Iop_Sub32Fx4,		"Sub32Fx4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_sub},
{Iop_Mul32Fx4,		"Mul32Fx4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_mul},
{Iop_Div32Fx4,		"Div32Fx4",	{Ity_V32, 	Ity_V32, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_div},
{Iop_Recip32Fx4,	"Recip32Fx4",	{Ity_V32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_recip},
{Iop_Sqrt32Fx4,		"Sqrt32Fx4",	{Ity_V32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_sqrt},
{Iop_RSqrt32Fx4,	"RSqrt32Fx4",	{Ity_V32, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V32,	mpf_rsqrt},

{Iop_Mul64Fx2,		"Mul64Fx2",	{Ity_V64, 	Ity_V64, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_mul},
{Iop_Div64Fx2,		"Div64Fx2",	{Ity_V64, 	Ity_V64, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_div},
{Iop_Recip64Fx2,	"Recip64Fx2",	{Ity_V64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_recip},
{Iop_Sqrt64Fx2,		"Sqrt64Fx2",	{Ity_V64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_sqrt},
{Iop_RSqrt64Fx2,	"RSqrt64Fx2",	{Ity_V64, 	Ity_INVALID, 	Ity_INVALID,	Ity_INVALID},	Ity_V64,	mpf_rsqrt},

};

/*

Instructions a ajouter :

Q-ops
		case Iop_MAddF64:
		case Iop_MAddF64r32:
		case Iop_MSubF64:
		case Iop_MSubF64r32:

Tri-ops
		case Iop_ScaleF64:
		case Iop_Yl2xF64:
		case Iop_Yl2xp1F64:
		case Iop_AtanF64:
		case Iop_PRemF64:
		case Iop_PRem1F64:
		case Iop_PRemC3210F64:
		case Iop_PRem1C3210F64:

Bin-ops
		case Iop_RoundF64toInt:
		case Iop_RoundF64toF32:
		case Iop_SinF64:
		case Iop_CosF64:
		case Iop_TanF64:
		case Iop_2xm1F64:

Un-ops
		case Iop_I32UtoFx4:
		case Iop_I32StoFx4:
		case Iop_QFtoI32Ux4_RZ:
		case Iop_QFtoI32Sx4_RZ:
		case Iop_RoundF32x4_RM:
		case Iop_RoundF32x4_RP:
		case Iop_RoundF32x4_RN:
		case Iop_RoundF32x4_RZ:
		case Iop_Est5FRSqrt:
		case Iop_RoundF64toF64_NEAREST:
		case Iop_RoundF64toF64_NegINF:
		case Iop_RoundF64toF64_PosINF:
		case Iop_RoundF64toF64_ZERO:
		case Iop_TruncF64asF32:
		case Iop_ReinterpF64asI64:
		case Iop_ReinterpI64asF64:
		case Iop_ReinterpI32asF32:
*/

static inline float as_float(unsigned long v)
{
	union {
		unsigned long in;
		float out;
	} u;
	u.in = v;
	return u.out;
}

static inline double as_double(unsigned long v1,unsigned long v2)
{
	union {
		unsigned long in[2];
		double out;
	} u;
	u.in[0] = v2;
	u.in[1] = v1;
	return u.out;
}

static inline unsigned long long as_i64(unsigned long v1, unsigned long v2)
{
	union {
		unsigned long in[2];
		unsigned long long out;
	} u;
	u.in[0] = v2;
	u.in[1] = v1;
	return u.out;
}

static inline unsigned int as_i32(unsigned long v)
{
	return v&0xffffffff;
}

static inline unsigned short as_i16(unsigned long v)
{
	return v&0xffff;
}

static void print(IRType type, unsigned long* v)
{
#ifdef DEBUG
	if (is_64bit)
	{
		switch(type)
		{
			case Ity_F64: 
			case Ity_V64:
				{
					float f = as_double(v[0],v[0]);
					print_float(f);
					break;
				}
			case Ity_F32: 
			case Ity_V32: 
				{
					float f = as_float(v[0]);
					print_float(f);
					break;
				}
			case Ity_I64:
				{
					VG_(printf)("%llx ",as_i64(v[0],v[0]));
					break;
				}
			case Ity_I32:
				{
					VG_(printf)("%x ",as_i32(v[0]));
					break;
				}
			case Ity_I16:
				{
					VG_(printf)("%x ",as_i16(v[0]));
					break;
				}

			default:
			case Ity_INVALID: panic(); break;
		}
		VG_(printf)("(%lx) ",v[0]);
	}
	else
	{
	}
#endif
}

static int double_is_invalid(unsigned long long x)
{
	int hx = (x & 0xffffffff00000000ULL) >> 32;

	// nan or infinite
	return ((hx & 0x7ff00000) == 0x7ff00000);
}

static int float_is_invalid(unsigned int x)
{
	// nan or infinite
	return ((x & 0x7fffffff) >= 0x7f800000);
}

static int convert_to_mpf(mpf_t result, IRType type, unsigned long* v, int* success)
{
	int size = 0;
	print(type,v);
	if (is_64bit)
		switch(type)
		{
			case Ity_F64: 
			case Ity_V64: 
				{
					// if this number is nan, bail out
					if (double_is_invalid(v[0]))
					{
						mpf_init_set_ui(result,1);
						*success = 0;
						return 1;
					}
					mpf_init_set_d (result, as_double(v[0],v[0])); 
					size = 1; 
					break;
				}
			case Ity_F32: 
			case Ity_V32: 
				{
					// if this number is nan, bail out
					if (float_is_invalid(v[0]))
					{
						mpf_init_set_ui(result,1);
						*success = 0;
						return 0;
					}
					mpf_init_set_d (result, as_float(v[0])); 
					size = 1; 
					break;
				}
			case Ity_I64: mpf_init_set_d (result, as_i64(v[0],v[0])   ); size = 1; break;
			case Ity_I32: mpf_init_set_si(result, as_i32(v[0])        ); size = 1; break;
			case Ity_I16: mpf_init_set_si(result, as_i16(v[0])        ); size = 1; break;

			default:
			case Ity_INVALID: panic(); break;
		}
	else
		switch(type)
		{
			case Ity_F64: 
			case Ity_V64: 
				{
					// if this number is nan, bail out
					unsigned long long val = v[0];
					val = (val<<32) | v[1];
					if (double_is_invalid(val))
					{
						mpf_init_set_ui(result,1);
						*success = 0;
						return 1;
					}
					mpf_init_set_d (result, as_double(v[0],v[1])); 
					size = 2; 
					break;
				}
			case Ity_F32: 
			case Ity_V32: 
				{
					// if this number is nan, bail out
					if (float_is_invalid(v[0]))
					{
						mpf_init_set_ui(result,1);
						*success = 0;
						return 0;
					}
					mpf_init_set_d (result, as_float(v[0])); 
					size = 1; 
					break;
				}
			case Ity_I64: mpf_init_set_d (result, as_i64(v[0],v[1])   ); size = 2; break;
			case Ity_I32: mpf_init_set_si(result, as_i32(v[0])        ); size = 1; break;
			case Ity_I16: mpf_init_set_si(result, as_i16(v[0])        ); size = 1; break;

			default:
			case Ity_INVALID: panic(); break;
		}

	return size;
}

typedef struct file_line_error
{
	mpf_t value;
	Char filename [FILE_LEN];
	Int line;
	Bool invalid_result;
	Bool invalid_param;
	struct file_line_error* next;
} file_line_error;

static file_line_error* error_list;

// the tracing callback
static VG_REGPARM(0) void trace_op(IROp op, file_line_error* err, unsigned long arg1, unsigned long arg2, unsigned long arg3, unsigned long arg4, unsigned long arg5, unsigned long arg6, unsigned long arg7, unsigned long arg8, unsigned long arg9)
{
	int i,j;
	unsigned long args[9] = {arg1,arg2,arg3,arg4,arg5,arg6,arg7,arg8,arg9};
	mpf_t args_gmp[4];
	mpf_t result;
	mpf_t fp_result;
	mpf_t errorval;
	int num_param = 0;
	int arg_index = 0;
	int size = 0;
	int ok = 1;

	// look for our op
	for(i=0;i<sizeof(fp_instr_list)/sizeof(*fp_instr_list);i++)
		if (fp_instr_list[i].op == op)
			break;

	if (i == sizeof(fp_instr_list)/sizeof(*fp_instr_list))
		panic();

	// count our params
	for(j=0;j<4;j++)
		if (fp_instr_list[i].t[j] != Ity_INVALID)
			num_param++;

#ifdef DEBUG
	if (err)
		VG_(printf)("%s:%d ",err->filename,err->line);
	VG_(printf)("%d-op callback op 0x%x %s params ",num_param,op,fp_instr_list[i].name);
#endif

	// convert our args into gmp format
	for(j=0;j<num_param;j++)
	{
		size = convert_to_mpf(args_gmp[j], fp_instr_list[i].t[j], &args[arg_index],&ok);
		arg_index +=size;
	}

	if (!ok)
	{
		if (err)
			err->invalid_param = True;
#ifdef DEBUG
		VG_(printf)("invalid params ");
#endif
		goto fail_invalid_param;
	}

#ifdef DEBUG
	VG_(printf)("result ");
#endif

	// convert the result
	size = convert_to_mpf(fp_result, fp_instr_list[i].ret, &args[arg_index],&ok);
	if (!ok)
	{
		if (err)
			err->invalid_result = True;
#ifdef DEBUG
		VG_(printf)("invalid result ");
#endif
		goto fail_invalid_result;
	}

	mpf_init(result);
	mpf_init(errorval);

	// execute the gmp equivalent
	if (fp_instr_list[i].gmp_exec)
	{
		mpf_t abs_errorval;

		mpf_init(abs_errorval);
		mpf_init(result);
		switch(num_param)
		{
			case 1:((GMP_1FUNC)fp_instr_list[i].gmp_exec)(result, args_gmp[0]);break;
			case 2:((GMP_2FUNC)fp_instr_list[i].gmp_exec)(result, args_gmp[0], args_gmp[1]);break;
			case 3:((GMP_3FUNC)fp_instr_list[i].gmp_exec)(result, args_gmp[0], args_gmp[1], args_gmp[2]);break;
			case 4:((GMP_4FUNC)fp_instr_list[i].gmp_exec)(result, args_gmp[0], args_gmp[1], args_gmp[2], args_gmp[3]);break;
			default:panic();break;
		}
		// compute and accumulate the difference
		mpf_sub(errorval, result, fp_result);
		mpf_abs(abs_errorval, errorval);
		if (err)
			mpf_add(err->value, err->value, abs_errorval);
#ifdef DEBUG
		signed long int exponent;
		float derror;

		derror = mpf_get_d_2exp(&exponent, errorval);
		VG_(printf)("  error ");
		print_double(derror);
		VG_(printf)("* 2^%ld",exponent);
#endif
	}


	mpf_clear(result);
	mpf_clear(errorval);

fail_invalid_result:
	mpf_clear(fp_result);
fail_invalid_param:
	for(j=0;j<num_param;j++)
		mpf_clear(args_gmp[j]);
#ifdef DEBUG
	VG_(printf)("\n");
#endif
}

// adds statements to readback and convert different data types
static int add_callback_read(IRExpr** result, IRSB* bb, IRType type, IRExpr* arg)
{
	int size = 0;
	if (is_64bit)
		switch(type)
		{
			case Ity_F64:
				{
					IRTemp t1 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I64);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_ReinterpF64asI64, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t1);
					size = 1;
					break;
				}
			case Ity_F32:
				{
					IRExpr* i32_data;
					IRTemp t1 = IRTemp_INVALID;
					IRTemp t2 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I32);
					t2 = newIRTemp(bb->tyenv, Ity_I64);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_ReinterpF32asI32, arg));
					addStmtToIRSB(bb, inst);
					i32_data = IRExpr_RdTmp(t1);
					inst = IRStmt_WrTmp(t2, IRExpr_Unop(Iop_32Uto64, i32_data));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t2);
					size = 1;
					break;
				}
			case Ity_V64:
			case Ity_V32:
				{
					IRTemp t1 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I64);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_V128to64, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t1);
					size = 1;
					break;
				}
			case Ity_I64:
				{
					IRTemp t = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t = newIRTemp(bb->tyenv, Ity_I64);
					inst = IRStmt_WrTmp(t, arg);
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t);
					size = 1;
					break;
				}
			case Ity_I32:
				{
					IRTemp t = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t = newIRTemp(bb->tyenv, Ity_I64);
					inst = IRStmt_WrTmp(t, IRExpr_Unop(Iop_32Uto64, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t);
					size = 1;
					break;
				}
			case Ity_I16:
				{
					IRTemp t = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t = newIRTemp(bb->tyenv, Ity_I64);
					inst = IRStmt_WrTmp(t, IRExpr_Unop(Iop_16Uto64, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t);
					size = 1;
					break;
				}
			default:
			case Ity_INVALID: panic(); break;
		}
	else // 32 bit
		switch(type)
		{
			case Ity_F64:
				{
					IRExpr* float_data;
					IRTemp t1 = IRTemp_INVALID;
					IRTemp t2 = IRTemp_INVALID;
					IRTemp t3 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I64);
					t2 = newIRTemp(bb->tyenv, Ity_I32);
					t3 = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_ReinterpF64asI64, arg));
					addStmtToIRSB(bb, inst);
					float_data = IRExpr_RdTmp(t1);
					inst = IRStmt_WrTmp(t2, IRExpr_Unop(Iop_64HIto32, float_data));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t2);
					inst = IRStmt_WrTmp(t3, IRExpr_Unop(Iop_64to32, float_data));
					addStmtToIRSB(bb, inst);
					result[1] = IRExpr_RdTmp(t3);
					size = 2;
					break;
				}
			case Ity_F32:
				{
					IRTemp t1 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_ReinterpF32asI32, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t1);
					size = 1;
					break;
				}
			case Ity_V64:
				{
					IRExpr* float_data;
					IRTemp t1 = IRTemp_INVALID;
					IRTemp t2 = IRTemp_INVALID;
					IRTemp t3 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I64);
					t2 = newIRTemp(bb->tyenv, Ity_I32);
					t3 = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_V128to64, arg));
					addStmtToIRSB(bb, inst);
					float_data = IRExpr_RdTmp(t1);
					inst = IRStmt_WrTmp(t2, IRExpr_Unop(Iop_64HIto32, float_data));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t2);
					inst = IRStmt_WrTmp(t3, IRExpr_Unop(Iop_64to32, float_data));
					addStmtToIRSB(bb, inst);
					result[1] = IRExpr_RdTmp(t3);
					size = 2;
					break;
				}
			case Ity_V32:
				{
					IRTemp t1 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t1 = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t1, IRExpr_Unop(Iop_V128to32, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t1);
					size = 1;
					break;
				}
			case Ity_I64:
				{
					IRTemp t2 = IRTemp_INVALID;
					IRTemp t3 = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t2 = newIRTemp(bb->tyenv, Ity_I32);
					t3 = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t2, IRExpr_Unop(Iop_64HIto32, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t2);
					inst = IRStmt_WrTmp(t3, IRExpr_Unop(Iop_64to32, arg));
					addStmtToIRSB(bb, inst);
					result[1] = IRExpr_RdTmp(t3);
					size = 2;
					break;
				}
			case Ity_I32:
				{
					IRTemp t = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t, arg);
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t);
					size = 1;
					break;
				}
			case Ity_I16:
				{
					IRTemp t = IRTemp_INVALID;
					IRStmt* inst = NULL;

					t = newIRTemp(bb->tyenv, Ity_I32);
					inst = IRStmt_WrTmp(t, IRExpr_Unop(Iop_16Uto32, arg));
					addStmtToIRSB(bb, inst);
					result[0] = IRExpr_RdTmp(t);
					size = 1;
					break;
				}
			default:
			case Ity_INVALID: panic(); break;
		}

	return size;
}

static int instrument_inst(IRStmt* st, IRSB* bb, IRExpr* e)
{
	IROp op; 
	IRExpr* arg[4];
	IRExpr* data[9];
	int i,j;
	int num_args = 0;
	IRExpr** argv;
	IRDirty* di;
	int data_index = 0;
	IRExpr* read_res;

	switch(e->tag)
	{
		case Iex_Unop:
			arg[0] = e->Iex.Unop.arg;
			arg[1] = NULL;
			arg[2] = NULL;
			arg[3] = NULL;
			op = e->Iex.Unop.op;
			num_args = 1;
			break;
		case Iex_Binop:
			arg[0] = e->Iex.Binop.arg1;
			arg[1] = e->Iex.Binop.arg2;
			arg[2] = NULL;
			arg[3] = NULL;
			op = e->Iex.Binop.op;
			num_args = 2;
			break;
		case Iex_Triop:
			arg[0] = e->Iex.Triop.arg1;
			arg[1] = e->Iex.Triop.arg2;
			arg[2] = e->Iex.Triop.arg3;
			arg[3] = NULL;
			op = e->Iex.Triop.op;
			num_args = 3;
			break;
		case Iex_Qop:
			arg[0] = e->Iex.Qop.arg1;
			arg[1] = e->Iex.Qop.arg2;
			arg[2] = e->Iex.Qop.arg3;
			arg[3] = e->Iex.Qop.arg4;
			op = e->Iex.Qop.op;
			num_args = 4;
			break;
		default:
			return 0;
	}

	for(j=0;j<9;j++)
		data[j] = NULL;

	for(i=0;i<sizeof(fp_instr_list)/sizeof(*fp_instr_list);i++)
		if (fp_instr_list[i].op == op)
			break;

	if (i==sizeof(fp_instr_list)/sizeof(*fp_instr_list))
		return 0;

	// grab all the args
	for(j=0;j<num_args;j++)
	{
		data_index += add_callback_read(&data[data_index], bb, fp_instr_list[i].t[j], arg[j]);
	}

	// grab the result
	read_res = IRExpr_RdTmp(st->Ist.WrTmp.tmp);
	data_index += add_callback_read(&data[data_index], bb, fp_instr_list[i].ret, read_res);

	// add the callback
	if (is_64bit)
	{
		if (data_index>6) panic();
		argv = mkIRExprVec_7( mkIRExpr_HWord(op),
				mkIRExpr_HWord(error_list),
				data[0], 
				data[1], 
				data[2], 
				data[3], 
				data[4] 
				);
	}
	else
	{
		if (data_index>9) panic();
		argv = mkIRExprVec_11( mkIRExpr_HWord(op),
				mkIRExpr_HWord(error_list),
				data[0], 
				data[1], 
				data[2], 
				data[3], 
				data[4], 
				data[5], 
				data[6], 
				data[7], 
				data[8]
				);
	}
	di = unsafeIRDirty_0_N(0,
			"trace_op",
			VG_(fnptr_to_fnentry)(trace_op),
			argv);
	addStmtToIRSB(bb, IRStmt_Dirty(di));

	return 1;
}

static void ft_post_clo_init(void)
{
	mpf_set_default_prec(1024);
	error_list = NULL;
}

static
IRSB* ft_instrument ( VgCallbackClosure* closure,
                      IRSB* bb_in,
                      VexGuestLayout* layout, 
                      VexGuestExtents* vge,
                      IRType gWordTy, IRType hWordTy )
{
	IRSB* bb_out;
	int i;

	if (gWordTy == Ity_I64)
		is_64bit = 1;
	else if (gWordTy == Ity_I32)
		is_64bit = 0;
	else
		panic();

	bb_out = deepCopyIRSBExceptStmts(bb_in);

	for (i = 0; i < bb_in->stmts_used; i++) {
		IRStmt* st = bb_in->stmts[i];
		IRExpr* e;

		addStmtToIRSB(bb_out, st);
		if (st->tag == Ist_WrTmp )
		{
			e = st->Ist.WrTmp.data;
			instrument_inst(st,bb_out,e);
		}
		else if (st->tag == Ist_IMark)
		{
			Bool found_dirname;
			Bool found_file_line;
			char file[FILE_LEN];
			char dir[FILE_LEN];
			file_line_error* er;
			er = (file_line_error*)VG_(malloc)("err",sizeof(*er));
			er->next = error_list;
			found_file_line = VG_(get_filename_linenum)(
					st->Ist.IMark.addr, 
					file, FILE_LEN,
					dir,  FILE_LEN, &found_dirname,
					&er->line
					);
			tl_assert(VG_(strlen)(dir) + VG_(strlen)(file) + 1 < FILE_LEN);
			VG_(strcat)(dir, "/");
			VG_(strcat)(dir, file);
			VG_(strcpy)(er->filename, dir);

			if (  ( error_list && (error_list->line == er->line) && (VG_(strncmp)(error_list->filename,er->filename,sizeof(er->filename)) == 0) ) ||
			      (!found_dirname) || (!found_file_line) )
			{
				VG_(free)(er);
			}
			else
			{
				mpf_init_set_ui(er->value, 0);
				er->invalid_result = False;
				er->invalid_param = False;
				error_list = er;
			}
		}
	}

	return bb_out;
}

static void ft_fini(Int exitcode)
{
	file_line_error* err = error_list;
	file_line_error* p;
	while(err)
	{
		signed long int exponent;
		double derror = mpf_get_d_2exp(&exponent, err->value);
		if ( (derror != 0.0) || (err->invalid_result) || (err->invalid_param) )
		{
			VG_(printf)("%s:%d ",err->filename,err->line);
			if (err->invalid_result)
			{
				VG_(printf)("invalid result ");
			}
			if (err->invalid_param)
			{
				VG_(printf)("invalid parameter ");
			}
			if (derror != 0.0)
			{
				VG_(printf)("error ");
				print_double(derror);
				VG_(printf)("* 2^%ld",exponent);
			}
			VG_(printf)("\n");
		}
		mpf_clear(err->value);
		p = err;
		err = err->next;
		VG_(free)(p);
	}
}

static void ft_pre_clo_init(void)
{
   VG_(details_name)            ("Fatgrind");
   VG_(details_version)         (NULL);
   VG_(details_description)     ("a floating point accuracy tracker");
   VG_(details_copyright_author)(
      "Copyright (C) 2009, and GNU GPL'd, by Stephane Marchesin.");
   VG_(details_bug_reports_to)  (VG_BUGS_TO);

   VG_(basic_tool_funcs)        (ft_post_clo_init,
                                 ft_instrument,
                                 ft_fini);

   /* No needs, no core events to track */
}

VG_DETERMINE_INTERFACE_VERSION(ft_pre_clo_init)