/*---------------------------------------------------------------------------+
 |  poly_l2.c                                                                |
 |                                                                           |
 | Compute the base 2 log of a FPU_REG, using a polynomial approximation.    |
 |                                                                           |
 | Copyright (C) 1992,1993                                                   |
 |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
 |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
 |                                                                           |
 |                                                                           |
 +---------------------------------------------------------------------------*/


#include "exception.h"
#include "reg_constant.h"
#include "fpu_emu.h"
#include "control_w.h"



#define	HIPOWER	9
static unsigned short const	lterms[HIPOWER][4] =
	{
	/* Ideal computation with these coeffs gives about
	   64.6 bit rel accuracy. */
	{ 0xe177, 0xb82f, 0x7652, 0x7154 },
	{ 0xee0f, 0xe80f, 0x2770, 0x7b1c },
	{ 0x0fc0, 0xbe87, 0xb143, 0x49dd },
	{ 0x78b9, 0xdadd, 0xec54, 0x34c2 },
	{ 0x003a, 0x5de9, 0x628b, 0x2909 },
	{ 0x5588, 0xed16, 0x4abf, 0x2193 },
	{ 0xb461, 0x85f7, 0x347a, 0x1c6a },
	{ 0x0975, 0x87b3, 0xd5bf, 0x1876 },
	{ 0xe85c, 0xcec9, 0x84e7, 0x187d }
	};




/*--- poly_l2() -------------------------------------------------------------+
 |   Base 2 logarithm by a polynomial approximation.                         |
 +---------------------------------------------------------------------------*/
void	poly_l2(FPU_REG const *arg, FPU_REG *result)
{
  short		  exponent;
  char		  zero;		/* flag for an Xx == 0 */
  unsigned short  bits, shift;
  unsigned long long       Xsq;
  FPU_REG	  accum, denom, num, Xx;


  exponent = arg->exp - EXP_BIAS;

  accum.tag = TW_Valid;	/* set the tags to Valid */

  if ( arg->sigh > (unsigned)0xb504f334 )
    {
      /* This is good enough for the computation of the polynomial
	 sum, but actually results in a loss of precision for
	 the computation of Xx. This will matter only if exponent
	 becomes zero. */
      exponent++;
      accum.sign = 1;	/* sign to negative */
      num.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
      reg_u_div(&CONST_1, arg, &num, FULL_PRECISION);
    }
  else
    {
      accum.sign = 0;	/* set the sign to positive */
      num.sigl = arg->sigl;		/* copy the mantissa */
      num.sigh = arg->sigh;
    }


  /* shift num left, lose the ms bit */
  num.sigh <<= 1;
  if ( num.sigl & 0x80000000 ) num.sigh |= 1;
  num.sigl <<= 1;

  denom.sigl = num.sigl;
  denom.sigh = num.sigh;
  poly_div4(&significand(&denom));
  denom.sigh += 0x80000000;			/* set the msb */
  Xx.exp = EXP_BIAS;  /* needed to prevent errors in div routine */
  reg_u_div(&num, &denom, &Xx, FULL_PRECISION);

  zero = !(Xx.sigh | Xx.sigl);
  
  mul64(&significand(&Xx), &significand(&Xx), &Xsq);
  poly_div16(&Xsq);

  accum.exp = -1;		/* exponent of accum */

  /* Do the basic fixed point polynomial evaluation */
  polynomial((unsigned *)&accum.sigl, (unsigned *)&Xsq, lterms, HIPOWER-1);

  if ( !exponent )
    {
      /* If the exponent is zero, then we would lose precision by
	 sticking to fixed point computation here */
      /* We need to re-compute Xx because of loss of precision. */
      FPU_REG   lXx;
      char	sign;
      
      sign = accum.sign;
      accum.sign = 0;

      /* make accum compatible and normalize */
      accum.exp = EXP_BIAS + accum.exp;
      normalize(&accum);

      if ( zero )
	{
	  reg_move(&CONST_Z, result);
	}
      else
	{
	  /* we need to re-compute lXx to better accuracy */
	  num.tag = TW_Valid;		/* set the tags to Vaild */
	  num.sign = 0;		/* set the sign to positive */
	  num.exp = EXP_BIAS - 1;
	  if ( sign )
	    {
	      /* The argument is of the form 1-x */
	      /* Use  1-1/(1-x) = x/(1-x) */
	      significand(&num) = - significand(arg);
	      normalize(&num);
	      reg_div(&num, arg, &num, FULL_PRECISION);
	    }
	  else
	    {
	      normalize(&num);
	    }

	  denom.tag = TW_Valid;	/* set the tags to Valid */
	  denom.sign = SIGN_POS;	/* set the sign to positive */
	  denom.exp = EXP_BIAS;
	  
	  reg_div(&num, &denom, &lXx, FULL_PRECISION);

	  reg_u_mul(&lXx, &accum, &accum, FULL_PRECISION);

	  reg_u_add(&lXx, &accum, result, FULL_PRECISION);
	  
	  normalize(result);
	}

      result->sign = sign;
      return;
    }

  mul64(&significand(&accum),
	&significand(&Xx), &significand(&accum));

  significand(&accum) += significand(&Xx);

  if ( Xx.sigh > accum.sigh )
    {
      /* There was an overflow */

      poly_div2(&significand(&accum));
      accum.sigh |= 0x80000000;
      accum.exp++;
    }

  /* When we add the exponent to the accum result later, we will
     require that their signs are the same. Here we ensure that
     this is so. */
  if ( exponent && ((exponent < 0) ^ (accum.sign)) )
    {
      /* signs are different */

      accum.sign = !accum.sign;

      /* An exceptional case is when accum is zero */
      if ( accum.sigl | accum.sigh )
	{
	  /* find 1-accum */
	  /* Shift to get exponent == 0 */
	  if ( accum.exp < 0 )
	    {
	      poly_div2(&significand(&accum));
	      accum.exp++;
	    }
	  /* Just negate, but throw away the sign */
	  significand(&accum) = - significand(&accum);
	  if ( exponent < 0 )
	    exponent++;
	  else
	    exponent--;
	}
    }

  shift = exponent >= 0 ? exponent : -exponent ;
  bits = 0;
  if ( shift )
    {
      if ( accum.exp )
	{
	  accum.exp++;
	  poly_div2(&significand(&accum));
	}
      while ( shift )
	{
	  poly_div2(&significand(&accum));
	  if ( shift & 1)
	    accum.sigh |= 0x80000000;
	  shift >>= 1;
	  bits++;
	}
    }

  /* Convert to 64 bit signed-compatible */
  accum.exp += bits + EXP_BIAS - 1;

  reg_move(&accum, result);
  normalize(result);

  return;
}


