[GitHub/mt8127/android_kernel_alcatel_ttab.git] / arch / x86 / math-emu / poly_tan.c

/*---------------------------------------------------------------------------+
 |  poly_tan.c                                                               |
 |                                                                           |
 | Compute the tan of a FPU_REG, using a polynomial approximation.           |
 |                                                                           |
 | Copyright (C) 1992,1993,1994,1997,1999                                    |
 |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
 |                       Australia.  E-mail   billm@melbpc.org.au            |
 |                                                                           |
 |                                                                           |
 +---------------------------------------------------------------------------*/

#include "exception.h"
#include "reg_constant.h"
#include "fpu_emu.h"
#include "fpu_system.h"
#include "control_w.h"
#include "poly.h"

#define	HiPOWERop	3	/* odd poly, positive terms */
static const unsigned long long oddplterm[HiPOWERop] = {
	0x0000000000000000LL,
	0x0051a1cf08fca228LL,
	0x0000000071284ff7LL
};

#define	HiPOWERon	2	/* odd poly, negative terms */
static const unsigned long long oddnegterm[HiPOWERon] = {
	0x1291a9a184244e80LL,
	0x0000583245819c21LL
};

#define	HiPOWERep	2	/* even poly, positive terms */
static const unsigned long long evenplterm[HiPOWERep] = {
	0x0e848884b539e888LL,
	0x00003c7f18b887daLL
};

#define	HiPOWERen	2	/* even poly, negative terms */
static const unsigned long long evennegterm[HiPOWERen] = {
	0xf1f0200fd51569ccLL,
	0x003afb46105c4432LL
};

static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;

/*--- poly_tan() ------------------------------------------------------------+
 |                                                                           |
 +---------------------------------------------------------------------------*/
