unit polynom_xsig;

interface

uses fpu_types;

procedure polynomial_Xsig (accum: PXsig; const x: pu64; const terms: parr64;
  const n: intg);

implementation

procedure polynomial_Xsig (accum: PXsig; const x: pu64; const terms: parr64;
  const n: intg);
  var
    i: intg;
    acc, Xprod: Xsig;
    lprod: u32;
    xlwr, xupr, prod: u64;
    overflowed: byte;
  begin

    xlwr := u32(x^);
    xupr := u32(x^ shr 32);

    acc.msw  := terms^[n] shr 32;
    acc.midw := terms^[n];
    acc.lsw  := 0;
    overflowed := 0;

    i := n - 1;
    while i >= 0 do
      begin
      (* Split the product into five parts to get a 16 byte result *)

      (* first word by first word *)
      prod := acc.msw * xupr;
      Xprod.midw := prod;
      Xprod.msw := prod shr 32;

      (* first word by second word *)
      prod  := acc.msw * xlwr;
      Xprod.lsw := prod;
      lprod := prod shr 32;
      Xprod.midw := Xprod.midw + lprod;
      if (lprod > Xprod.midw) then
        Inc(Xprod.msw);

      (* second word by first word *)
      prod := acc.midw * xupr;
      Xprod.lsw := Xprod.lsw + prod;
      if (u32(prod) > Xprod.lsw) then
        begin
        Inc(Xprod.midw);
        if (Xprod.midw = 0) then
          Inc(Xprod.msw);
        end;
      lprod := prod shr 32;
      Xprod.midw := Xprod.midw + lprod;
      if (lprod > Xprod.midw) then
        Inc(Xprod.msw);

      (* second word by second word *)
      prod  := acc.midw * xlwr;
      lprod := prod shr 32;
      Xprod.lsw := Xprod.lsw + lprod;
      if (lprod > Xprod.lsw) then
        begin
        Inc(Xprod.midw);
        if (Xprod.midw = 0) then
          Inc(Xprod.msw);
        end;

      (* third word by first word *)
      prod  := acc.lsw * xupr;
      lprod := prod shr 32;
      Xprod.lsw := Xprod.lsw + lprod;
      if (lprod > Xprod.lsw) then
        begin
        Inc(Xprod.midw);
        if (Xprod.midw = 0) then
          Inc(Xprod.msw);
        end;

      if (overflowed) <> 0 then
        begin
        Xprod.midw := Xprod.midw + xlwr;
        if (u32(xlwr) > Xprod.midw) then
          Inc(Xprod.msw);
        Xprod.msw  := Xprod.msw + xupr;
        overflowed := 0;    (* We don't check this addition for overflow *)
        Dec(i);
        end;

      acc.lsw  := Xprod.lsw;
      acc.midw := u32(terms^[i]) + Xprod.midw;
      acc.msw  := (terms^[i] shr 32) + Xprod.msw;
      if (Xprod.msw > acc.msw) then
        overflowed := 1;
      if (u32(terms^[i]) > acc.midw) then
        begin
        Inc(acc.msw);
        if (acc.msw = 0) then
          overflowed := 1;
        end;
      end;

    (* We don't check the addition to accum for overflow *)
    accum^.lsw := accum^.lsw + acc.lsw;
    if (acc.lsw > accum^.lsw) then
      begin
      Inc(accum^.midw);
      if (accum^.midw = 0) then
        Inc(accum^.msw);
      end;
    accum^.midw := accum^.midw + acc.midw;
    if (acc.midw > accum^.midw) then
      begin
      Inc(accum^.msw);
      end;
    accum^.msw := accum^.msw + acc.msw;
  end;


end.
