unit poly_sin_unit;

interface

uses fpu_types;

const
  N_COEFF_P = 4 - 1;

const
  N_COEFF_N = 4 - 1;

const
  N_COEFF_PH = 4 - 1;

const
  N_COEFF_NH = 4 - 1;

const
  pos_terms_l: array[0..N_COEFF_P] of u64 =
    (
    $aaaaaaaaaaaaaaab,
    $00d00d00d00cf906,
    $000006b99159a8bb,
    $000000000d7392e6);

const
  neg_terms_l: array[0..N_COEFF_N] of u64 =
    (
    $2222222222222167,
    $0002e3bc74aab624,
    $0000000b09229062,
    $00000000000c7973);

const
  pos_terms_h: array[0..N_COEFF_PH] of u64 =
    (
    $0000000000000000,
    $05b05b05b05b0406,
    $000049f93edd91a9,
    $00000000c9c9ed62);

const
  neg_terms_h: array[0..N_COEFF_NH] of u64 =
    (
    $aaaaaaaaaaaaaa98,
    $001a01a01a019064,
    $0000008f76c68a77,
    $0000000000d58f5e);

procedure poly_sine (st0_ptr: PFPU_REG);
procedure poly_cos (st0_ptr: PFPU_REG);

implementation

uses fpu_system, mul_Xsig, Service_fpu, shr_Xsig_unit, polynom_xsig,
  round_Xsig_unit, errors;

procedure poly_sine (st0_ptr: PFPU_REG);
  var
    exponent_, echange: intg;
    accumulator, argSqrd, argTo4: Xsig;
    fix_up, adj: s32;
    fixed_arg: u64;
    res: FPU_REG;
  begin

    exponent_ := exponent(st0_ptr);

    accumulator.lsw  := 0;
    accumulator.midw := 0;
    accumulator.msw  := 0;

    (* Split into two ranges, for arguments below and above 1.0 *)
    (* The boundary between upper and lower is approx 0.88309101259 *)
    if ((exponent_ < -1) or ((exponent_ = -1) and (st0_ptr^.sigh <= $e21240aa))) then
      begin
      (* The argument is <= 0.883091012ù59 *)

      argSqrd.msw  := st0_ptr^.sigh;
      argSqrd.midw := st0_ptr^.sigl;
      argSqrd.lsw  := 0;
      mul64_Xsig(@argSqrd, significand(st0_ptr));
      shr_Xsig(@argSqrd, 2 * (-1 - exponent_));
      argTo4.msw  := argSqrd.msw;
      argTo4.midw := argSqrd.midw;
      argTo4.lsw  := argSqrd.lsw;
      mul_Xsig_Xsig(@argTo4, @argTo4);

      polynomial_Xsig(@accumulator, @pu64(argTo4.msw)^, @neg_terms_l,
        N_COEFF_N - 1);
      mul_Xsig_Xsig(@accumulator, @argSqrd);
      negate_Xsig(@accumulator);

      polynomial_Xsig(@accumulator, @pu64(argTo4.msw)^, @pos_terms_l,
        N_COEFF_P - 1);

      shr_Xsig(@accumulator, 2);    (* Divide by four *)
      accumulator.msw := accumulator.msw or $8000;  (* Add 1.0 *)

      mul64_Xsig(@accumulator, significand(st0_ptr));
      mul64_Xsig(@accumulator, significand(st0_ptr));
      mul64_Xsig(@accumulator, significand(st0_ptr));

      (* Divide by four, FPU_REG compatible, etc *)
      exponent_ := 3 * exponent_;

      (* The minimum exponent_ difference is 3 *)
      shr_Xsig(@accumulator, exponent(st0_ptr) - exponent_);

      negate_Xsig(@accumulator);
      pu64(accumulator.msw)^ := pu64(accumulator.msw)^ + significand(st0_ptr)^;

      echange := round_Xsig(@accumulator);

      setexponentpos(@res, exponent(st0_ptr) + echange);
      end
    else
      begin
      (* The argument is > 0.88309101259 *)
      (* We use sin(st(0)) := cos(pi/2-st(0)) *)

      fixed_arg := significand(st0_ptr)^;

      if (exponent_ = 0) then
        begin
        (* The argument is >= 1.0 *)

        (* Put the binary point at the left. *)
        fixed_arg := fixed_arg shl 1;
        end;
      (* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 *)
      fixed_arg := $921fb54442d18469 - fixed_arg;
      (* There is a special case which arises due to rounding, to fix here. *)
      if (fixed_arg = $ffffffffffffffff) then
        fixed_arg := 0;

      XSIG_LL(argSqrd)^ := fixed_arg;
      argSqrd.lsw := 0;
      mul64_Xsig(@argSqrd, @fixed_arg);

      XSIG_LL(argTo4)^ := XSIG_LL(argSqrd)^;
      argTo4.lsw := argSqrd.lsw;
      mul_Xsig_Xsig(@argTo4, @argTo4);

      polynomial_Xsig(@accumulator, XSIG_LL(argTo4), @neg_terms_h,
        N_COEFF_NH - 1);
      mul_Xsig_Xsig(@accumulator, @argSqrd);
      negate_Xsig(@accumulator);

      polynomial_Xsig(@accumulator, XSIG_LL(argTo4), @pos_terms_h,
        N_COEFF_PH - 1);
      negate_Xsig(@accumulator);

      mul64_Xsig(@accumulator, @fixed_arg);
      mul64_Xsig(@accumulator, @fixed_arg);

      shr_Xsig(@accumulator, 3);
      negate_Xsig(@accumulator);

      add_Xsig_Xsig(@accumulator, @argSqrd);

      shr_Xsig(@accumulator, 1);

      accumulator.lsw := accumulator.lsw or 1;
      (* A zero accumulator here would cause problems *)
      negate_Xsig(@accumulator);

      (* The basic computation is complete. Now fix the answer to
   compensate for the error due to the approximation used for
   pi/2
   *)

      (* This has an exponent_ of -65 *)
      fix_up := $898cc517;
      (* The fix-up needs to be improved for larger args *)
      if (argSqrd.msw and $ffc00000) <> 0 then
        begin
        (* Get about 32 bit precision in these: *)
        fix_up := fix_up - Trunc(mul_32_32($898cc517, argSqrd.msw) / 6);
        end;
      fix_up := mul_32_32(fix_up, LL_MSW(@fixed_arg)^);

      adj := accumulator.lsw;    (* temp save *)
      accumulator.lsw := accumulator.lsw - fix_up;
      if (accumulator.lsw > adj) then
        Dec(XSIG_LL(accumulator)^);

      echange := round_Xsig(@accumulator);

      setexponentpos(@res, echange - 1);
      end;

    significand(@res)^ := XSIG_LL(accumulator)^;
    setsign(@res, getsign(st0_ptr));
    FPU_copy_to_reg0(@res, TAG_Valid);

