unit poly_atan_unit;

interface

uses fpu_types;

const
  HIPOWERon = 6;  // odd poly, negative terms

const
  oddnegterms: array[0..HIPOWERon - 1] of u64 =
    (
    $0000000000000000,  //Dummy (not for - 1.0)
    $015328437f756467,
    $0005dda27b73dec6,
    $0000226bf2bfb91a,
    $000000ccc439c5f7,
    $0000000355438407);

const
  HIPOWERop = 6;  // odd poly, positive terms

const
  oddplterms: array[0..HIPOWERop - 1] of u64 = (
    //  $aaaaaaaaaaaaaaab,  transferred to fixedpterm[] */
    $0db55a71875c9ac2,
    $0029fce2d67880b0,
    $0000dfd3908b4596,
    $00000550fd61dab4,
    $0000001c9422b3f9,
    $000000003e3301e1);

const
  denomterm: u64 = $ebd9b842c5c53a0e;

const
  fixedpterm: Xsig = (lsw: $aaaaaaaa; midw: $aaaaaaaa; msw: $aaaaaaaa);

const
  pi_signif: Xsig = (lsw: $c90fdaa2; midw: $2168c234; msw: $c4c6628b);

procedure poly_atan (st0_ptr: PFPU_REG; st0_tag: u_char; st1_ptr: PFPU_REG;
  st1_tag: u_char);

implementation

uses Service_fpu, reg_convert, fpu_system, div_xsig_unit, round_Xsig_unit,
  shr_Xsig_unit, mul_Xsig, poly_2xm1_unit, polynom_xsig, reg_round, errors;

