(380 . 11)(380 . 10)
               -- Multiply, give bottom and top halves
            when '*' =>
               Want(2);
               -- Ch9: Comba's algorithm
               FZ_Mul_Comba(X     => Stack(SP - 1),
                            Y     => Stack(SP),
                            XY_Lo => Stack(SP - 1),
                            XY_Hi => Stack(SP));
               FZ_Mult(X     => Stack(SP - 1),
                       Y     => Stack(SP),
                       XY_Lo => Stack(SP - 1),
                       XY_Hi => Stack(SP));
               
               -- Modular Multiplication
            when 'M' =>

(19 . 9)(19 . 52)

with Word_Ops; use Word_Ops;


-- Fundamental Arithmetic operators on FZ:
package body FZ_Arith is
   
   -- Destructive Add: X := X + Y; Overflow := Carry; optional OF_In
   procedure FZ_Add_D(X          : in out FZ;
                      Y          : in     FZ;
                      Overflow   : out    WBool;
                      OF_In      : in     WBool := 0) is
      Carry : WBool := OF_In;
   begin
      for i in 0 .. Word_Index(X'Length - 1) loop
         declare
            A : constant Word := X(X'First + i);
            B : constant Word := Y(Y'First + i);
            S : constant Word := A + B + Carry;
         begin
            X(X'First + i) := S;
            Carry          := W_Carry(A, B, S);
         end;
      end loop;
      Overflow := Carry;
   end FZ_Add_D;
   pragma Inline_Always(FZ_Add_D);
   
   
   -- Destructive Add: X := X + W; Overflow := Carry
   procedure FZ_Add_D_W(X        : in out FZ;
                        W        : in     Word;
                        Overflow : out    WBool) is
      Carry : Word := W;
   begin
      for i in X'Range loop
         declare
            A : constant Word := X(I);
            S : constant Word := A + Carry;
         begin
            X(i)  := S;
            Carry := W_Carry(A, 0, S);
         end;
      end loop;
      Overflow := Carry;
   end FZ_Add_D_W;
   pragma Inline_Always(FZ_Add_D_W);

   
   -- Sum := X + Y; Overflow := Carry
   procedure FZ_Add(X          : in  FZ;
                    Y          : in  FZ;
(29 . 13)(72 . 13)
                    Overflow   : out WBool) is
      Carry : WBool := 0;
   begin
      for i in X'Range loop
      for i in 0 .. Word_Index(X'Length - 1) loop
         declare
            A : constant Word := X(I);
            B : constant Word := Y(I);
            A : constant Word := X(X'First + i);
            B : constant Word := Y(Y'First + i);
            S : constant Word := A + B + Carry;
         begin
            Sum(i) := S;
            Sum(Sum'First + i) := S;
            Carry  := W_Carry(A, B, S);
         end;
      end loop;
(103 . 4)(146 . 49)
   end FZ_Sub;
   pragma Inline_Always(FZ_Sub);
   
   
   -- Destructive: If Cond is 1, NotN := ~N; otherwise NotN := N.
   procedure FZ_Not_Cond_D(N    : in out FZ;
                           Cond : in     WBool)is
      
      -- The inversion mask
      Inv : constant Word := 0 - Cond;
      
   begin
      
      for i in N'Range loop
         
         -- Invert (or, if Cond is 0, do nothing)
         N(i) := N(i) xor Inv;
         
      end loop;
      
   end FZ_Not_Cond_D;
   pragma Inline_Always(FZ_Not_Cond_D);
   
   
   -- Subtractor that gets absolute value if underflowed, in const. time
   procedure FZ_Sub_Abs(X          : in FZ;
                        Y          : in FZ;
                        Difference : out FZ;
                        Underflow  : out WBool) is
      
      O : Word := 0;
      pragma Unreferenced(O);
      
   begin
      
      -- First, we subtract normally
      FZ_Sub(X, Y, Difference, Underflow);
      
      -- If borrow - negate,
      FZ_Not_Cond_D(Difference, Underflow);
      
      -- ... and also increment.
      FZ_Add_D_W(Difference, Underflow, O);
      
   end FZ_Sub_Abs;
   pragma Inline_Always(FZ_Sub_Abs);
   
   
end FZ_Arith;

(20 . 11)(20 . 24)
with Words;   use Words;
with FZ_Type; use FZ_Type;


-- Fundamental Arithmetic operators on FZ:
package FZ_Arith is
   
   pragma Pure;
   
   -- Destructive Add: X := X + Y; Overflow := Carry; optional OF_In
   procedure FZ_Add_D(X          : in out FZ;
                      Y          : in     FZ;
                      Overflow   : out    WBool;
                      OF_In      : in     WBool := 0);
   pragma Precondition(X'Length = Y'Length);
   
      -- Destructive Add: X := X + W; Overflow := Carry
   procedure FZ_Add_D_W(X        : in out FZ;
                        W        : in     Word;
                        Overflow : out    WBool);
   
   -- Sum := X + Y; Overflow := Carry
   procedure FZ_Add(X          : in    FZ;
                    Y          : in    FZ;
(55 . 4)(68 . 15)
                    Underflow  : out WBool);
   pragma Precondition(X'Length = Y'Length and X'Length = Difference'Length);
   
   -- Destructive: If Cond is 1, NotN := ~N; otherwise NotN := N.
   procedure FZ_Not_Cond_D(N    : in out FZ;
                           Cond : in     WBool);
   
   -- Subtractor that gets absolute value if underflowed, in const. time
   procedure FZ_Sub_Abs(X          : in FZ;
                        Y          : in FZ;
                        Difference : out FZ;
                        Underflow  : out WBool);
   pragma Precondition(X'Length = Y'Length and X'Length = Difference'Length);
   
end FZ_Arith;

(45 . 7)(45 . 7)
   begin
      
      -- XY_Lo:XY_Hi := X * Y
      FZ_Mul_Comba(X, Y, XY_Lo, XY_Hi);
      FZ_Mult(X, Y, XY_Lo, XY_Hi);
      
      -- Product := XY mod M
      FZ_Mod(XY, Modulus, Product);

(20 . 15)(20 . 15)
with Words;    use Words;
with Word_Ops; use Word_Ops;
with W_Mul;    use W_Mul;
with FZ_Arith; use FZ_Arith;


package body FZ_Mul is
   
   -- Comba's multiplier.
      
   -- Comba's multiplier. (CAUTION: UNBUFFERED)
   procedure FZ_Mul_Comba(X     : in  FZ;
                          Y     : in  FZ;
                          XY_Lo : out FZ;
                          XY_Hi : out FZ) is
                          XY    : out FZ) is
      
      -- Words in each multiplicand
      L : constant Word_Index := X'Length;
(39 . 11)(39 . 8)
      -- 3-word Accumulator
      A2, A1, A0 : Word := 0;
      
      -- Register holding Product; indexed from zero
      XY : FZ(0 .. LP - 1);
      
      -- Type for referring to a column of XY
      subtype ColXY is Word_Index range XY'Range;
      subtype ColXY is Word_Index range 0 .. LP - 1;
      
      -- Compute the Nth (indexed from zero) column of the Product
      procedure Col(N : in ColXY; U : in ColXY; V : in ColXY) is
(86 . 7)(83 . 7)
         end loop;
         
         -- We now have the Nth (indexed from zero) word of XY
         XY(N) := A0;
         XY(XY'First + N) := A0;
         
         -- Right-Shift the Accumulator by one Word width:
         A0 := A1;
(115 . 11)(112 . 163)
      -- The very last Word of the Product:
      XY(XY'Last) := A0;
      
      -- Output the Product's lower and upper FZs:
      XY_Lo := XY(0 .. L - 1);
      XY_Hi := XY(L .. XY'Last);
      
   end FZ_Mul_Comba;
   pragma Inline_Always(FZ_Mul_Comba);
   
   
   -- Karatsuba's Multiplier. (CAUTION: UNBUFFERED)
   procedure Mul_Karatsuba(X  : in  FZ;
                           Y  : in  FZ;
                           XY : out FZ) is
      
      -- L is the wordness of a multiplicand. Guaranteed to be a power of two.
      L : constant Word_Count := X'Length;
      
      -- An 'LSeg' is the same length as either multiplicand.
      subtype LSeg is FZ(1 .. L);
      
      -- K is HALF of the length of a multiplicand.
      K : constant Word_Index := L / 2;
      
      -- A 'KSeg' is the same length as HALF of a multiplicand.
      subtype KSeg is FZ(1 .. K);
      
      -- The three L-sized variables of the product equation, i.e.:
      -- XY = LL + 2^(K*Bitness)(LL + HH + (-1^DD_Sub)*DD) + 2^(2*K*Bitness)HH
      LL, DD, HH : LSeg;
      
      -- K-sized terms of Dx * Dy = DD
      Dx, Dy     : KSeg;  -- Dx = abs(XLo - XHi) , Dy = abs(YLo - YHi)
      
      -- Subtraction borrows, signs of (XL - XH) and (YL - YH),
      Cx, Cy     : WBool; -- so that we can calculate (-1^DD_Sub)
      
      -- Bottom and Top K-sized halves of the multiplicand X.
      XLo        : KSeg  renames    X(  X'First       ..   X'Last - K );
      XHi        : KSeg  renames    X(  X'First + K   ..   X'Last     );
      
      -- Bottom and Top K-sized halves of the multiplicand Y.
      YLo        : KSeg  renames    Y(  Y'First       ..   Y'Last - K );
      YHi        : KSeg  renames    Y(  Y'First + K   ..   Y'Last     );
      
      -- L-sized middle segment of the product XY (+/- K from the midpoint).
      XYMid      : LSeg  renames   XY( XY'First + K   ..  XY'Last - K );
      
      -- Bottom and Top L-sized halves of the product XY.
      XYLo       : LSeg  renames   XY( XY'First       ..  XY'Last - L );
      XYHi       : LSeg  renames   XY( XY'First + L   ..  XY'Last     );
      
      -- Topmost K-sized quarter segment of the product XY, or 'tail'
      XYHiHi     : KSeg  renames XYHi( XYHi'First + K .. XYHi'Last    );
      
      -- Whether the DD term is being subtracted.
      DD_Sub     : WBool;
      
      -- Carry from individual term additions.
      C          : WBool;
      
      -- Tail-Carry accumulator, for the final ripple
      TC         : Word;
      
   begin
      
      -- Recurse: LL := XL * YL
      FZ_Multiply(XLo, YLo, LL);
      
      -- Recurse: HH := XH * YH
      FZ_Multiply(XHi, YHi, HH);
      
      -- Dx := |XL - XH| , Cx := Borrow (i.e. 1 iff XL < XH)
      FZ_Sub_Abs(X => XLo, Y => XHi, Difference => Dx, Underflow => Cx);
      
      -- Dy := |YL - YH| , Cy := Borrow (i.e. 1 iff YL < YH)
      FZ_Sub_Abs(X => YLo, Y => YHi, Difference => Dy, Underflow => Cy);
      
      -- Recurse: DD := Dx * Dy
      FZ_Multiply(Dx, Dy, DD);
      
      -- Whether (XL - XH)(YL - YH) is positive, and so DD must be subtracted:
      DD_Sub := 1 - (Cx xor Cy);
      
      -- XY := LL + 2^(2 * K * Bitness) * HH
      XYLo := LL;
      XYHi := HH;
      
      -- XY += 2^(K * Bitness) * HH, but carry goes in Tail Carry accum.
      FZ_Add_D(X => XYMid, Y => HH, Overflow => TC);
      
      -- XY += 2^(K * Bitness) * LL, ...
      FZ_Add_D(X => XYMid, Y => LL, Overflow => C);
      
      -- ... but the carry goes into the Tail Carry accumulator.
      TC := TC + C;
      
      -- XY += 2^(K * Bitness) * (-1^DD_Sub) * DD
      FZ_Not_Cond_D(N => DD, Cond => DD_Sub); -- invert DD if 2s-complementing
      FZ_Add_D(OF_In    => DD_Sub,            -- ... and then must increment
               X        => XYMid,
               Y        => DD,
               Overflow => C); -- carry will go in Tail Carry accumulator
      
      -- Compute the final Tail Carry for the ripple
      TC := TC + C - DD_Sub;
      
      -- Barring a cosmic ray, 0 <= TC <= 2 .
      pragma Assert(TC <= 2);
      
      -- Ripple the Tail Carry into the tail.
      FZ_Add_D_W(X => XYHiHi, W => TC, Overflow => C);
      
      -- Barring a cosmic ray, the tail ripple will NOT overflow.
      pragma Assert(C = 0);

   end Mul_Karatsuba;
   -- CAUTION: Inlining prohibited for Mul_Karatsuba !
   
   
   -- Multiplier. (CAUTION: UNBUFFERED)
   procedure FZ_Multiply(X     : in  FZ;
                         Y     : in  FZ;
                         XY    : out FZ) is
      
      -- The length of either multiplicand
      L : constant Word_Count := X'Length;
      
   begin
      
      if L <= Karatsuba_Thresh then
         
         -- Base case:
         FZ_Mul_Comba(X, Y, XY);
         
      else
         
         -- Recursive case:
         Mul_Karatsuba(X, Y, XY);
         
      end if;
      
   end FZ_Multiply;
   pragma Inline_Always(FZ_Multiply);
   
   
   -- Multiplier. Preserves the inputs.
   procedure FZ_Mult(X     : in  FZ;
                     Y     : in  FZ;
                     XY_Lo : out FZ;
                     XY_Hi : out FZ) is
      
      -- Product buffer.
      P : FZ(1 .. 2 * X'Length);
      
   begin
      
      FZ_Multiply(X, Y, P);
      
      XY_Lo := P(P'First            .. P'First + X'Length - 1);
      XY_Hi := P(P'First + X'Length .. P'Last);
      
   end FZ_Mult;
   pragma Inline_Always(FZ_Mult);
   
end FZ_Mul;

(24 . 11)(24 . 37)
   
   pragma Pure;
   
   -- Comba's multiplier.
   -- Karatsuba Threshhold - at or below this many words, we use Comba mult.
   Karatsuba_Thresh : constant Indices := 8;
   
   -- Multiply. (CAUTION: UNBUFFERED)
   procedure FZ_Multiply(X     : in  FZ;
                         Y     : in  FZ;
                         XY    : out FZ);
   pragma Precondition(X'Length = Y'Length and
                         XY'Length = (X'Length + Y'Length));
   
   -- Comba's multiplier. (CAUTION: UNBUFFERED)
   procedure FZ_Mul_Comba(X     : in  FZ;
                          Y     : in  FZ;
                          XY_Lo : out FZ;
                          XY_Hi : out FZ);
                          XY    : out FZ);
   pragma Precondition(X'Length = Y'Length and
                         XY'Length = (X'Length + Y'Length));
   
   -- Karatsuba's Multiplier. (CAUTION: UNBUFFERED)
   procedure Mul_Karatsuba(X  : in  FZ;
                           Y  : in  FZ;
                           XY : out FZ);
   pragma Precondition(X'Length = Y'Length and
                         XY'Length = (X'Length + Y'Length) and
                         X'Length mod 2 = 0);
   -- CAUTION: Inlining prohibited for Mul_Karatsuba !
   
   -- Multiplier. Preserves the inputs.
   procedure FZ_Mult(X     : in  FZ;
                     Y     : in  FZ;
                     XY_Lo : out FZ;
                     XY_Hi : out FZ);
   pragma Precondition(X'Length = Y'Length and
                         XY_Lo'Length = XY_Hi'Length and
                         XY_Lo'Length = ((X'Length + Y'Length) / 2));

(22 . 8)(22 . 16)

package body W_Mul is
   
   function Mul_HalfWord_Iron(X : in HalfWord;
                              Y : in HalfWord) return Word is
   begin
      return X * Y;
   end Mul_HalfWord_Iron;
   pragma Inline_Always(Mul_HalfWord_Iron);
   
   
   -- Multiply half-words X and Y, producing a Word-sized product
   function Mul_HalfWord(X : in HalfWord; Y : in HalfWord) return Word is
   function Mul_HalfWord_Soft(X : in HalfWord; Y : in HalfWord) return Word is
      
      -- X-Slide
      XS : Word := X;
(68 . 8)(76 . 8)
      -- Return the Product
      return XY;
      
   end Mul_HalfWord;
   pragma Inline_Always(Mul_HalfWord);
   end Mul_HalfWord_Soft;
   pragma Inline_Always(Mul_HalfWord_Soft);
   
   
   -- Get the bottom half of a Word
(107 . 16)(115 . 16)
      YH : constant HalfWord := TopHW(Y);
      
      -- XL * YL
      LL : constant Word := Mul_HalfWord(XL, YL);
      LL : constant Word := Mul_HalfWord_Iron(XL, YL);
      
      -- XL * YH
      LH : constant Word := Mul_HalfWord(XL, YH);
      LH : constant Word := Mul_HalfWord_Iron(XL, YH);
      
      -- XH * YL
      HL : constant Word := Mul_HalfWord(XH, YL);
      HL : constant Word := Mul_HalfWord_Iron(XH, YL);
      
      -- XH * YH
      HH : constant Word := Mul_HalfWord(XH, YH);
      HH : constant Word := Mul_HalfWord_Iron(XH, YH);
      
      -- Carry
      CL : constant Word := TopHW(TopHW(LL) + BottomHW(LH) + BottomHW(HL));

(31 . 8)(31 . 11)
   -- The number of bytes in a Half-Word
   HalfByteness : constant Positive := Byteness / 2;
   
   -- Multiply half-words X and Y, producing a Word-sized product
   function Mul_HalfWord(X : in HalfWord; Y : in HalfWord) return Word;
   -- Multiply half-words X and Y, producing a Word-sized product (Iron)
   function Mul_HalfWord_Iron(X : in HalfWord; Y : in HalfWord) return Word;
   
   -- Multiply half-words X and Y, producing a Word-sized product (Egyptian)
   function Mul_HalfWord_Soft(X : in HalfWord; Y : in HalfWord) return Word;
   
   -- Get the bottom half of a Word
   function BottomHW(W : in Word) return HalfWord;
(40 . 7)(43 . 7)
   -- Get the top half of a Word
   function TopHW(W : in Word) return HalfWord;
   
   -- Carry out X*Y mult, return lower word XY_LW and upper word XY_HW.
   -- Carry out X*Y mult, return lower word XY_LW and upper word XY_HW (Iron)
   procedure Mul_Word(X : in Word; Y : in Word; XY_LW : out Word; XY_HW : out Word);
   
end W_Mul;