{$ifdef PARANOID}
  if ( (exponent_(@res) >= 0)
      @ and (significand(@res) > BX_CONST64($800000000000)) )
    begin
      EXCEPTION(EX_INTERNAL|$150);
    end;
{$endif}(* PARANOID *)

  end;


(*--- poly_cos() ------------------------------------------------------------+
or                                                                          |
 +---------------------------------------------------------------------------*)
procedure poly_cos (st0_ptr: PFPU_REG);
  var
    res: FPU_REG;
    exponent_, exp2, echange: s32;
    accumulator, argSqrd, fix_up, argTo4: Xsig;
    fixed_arg: u64;
  begin

{$ifdef PARANOID}
  if ( (exponent_(st0_ptr) > 0)
      or ((exponent_(st0_ptr) = 0)
	  @ and (significand(st0_ptr) > BX_CONST64($c90fdaa22168c234))) )
    begin
      EXCEPTION(EX_Invalid);
      FPU_copy_to_reg0(@CONST_QNaN, TAG_Special);
      exit;
    end;
{$endif}(* PARANOID *)

    exponent_ := exponent(st0_ptr);

    accumulator.lsw  := 0;
    accumulator.midw := 0;
    accumulator.msw  := 0;

    if ((exponent_ < -1) or ((exponent_ = -1) and (st0_ptr^.sigh <= $b00d6f54))) then
      begin
      (* arg is < 0.687705 *)

      argSqrd.msw  := st0_ptr^.sigh;
      argSqrd.midw := st0_ptr^.sigl;
      argSqrd.lsw  := 0;
      mul64_Xsig(@argSqrd, significand(st0_ptr));

      if (exponent_ < -1) then
        begin
        (* shift the argument right by the required places *)
        shr_Xsig(@argSqrd, 2 * (-1 - exponent_));
        end;

      argTo4.msw  := argSqrd.msw;
      argTo4.midw := argSqrd.midw;
      argTo4.lsw  := argSqrd.lsw;
      mul_Xsig_Xsig(@argTo4, @argTo4);

      polynomial_Xsig(@accumulator, XSIG_LL(argTo4), @neg_terms_h,
        N_COEFF_NH - 1);
      mul_Xsig_Xsig(@accumulator, @argSqrd);
      negate_Xsig(@accumulator);

      polynomial_Xsig(@accumulator, XSIG_LL(argTo4), @pos_terms_h,
        N_COEFF_PH - 1);
      negate_Xsig(@accumulator);

      mul64_Xsig(@accumulator, significand(st0_ptr));
      mul64_Xsig(@accumulator, significand(st0_ptr));
      shr_Xsig(@accumulator, -2 * (1 + exponent_));

      shr_Xsig(@accumulator, 3);
      negate_Xsig(@accumulator);

      add_Xsig_Xsig(@accumulator, @argSqrd);

      shr_Xsig(@accumulator, 1);

      (* It doesn't matter if accumulator is all zero here, the
   following code will work ok *)
      negate_Xsig(@accumulator);

      if (accumulator.lsw and $8000) <> 0 then
        Inc(XSIG_LL(accumulator)^);
      if (accumulator.msw = 0) then
        begin
        (* The res is 1.0 *)
        FPU_copy_to_reg0(@CONST_1, TAG_Valid);
        exit;
        end
      else
        begin
        significand(@res)^ := XSIG_LL(accumulator)^;

        (* will be a valid positive nr with expon := -1 *)
        setexponentpos(@res, -1);
        end;
      end
    else
      begin
      fixed_arg := significand(st0_ptr)^;

      if (exponent_ = 0) then
        begin
        (* The argument is >= 1.0 *)

        (* Put the binary point at the left. *)
        fixed_arg := fixed_arg shl 1;
        end;
      (* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 *)
      fixed_arg := $921fb54442d18469 - fixed_arg;
      (* There is a special case which arises due to rounding, to fix here. *)
      if (fixed_arg = $ffffffffffffffff) then
        fixed_arg := 0;

      exponent_ := -1;
      exp2 := -1;

      (* A shift is needed here only for a narrow range of arguments,
   i.e. for fixed_arg approx 2^-32, but we pick up more... *)
      if ((LL_MSW(@fixed_arg)^ and $ffff0000) = 0) then
        begin
        fixed_arg := fixed_arg shl 16;
        exponent_ := exponent_ - 16;
        exp2 := exp2 - 16;
        end;

      XSIG_LL(argSqrd)^ := fixed_arg;
      argSqrd.lsw := 0;
      mul64_Xsig(@argSqrd, @fixed_arg);

      if (exponent_ < -1) then
        begin
        (* shift the argument right by the required places *)
        shr_Xsig(@argSqrd, 2 * (-1 - exponent_));
        end;

      argTo4.msw  := argSqrd.msw;
      argTo4.midw := argSqrd.midw;
      argTo4.lsw  := argSqrd.lsw;
      mul_Xsig_Xsig(@argTo4, @argTo4);

      polynomial_Xsig(@accumulator, XSIG_LL(argTo4), @neg_terms_l,
        N_COEFF_N - 1);
      mul_Xsig_Xsig(@accumulator, @argSqrd);
      negate_Xsig(@accumulator);

      polynomial_Xsig(@accumulator, XSIG_LL(argTo4), @pos_terms_l,
        N_COEFF_P - 1);

      shr_Xsig(@accumulator, 2);    (* Divide by four *)
      accumulator.msw := accumulator.msw or $8000;  (* Add 1.0 *)

      mul64_Xsig(@accumulator, @fixed_arg);
      mul64_Xsig(@accumulator, @fixed_arg);
      mul64_Xsig(@accumulator, @fixed_arg);

      (* Divide by four, FPU_REG compatible, etc *)
      exponent_ := 3 * exponent_;

      (* The minimum exponent_ difference is 3 *)
      shr_Xsig(@accumulator, exp2 - exponent_);

      negate_Xsig(@accumulator);
      XSIG_LL(accumulator)^ := XSIG_LL(accumulator)^ + fixed_arg;

      (* The basic computation is complete. Now fix the answer to
   compensate for the error due to the approximation used for
   pi/2
   *)

      (* This has an exponent_ of -65 *)
      XSIG_LL(fix_up)^ := $898cc51701b839a2;
      fix_up.lsw := 0;

      (* The fix-up needs to be improved for larger args *)
      if (argSqrd.msw and $ffc00000) <> 0 then
        begin
        (* Get about 32 bit precision in these: *)
        fix_up.msw := Trunc(fix_up.msw - mul_32_32($898cc517, argSqrd.msw) / 2);
        fix_up.msw := Trunc(fix_up.msw + mul_32_32($898cc517, argTo4.msw) / 24);
        end;

      exp2 := exp2 + norm_Xsig(@accumulator);
      shr_Xsig(@accumulator, 1); (* Prevent overflow *)
      Inc(exp2);
      shr_Xsig(@fix_up, 65 + exp2);

      add_Xsig_Xsig(@accumulator, @fix_up);

      echange := round_Xsig(@accumulator);

      setexponentpos(@res, exp2 + echange);
      significand(@res)^ := XSIG_LL(accumulator)^;
      end;

    FPU_copy_to_reg0(@res, TAG_Valid);

{$ifdef PARANOID}
  if ( (exponent_(@res) >= 0)
      @ and (significand(@res) > BX_CONST64($800000000000)) )
    begin
      EXCEPTION(EX_INTERNAL|$151);
    end;
{$endif}(* PARANOID *)

  end;

end.