void poly_tan(FPU_REG *st0_ptr)
{
	long int exponent;
	int invert;
	Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
	    argSignif, fix_up;
	unsigned long adj;

	exponent = exponent(st0_ptr);

#ifdef PARANOID
	if (signnegative(st0_ptr)) {	/* Can't hack a number < 0.0 */
		arith_invalid(0);
		return;
	}			/* Need a positive number */
#endif /* PARANOID */

	/* Split the problem into two domains, smaller and larger than pi/4 */
	if ((exponent == 0)
	    || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) {
		/* The argument is greater than (approx) pi/4 */
		invert = 1;
		accum.lsw = 0;
		XSIG_LL(accum) = significand(st0_ptr);

		if (exponent == 0) {
			/* The argument is >= 1.0 */
			/* Put the binary point at the left. */
			XSIG_LL(accum) <<= 1;
		}
		/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
		XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
		/* This is a special case which arises due to rounding. */
		if (XSIG_LL(accum) == 0xffffffffffffffffLL) {
			FPU_settag0(TAG_Valid);
			significand(st0_ptr) = 0x8a51e04daabda360LL;
			setexponent16(st0_ptr,
				      (0x41 + EXTENDED_Ebias) | SIGN_Negative);
			return;
		}

		argSignif.lsw = accum.lsw;
		XSIG_LL(argSignif) = XSIG_LL(accum);
		exponent = -1 + norm_Xsig(&argSignif);
	} else {
		invert = 0;
		argSignif.lsw = 0;
		XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);

		if (exponent < -1) {
			/* shift the argument right by the required places */
			if (FPU_shrx(&XSIG_LL(accum), -1 - exponent) >=
			    0x80000000U)
				XSIG_LL(accum)++;	/* round up */
		}
	}

	XSIG_LL(argSq) = XSIG_LL(accum);
	argSq.lsw = accum.lsw;
	mul_Xsig_Xsig(&argSq, &argSq);
	XSIG_LL(argSqSq) = XSIG_LL(argSq);
	argSqSq.lsw = argSq.lsw;
	mul_Xsig_Xsig(&argSqSq, &argSqSq);

	/* Compute the negative terms for the numerator polynomial */
	accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
	polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm,
			HiPOWERon - 1);
	mul_Xsig_Xsig(&accumulatoro, &argSq);
	negate_Xsig(&accumulatoro);
	/* Add the positive terms */
	polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm,
			HiPOWERop - 1);

	/* Compute the positive terms for the denominator polynomial */
	accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
	polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm,
			HiPOWERep - 1);
	mul_Xsig_Xsig(&accumulatore, &argSq);
	negate_Xsig(&accumulatore);
	/* Add the negative terms */
	polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm,
			HiPOWERen - 1);
	/* Multiply by arg^2 */
	mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
	mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
	/* de-normalize and divide by 2 */
	shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1);
	negate_Xsig(&accumulatore);	/* This does 1 - accumulator */

	/* Now find the ratio. */
	if (accumulatore.msw == 0) {
		/* accumulatoro must contain 1.0 here, (actually, 0) but it
		   really doesn't matter what value we use because it will
		   have negligible effect in later calculations
		 */
		XSIG_LL(accum) = 0x8000000000000000LL;
		accum.lsw = 0;
	} else {
		div_Xsig(&accumulatoro, &accumulatore, &accum);
	}

	/* Multiply by 1/3 * arg^3 */
	mul64_Xsig(&accum, &XSIG_LL(argSignif));
	mul64_Xsig(&accum, &XSIG_LL(argSignif));
	mul64_Xsig(&accum, &XSIG_LL(argSignif));
	mul64_Xsig(&accum, &twothirds);
	shr_Xsig(&accum, -2 * (exponent + 1));

	/* tan(arg) = arg + accum */
	add_two_Xsig(&accum, &argSignif, &exponent);

	if (invert) {
		/* We now have the value of tan(pi_2 - arg) where pi_2 is an
		   approximation for pi/2
		 */
		/* The next step is to fix the answer to compensate for the
		   error due to the approximation used for pi/2
		 */

		/* This is (approx) delta, the error in our approx for pi/2
		   (see above). It has an exponent of -65
		 */
		XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
		fix_up.lsw = 0;

		if (exponent == 0)
			adj = 0xffffffff;	/* We want approx 1.0 here, but
						   this is close enough. */
		else if (exponent > -30) {
			adj = accum.msw >> -(exponent + 1);	/* tan */
			adj = mul_32_32(adj, adj);	/* tan^2 */
		} else
			adj = 0;
		adj = mul_32_32(0x898cc517, adj);	/* delta * tan^2 */

		fix_up.msw += adj;
		if (!(fix_up.msw & 0x80000000)) {	/* did fix_up overflow ? */
			/* Yes, we need to add an msb */
			shr_Xsig(&fix_up, 1);
			fix_up.msw |= 0x80000000;
			shr_Xsig(&fix_up, 64 + exponent);
		} else
			shr_Xsig(&fix_up, 65 + exponent);

		add_two_Xsig(&accum, &fix_up, &exponent);

		/* accum now contains tan(pi/2 - arg).
		   Use tan(arg) = 1.0 / tan(pi/2 - arg)
		 */
		accumulatoro.lsw = accumulatoro.midw = 0;
		accumulatoro.msw = 0x80000000;
		div_Xsig(&accumulatoro, &accum, &accum);
		exponent = -exponent - 1;
	}

	/* Transfer the result */
	round_Xsig(&accum);
	FPU_settag0(TAG_Valid);
	significand(st0_ptr) = XSIG_LL(accum);
	setexponent16(st0_ptr, exponent + EXTENDED_Ebias);	/* Result is positive. */

}
Commit	Line	Data
1da177e4 LT	1	/*---------------------------------------------------------------------------+
	2	\| poly_tan.c \|
	3	\| \|
	4	\| Compute the tan of a FPU_REG, using a polynomial approximation. \|
	5	\| \|
	6	\| Copyright (C) 1992,1993,1994,1997,1999 \|
	7	\| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, \|
	8	\| Australia. E-mail billm@melbpc.org.au \|
	9	\| \|
	10	\| \|
	11	+---------------------------------------------------------------------------*/
	12
	13	#include "exception.h"
	14	#include "reg_constant.h"
	15	#include "fpu_emu.h"
	16	#include "fpu_system.h"
	17	#include "control_w.h"
	18	#include "poly.h"
	19
