ncx.c

一个用来实现偏微分方程中网格的计算库

	ix_int xx = (ix_int)(*ip);
	put_ix_int(xp, &xx);
#   if IX_INT_MAX < SHORT_MAX
	if(*ip > X_INT_MAX || *ip < X_INT_MIN)
		return NC_ERANGE;
#   endif
	return ENOERR;
#endif
}

int
ncx_put_int_int(void *xp, const int *ip)
{
#if SIZEOF_IX_INT == SIZEOF_INT && IX_INT_MAX == INT_MAX
	put_ix_int(xp, (ix_int *)ip);
	return ENOERR;
#else
	ix_int xx = (ix_int)(*ip);
	put_ix_int(xp, &xx);
#   if IX_INT_MAX < INT_MAX
	if(*ip > X_INT_MAX || *ip < X_INT_MIN)
		return NC_ERANGE;
#   endif
	return ENOERR;
#endif
}

int
ncx_put_int_long(void *xp, const long *ip)
{
#if SIZEOF_IX_INT == SIZEOF_LONG && IX_INT_MAX == LONG_MAX
	put_ix_int(xp, (ix_int *)ip);
	return ENOERR;
#else
	ix_int xx = (ix_int)(*ip);
	put_ix_int(xp, &xx);
#   if IX_INT_MAX < LONG_MAX
	if(*ip > X_INT_MAX || *ip < X_INT_MIN)
		return NC_ERANGE;
#   endif
	return ENOERR;
#endif
}

