unit poly_tan_unit;

interface

uses fpu_types;

const
  HiPOWERop = 3;

const
  oddplterm: array [0..HiPOWERop - 1] of u64 = (
    $0000000000000000,
    $0051a1cf08fca228,
    $0000000071284ff7);

const
  HiPOWERon = 2;  { odd poly, negative terms }

const
  oddnegterm: array[0..HiPOWERon - 1] of u64 = (
    $1291a9a184244e80,
    $0000583245819c21);

const
  HiPOWERep = 2;  { even poly, positive terms }

const
  evenplterm: array[0..HiPOWERep - 1] of u64 = (
    $0e848884b539e888,
    $00003c7f18b887da);

const
  HiPOWERen = 2;  { even poly, negative terms }

const
  evennegterm: array[0..HiPOWERen - 1] of u64 = (
    $f1f0200fd51569cc,
    $003afb46105c4432);

const
  twothirds: u64 = $aaaaaaaaaaaaaaab;

procedure poly_tan (st0_ptr: PFPU_REG; invert: intg);

implementation

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

procedure poly_tan (st0_ptr: PFPU_REG; invert: intg);
  var
    exponent_: s32;
    argSq, argSqSq, accumulatoro, accumulatore, accum, argSignif: Xsig;
  begin

    exponent_ := exponent(st0_ptr);


{$ifdef PARANOID}
  if ( signnegative(st0_ptr) )<>0 then	(* Can't hack a number < 0.0 *)
    begin arith_invalid(0); exit; end;  (* Need a positive number *)
{$endif}(* PARANOID *)

    if ((exponent_ >= 0) or ((exponent_ = -1) and (st0_ptr^.sigh > $c90fdaa2))) then
      begin
      fpu_EXCEPTION($250);
      end
    else
      begin
      argSignif.lsw := 0;
      pu64(accum.msw)^ := significand(st0_ptr)^;
      pu64(accum.msw)^ := pu64(accum.msw)^;
      ;

      if (exponent_ < -1) then
        begin
        (* shift the argument right by the required places *)
        if (FPU_shrx(@pu64(accum.msw)^, -1 - exponent_) >= $8000) then
          Inc(pu64(accum.msw)^);  (* round up *)
        end;
      end;

    pu64(argSq.msw)^ := pu64(accum.msw)^;
    argSq.lsw := accum.lsw;
    mul_Xsig_Xsig(@argSq, @argSq);
    pu64(argSqSq.msw)^ := pu64(argSq.msw)^;
    argSqSq.lsw := argSq.lsw;
    mul_Xsig_Xsig(@argSqSq, @argSqSq);

    (* Compute the negative terms for the numerator polynomial *)
    accumulatoro.msw  := 0;
    accumulatoro.midw := 0;
    accumulatoro.lsw  := 0;
    polynomial_Xsig(@accumulatoro, @pu64(argSqSq.msw)^, @oddnegterm, HiPOWERon - 1);
    mul_Xsig_Xsig(@accumulatoro, @argSq);
    negate_Xsig(@accumulatoro);
    (* Add the positive terms *)
    polynomial_Xsig(@accumulatoro, @pu64(argSqSq.msw)^, @oddplterm, HiPOWERop - 1);


    (* Compute the positive terms for the denominator polynomial *)
    accumulatore.msw  := 0;
    accumulatore.midw := 0;
    accumulatore.lsw  := 0;
    polynomial_Xsig(@accumulatore, @pu64(argSqSq.msw)^, @evenplterm, HiPOWERep - 1);
    mul_Xsig_Xsig(@accumulatore, @argSq);
    negate_Xsig(@accumulatore);
    (* Add the negative terms *)
    polynomial_Xsig(@accumulatore, @pu64(argSqSq.msw)^, @evennegterm, HiPOWERen - 1);
    (* Multiply by arg^2 *)
    mul64_Xsig(@accumulatore, @pu64(argSignif.msw)^);
    mul64_Xsig(@accumulatore, @pu64(argSignif.msw)^);
    (* de-normalize and divide by 2 *)
    shr_Xsig(@accumulatore, -2 * (1 + exponent_) + 1);
    negate_Xsig(@accumulatore);      (* This does 1 - accumulator *)

    (* Now find the ratio. *)
    if (accumulatore.msw = 0) then
      begin
      (* accumulatoro must contain 1.0 here, (actually, 0) but it
   really doesn't matter what value we use because it will
   have negligible effect in later calculations
   *)
      pu64(accum.msw)^ := $800000000000;
      accum.lsw := 0;
      end
    else
      begin
      div_Xsig(@accumulatoro, @accumulatore, @accum);
      end;

    (* Multiply by 1/3 * arg^3 *)
    mul64_Xsig(@accum, @pu64(argSignif.msw)^);
    mul64_Xsig(@accum, @pu64(argSignif.msw)^);
    mul64_Xsig(@accum, @pu64(argSignif.msw)^);
    mul64_Xsig(@accum, @twothirds);
    shr_Xsig(@accum, -2 * (exponent_ + 1));


    (* tan(arg) := arg + accum *)
    add_two_Xsig(@accum, @argSignif, @exponent);

    if (invert) <> 0 then
      begin
      (* accum now contains tan(pi/2 - arg).
   Use tan(arg) := 1.0 / tan(pi/2 - arg)
   *)
      accumulatoro.lsw  := 0;
      accumulatoro.midw := 0;
      accumulatoro.msw  := $8000;
      div_Xsig(@accumulatoro, @accum, @accum);
      exponent_ := -exponent_;
      end;


    (* Transfer the result *)
    exponent_ := exponent_ + pu64(accum.msw)^;
    FPU_settag0(TAG_Valid);
    significand(st0_ptr)^ := pu64(accum.msw)^;
    setexponent16(st0_ptr, exponent_ + EXTENDED_Ebias);  (* Result is positive. *)

  end;


end.