procedure poly_atan (st0_ptr: PFPU_REG; st0_tag: u_char; st1_ptr: PFPU_REG;
  st1_tag: u_char);
  var
    transformed, inverted, sign1, sign2: u_char;
    exponent_: s32;
    dummy_exp: s32;
    accumulator, Numer, Denom, accumulatore, argSignif, argSq, argSqSq: Xsig;
    tag: u_char;
  begin

    sign1 := getsign(st0_ptr);
    sign2 := getsign(st1_ptr);
    if (st0_tag = TAG_Valid) then
      begin
      exponent_ := exponent(st0_ptr);
      end
    else
      begin
      (* This gives non-compatible stack contents... *)
      FPU_to_exp16(st0_ptr, st0_ptr);
      exponent_ := exponent16(st0_ptr);
      end;
    if (st1_tag = TAG_Valid) then
      begin
      exponent_ := exponent_ - exponent(st1_ptr);
      end
    else
      begin
      (* This gives non-compatible stack contents... *)
      FPU_to_exp16(st1_ptr, st1_ptr);
      exponent_ := exponent_ - exponent16(st1_ptr);
      end;

    if ((exponent_ < 0) or ((exponent_ = 0) and
      ((st0_ptr^.sigh < st1_ptr^.sigh) or ((st0_ptr^.sigh = st1_ptr^.sigh) and
      (st0_ptr^.sigl < st1_ptr^.sigl))))) then
      begin
      inverted  := 1;
      Numer.lsw := 0;
      Denom.lsw := 0;
      XSIG_LL(Numer)^ := significand(st0_ptr)^;
      XSIG_LL(Denom)^ := significand(st1_ptr)^;
      end
    else
      begin
      inverted  := 0;
      exponent_ := -exponent_;
      Numer.lsw := 0;
      Denom.lsw := 0;
      XSIG_LL(Numer)^ := significand(st1_ptr)^;
      XSIG_LL(Denom)^ := significand(st0_ptr)^;
      end;
    div_Xsig(@Numer, @Denom, @argSignif);
    exponent_ := exponent_ + norm_Xsig(@argSignif);

    if ((exponent_ >= -1) or ((exponent_ = -2) and
      (argSignif.msw > $d413ccd0))) then
      begin
      (* The argument is greater than sqrt(2)-1 (:=0.414213562...) *)
      (* Convert the argument by an identity for atan *)
      transformed := 1;

      if (exponent_ >= 0) then
        begin
{$ifdef PARANOID}
	  if ( ( (exponent_ = 0) and
		 (argSignif.lsw = 0) @ and (argSignif.midw = 0) @@
		 (argSignif.msw = $80000000) ) )= false then
	    begin
	      EXCEPTION(EX_INTERNAL|$104);  (* There must be a logic error *)
	      exit;
	    end;
{$endif}(* PARANOID *)
        argSignif.msw := 0;   (* Make the transformed arg ^. 0.0 *)
        end
      else
        begin
        Denom.lsw := argSignif.lsw;
        Numer.lsw := Denom.lsw;
        XSIG_LL(Denom)^ := XSIG_LL(argSignif)^;
        XSIG_LL(Numer)^ := XSIG_LL(Denom)^;

        if (exponent_ < -1) then
          shr_Xsig(@Numer, -1 - exponent_);
        negate_Xsig(@Numer);

        shr_Xsig(@Denom, -exponent_);
        Denom.msw := Denom.msw or $80000000;

        div_Xsig(@Numer, @Denom, @argSignif);

        exponent_ := -1 + norm_Xsig(@argSignif);
        end;
      end
    else
      begin
      transformed := 0;
      end;

    argSq.lsw  := argSignif.lsw;
    argSq.midw := argSignif.midw;
    argSq.msw  := argSignif.msw;
    mul_Xsig_Xsig(@argSq, @argSq);

    argSqSq.lsw  := argSq.lsw;
    argSqSq.midw := argSq.midw;
    argSqSq.msw  := argSq.msw;
    mul_Xsig_Xsig(@argSqSq, @argSqSq);

    accumulatore.lsw := argSq.lsw;
    XSIG_LL(accumulatore)^ := XSIG_LL(argSq)^;

    shr_Xsig(@argSq, 2 * (-1 - exponent_ - 1));
    shr_Xsig(@argSqSq, 4 * (-1 - exponent_ - 1));

  (* Now have argSq etc with binary point at the left
     .1xxxxxxxx *)

    (* Do the basic fixed point polynomial evaluation *)
    accumulator.msw  := 0;
    accumulator.midw := 0;
    accumulator.lsw  := 0;
    polynomial_Xsig(@accumulator, XSIG_LL(argSqSq), @oddplterms, HIPOWERop - 1);
    mul64_Xsig(@accumulator, XSIG_LL(argSq));
    negate_Xsig(@accumulator);
    polynomial_Xsig(@accumulator, XSIG_LL(argSqSq), @oddnegterms, HIPOWERon - 1);
    negate_Xsig(@accumulator);
    add_two_Xsig(@accumulator, @fixedpterm, @dummy_exp);

    mul64_Xsig(@accumulatore, @denomterm);
    shr_Xsig(@accumulatore, 1 + 2 * (-1 - exponent_));
    accumulatore.msw := accumulatore.msw or $80000000;

    div_Xsig(@accumulator, @accumulatore, @accumulator);

    mul_Xsig_Xsig(@accumulator, @argSignif);
    mul_Xsig_Xsig(@accumulator, @argSq);

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

    if (transformed) <> 0 then
      begin
      (* compute pi/4 - accumulator *)
      shr_Xsig(@accumulator, -1 - exponent_);
      negate_Xsig(@accumulator);
      add_Xsig_Xsig(@accumulator, @pi_signif);
      exponent_ := -1;
      end;

    if (inverted) <> 0 then
      begin
      (* compute pi/2 - accumulator *)
      shr_Xsig(@accumulator, -exponent_);
      negate_Xsig(@accumulator);
      add_Xsig_Xsig(@accumulator, @pi_signif);
      exponent_ := 0;
      end;

    if (sign1) <> 0 then
      begin
      (* compute pi - accumulator *)
      shr_Xsig(@accumulator, 1 - exponent_);
      negate_Xsig(@accumulator);
      add_Xsig_Xsig(@accumulator, @pi_signif);
      exponent_ := 1;
      end;

    exponent_ := exponent_ + round_Xsig(@accumulator);

    significand(st1_ptr)^ := XSIG_LL(accumulator)^;
    setexponent16(st1_ptr, exponent_);

    tag := FPU_round(st1_ptr, 1, 0, FULL_PRECISION, sign2);
    FPU_settagi(1, tag);


    set_precision_flag_up();  (* We do not really know if up or down,
             use this as the default. *)

  end;

end.