1da177e4	20	#define HiPOWERop 3 /* odd poly, positive terms */
3d0d14f9 IM	21	static const unsigned long long oddplterm[HiPOWERop] = {
	22	0x0000000000000000LL,
	23	0x0051a1cf08fca228LL,
	24	0x0000000071284ff7LL
1da177e4 LT	25	};
	26
	27	#define HiPOWERon 2 /* odd poly, negative terms */
3d0d14f9 IM	28	static const unsigned long long oddnegterm[HiPOWERon] = {
	29	0x1291a9a184244e80LL,
	30	0x0000583245819c21LL
1da177e4 LT	31	};
	32
	33	#define HiPOWERep 2 /* even poly, positive terms */
3d0d14f9 IM	34	static const unsigned long long evenplterm[HiPOWERep] = {
	35	0x0e848884b539e888LL,
	36	0x00003c7f18b887daLL
1da177e4 LT	37	};
	38
	39	#define HiPOWERen 2 /* even poly, negative terms */
3d0d14f9 IM	40	static const unsigned long long evennegterm[HiPOWERen] = {
	41	0xf1f0200fd51569ccLL,
	42	0x003afb46105c4432LL
1da177e4 LT	43	};
	44
	45	static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
	46
1da177e4 LT	47	/*--- poly_tan() ------------------------------------------------------------+
	48	\| \|
	49	+---------------------------------------------------------------------------*/
