/arch/x86/math-emu/poly_l2.c

  1/*---------------------------------------------------------------------------+
  2 |  poly_l2.c                                                                |
  3 |                                                                           |
  4 | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
  5 |                                                                           |
  6 | Copyright (C) 1992,1993,1994,1997                                         |
  7 |                  W. Metzenthen, 22 Parker St, Ormond, Vic 3163, Australia |
  8 |                  E-mail   billm@suburbia.net                              |
  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
 20static void log2_kernel(FPU_REG const *arg, u_char argsign,
 21			Xsig * accum_result, long int *expon);
 22
 23/*--- poly_l2() -------------------------------------------------------------+
 24 |   Base 2 logarithm by a polynomial approximation.                         |
 25 +---------------------------------------------------------------------------*/
 26void poly_l2(FPU_REG *st0_ptr, FPU_REG *st1_ptr, u_char st1_sign)
 27{
 28	long int exponent, expon, expon_expon;
 29	Xsig accumulator, expon_accum, yaccum;
 30	u_char sign, argsign;
 31	FPU_REG x;
 32	int tag;
 33
 34	exponent = exponent16(st0_ptr);
 35
 36	/* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) */
 37	if (st0_ptr->sigh > (unsigned)0xb504f334) {
 38		/* Treat as  sqrt(2)/2 < st0_ptr < 1 */
 39		significand(&x) = -significand(st0_ptr);
 40		setexponent16(&x, -1);
 41		exponent++;
 42		argsign = SIGN_NEG;
 43	} else {
 44		/* Treat as  1 <= st0_ptr < sqrt(2) */
 45		x.sigh = st0_ptr->sigh - 0x80000000;
 46		x.sigl = st0_ptr->sigl;
 47		setexponent16(&x, 0);
 48		argsign = SIGN_POS;
 49	}
 50	tag = FPU_normalize_nuo(&x);
 51
 52	if (tag == TAG_Zero) {
 53		expon = 0;
 54		accumulator.msw = accumulator.midw = accumulator.lsw = 0;
 55	} else {
 56		log2_kernel(&x, argsign, &accumulator, &expon);
 57	}
 58
 59	if (exponent < 0) {
 60		sign = SIGN_NEG;
 61		exponent = -exponent;
 62	} else
 63		sign = SIGN_POS;
 64	expon_accum.msw = exponent;
 65	expon_accum.midw = expon_accum.lsw = 0;
 66	if (exponent) {
 67		expon_expon = 31 + norm_Xsig(&expon_accum);
 68		shr_Xsig(&accumulator, expon_expon - expon);
 69
 70		if (sign ^ argsign)
 71			negate_Xsig(&accumulator);
 72		add_Xsig_Xsig(&accumulator, &expon_accum);
 73	} else {
 74		expon_expon = expon;
 75		sign = argsign;
 76	}
 77
 78	yaccum.lsw = 0;
 79	XSIG_LL(yaccum) = significand(st1_ptr);
 80	mul_Xsig_Xsig(&accumulator, &yaccum);
 81
 82	expon_expon += round_Xsig(&accumulator);
 83
 84	if (accumulator.msw == 0) {
 85		FPU_copy_to_reg1(&CONST_Z, TAG_Zero);
 86		return;
 87	}
 88
 89	significand(st1_ptr) = XSIG_LL(accumulator);
 90	setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);
 91
 92	tag = FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign ^ st1_sign);
 93	FPU_settagi(1, tag);
 94
 95	set_precision_flag_up();	/* 80486 appears to always do this */
 96
 97	return;
 98
 99}
