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