MODULE MPInt;

IMPORT MPZ, Math, Text, Scan, Fmt, Lex, FloatMode;

CONST MP_NEG = MPZ.MP_NEG;
      MP_ZERO= MPZ.MP_ZERO;
      MP_POS = MPZ.MP_POS;

TYPE NumDgt = MPZ.NumDgt;
     
(* ----------- Internal Procedures -------------- *)

PROCEDURE Slice(y,x: T; p,s: SizeT) =
  (* Copies to "y" the "s" digits of "x" starting at position "p", 
     that is, "y[yp .. yp+s-1] := x[xp+p ... xp+p+s-1], where "xp" 
     and "yp" are respectively the first position of "x" and "y". If 
     if "yp+p+s-1" is out of range then copies only the existent digits *)
  BEGIN
    WITH 
      xp = MPZ.FirstPosition(x^),
      yp = MPZ.FirstPosition(y^),
      ipx = xp+p-1
    DO
      IF ipx+s-1 > LAST(x^)  THEN s:=LAST(x^)-ipx+1 END;     
      SUBARRAY(y^,yp,s) := SUBARRAY(x^,ipx,s);
    END
  END Slice;

PROCEDURE Assign(y,x: T; s: SizeT) =
  (* Copies "s" digits from "x" to "y". *)
  BEGIN

    Slice(y,x,1,s);
    
    MPZ.SetSign(y^,MPZ.GetSign(x^));
  END Assign;

