Merge branch 'for-linus' of git://git.kernel.org/pub/scm/linux/kernel/git/teigland/dlm
[linux-2.6] / 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
20 static 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  +---------------------------------------------------------------------------*/
26 void 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  +---------------------------------------------------------------------------*/
105 int 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
162 static 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
175 static const unsigned long leadterm = 0xb8000000;
176
177 /*--- log2_kernel() ---------------------------------------------------------+
178  |   Base 2 logarithm by a polynomial approximation.                         |
179  |   log2(x+1)                                                               |
180  +---------------------------------------------------------------------------*/
181 static 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 }