/*--- poly_l2p1() -----------------------------------------------------------+
 |   Base 2 logarithm by a polynomial approximation.                         |
 |   log2(x+1)                                                               |
 +---------------------------------------------------------------------------*/
int	poly_l2p1(FPU_REG const *arg, FPU_REG *result)
{
  char		sign = 0;
  unsigned long long     Xsq;
  FPU_REG      	arg_pl1, denom, accum, local_arg, poly_arg;


  sign = arg->sign;

  reg_add(arg, &CONST_1, &arg_pl1, FULL_PRECISION);

  if ( (arg_pl1.sign) | (arg_pl1.tag) )
    {			/* We need a valid positive number! */
      return 1;
    }

  reg_add(&CONST_1, &arg_pl1, &denom, FULL_PRECISION);
  reg_div(arg, &denom, &local_arg, FULL_PRECISION);
  local_arg.sign = 0;	/* Make the sign positive */

  /* Now we need to check that  |local_arg| is less than
     3-2*sqrt(2) = 0.17157.. = .0xafb0ccc0 * 2^-2 */

  if ( local_arg.exp >= EXP_BIAS - 3 )
    {
      if ( (local_arg.exp > EXP_BIAS - 3) ||
	  (local_arg.sigh > (unsigned)0xafb0ccc0) )
	{
	  /* The argument is large */
	  poly_l2(&arg_pl1, result); return 0;
	}
    }

  /* Make a copy of local_arg */
  reg_move(&local_arg, &poly_arg);

  /* Get poly_arg bits aligned as required */
  shrx((unsigned *)&(poly_arg.sigl), -(poly_arg.exp - EXP_BIAS + 3));

  mul64(&significand(&poly_arg), &significand(&poly_arg), &Xsq);
  poly_div16(&Xsq);

  /* Do the basic fixed point polynomial evaluation */
  polynomial(&(accum.sigl), (unsigned *)&Xsq, lterms, HIPOWER-1);

  accum.tag = TW_Valid;	/* set the tags to Valid */
  accum.sign = SIGN_POS;	/* and make accum positive */

  /* make accum compatible and normalize */
  accum.exp = EXP_BIAS - 1;
  normalize(&accum);

  reg_u_mul(&local_arg, &accum, &accum, FULL_PRECISION);

  reg_u_add(&local_arg, &accum, result, FULL_PRECISION);

  /* Multiply the result by 2 */
  result->exp++;

  result->sign = sign;
  
  return 0;
}