{$ifdef LINUX}
{$include ../defines.pas}
{$else}
{$include ..\defines.pas}
{$endif}
unit poly_l2_unit;

interface

uses fpu_types;

procedure log2_kernel(const arg:PFPU_REG; argsign:u_char; accum_result:PXsig;
			expon:ps32);
procedure	poly_l2(st0_ptr:PFPU_REG; st1_ptr:PFPU_REG; st1_sign:u_char);         

function poly_l2p1(sign0:u_char; sign1:u_char;
		  st0_ptr:PFPU_REG; st1_ptr:PFPU_REG; dest:PFPU_REG):intg;
{$include fpu_const.pas}
implementation

uses Service_fpu, CONFIG, reg_norm, round_Xsig_unit, shr_Xsig_unit,
  fpu_system, mul_Xsig, reg_round, errors, div_xsig_unit, polynom_xsig;

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

  exponent_ := exponent16(st0_ptr);

  (* From st0_ptr, make a number > sqrt(2)/2 and < sqrt(2) *)
  if ( st0_ptr^.sigh > unsigned($b504f334) ) then
    begin
      (* Treat as  sqrt(2)/2 < st0_ptr < 1 *)
      significand(@x)^ := - significand(st0_ptr)^;
      setexponent16(@x, -1);
      inc(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 ) then
    begin
      expon := 0;
      accumulator.msw := 0;
      accumulator.midw := 0;
      accumulator.lsw := 0;
    end
  else
    begin
      log2_kernel(@x, argsign, @accumulator, @expon);
    end;

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

      if ( sign or argsign )<>0 then
	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 ) then
    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)                                                               |
 +---------------------------------------------------------------------------*)
function poly_l2p1(sign0:u_char; sign1:u_char;
		  st0_ptr:PFPU_REG; st1_ptr:PFPU_REG; dest:PFPU_REG):intg;
var
  tag:u_char;
  exponent_:s32;
  accumulator, yaccum:Xsig;
begin

  if ( exponent16(st0_ptr) < 0 ) then
    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 ) then
         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 ) then
	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 ) then
	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 ) then
      begin
	    result := 1;
       exit;
      end;
{$endif} (* PECULIAR_486 *)
	end;

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

  if ( exponent(dest) <= EXP_UNDER ) then
    fpu_EXCEPTION(EX_Underflow);

  result := 0;

end;


const	HIPOWER =	10;
const logterms:array[0..HIPOWER-1] of u64 = (
  $2a8eca5705fc2ef0,
  $f6384ee1d01febce,
  $093bb62877cdf642,
  $006985d8a9ec439b,
  $0005212c4f55a9c8,
  $00004326a16927f0,
  $0000038d1d80a0e7,
  $0000003141cc80c6,
  $00000002b1668c9f,
  $000000002c7a46aa
);

const leadterm:u32 = $b8000000;


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

  exponent_ := exponent16(arg);
  Numer.lsw := 0;
  Denom.lsw := 0;
  XSIG_LL(Denom)^ := significand(arg)^;
  XSIG_LL(Numer)^ := XSIG_LL(Denom)^;
  if ( argsign = SIGN_POS ) then
    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 )<> 0 then
	begin
	  div_Xsig(@Numer, @Denom, @argSignif);
	  inc(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 )<>0 then
    inc(Xsq);

  accumulator.msw := 0;
  accumulator.midw := 0;
  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.
