PageRenderTime 18ms CodeModel.GetById 1ms app.highlight 14ms RepoModel.GetById 1ms app.codeStats 0ms

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

https://bitbucket.org/evzijst/gittest
C | 272 lines | 192 code | 48 blank | 32 comment | 24 complexity | 1240a719feb33849a3763e6ca2751d4d MD5 | raw file
  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
 22static 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 +---------------------------------------------------------------------------*/
 29void	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 +---------------------------------------------------------------------------*/
119int	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
182static 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
196static const unsigned long leadterm = 0xb8000000;
197
198
199/*--- log2_kernel() ---------------------------------------------------------+
200 |   Base 2 logarithm by a polynomial approximation.                         |
201 |   log2(x+1)                                                               |
202 +---------------------------------------------------------------------------*/
203static 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}