unit poly_l2;

interface

implementation

procedure	poly_l2(st0_ptr:PFPU_REG; st1_ptr:PFPU_REG; st1_sign:u_char);
begin
  s32	       exponent, expon, expon_expon;
  Xsig         accumulator, expon_accum, yaccum;
  u_char       sign, argsign;
  FPU_REG      x;
  int          tag;

  exponent := exponent16(st0_ptr);

  (* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) *)
  if ( st0_ptr^.sigh > (unsigned)$b504f334 )
    begin
      (* Treat as  sqrt(2)/2 < st0_ptr < 1 *)
      significand(@x) := - significand(st0_ptr);
      setexponent16(@x, -1);
      exponent++;
      argsign := SIGN_NEG;
    end;
  else
    begin
      (* Treat as  1 <= st0_ptr < sqrt(2) *)
      x.sigh := st0_ptr^.sigh - $80000000;
      x.sigl := st0_ptr^.sigl;
      setexponent16(@x, 0);
      argsign := SIGN_POS;
    end;
  tag := FPU_normalize_nuo(@x, 0);

  if ( tag = TAG_Zero )
    begin
      expon := 0;
      accumulator.msw := accumulator.midw := accumulator.lsw := 0;
    end;
  else
    begin
      log2_kernel(@x, argsign, @accumulator, @expon);
    end;

  if ( exponent < 0 )
    begin
      sign := SIGN_NEG;
      exponent := -exponent;
    end;
  else
    sign := SIGN_POS;
  expon_accum.msw := exponent; expon_accum.midw := expon_accum.lsw := 0;
  if ( exponent )
    begin
      expon_expon := 31 + norm_Xsig(@expon_accum);
      shr_Xsig(@accumulator, expon_expon - expon);

      if ( sign or argsign )
	negate_Xsig(@accumulator);
      add_Xsig_Xsig(@accumulator, @expon_accum);
    end;
  else
    begin
      expon_expon := expon;
      sign := argsign;
    end;

  yaccum.lsw := 0; XSIG_LL(yaccum) := significand(st1_ptr);
  mul_Xsig_Xsig(@accumulator, @yaccum);

  expon_expon :=  expon_expon + round_Xsig(@accumulator);

  if ( accumulator.msw = 0 )
    begin
      FPU_copy_to_reg1(@CONST_Z, TAG_Zero);
      exit;
    end;

  significand(st1_ptr) := XSIG_LL(accumulator);
  setexponent16(st1_ptr, expon_expon + exponent16(st1_ptr) + 1);

  tag := FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign or st1_sign);
  FPU_settagi(1, tag);

  set_precision_flag_up();  (* 80486 appears to always do this *)

  exit;

end;


(*--- poly_l2p1() -----------------------------------------------------------+
or  Base 2 logarithm by a polynomial approximation.                         |
or  log2(x+1)                                                               |
 +---------------------------------------------------------------------------*)
int	poly_l2p1(u_char sign0, u_char sign1,
		  FPU_REG *st0_ptr, FPU_REG *st1_ptr, FPU_REG *dest)