PROCEDURE DivRem(n,d : T; VAR q,r : T) = 
  (* computes "q" and "r" such that "n=qd+r" AND "0 <= |r| < |q|" *)
  VAR nsize,dsize,qsize,rsize : SizeT;
      nsign,dsign,qsign       : Sign;
      norm_steps              : NumDgt;
  BEGIN
    <* ASSERT NOT IsZero(d) *>

    dsign:=MPZ.GetSign(d^);
    nsign:=MPZ.GetSign(n^);
    
    CASE CompAbs(n,d) OF 
    | -1 => BEGIN (* |n| < |d| *)
              q:=Create(1);   (* creates q=0 *)
              r:=Copy(n);
            END;
    |  0 => BEGIN  (* |n| = |d| *)
              q:=FromInteger(1);
              MPZ.SetSign(q^,nsign*dsign);
              r:=Create(1);   (* creates r=0 *)
            END;
    | +1 => BEGIN  (* |n| > |d| *)
              nsize:=MPZ.NumDigits(n^);
              dsize:=MPZ.NumDigits(d^);

           (* Ensure space is enough for quotient and remainder.  We
              need space for an extra limb in the remainder, because
              it's up-shifted (normalized) below.  *)
              
              rsize:=nsize+1;
              qsize:=rsize-dsize;

              qsign:=nsign*dsign;
              q:=Create(qsize);

              IF dsize = 1 THEN 
                 (* since the denominator is an one digit number, call "divide" *)
                 r:=Create(dsize);               
                 EVAL MPZ.DivRem(q^,r^,n^,d^,nsize,dsize);
              ELSE
                (* Normalizes the denominator, i.e. shifts it to left 
                   until its most significant bit is 1. Also, shift the 
                   numerator the same number of steps (to keep the 
                   quotient the invariant).  *)

                norm_steps := MPZ.CountLeadingZerosInDigit(MPZ.GetDigit(d^,dsize));
                 
                IF norm_steps # 0 THEN 
                  (* Shift up the denominator setting the most
	             significant bit of the most significant limb.
	             Use temporary storage not to clobber the original
	             contents of the denominator.  *)

                  WITH
                    t = Create(dsize)
                  DO
                    EVAL MPZ.LShift(t^,d^,dsize,norm_steps);
                    MPZ.SetSign(t^,dsign); 
                    d:=t;
                  END;

                  (* Shift up the numerator, possibly introducing a
	             new most significant word.  Move the shifted
	             numerator in the remainder meanwhile.  *)
               
                  r := Create(nsize);
                  WITH 
                    dgt = MPZ.LShift(r^,n^,nsize,norm_steps)
                  DO 
                    MPZ.SetSign(r^,nsign);
                    IF dgt # 0 THEN 
                      WITH 
                        newr = Create(nsize+1)
                      DO 
                        Assign(newr,r,nsize);
                        MPZ.SetDigit(newr^,nsize+1,dgt);
                        r:=newr
                      END;
                    END;
                  END;  
                ELSE
                   (* The denominator is already normalized, as
	              required.  Move the numerator to the remainder. *)
                   r:=Copy(n)
                END;

                rsize:=MPZ.NumDigits(r^);
                qsize:=rsize-dsize;
                WITH
                  raux = Create(rsize),
                  qdgt = MPZ.DivRem(q^,raux^,r^,d^,rsize,dsize)
                DO
                  MPZ.SetSign(q^,qsign);
                  (* Includes "qdgt" in "q" *)
                  IF qdgt # 0 THEN 
                    WITH
                      newq = Create(qsize+1)
                    DO
                      Assign(newq,q,qsize);
                      MPZ.SetDigit(newq^,qsize+1,qdgt);
                      q:=newq
                    END 
                  END;

                  (* gets the piece of "raux" corresponding to the remainder *)
                  r:=Create(dsize);
                  Assign(r,raux,dsize)
                END; 

                
                (* Desnormalizes the remainder "r" *)
                
                IF norm_steps # 0 THEN 
                  rsize:=MPZ.NumDigits(r^);

                  WITH
                    raux = Create(rsize)
                  DO
                    EVAL MPZ.RShift(raux^,r^,rsize,norm_steps);
                    MPZ.SetSign(raux^,nsign);
                  
                    (* Eliminates leading zeros introduced by "rshift" *)
                    r:=Copy(raux)
                  END;
                END;
              END; 
              MPZ.SetSign(q^,qsign);
              IF MPZ.Num_IsZero(r^) THEN 
                MPZ.SetSign(r^,MP_ZERO)
              ELSE 
                MPZ.SetSign(r^,nsign);
              END
            END   
    END
  END DivRem;

PROCEDURE NormalizeToGCD(x: T; VAR x_zero_limbs, x_zero_bits: NumDgt) : T =
  (* puts "x" on the format required by "mpn_gcd" that is an odd number;
     returns "x" in this format and puts in "x_zero_limbs" the number of 
     limbs of "x" that are zero and in "x_zero_bits" the number of zero-bits 
     in the least significant limb different of zero. So, the number returned 
     corresponds to "x" right shifted "x_zero_limbs * SIZEOF(limb) + x_zero_bits" 
     positions *)
  VAR r : T;
  BEGIN
    x_zero_limbs := MPZ.CountLeastSignificantLimbZero(x^);
    WITH 
      xp = MPZ.FirstPosition(x^),
      xl = MPZ.GetSize(x^),
      rl = xl-x_zero_limbs,
      msl = MPZ.GetDigit(x^,xp + x_zero_limbs)
    DO 
      x_zero_bits := MPZ.CountTrailingZerosInDigit(msl);

      (* copies "x" to "r" excluding the trailing zero limbs *)
      IF x_zero_limbs > 0 THEN 
        r := Create(rl);
        Slice(r,x,xp+x_zero_limbs,xl);  
        MPZ.SetSign(r^,MPZ.GetSign(x^));
      ELSE 
        r := Copy(x)
      END;

      (* the first limb of "r" is "msl", that is a non zero limb;  *)
      IF x_zero_bits > 0 THEN 
        WITH
          raux = Create(rl)
        DO
          EVAL MPZ.RShift(raux^,r^,rl,x_zero_bits);
          MPZ.SetSign(raux^,MPZ.GetSign(r^));

          r:=Copy(raux);   (* to eliminate the most significant limb if it's zero *)
        END;
      END;
      RETURN r
    END;
  END NormalizeToGCD;

(* ----------- Procedures defined in the interface -------------- *)

PROCEDURE NumDigits(READONLY x : T) : SizeT = 
  BEGIN
    RETURN MPZ.NumDigits(x^)
  END NumDigits;

PROCEDURE IsZero(READONLY x: T) : BOOLEAN = 
  BEGIN
    RETURN MPZ.IsZero(x^)
  END IsZero;

PROCEDURE GetSign(READONLY x : T) : Sign = 
  BEGIN
    RETURN MPZ.GetSign(x^)
  END GetSign;

PROCEDURE Create(s : SizeT) : T =
  VAR z : T;
  BEGIN
    z := NEW(T,s+1);
    MPZ.SetSign(z^,MP_ZERO);
    RETURN z;
  END Create;

PROCEDURE Copy(x : T) : T =
  BEGIN
    WITH
      xl = MPZ.NumDigits(x^),
      y  = Create(xl)
    DO
      Assign(y,x,xl);
      IF MPZ.Num_IsZero(y^) THEN 
        MPZ.SetSign(y^,MP_ZERO)
      ELSE 
        MPZ.SetSign(y^,MPZ.GetSign(x^))
      END;
      RETURN y
    END;          
  END Copy;

PROCEDURE Neg(x : T) : T =
  BEGIN
    WITH y = Copy(x) DO
      MPZ.SetSign(y^,-MPZ.GetSign(x^));
      RETURN y
    END
  END Neg;

PROCEDURE Abs(x : T) : T =
  BEGIN
    WITH y = Copy(x) DO 
      MPZ.SetSign(y^,ABS(MPZ.GetSign(x^)));
      RETURN y
    END
  END Abs;
  
PROCEDURE Comp(x,y : T) : Sign =
  BEGIN
    WITH
      sx = MPZ.GetSign(x^),
      sy = MPZ.GetSign(y^),
      d  = sx-sy
    DO
      IF d # 0 THEN 
        RETURN d DIV ABS(d)
      ELSE
        WITH
          xsize = MPZ.NumDigits(x^),
          ysize = MPZ.NumDigits(y^),
          dsize = xsize - ysize
        DO
          IF xsize # ysize THEN
            RETURN dsize DIV ABS(dsize)
          ELSE
            RETURN MPZ.Comp(x^,y^,ysize)
          END
        END
      END
    END
  END Comp;

PROCEDURE CompAbs(x,y : T) : Sign =
  BEGIN
    RETURN Comp(Abs(x),Abs(y))
  END CompAbs;
  
PROCEDURE Add(x,y : T) : T =
  VAR z : T;
  BEGIN
    WITH
      xsign = MPZ.GetSign(x^),
      ysign = MPZ.GetSign(y^),
      xsize = MPZ.NumDigits(x^),
      ysize = MPZ.NumDigits(y^)
    DO
      IF IsZero(x) THEN 
        RETURN Copy(y)
      ELSIF IsZero(y) THEN 
        RETURN Copy(x)
      ELSIF xsign = ysign THEN (* same signals; so adds magnitudes *)
        IF Comp(x,y) >= 0 THEN (* x >=  y  *) 
          z:=Create(xsize+1);
          MPZ.Add(z^,x^,y^,xsize,ysize)
        ELSE (* x < y *)
          z:=Create(ysize+1);
          MPZ.Add(z^,y^,x^,ysize,xsize)
        END;
        MPZ.SetSign(z^,xsign)
      ELSE (* different signals; so subtracts magnitudes *) 
        CASE CompAbs(x,y) OF
        |  0 => BEGIN (* x = -y => z = 0 *)
                  z:=Create(1);
                END;
        | +1 => BEGIN (* |x| > |y| *)
                  z:=Create(xsize);
                  MPZ.Sub(z^,x^,y^,xsize,ysize);
                  MPZ.SetSign(z^,xsign);
                END;
        | -1 => BEGIN (* |x| < |y| *)
                  z:=Create(ysize);
                  MPZ.Sub(z^,y^,x^,ysize,xsize);
                  MPZ.SetSign(z^,ysign)
                END;
        ELSE
        END
      END
    END;
    RETURN z 
  END Add;
  
PROCEDURE Sub(x,y : T) : T =
  BEGIN
    IF MPZ.IsZero(y^) THEN 
      RETURN Copy(x)
    ELSIF MPZ.IsZero(x^) THEN 
      RETURN Copy(Neg(y))
    ELSE
      RETURN Add(x,Neg(y));
    END
  END Sub;

PROCEDURE Mult(x,y : T) : T =
  VAR z : T;
  BEGIN
    WITH
      xsize=MPZ.NumDigits(x^),
      ysize=MPZ.NumDigits(y^)
    DO
      z:=Create(xsize+ysize);
      IF xsize >= ysize THEN 
        MPZ.Mult(z^,x^,y^,xsize,ysize);
      ELSE
        MPZ.Mult(z^,y^,x^,ysize,xsize);
      END
    END; 
    RETURN z
  END Mult;


PROCEDURE Div(x,y : T) : T = 
  VAR q,r : T;
  BEGIN
    <* ASSERT NOT MPZ.IsZero(y^) *> 
    DivRem(x,y,q,r);
    RETURN q
  END Div;
  
PROCEDURE Mod(x,y : T) : T = 
  VAR q,r : T;
  BEGIN
    DivRem(x,y,q,r);
    RETURN r
  END Mod;

PROCEDURE GCD(x,y : T) : T = 
  VAR x_zero_limbs, 
      y_zero_limbs, 
      g_zero_limbs : NumDgt; (* number of least significant zero-limbs *)
      x_zero_bits, 
      y_zero_bits,  
      g_zero_bits  :  NumDgt; (* number of zero-bits in the least significant limb*)
      xsize, ysize, gsize : SizeT;
      z : T;
      
  BEGIN
    IF IsZero(x) THEN        (* gcd(0,y) = |y| *) 
      IF IsZero(y) THEN 
        RETURN DecDig[1];    (* gcd(0,0) = 1 *)
      ELSE 
        RETURN Abs(y); 
      END
    ELSIF IsZero(y) THEN     (* gcd(x,0) = |x| *)
      RETURN Abs(x);
    ELSE
      xsize := MPZ.NumDigits(x^);
      ysize := MPZ.NumDigits(y^);
      IF (xsize = 1) OR (ysize = 1) THEN 
        (* since one of the operands has size 1 then use the routine for 
           one digit operand *)

        z := Create(1);
        EVAL MPZ.GCD(z^,x^,y^,xsize,ysize);
        RETURN z
      ELSE 
        WITH 
          u = NormalizeToGCD(x,x_zero_limbs,x_zero_bits),
          v = NormalizeToGCD(y,y_zero_limbs,y_zero_bits),
          usize = MPZ.NumDigits(u^),
          vsize = MPZ.NumDigits(v^)
        DO
          (* sets "g_zero_limbs = MIN(x_zero_limbs,y_zero_limbs)" and 
             if "x_zero_limbs = y_zero_limbs" then sets "g_zero_bits" to 
             MIN(x_zero_bits,y_zero_bits) *)

          IF x_zero_limbs < y_zero_limbs THEN 
            g_zero_limbs := x_zero_limbs;
            g_zero_bits := x_zero_bits
          ELSIF x_zero_limbs > y_zero_limbs THEN 
            g_zero_limbs := y_zero_limbs;
            g_zero_bits := y_zero_bits
          ELSE 
            g_zero_limbs := x_zero_limbs;
            IF x_zero_bits <= y_zero_bits THEN 
              g_zero_bits := x_zero_bits
            ELSE 
              g_zero_bits := y_zero_bits
            END
          END;

          (* "mpn.gcd" requires that the first argument must NOT have 
              more bits then the second  *)
          
          IF (usize < vsize) OR (usize = vsize AND CompAbs(u,v) < 0) THEN 
            z:=Create(usize);            
            gsize:=MPZ.GCD(z^,u^,v^,usize,vsize)
          ELSE 
            z:=Create(vsize);            
            gsize:=MPZ.GCD(z^,v^,u^,vsize,usize)
          END;
          MPZ.SetSign(z^,MP_POS);
          
          (* Desnormalize the result, that is, shift to left the most significant 
             non zero limb "g_zero_bits" positions and inserts the zeros 
             corresponding to trailing zeros removed from the operands *)

          IF g_zero_bits # 0 THEN 
            WITH 
              zaux=Create(gsize),
              dgt = MPZ.LShift(zaux^,z^,gsize,g_zero_bits)
            DO 
              MPZ.SetSign(zaux^,MPZ.GetSign(z^));
              IF dgt # 0 THEN 
                gsize:=gsize+1;
                WITH 
                  zaux2 = Create(gsize)
                DO 
                  Assign(zaux2,zaux,gsize-1);
                  MPZ.SetDigit(zaux2^,gsize,dgt);
                  z:=Copy(zaux2)
                END;
              ELSE 
                z:=Copy(zaux)
              END
            END
          END;  

          IF g_zero_limbs # 0 THEN 
            gsize:=gsize+g_zero_limbs;
            WITH 
              zaux = Create(gsize),
              zsize = MPZ.NumDigits(z^)
            DO 
              MPZ.FillWithZero(zaux^,1,g_zero_limbs);
              FOR i:=1 TO zsize DO 
                zaux^[i+g_zero_limbs] := z^[i]
              END;
              MPZ.SetSign(zaux^,MPZ.GetSign(z^));
              z:=Copy(zaux)
            END;
          END;
          RETURN z
        END;
      END
    END
  END GCD;

PROCEDURE Sqrt(x : T; VAR r: T) : T = 
  VAR xsign : Sign;
  BEGIN
    xsign:= MPZ.GetSign(x^);
    <* ASSERT xsign # MP_NEG *>
    
    IF xsign = MP_ZERO THEN 
      r:=Create(1);
      RETURN Create(1)  
    ELSE 
      WITH
        xsize = MPZ.NumDigits(x^),
        srsize = (xsize DIV 2) + (xsize MOD 2),
        sr  = Create(srsize),
        aux = Create(xsize),
        rsize = MPZ.Sqrt(sr^,aux^,x^,xsize)       
      DO
        MPZ.SetSign(sr^,MP_POS);  (* Since z > 0 *)
        IF rsize > 0 THEN 
          MPZ.SetSign(aux^,MP_POS);
          r:= Create(rsize);              
          Assign(r,aux,rsize);
        ELSE
          r:=Create(1)    (* makes r=0 *)
        END;
        RETURN sr
      END
    END
  END Sqrt;

PROCEDURE IsPerfectSquare(x: T) : BOOLEAN =
  BEGIN
    RETURN MPZ.IsPerfectSquare(x^)
  END IsPerfectSquare;

PROCEDURE ToString(x : T; base : BaseT:=10) : TEXT = 
  BEGIN
    WITH aux = Copy(x) DO 
      RETURN MPZ.ToString(aux^,base)
    END
  END ToString;

PROCEDURE FromString(s : TEXT; base : BaseT:=10) : T =
  VAR z                  : T;
      zl,i,slen, strlen  : CARDINAL;
      str                : TEXT;
      sign               : Sign;
  BEGIN
    slen := Text.Length(s);
    (* Throw away any initial space *)
    i:=0;
    WHILE (i<slen) AND (Text.GetChar(s,i)=' ') DO i:=i+1 END;

    (* Set up sign, base, initialize result *)
    IF Text.GetChar(s,i)='-' THEN 
      sign:=MP_NEG; 
      i:=i+1
    ELSE
      sign:=MP_POS;
      IF Text.GetChar(s,i)='+' THEN i:=i+1 END
    END;

    (* copies to str the effective part of s   *)
    str:=Text.Sub(s,i,slen);
    strlen:=Text.Length(str);

    (* computes the size of MPNat must have to store str *)
    WITH
      br = FLOAT (base,LONGREAL),
      nr = FLOAT (strlen,LONGREAL),
      dgr= FLOAT (MPZ.Dgt_Size,LONGREAL),
      zr = nr * Math.log(br) / (Math.log(2.0d0) * dgr) + 1.0d0
    DO
      zl:=CEILING(zr)
    END;

    z:=Create(zl);

    MPZ.FromString(z^, str, strlen, base);
    
    (* Set sign of result  *)
    IF MPZ.Num_IsZero(z^)  THEN sign:=MP_ZERO END;
    MPZ.SetSign(z^,sign);
    
    RETURN Copy(z)   (* to normalize the "z" *)
  END FromString;

PROCEDURE FromInteger(i : INTEGER) : T =
  BEGIN
    IF i = 0 THEN 
      RETURN Create(1)
    ELSE 
      WITH
        num = ABS(i),
        sign = i DIV num,
        z = Create(1)
      DO
        z[MPZ.FirstPosition(z^)]:=num;
        MPZ.SetSign(z^,sign);
        RETURN z
      END
    END;
  END FromInteger;

PROCEDURE ToInteger(z : T) : INTEGER =
  BEGIN
    IF MPZ.NumDigits(z^) > 1 THEN 
      RETURN LAST(INTEGER)
    END;
    
    IF MPZ.GetSign(z^) = MP_NEG THEN 
      RETURN -z[MPZ.FirstPosition(z^)]
    ELSE
      RETURN z[MPZ.FirstPosition(z^)]
    END
  END ToInteger;

PROCEDURE Float(z : T) : LONGREAL =
<* FATAL FloatMode.Trap, Lex.Error *> 
  BEGIN
    WITH
      strz = ToString(z)
    DO
      RETURN Scan.LongReal(strz)
    END 
  END Float;

PROCEDURE Trunc(x : LONGREAL) : T =
  CONST Min = FLOAT(FIRST(INTEGER),LONGREAL);
        Max = FLOAT(LAST(INTEGER),LONGREAL);
  VAR strx : TEXT;
  BEGIN
    IF Min <= x AND x <= Max THEN 
      RETURN FromInteger(TRUNC(x))
    ELSE
      strx := Fmt.LongReal(x,Fmt.Style.Fix,0);
      RETURN FromString(strx)
    END
  END Trunc;

PROCEDURE Floor(x : LONGREAL) : T =
  BEGIN
    WITH 
      xi = Trunc(x)
    DO 
      IF x >= 0.0d0 THEN 
        RETURN xi
      ELSE
        RETURN Sub(xi,FromInteger(1))
      END
    END 
  END Floor;

PROCEDURE Ceiling(x : LONGREAL) : T =
  BEGIN
    WITH 
      xi = Trunc(x)
    DO 
      IF x <= 0.0d0 THEN 
        RETURN xi
      ELSE
        RETURN Add(xi,FromInteger(1))
      END
    END 
  END Ceiling;

PROCEDURE Round(x : LONGREAL) : T =
  BEGIN
    WITH 
      xi = Floor(x),
      d  = x - Float(xi)
    DO 
      IF d > 0.5d0 THEN 
        RETURN Add(xi,FromInteger(1))
      ELSE
        RETURN xi
      END
    END 
  END Round;

PROCEDURE FromLongReal(x : LONGREAL; eps : LONGREAL) : Pair =
  VAR inf   := Pair{n:=FromInteger(1),m:=FromInteger(0)};  
      sup   := Pair{n:=FromInteger(1),m:=FromInteger(1)};
      med   : Pair;
      fmed  : LONGREAL;
      xi    : T;
      xf    : LONGREAL;
      ok    : BOOLEAN := FALSE;
      d     : LONGREAL;
  BEGIN
    xi := Round(x);
    xf := ABS(x - Float(xi));
    IF (xf < eps) THEN 
      med.n:=FromInteger(1);
      med.m:=xi 
    ELSE
      xi := Trunc(x); 
      xf := ABS(x - Float(xi));
      IF x < 0.0d0 THEN xi:=Neg(xi) END;
      REPEAT
        med.n:=Add(inf.n,sup.n);
        med.m:=Add(inf.m,sup.m);
        fmed:=Float(med.m)/Float(med.n);
        d:=xf - fmed;
        IF (ABS(d) < eps) THEN 
          ok:=TRUE
        ELSIF d < 0.0d0 THEN
             sup:=med
           ELSE 
             inf:=med
        END;
      UNTIL ok; 
      med.m:=Add(med.m,Mult(xi,med.n));
      IF x < 0.0d0 THEN med.m:=Neg(med.m) END;
    END;
    RETURN med
  END FromLongReal;

PROCEDURE Random(l : INTEGER) : T =
  BEGIN
    IF l = 0 THEN 
      RETURN Create(1)
    ELSE
      WITH
        zl = ABS(l),
        zs = l DIV zl,
        z  = Create(zl)
      DO
        MPZ.Random(z^,zl,zs);
        RETURN z
      END;
    END 
  END Random;

BEGIN
  FOR i:=0 TO 9 DO 
    DecDig[i]:=FromInteger(i);
  END;
END MPInt.