e8d591dc	50	void poly_tan(FPU_REG *st0_ptr)
1da177e4	51	{
3d0d14f9 IM	52	long int exponent;
	53	int invert;
	54	Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
	55	argSignif, fix_up;
	56	unsigned long adj;
1da177e4	57
3d0d14f9	58	exponent = exponent(st0_ptr);
1da177e4 LT	59
1da177e4 LT	60	#ifdef PARANOID
3d0d14f9 IM	61	if (signnegative(st0_ptr)) { /* Can't hack a number < 0.0 */
	62	arith_invalid(0);
	63	return;
	64	} /* Need a positive number */
1da177e4 LT	65	#endif /* PARANOID */
1da177e4 LT	66
3d0d14f9 IM	67	/* Split the problem into two domains, smaller and larger than pi/4 */
	68	if ((exponent == 0)
	69	\|\| ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) {
	70	/* The argument is greater than (approx) pi/4 */
	71	invert = 1;
	72	accum.lsw = 0;
	73	XSIG_LL(accum) = significand(st0_ptr);
	74
	75	if (exponent == 0) {
	76	/* The argument is >= 1.0 */
	77	/* Put the binary point at the left. */
	78	XSIG_LL(accum) <<= 1;
	79	}
	80	/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
	81	XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
	82	/* This is a special case which arises due to rounding. */
	83	if (XSIG_LL(accum) == 0xffffffffffffffffLL) {
	84	FPU_settag0(TAG_Valid);
	85	significand(st0_ptr) = 0x8a51e04daabda360LL;
	86	setexponent16(st0_ptr,
	87	(0x41 + EXTENDED_Ebias) \| SIGN_Negative);
	88	return;
	89	}
	90
	91	argSignif.lsw = accum.lsw;
	92	XSIG_LL(argSignif) = XSIG_LL(accum);
	93	exponent = -1 + norm_Xsig(&argSignif);
	94	} else {
	95	invert = 0;
	96	argSignif.lsw = 0;
	97	XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr);
	98
	99	if (exponent < -1) {
	100	/* shift the argument right by the required places */
	101	if (FPU_shrx(&XSIG_LL(accum), -1 - exponent) >=
	102	0x80000000U)
	103	XSIG_LL(accum)++; /* round up */
	104	}
1da177e4 LT	105	}
1da177e4 LT	106
3d0d14f9 IM	107	XSIG_LL(argSq) = XSIG_LL(accum);
	108	argSq.lsw = accum.lsw;
	109	mul_Xsig_Xsig(&argSq, &argSq);
	110	XSIG_LL(argSqSq) = XSIG_LL(argSq);
	111	argSqSq.lsw = argSq.lsw;
	112	mul_Xsig_Xsig(&argSqSq, &argSqSq);
	113
	114	/* Compute the negative terms for the numerator polynomial */
	115	accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
	116	polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm,
	117	HiPOWERon - 1);
	118	mul_Xsig_Xsig(&accumulatoro, &argSq);
	119	negate_Xsig(&accumulatoro);
	120	/* Add the positive terms */
	121	polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm,
	122	HiPOWERop - 1);
	123
	124	/* Compute the positive terms for the denominator polynomial */
	125	accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
	126	polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm,
	127	HiPOWERep - 1);
	128	mul_Xsig_Xsig(&accumulatore, &argSq);
	129	negate_Xsig(&accumulatore);
	130	/* Add the negative terms */
	131	polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm,
	132	HiPOWERen - 1);
	133	/* Multiply by arg^2 */
	134	mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
	135	mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
	136	/* de-normalize and divide by 2 */
	137	shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1);
	138	negate_Xsig(&accumulatore); /* This does 1 - accumulator */
	139
	140	/* Now find the ratio. */
	141	if (accumulatore.msw == 0) {
	142	/* accumulatoro must contain 1.0 here, (actually, 0) but it
	143	really doesn't matter what value we use because it will
	144	have negligible effect in later calculations
	145	*/
	146	XSIG_LL(accum) = 0x8000000000000000LL;
	147	accum.lsw = 0;
	148	} else {
	149	div_Xsig(&accumulatoro, &accumulatore, &accum);
1da177e4	150	}
3d0d14f9 IM	151
	152	/* Multiply by 1/3 * arg^3 */
	153	mul64_Xsig(&accum, &XSIG_LL(argSignif));
	154	mul64_Xsig(&accum, &XSIG_LL(argSignif));
	155	mul64_Xsig(&accum, &XSIG_LL(argSignif));
	156	mul64_Xsig(&accum, &twothirds);
	157	shr_Xsig(&accum, -2 * (exponent + 1));
	158
	159	/* tan(arg) = arg + accum */
	160	add_two_Xsig(&accum, &argSignif, &exponent);
	161
	162	if (invert) {
	163	/* We now have the value of tan(pi_2 - arg) where pi_2 is an
	164	approximation for pi/2
	165	*/
	166	/* The next step is to fix the answer to compensate for the
	167	error due to the approximation used for pi/2
	168	*/
	169
	170	/* This is (approx) delta, the error in our approx for pi/2
	171	(see above). It has an exponent of -65
	172	*/
	173	XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
	174	fix_up.lsw = 0;
	175
	176	if (exponent == 0)
	177	adj = 0xffffffff; /* We want approx 1.0 here, but
	178	this is close enough. */
	179	else if (exponent > -30) {
	180	adj = accum.msw >> -(exponent + 1); /* tan */
	181	adj = mul_32_32(adj, adj); /* tan^2 */
	182	} else
	183	adj = 0;
	184	adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */
	185
	186	fix_up.msw += adj;
	187	if (!(fix_up.msw & 0x80000000)) { /* did fix_up overflow ? */
	188	/* Yes, we need to add an msb */
	189	shr_Xsig(&fix_up, 1);
	190	fix_up.msw \|= 0x80000000;
	191	shr_Xsig(&fix_up, 64 + exponent);
	192	} else
	193	shr_Xsig(&fix_up, 65 + exponent);
	194
	195	add_two_Xsig(&accum, &fix_up, &exponent);
	196
	197	/* accum now contains tan(pi/2 - arg).
	198	Use tan(arg) = 1.0 / tan(pi/2 - arg)
	199	*/
	200	accumulatoro.lsw = accumulatoro.midw = 0;
	201	accumulatoro.msw = 0x80000000;
	202	div_Xsig(&accumulatoro, &accum, &accum);
	203	exponent = -exponent - 1;
1da177e4	204	}
3d0d14f9 IM	205
	206	/* Transfer the result */
	207	round_Xsig(&accum);
	208	FPU_settag0(TAG_Valid);
	209	significand(st0_ptr) = XSIG_LL(accum);
	210	setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */
1da177e4 LT	211
1da177e4 LT	212	}