begin
  u_char       	tag;
  s32        	exponent;
  Xsig         	accumulator, yaccum;

  if ( exponent16(st0_ptr) < 0 )
    begin
      log2_kernel(st0_ptr, sign0, @accumulator, @exponent);

      yaccum.lsw := 0;
      XSIG_LL(yaccum) := significand(st1_ptr);
      mul_Xsig_Xsig(@accumulator, @yaccum);

      exponent :=      exponent + round_Xsig(@accumulator);

      exponent :=      exponent + exponent16(st1_ptr) + 1;
      if ( exponent < EXP_WAY_UNDER ) exponent := EXP_WAY_UNDER;

      significand(dest) := XSIG_LL(accumulator);
      setexponent16(dest, exponent);

      tag := FPU_round(dest, 1, 0, FULL_PRECISION, sign0 or sign1);
      FPU_settagi(1, tag);

      if ( tag = TAG_Valid )
	set_precision_flag_up();   (* 80486 appears to always do this *)
    end;
  else
    begin
      (* The magnitude of st0_ptr is far too large. *)

      if ( sign0 !:= SIGN_POS )
	begin
	  (* Trying to get the log of a negative number. *)
{$ifdef} PECULIAR_486   (* Stupid 80486 doesn't worry about log(negative). *)
	  changesign(st1_ptr);
{$else}
	  if ( arith_invalid(1) < 0 )
	    return 1;
{$endif} (* PECULIAR_486 *)
	end;

      (* 80486 appears to do this *)
      if ( sign0 = SIGN_NEG )
	set_precision_flag_down();
      else
	set_precision_flag_up();
    end;

  if ( exponent(dest) <= EXP_UNDER )
    EXCEPTION(EX_Underflow);

  return 0;

end;




#undef HIPOWER
#define	HIPOWER	10
static const u64 logterms[HIPOWER] :=
begin
  BX_CONST64($2a8eca5705fc2ef0),
  BX_CONST64($f6384ee1d01febce),
  BX_CONST64($093bb62877cdf642),
  BX_CONST64($006985d8a9ec439b),
  BX_CONST64($0005212c4f55a9c8),
  BX_CONST64($00004326a16927f0),
  BX_CONST64($0000038d1d80a0e7),
  BX_CONST64($0000003141cc80c6),
  BX_CONST64($00000002b1668c9f),
  BX_CONST64($000000002c7a46aa)
end;;

static const u32 leadterm := $b8000000;


(*--- log2_kernel() ---------------------------------------------------------+
or  Base 2 logarithm by a polynomial approximation.                         |
or  log2(x+1)                                                               |
 +---------------------------------------------------------------------------*)
static procedure log2_kernel(FPU_REG const *arg, u_char argsign, Xsig *accum_result,
			s32 *expon)
begin
  s32    exponent, adj;
  u64    Xsq;
  Xsig   accumulator, Numer, Denom, argSignif, arg_signif;

  exponent := exponent16(arg);
  Numer.lsw := Denom.lsw := 0;
  XSIG_LL(Numer) := XSIG_LL(Denom) := significand(arg);
  if ( argsign = SIGN_POS )
    begin
      shr_Xsig(@Denom, 2 - (1 + exponent));
      Denom.msw :=      Denom.msw or $80000000;
      div_Xsig(@Numer, @Denom, @argSignif);
    end;
  else
    begin
      shr_Xsig(@Denom, 1 - (1 + exponent));
      negate_Xsig(@Denom);
      if ( Denom.msw  and $80000000 )
	begin
	  div_Xsig(@Numer, @Denom, @argSignif);
	  exponent ++;
	end;
      else
	begin
	  (* Denom must be 1.0 *)
	  argSignif.lsw := Numer.lsw; argSignif.midw := Numer.midw;
	  argSignif.msw := Numer.msw;
	end;
    end;

#ifndef PECULIAR_486
  (* Should check here that  |local_arg|  is within the valid range *)
  if ( exponent >= -2 )
    begin
      if ( (exponent > -2) or
	  (argSignif.msw > (unsigned)$afb0ccc0) )
	begin
	  (* The argument is too large *)
	end;
    end;
{$endif} (* PECULIAR_486 *)

  arg_signif.lsw := argSignif.lsw; XSIG_LL(arg_signif) := XSIG_LL(argSignif);
  adj := norm_Xsig(@argSignif);
  accumulator.lsw := argSignif.lsw; XSIG_LL(accumulator) := XSIG_LL(argSignif);
  mul_Xsig_Xsig(@accumulator, @accumulator);
  shr_Xsig(@accumulator, 2*(-1 - (1 + exponent + adj)));
  Xsq := XSIG_LL(accumulator);
  if ( accumulator.lsw  and $80000000 )
    Xsq++;

  accumulator.msw := accumulator.midw := accumulator.lsw := 0;
  (* Do the basic fixed point polynomial evaluation *)
  polynomial_Xsig(@accumulator, @Xsq, logterms, HIPOWER-1);

  mul_Xsig_Xsig(@accumulator, @argSignif);
  shr_Xsig(@accumulator, 6 - adj);

  mul32_Xsig(@arg_signif, leadterm);
  add_two_Xsig(@accumulator, @arg_signif, @exponent);

  *expon := exponent + 1;
  accum_result^.lsw := accumulator.lsw;
  accum_result^.midw := accumulator.midw;
  accum_result^.msw := accumulator.msw;

end;

end.
