poly_tan.c

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