int
ncx_put_int_float(void *xp, const float *ip)
{
	ix_int xx = (ix_int)(*ip);
	put_ix_int(xp, &xx);
	if(*ip > (double)X_INT_MAX || *ip < (double)X_INT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_put_int_double(void *xp, const double *ip)
{
	ix_int xx = (ix_int)(*ip);
	put_ix_int(xp, &xx);
	if(*ip > X_INT_MAX || *ip < X_INT_MIN)
		return NC_ERANGE;
	return ENOERR;
}
 

/* x_float */

#if X_SIZEOF_FLOAT == SIZEOF_FLOAT && !defined(NO_IEEE_FLOAT)

static void
get_ix_float(const void *xp, float *ip)
{
#ifdef WORDS_BIGENDIAN
	(void) memcpy(ip, xp, sizeof(float));
#else
	swap4b(ip, xp);
#endif
}

static void
put_ix_float(void *xp, const float *ip)
{
#ifdef WORDS_BIGENDIAN
	(void) memcpy(xp, ip, X_SIZEOF_FLOAT);
#else
	swap4b(xp, ip);
#endif
}

#elif vax

/* What IEEE single precision floating point looks like on a Vax */
struct	ieee_single {
	unsigned int	exp_hi       : 7;
	unsigned int	sign         : 1;
	unsigned int 	mant_hi      : 7;
	unsigned int	exp_lo       : 1;
	unsigned int	mant_lo_hi   : 8;
	unsigned int	mant_lo_lo   : 8;
};

/* Vax single precision floating point */
struct	vax_single {
	unsigned int	mantissa1 : 7;
	unsigned int	exp       : 8;
	unsigned int	sign      : 1;
	unsigned int	mantissa2 : 16;
};

#define VAX_SNG_BIAS	0x81
#define IEEE_SNG_BIAS	0x7f

static struct sgl_limits {
	struct vax_single s;
	struct ieee_single ieee;
} max = {
	{ 0x7f, 0xff, 0x0, 0xffff },	/* Max Vax */
	{ 0x7f, 0x0, 0x0, 0x1, 0x0, 0x0 }		/* Max IEEE */
};
static struct sgl_limits min = {
	{ 0x0, 0x0, 0x0, 0x0 },	/* Min Vax */
	{ 0x0, 0x0, 0x0, 0x0, 0x0, 0x0 }		/* Min IEEE */
};

static void
get_ix_float(const void *xp, float *ip)
{
		struct vax_single *const vsp = (struct vax_single *) ip;
		const struct ieee_single *const isp =
			 (const struct ieee_single *) xp;
		unsigned exp = isp->exp_hi << 1 | isp->exp_lo;

		switch(exp) {
		case 0 :
			/* ieee subnormal */
			if(isp->mant_hi == min.ieee.mant_hi
				&& isp->mant_lo_hi == min.ieee.mant_lo_hi
				&& isp->mant_lo_lo == min.ieee.mant_lo_lo)
			{
				*vsp = min.s;
			}
			else
			{
				unsigned mantissa = (isp->mant_hi << 16)
					 | isp->mant_lo_hi << 8
					 | isp->mant_lo_lo;
				unsigned tmp = mantissa >> 20;
				if(tmp >= 4) {
					vsp->exp = 2;
				} else if (tmp >= 2) {
					vsp->exp = 1;
				} else {
					*vsp = min.s;
					break;
				} /* else */
				tmp = mantissa - (1 << (20 + vsp->exp ));
				tmp <<= 3 - vsp->exp;
				vsp->mantissa2 = tmp;
				vsp->mantissa1 = (tmp >> 16);
			}
			break;
		case 0xfe :
		case 0xff :
			*vsp = max.s;
			break;
		default :
			vsp->exp = exp - IEEE_SNG_BIAS + VAX_SNG_BIAS;
			vsp->mantissa2 = isp->mant_lo_hi << 8 | isp->mant_lo_lo;
			vsp->mantissa1 = isp->mant_hi;
		}

		vsp->sign = isp->sign;

}


static void
put_ix_float(void *xp, const float *ip)
{
		const struct vax_single *const vsp =
			 (const struct vax_single *)ip;
		struct ieee_single *const isp = (struct ieee_single *) xp;

		switch(vsp->exp){
		case 0 :
			/* all vax float with zero exponent map to zero */
			*isp = min.ieee;
			break;
		case 2 :
		case 1 :
		{
			/* These will map to subnormals */
			unsigned mantissa = (vsp->mantissa1 << 16)
					 | vsp->mantissa2;
			mantissa >>= 3 - vsp->exp;
			mantissa += (1 << (20 + vsp->exp));
			isp->mant_lo_lo = mantissa;
			isp->mant_lo_hi = mantissa >> 8;
			isp->mant_hi = mantissa >> 16;
			isp->exp_lo = 0;
			isp->exp_hi = 0;
		}
			break;
		case 0xff : /* max.s.exp */
			if( vsp->mantissa2 == max.s.mantissa2
				&& vsp->mantissa1 == max.s.mantissa1)
			{
				/* map largest vax float to ieee infinity */
				*isp = max.ieee;
				break;
			} /* else, fall thru */
		default :
		{
			unsigned exp = vsp->exp - VAX_SNG_BIAS + IEEE_SNG_BIAS;
			isp->exp_hi = exp >> 1;
			isp->exp_lo = exp;
			isp->mant_lo_lo = vsp->mantissa2;
			isp->mant_lo_hi = vsp->mantissa2 >> 8;
			isp->mant_hi = vsp->mantissa1;
		}
		}

		isp->sign = vsp->sign;

}

	/* vax */
#elif defined(_CRAY) && !defined(__crayx1)

/*
 * Return the number of bytes until the next "word" boundary
 * N.B. This is based on the very wierd YMP address structure,
 * which puts the address within a word in the leftmost 3 bits
 * of the address.
 */
static size_t
word_align(const void *vp)
{
	const size_t rem = ((size_t)vp >> (64 - 3)) & 0x7;
	return (rem != 0);
}

struct ieee_single_hi {
	unsigned int	sign	: 1;
	unsigned int	 exp	: 8;
	unsigned int	mant	:23;
	unsigned int	pad	:32;
};
typedef struct ieee_single_hi ieee_single_hi;

struct ieee_single_lo {
	unsigned int	pad	:32;
	unsigned int	sign	: 1;
	unsigned int	 exp	: 8;
	unsigned int	mant	:23;
};
typedef struct ieee_single_lo ieee_single_lo;

static const int ieee_single_bias = 0x7f;

struct ieee_double {
	unsigned int	sign	: 1;
	unsigned int	 exp	:11;
	unsigned int	mant	:52;
};
typedef struct ieee_double ieee_double;

static const int ieee_double_bias = 0x3ff;

#if defined(NO_IEEE_FLOAT)

struct cray_single {
	unsigned int	sign	: 1;
	unsigned int	 exp	:15;
	unsigned int	mant	:48;
};
typedef struct cray_single cray_single;

static const int cs_ieis_bias = 0x4000 - 0x7f;

static const int cs_id_bias = 0x4000 - 0x3ff;


static void
get_ix_float(const void *xp, float *ip)
{

	if(word_align(xp) == 0)
	{
		const ieee_single_hi *isp = (const ieee_single_hi *) xp;
		cray_single *csp = (cray_single *) ip;

		if(isp->exp == 0)
		{
			/* ieee subnormal */
			*ip = (double)isp->mant;
			if(isp->mant != 0)
			{
				csp->exp -= (ieee_single_bias + 22);
			}
		}
		else
		{
			csp->exp  = isp->exp + cs_ieis_bias + 1;
			csp->mant = isp->mant << (48 - 1 - 23);
			csp->mant |= (1 << (48 - 1));
		}
		csp->sign = isp->sign;


	}
	else
	{
		const ieee_single_lo *isp = (const ieee_single_lo *) xp;
		cray_single *csp = (cray_single *) ip;

		if(isp->exp == 0)
		{
			/* ieee subnormal */
			*ip = (double)isp->mant;
			if(isp->mant != 0)
			{
				csp->exp -= (ieee_single_bias + 22);
			}
		}
		else
		{
			csp->exp  = isp->exp + cs_ieis_bias + 1;
			csp->mant = isp->mant << (48 - 1 - 23);
			csp->mant |= (1 << (48 - 1));
		}
		csp->sign = isp->sign;


	}
}

static void
put_ix_float(void *xp, const float *ip)
{
	if(word_align(xp) == 0)
	{
		ieee_single_hi *isp = (ieee_single_hi*)xp;
	const cray_single *csp = (const cray_single *) ip;
	int ieee_exp = csp->exp - cs_ieis_bias -1;

	isp->sign = csp->sign;

	if(ieee_exp >= 0xff)
	{
		/* NC_ERANGE => ieee Inf */
		isp->exp = 0xff;
		isp->mant = 0x0;
	}
	else if(ieee_exp > 0)
	{
		/* normal ieee representation */
		isp->exp  = ieee_exp;
		/* assumes cray rep is in normal form */
		assert(csp->mant & 0x800000000000);
		isp->mant = (((csp->mant << 1) &
				0xffffffffffff) >> (48 - 23));
	}
	else if(ieee_exp > -23)
	{
		/* ieee subnormal, right  */
		const int rshift = (48 - 23 - ieee_exp);

		isp->mant = csp->mant >> rshift;

#if 0
		if(csp->mant & (1 << (rshift -1)))
		{
			/* round up */
			isp->mant++;
		}
#endif

		isp->exp  = 0;
	}
	else
	{
		/* smaller than ieee can represent */
		isp->exp = 0;
		isp->mant = 0;
	}

	}
	else
	{
		ieee_single_lo *isp = (ieee_single_lo*)xp;
	const cray_single *csp = (const cray_single *) ip;
	int ieee_exp = csp->exp - cs_ieis_bias -1;

	isp->sign = csp->sign;

	if(ieee_exp >= 0xff)
	{
		/* NC_ERANGE => ieee Inf */
		isp->exp = 0xff;
		isp->mant = 0x0;
	}
	else if(ieee_exp > 0)
	{
		/* normal ieee representation */
		isp->exp  = ieee_exp;
		/* assumes cray rep is in normal form */
		assert(csp->mant & 0x800000000000);
		isp->mant = (((csp->mant << 1) &
				0xffffffffffff) >> (48 - 23));
	}
	else if(ieee_exp > -23)
	{
		/* ieee subnormal, right  */
		const int rshift = (48 - 23 - ieee_exp);

		isp->mant = csp->mant >> rshift;

#if 0
		if(csp->mant & (1 << (rshift -1)))
		{
			/* round up */
			isp->mant++;
		}
#endif

		isp->exp  = 0;
	}
	else
	{
		/* smaller than ieee can represent */
		isp->exp = 0;
		isp->mant = 0;
	}

	}
}

#else
	/* IEEE Cray with only doubles */
static void
get_ix_float(const void *xp, float *ip)
{

	ieee_double *idp = (ieee_double *) ip;

	if(word_align(xp) == 0)
	{
		const ieee_single_hi *isp = (const ieee_single_hi *) xp;
		if(isp->exp == 0 && isp->mant == 0)
		{
			idp->exp = 0;
			idp->mant = 0;
		}
		else
		{
			idp->exp = isp->exp + (ieee_double_bias - ieee_single_bias);
			idp->mant = isp->mant << (52 - 23);
		}
		idp->sign = isp->sign;
	}
	else
	{
		const ieee_single_lo *isp = (const ieee_single_lo *) xp;
		if(isp->exp == 0 && isp->mant == 0)
		{
			idp->exp = 0;
			idp->mant = 0;
		}
		else
		{
			idp->exp = isp->exp + (ieee_double_bias - ieee_single_bias);
			idp->mant = isp->mant << (52 - 23);
		}
		idp->sign = isp->sign;
	}
}

static void
put_ix_float(void *xp, const float *ip)
{
	const ieee_double *idp = (const ieee_double *) ip;
	if(word_align(xp) == 0)
	{
		ieee_single_hi *isp = (ieee_single_hi*)xp;
		if(idp->exp > (ieee_double_bias - ieee_single_bias))
			isp->exp = idp->exp - (ieee_double_bias - ieee_single_bias);
		else
			isp->exp = 0;
		isp->mant = idp->mant >> (52 - 23);
		isp->sign = idp->sign;
	}
	else
	{
		ieee_single_lo *isp = (ieee_single_lo*)xp;
		if(idp->exp > (ieee_double_bias - ieee_single_bias))
			isp->exp = idp->exp - (ieee_double_bias - ieee_single_bias);
		else
			isp->exp = 0;
		isp->mant = idp->mant >> (52 - 23);
		isp->sign = idp->sign;
	}
}
#endif

#elif _SX && _FLOAT2
static void
get_ix_float(const void *xp, float *ip)
{
	const int ncnv = ie3_fl2(xp, ip, 4, 8, 1);
}

static void
put_ix_float(void *xp, const float *ip)
{
	const int ncnv = fl2_ie3(ip, xp, 8, 4, 1);
}
#else
#error "ix_float implementation"
#endif


int
ncx_get_float_schar(const void *xp, schar *ip)
{
	float xx;
	get_ix_float(xp, &xx);
	*ip = (schar) xx;
	if(xx > SCHAR_MAX || xx < SCHAR_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_float_uchar(const void *xp, uchar *ip)
{
	float xx;
	get_ix_float(xp, &xx);
	*ip = (uchar) xx;
	if(xx > UCHAR_MAX || xx < 0)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_float_short(const void *xp, short *ip)
{
	float xx;
	get_ix_float(xp, &xx);
	*ip = (short) xx;
	if(xx > SHORT_MAX || xx < SHORT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_float_int(const void *xp, int *ip)
{
	float xx;
	get_ix_float(xp, &xx);
	*ip = (int) xx;
	if(xx > (double)INT_MAX || xx < (double)INT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_float_long(const void *xp, long *ip)
{
	float xx;
	get_ix_float(xp, &xx);
	*ip = (long) xx;
	if(xx > LONG_MAX || xx < LONG_MIN)
		return NC_ERANGE;
	return ENOERR;
}

int
ncx_get_float_float(const void *xp, float *ip)
{
	/* TODO */
	get_ix_float(xp, ip);
	return ENOERR;
}

int
ncx_get_float_double(const void *xp, double *ip)
{
	/* TODO */
	float xx;
	get_ix_float(xp, &xx);
	*ip = xx;
	return ENOERR;
}


int
ncx_put_float_schar(void *xp, const schar *ip)
{
	float xx = (float) *ip;
	put_ix_float(xp, &xx);
	return ENOERR;
}

int
ncx_put_float_uchar(void *xp, const uchar *ip)
{
	float xx = (float) *ip;
	put_ix_float(xp, &xx);
	return ENOERR;
}

int
ncx_put_float_short(void *xp, const short *ip)
{
	float xx = (float) *ip;
	put_ix_float(xp, &xx);
#if 0	/* TODO: figure this out */
	if((float)(*ip) > X_FLOAT_MAX || (float)(*ip) < X_FLOAT_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_float_int(void *xp, const int *ip)
{
	float xx = (float) *ip;
	put_ix_float(xp, &xx);
#if 1	/* TODO: figure this out */
	if((float)(*ip) > X_FLOAT_MAX || (float)(*ip) < X_FLOAT_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_float_long(void *xp, const long *ip)
{
	float xx = (float) *ip;
	put_ix_float(xp, &xx);
#if 1	/* TODO: figure this out */
	if((float)(*ip) > X_FLOAT_MAX || (float)(*ip) < X_FLOAT_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_float_float(void *xp, const float *ip)
{
	put_ix_float(xp, ip);
#ifdef NO_IEEE_FLOAT
	if(*ip > X_FLOAT_MAX || *ip < X_FLOAT_MIN)
		return NC_ERANGE;
#endif
	return ENOERR;
}

int
ncx_put_float_double(void *xp, const double *ip)
{
	float xx = (float) *ip;
	put_ix_float(xp, &xx);
	if(*ip > X_FLOAT_MAX || *ip < X_FLOAT_MIN)
		return NC_ERANGE;
	return ENOERR;
}

/* x_double */

#if X_SIZEOF_DOUBLE == SIZEOF_DOUBLE  && !defined(NO_IEEE_FLOAT)

static void
get_ix_double(const void *xp, double *ip)
{
#ifdef WORDS_BIGENDIAN
	(void) memcpy(ip, xp, sizeof(double));
#else
	swap8b(ip, xp);
#endif
}

static void
put_ix_double(void *xp, const double *ip)
{
#ifdef WORDS_BIGENDIAN
	(void) memcpy(xp, ip, X_SIZEOF_DOUBLE);
#else
	swap8b(xp, ip);
#endif
}

#elif vax