100
101/*--- poly_l2p1() -----------------------------------------------------------+
102 |   Base 2 logarithm by a polynomial approximation.                         |
103 |   log2(x+1)                                                               |
104 +---------------------------------------------------------------------------*/
105int poly_l2p1(u_char sign0, u_char sign1,
106	      FPU_REG * st0_ptr, FPU_REG * st1_ptr, FPU_REG * dest)
107{
108	u_char tag;
109	long int exponent;
110	Xsig accumulator, yaccum;
111
112	if (exponent16(st0_ptr) < 0) {
113		log2_kernel(st0_ptr, sign0, &accumulator, &exponent);
114
115		yaccum.lsw = 0;
116		XSIG_LL(yaccum) = significand(st1_ptr);
117		mul_Xsig_Xsig(&accumulator, &yaccum);
118
119		exponent += round_Xsig(&accumulator);
120
121		exponent += exponent16(st1_ptr) + 1;
122		if (exponent < EXP_WAY_UNDER)
123			exponent = EXP_WAY_UNDER;
124
125		significand(dest) = XSIG_LL(accumulator);
126		setexponent16(dest, exponent);
127
128		tag = FPU_round(dest, 1, 0, FULL_PRECISION, sign0 ^ sign1);
129		FPU_settagi(1, tag);
130
131		if (tag == TAG_Valid)
132			set_precision_flag_up();	/* 80486 appears to always do this */
133	} else {
134		/* The magnitude of st0_ptr is far too large. */
135
136		if (sign0 != SIGN_POS) {
137			/* Trying to get the log of a negative number. */
138#ifdef PECULIAR_486		/* Stupid 80486 doesn't worry about log(negative). */
139			changesign(st1_ptr);
140#else
141			if (arith_invalid(1) < 0)
142				return 1;
143#endif /* PECULIAR_486 */
144		}
145
146		/* 80486 appears to do this */
147		if (sign0 == SIGN_NEG)
148			set_precision_flag_down();
149		else
150			set_precision_flag_up();
151	}
152
153	if (exponent(dest) <= EXP_UNDER)
154		EXCEPTION(EX_Underflow);
155
156	return 0;
157
158}
159
160#undef HIPOWER
161#define	HIPOWER	10
162static const unsigned long long logterms[HIPOWER] = {
163	0x2a8eca5705fc2ef0LL,
164	0xf6384ee1d01febceLL,
165	0x093bb62877cdf642LL,
166	0x006985d8a9ec439bLL,
167	0x0005212c4f55a9c8LL,
168	0x00004326a16927f0LL,
169	0x0000038d1d80a0e7LL,
170	0x0000003141cc80c6LL,
171	0x00000002b1668c9fLL,
172	0x000000002c7a46aaLL
173};
174
175static const unsigned long leadterm = 0xb8000000;
176
177/*--- log2_kernel() ---------------------------------------------------------+
178 |   Base 2 logarithm by a polynomial approximation.                         |
179 |   log2(x+1)                                                               |
180 +---------------------------------------------------------------------------*/
181static void log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
182			long int *expon)
183{
184	long int exponent, adj;
185	unsigned long long Xsq;
186	Xsig accumulator, Numer, Denom, argSignif, arg_signif;
187
188	exponent = exponent16(arg);
189	Numer.lsw = Denom.lsw = 0;
190	XSIG_LL(Numer) = XSIG_LL(Denom) = significand(arg);
191	if (argsign == SIGN_POS) {
192		shr_Xsig(&Denom, 2 - (1 + exponent));
193		Denom.msw |= 0x80000000;
194		div_Xsig(&Numer, &Denom, &argSignif);
195	} else {
196		shr_Xsig(&Denom, 1 - (1 + exponent));
197		negate_Xsig(&Denom);
198		if (Denom.msw & 0x80000000) {
199			div_Xsig(&Numer, &Denom, &argSignif);
200			exponent++;
201		} else {
202			/* Denom must be 1.0 */
203			argSignif.lsw = Numer.lsw;
204			argSignif.midw = Numer.midw;
205			argSignif.msw = Numer.msw;
206		}
207	}
208
209#ifndef PECULIAR_486
210	/* Should check here that  |local_arg|  is within the valid range */
211	if (exponent >= -2) {
212		if ((exponent > -2) || (argSignif.msw > (unsigned)0xafb0ccc0)) {
213			/* The argument is too large */
214		}
215	}
216#endif /* PECULIAR_486 */
217
218	arg_signif.lsw = argSignif.lsw;
219	XSIG_LL(arg_signif) = XSIG_LL(argSignif);
220	adj = norm_Xsig(&argSignif);
221	accumulator.lsw = argSignif.lsw;
222	XSIG_LL(accumulator) = XSIG_LL(argSignif);
223	mul_Xsig_Xsig(&accumulator, &accumulator);
224	shr_Xsig(&accumulator, 2 * (-1 - (1 + exponent + adj)));
225	Xsq = XSIG_LL(accumulator);
226	if (accumulator.lsw & 0x80000000)
227		Xsq++;
228
229	accumulator.msw = accumulator.midw = accumulator.lsw = 0;
230	/* Do the basic fixed point polynomial evaluation */
231	polynomial_Xsig(&accumulator, &Xsq, logterms, HIPOWER - 1);
232
233	mul_Xsig_Xsig(&accumulator, &argSignif);
234	shr_Xsig(&accumulator, 6 - adj);
235
236	mul32_Xsig(&arg_signif, leadterm);
237	add_two_Xsig(&accumulator, &arg_signif, &exponent);
238
239	*expon = exponent + 1;
240	accum_result->lsw = accumulator.lsw;
241	accum_result->midw = accumulator.midw;
242	accum_result->msw = accumulator.msw;
243
244}