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