UNSAFE MODULE MPInt;

IMPORT Word, Ctypes, MPNat, Math, Text, Scan, Fmt, Lex, FloatMode;
IMPORT Wr,Thread, Stdio, OldFmt;

CONST MP_NEG = -1;
      MP_ZERO= 0;
      MP_POS = 1;

TYPE DigitT = Word.T;

(* REVEAL T = BRANDED REF ARRAY OF MPNat.mp_limbT;   *)

REVEAL T = BRANDED REF ARRAY OF DigitT;
          (* the first position in this array will be used to define 
             the signal of the number, that is:
                  z[0]=MP_ZERO=> Number is zero
                  z[0]=MP_NEG => Number is negative
                  z[0]=MP_POS => Number is positive *)



(* ----------- Internal Procedures -------------- *)

PROCEDURE FirstPosition(z : T) : CARDINAL =
  BEGIN
    RETURN FIRST(z^)+1
  END FirstPosition;

PROCEDURE IntToNat(z : T) : MPNat.T =
  BEGIN
    RETURN LOOPHOLE(ADR(z[FirstPosition(z)]),MPNat.T)
  END IntToNat;

PROCEDURE GetSize(z : T) : CARDINAL =
  BEGIN
    RETURN NUMBER(z^)-1
  END GetSize;

PROCEDURE GetDigit(z : T; i: CARDINAL) : DigitT =
  BEGIN
    RETURN z[FirstPosition(z)+i-1]
  END GetDigit;

PROCEDURE SetDigit(z : T; i: CARDINAL; v: DigitT) =
  BEGIN
    IF i <= GetSize(z) THEN 
      z[FirstPosition(z)+i-1]:=v
    END 
  END SetDigit;

PROCEDURE SetSign(z : T; s: Sign) =
  (* sets the sign of the number represented by z. *)
  BEGIN
(*
Wr.PutText(Stdio.stdout,"EM MPInt.SetSign --- FIRST(z^) = "&Fmt.Int(FIRST(z^))&" z = "&OldFmt.Ref(z)&"\n"); 
Wr.Flush(Stdio.stdout);
*)
    z[FIRST(z^)] := s+1
  END SetSign;

PROCEDURE Assign(x,z: T; s: SizeT) =
  (* Copies the "s" digits from "z" to "x". *)
  BEGIN
    WITH zp = FirstPosition(z) DO

(*
Wr.PutText(Stdio.stdout," zp = "&Fmt.Int(zp)&"   LAST(z) = "&Fmt.Int(LAST(z^))&"\n");
Wr.Flush(Stdio.stdout);

FOR i:=0 TO LAST(z^) DO 
Wr.PutText(Stdio.stdout," z["&Fmt.Int(i)&"] = "&Fmt.Unsigned(z^[i])&"\n");
END;
Wr.Flush(Stdio.stdout);
*)


      SUBARRAY(x^,zp,s):=SUBARRAY(z^,zp,s)
    END;
    SetSign(x,GetSign(z));
  END Assign;

PROCEDURE Num_IsZero(num: T): BOOLEAN =
  VAR b : BOOLEAN:=TRUE;
      i : CARDINAL;
  BEGIN
    WITH
      nl= LAST(num^),
      np= FirstPosition(num)
    DO
      i:=nl;
      WHILE (i>=np) AND b DO 
        IF num[i] # 0 THEN b:=FALSE END;
        i:=i-1
      END;
      RETURN b
    END
  END Num_IsZero;

PROCEDURE DigitValueInBase(c: CHAR; base: BaseT:=10): CHAR =
(* Converts the char "c" to the char format used by MPN *)
  VAR rc : CHAR;
  BEGIN
    IF (c>='0') AND (c<='9') THEN 
      rc:=VAL(ORD(c)-ORD('0'),CHAR);
    ELSIF (c>='a') AND (c<='z') THEN 
      rc:=VAL(ORD(c)-ORD('a')+10,CHAR);
    ELSIF (c>='A') AND (c<='Z') THEN 
      rc:=VAL(ORD(c)-ORD('A')+10,CHAR);
    ELSE
      RETURN '\n'
    END;

    IF ORD(rc) < base THEN 
      RETURN rc
    ELSE
      RETURN '\n'
    END
  END DigitValueInBase;

PROCEDURE CountLeadingZerosInDigit(d : DigitT) : Ctypes.unsigned_long_int =
  VAR numzeros : Ctypes.unsigned_long_int;
      tp       : DigitT;
  BEGIN
    WITH
      (* Rotate(x,n) returns "i" bit equals "(i-n) MOD Word.Size" *) 
      msb = Word.Rotate(1,Word.Size-1)
    DO
      numzeros:=0;
      tp := d;
      WHILE (Word.And(tp,msb) = 0) AND (numzeros < Word.Size)  DO 
        tp:=Word.LeftShift(tp,1);
        numzeros:=numzeros+1
      END;
    END;
    RETURN numzeros
  END CountLeadingZerosInDigit;
  

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

PROCEDURE GetSign(z : T) : Sign =
  BEGIN
    RETURN z[FIRST(z^)]-1
  END GetSign;
  
PROCEDURE IsZero(z : T) : BOOLEAN =
  BEGIN
    RETURN GetSign(z) = MP_ZERO;
  END IsZero;

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

PROCEDURE NumDigits(z : T) : SizeT =
  (* Returns the number of digits IN "z", not considering leading 
     zeros. *)
  VAR s : CARDINAL; 
  BEGIN 
    IF IsZero(z) THEN 
      RETURN 1
    ELSE
      WITH
        zp=FirstPosition(z),
        zl=GetSize(z)
      DO
        s:=zl;
        WHILE z[s]=0 AND s > zp  DO s:=s-1 END;
        RETURN s
      END
    END
  END NumDigits;

PROCEDURE Copy(z : T) : T =
  BEGIN
    WITH
      zl = NumDigits(z),
      y  = Create(zl)
    DO
      Assign(y,z,zl);
      IF NumDigits(y)=1 AND GetDigit(y,1) = 0 THEN 
        SetSign(y,MP_ZERO)
      END;
      RETURN y
    END;          
  END Copy;

PROCEDURE Neg(z : T) : T =
  BEGIN
    WITH y = Copy(z) DO
      SetSign(y,-GetSign(z));
      RETURN y
    END
  END Neg;

PROCEDURE Abs(z : T) : T =
  BEGIN
    WITH y = Copy(z) DO 
      SetSign(y,ABS(GetSign(z)));
      RETURN y
    END
  END Abs;
  
PROCEDURE Comp(y,z : T) : Sign =
  BEGIN
    WITH
      sy = GetSign(y),
      sz = GetSign(z),
      d  = sy-sz
    DO
      IF d # 0 THEN 
        RETURN d DIV ABS(d)
      ELSE
        WITH
          ysize = NumDigits(y),
          zsize = NumDigits(z),
          dsize = ysize - zsize
        DO
          IF ysize # zsize THEN
            RETURN dsize DIV ABS(dsize)
          ELSE
            RETURN MPNat.Cmp(IntToNat(y),IntToNat(z),ysize)
          END
        END
      END
    END
  END Comp;

PROCEDURE CompAbs(y,z : T) : Sign =
  BEGIN
    RETURN Comp(Abs(y),Abs(z))
  END CompAbs;
  
PROCEDURE Add(y,z : T) : T =
  VAR n : T;
      c : DigitT;
  BEGIN
    WITH
      sy = GetSign(y),
      sz = GetSign(z),
      yp = IntToNat(y),
      zp = IntToNat(z),
      yl = NumDigits(y),
      zl = NumDigits(z)
    DO
      IF sy = sz THEN 
        (* Adds magnitudes if signs are the same  *)
        CASE Comp(y,z) OF 
        | +1 => BEGIN (* y >  z  *) 
                  n:=Create(yl+1);
                  c:=MPNat.Add(IntToNat(n),yp,yl,zp,zl);
                  SetDigit(n,yl+1,c)
                END;
        |  0 => BEGIN (* y = z. To improve speed, we use "lshift" *)
                  n:=Create(yl+1);
                  c:=MPNat.LShift(IntToNat(n),yp,yl,1) 
                END;
        | -1 => BEGIN (* y < z *)
                  n:=Create(zl+1);
                  c:=MPNat.Add(IntToNat(n),zp,zl,yp,yl);
                  SetDigit(n,zl+1,c)
                END
        END;
        SetSign(n,sy)
      ELSE (* Subtract magnitude if signs are different *)
        CASE CompAbs(y,z) OF
        |  0 => BEGIN (* y = -z => n = 0 *)
                  n:=Create(1);
                END;
        | +1 => BEGIN (* |y| > |z| *)
                  n:=Create(yl);
                  c:=MPNat.Sub(IntToNat(n),yp,yl,zp,zl);
                  SetSign(n,sy);
                END;
        | -1 => BEGIN (* |y| < |z| *)
                  n:=Create(zl);
                  c:=MPNat.Sub(IntToNat(n),zp,zl,yp,yl);
                  SetSign(n,sz)
                END;
        ELSE
        END
      END
    END;
    RETURN n
  END Add;
  
PROCEDURE Sub(y,z : T) : T =
  BEGIN
    IF IsZero(z) THEN 
      RETURN Copy(y)
    ELSE
      RETURN Add(y,Neg(z));
    END
  END Sub;

PROCEDURE Mult(y,z : T) : T =
  VAR n : T;
      c : DigitT;
  BEGIN
    WITH
      yl=NumDigits(y),
      zl=NumDigits(z)
    DO
      n:=Create(yl+zl);
      IF yl >= zl THEN 
        c:=MPNat.Mul(IntToNat(n),IntToNat(y),yl,IntToNat(z),zl);
      ELSE
        c:=MPNat.Mul(IntToNat(n),IntToNat(z),zl,IntToNat(y),yl);
      END
    END; 
    SetSign(n,GetSign(y)*GetSign(z));
    RETURN n
  END Mult;

PROCEDURE Divide(n,d : T; VAR q,r : T) : Sign RAISES{Thread.Alerted, Wr.Failure} =
  VAR nsize,dsize,qsize,rsize : SizeT;
      nsign,dsign,qsign       : Sign;
      np,dp,qp,rp             : MPNat.T;
      norm_steps              : Ctypes.unsigned_long_int;
  BEGIN
    dsign:=GetSign(d);

    IF dsign = MP_ZERO THEN 
      Wr.PutText(Stdio.stdout,"******* Division by ZERO ******\n");
      Wr.Flush(Stdio.stdout);
      WITH c=0 DO 
        EVAL 1 DIV c  (* to abort the program *)
      END
    END;

    nsign:=GetSign(n);
    
    CASE CompAbs(n,d) OF 
    | -1 => BEGIN (* |n| < |d| *)
              q:=Create(1);   (* creates q=0 *)
              r:=Copy(n);
              RETURN MP_ZERO;
            END;
    |  0 => BEGIN  (* |n| = |d| *)
              q:=FromInteger(1);
              SetSign(q,nsign*dsign);
              r:=Create(1);   (* creates r=0 *)
              RETURN MP_POS
            END;
    | +1 => BEGIN  (* |n| > |d| *)
              nsize:=NumDigits(n);
              dsize:=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);
              np:=IntToNat(n);
              qp:=IntToNat(q);

              IF dsize = 1 THEN 
                 WITH
                   ddgt = GetDigit(d,dsize),
                   rdgt = MPNat.DivMod_1(qp,np,nsize,ddgt)
                 DO
                   r:=Create(dsize);               
                   SetDigit(r,dsize,rdgt);
                 END
              ELSE
                dp:=IntToNat(d);
                WITH
                  ddgt = GetDigit(d,dsize)
                DO
                  norm_steps := CountLeadingZerosInDigit(ddgt);
                END;
                 
                (* Normalize the denominator, i.e. make its most
                   significant bit set by shifting it "norm_steps"
                   bits to the left.  Also shift the numerator the
                   same number of steps (to keep the quotient the
                   same!).  *)

                IF norm_steps # 0 THEN 
                  (* Shift up the denominator setting the most
	             significant bit OF the most significant word.
	             Use temporary storage not to clobber the original
	             contents of the denominator.  *)

                  WITH
                    t = Create(dsize),
                    tp = IntToNat(t),
                    dgt = MPNat.LShift(tp,dp,dsize,norm_steps)  
                  DO
                    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);
                  rp:=IntToNat(r);
                  WITH 
                    dgt = MPNat.LShift(rp,np,nsize,norm_steps)
                  DO 
                    SetSign(r,nsign);
                    IF dgt # 0 THEN 
                      WITH 
                        newr = Create(nsize+1)
                      DO 
                        Assign(newr,r,nsize);
                        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;

                rp:=IntToNat(r);
                dp:=IntToNat(d);
                rsize:=NumDigits(r);
                qsize:=rsize-dsize;
                WITH
                  qdgt = MPNat.DivMod(qp,0,rp,rsize,dp,dsize)
                DO
                  SetSign(q,qsign);
                  (* Includes "qdgt" in "q" *)
                  IF qdgt # 0 THEN 
                    WITH
                      newq = Create(qsize+1)
                    DO
                      Assign(newq,q,qsize);
                      SetDigit(newq,qsize+1,qdgt);
                      q:=newq
                    END 
                  END;
                END; 
                (* Desnormalizes the remainder "r" *)
                WITH
                  tr = Create(dsize)
                DO
                  Assign(tr,r,dsize);
                  (* Copies the value of "tr" to "r" eliminating
                     leading zeros. *)
                  r:=Copy(tr);
                END;

                IF (norm_steps # 0) AND (NOT IsZero(r)) THEN 
                  rp:=IntToNat(r);
                  rsize:=NumDigits(r);
                  WITH
                    rdgt = MPNat.RShift(rp,rp,rsize,norm_steps)
                  DO
                    (* Eliminates leading zeros introduced by "rshift" *)
                    r:=Copy(r)
                  END;
                END;
              END; 
              SetSign(q,qsign);
              IF Num_IsZero(r) THEN 
                SetSign(r,MP_ZERO)
              ELSE 
                SetSign(r,nsign);
              END;
              RETURN qsign           
            END;   
    END
  END Divide;

PROCEDURE Div(y,z : T) : T RAISES{Thread.Alerted, Wr.Failure} =
  VAR q,r : T;
  BEGIN
    WITH s = Divide(y,z,q,r) DO 
      RETURN q
    END
  END Div;
  
PROCEDURE Mod(y,z : T) : T RAISES{Thread.Alerted, Wr.Failure} =
  VAR q,r : T;
  BEGIN
    WITH s = Divide(y,z,q,r) DO 
      RETURN r
    END
  END Mod;

PROCEDURE GCD(y,z : T) : T =
  VAR u,v,r,raux,foo : T;
      dgt   : DigitT;
      rl    : SizeT;
  BEGIN
    IF IsZero(y) THEN (* GCD(0,z) = z *) 
      RETURN Abs(z); 
    ELSIF IsZero(z) THEN (* GCD(y,0) = y *)
      RETURN Abs(y);
    ELSE
      (* Makes copy of y and z since mpn_gcd destroys its operands  *)
      u:=Copy(y);
      v:=Copy(z);
      
      WITH
        ul = NumDigits(u),
        vl = NumDigits(v),
        up = IntToNat(u),
        vp = IntToNat(v)
      DO
        IF ul = 1 THEN 
          r:=Create(ul);
          dgt := MPNat.GCD_1(vp,vl,GetDigit(u,ul));
          SetDigit(r,ul,dgt)
        ELSIF vl = 1  THEN
          r:=Create(vl);
          dgt := MPNat.GCD_1(up,ul,GetDigit(v,vl));
          SetDigit(r,vl,dgt)
        ELSE
          IF (ul > vl) OR (ul = vl AND CompAbs(u,z) >= 0)  THEN 
            raux:=Create(vl);   

            foo := Create(ul);

            rl:=MPNat.GCDExt(IntToNat(raux),IntToNat(foo),up,ul,vp,vl);  

(*
r:=FromString("777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777777");

raux:=FromString("3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333");

Wr.PutText(Stdio.stdout,"EM MPInt.GCD --- antes de MPNat.GCDExt raux = "&ToString(raux)&"\n");
Wr.Flush(Stdio.stdout);

Wr.PutText(Stdio.stdout," Antes raux[0] = "&Fmt.Int(FirstPosition(raux))&"   LAST(raux) = "&Fmt.Int(LAST(raux^))&"\n");
Wr.Flush(Stdio.stdout);


WITH 
(*
rpa0 = ADR(raux[0]),
vrp0 = LOOPHOLE(rpa0,MPNat.T),
rpa1 = ADR(raux[1]),
vrp1 = LOOPHOLE(rpa1,MPNat.T),
rpa2 = ADR(raux[2]),
vrp2 = LOOPHOLE(rpa2,MPNat.T),
rpa3 = ADR(raux[3]),
vrp3 = LOOPHOLE(rpa3,MPNat.T),
rpa4 = ADR(raux[4]),
vrp4 = LOOPHOLE(rpa4,MPNat.T)
*)
rpa0 = ADR(r[0]),
vrp0 = LOOPHOLE(rpa0,MPNat.T),
rpa1 = ADR(r[1]),
vrp1 = LOOPHOLE(rpa1,MPNat.T),
rpa2 = ADR(r[2]),
vrp2 = LOOPHOLE(rpa2,MPNat.T),
rpa3 = ADR(r[3]),
vrp3 = LOOPHOLE(rpa3,MPNat.T),
rpa4 = ADR(r[4]),
vrp4 = LOOPHOLE(rpa4,MPNat.T)

DO 
Wr.PutText(Stdio.stdout,"  vrp0 = "&OldFmt.Addr(vrp0)&"  vrp1 = "&OldFmt.Addr(vrp1)&"  vrp2 = "&OldFmt.Addr(vrp2)&"  vrp3 = "&OldFmt.Addr(vrp3)&"  vrp4 = "&OldFmt.Addr(vrp4)&"\n");

Wr.PutText(Stdio.stdout,"  raux[0] = "&Fmt.Unsigned(vrp0^)&"  raux[1] = "&Fmt.Unsigned(vrp1^)&"  raux[2] = "&Fmt.Unsigned(vrp2^)&"  raux[3] = "&Fmt.Unsigned(vrp3^)&"  raux[4] = "&Fmt.Unsigned(vrp4^)&"\n");
Wr.Flush(Stdio.stdout);

WITH rp = LOOPHOLE(ADR(raux[1]),MPNat.T) DO 

(*
Wr.PutText(Stdio.stdout," Antes raux = "&OldFmt.Ref(raux)&"  rp = "&OldFmt.Addr(rp)&"  up = "&OldFmt.Addr(up)&"  vp = "&OldFmt.Addr(vp)&"\n");
Wr.Flush(Stdio.stdout);
*)

rl:=MPNat.GCDExt(rp,IntToNat(foo),up,ul,vp,vl); 

(*            rl:=MPNat.GCDExt(IntToNat(raux),IntToNat(foo),up,ul,vp,vl);   *)

(*
Wr.PutText(Stdio.stdout," raux[0] = "&Fmt.Int(FirstPosition(raux))&"   LAST(raux) = "&Fmt.Int(LAST(raux^))&"\n");
Wr.Flush(Stdio.stdout);
*)

END;

(*
Wr.PutText(Stdio.stdout,"EM MPInt.GCD --- depois de MPNat.GCDExt raux = "&ToString(r)&"\n");
Wr.Flush(Stdio.stdout);
*)

Wr.PutText(Stdio.stdout,"  vrp0 = "&OldFmt.Addr(vrp0)&"  vrp1 = "&OldFmt.Addr(vrp1)&"  vrp2 = "&OldFmt.Addr(vrp2)&"  vrp3 = "&OldFmt.Addr(vrp3)&"  vrp4 = "&OldFmt.Addr(vrp4)&"\n");

Wr.PutText(Stdio.stdout,"  raux[0] = "&Fmt.Unsigned(vrp0^)&"  raux[1] = "&Fmt.Unsigned(vrp1^)&"  raux[2] = "&Fmt.Unsigned(vrp2^)&"  raux[3] = "&Fmt.Unsigned(vrp3^)&"  raux[4] = "&Fmt.Unsigned(vrp4^)&"\n");
Wr.Flush(Stdio.stdout);

END;
*)

(*
VAR rpa,fpa,upa,vpa : MPNat.T;
BEGIN
FOR i:=0 TO 55 DO 
rpa := LOOPHOLE(ADR(raux[i]),MPNat.T);
IF i <= ul THEN 
upa := LOOPHOLE(ADR(u[i]),MPNat.T);
END;
IF i <= vl THEN 
fpa := LOOPHOLE(ADR(foo[i]),MPNat.T);
vpa := LOOPHOLE(ADR(v[i]),MPNat.T);
END;
Wr.PutText(Stdio.stdout,Fmt.Int(i)&"  rpa = "&OldFmt.Addr(rpa)&"  rpa[i] = "&Fmt.Unsigned(raux^[i])&"  fpa = "&OldFmt.Addr(fpa)&"  upa = "&OldFmt.Addr(upa)&"  vpa = "&OldFmt.Addr(vpa)&"\n");
END;
END;
Wr.Flush(Stdio.stdout);

*)



          ELSE
            raux:=Create(ul);
            foo := Create(ul);

Wr.PutText(Stdio.stdout,"EM MPInt.GCD ul < vl --- antes de MPNat.GCDExt raux = "&OldFmt.Ref(raux)&"\n");
Wr.Flush(Stdio.stdout);
            
            rl:=MPNat.GCDExt(IntToNat(raux),IntToNat(foo),vp,vl,up,ul)
          END;
          
Wr.PutText(Stdio.stdout,"EM MPInt.GCD --- voltei de MPNat.GCDExt com rl = "&Fmt.Int(rl)&" raux = "&OldFmt.Ref(raux)&"\n");
Wr.Flush(Stdio.stdout);

          r:=Create(rl);

Wr.PutText(Stdio.stdout,"EM MPInt.GCD --- depois de Create(rl) \n");
Wr.Flush(Stdio.stdout);
          
          Assign(r,raux,rl);

Wr.PutText(Stdio.stdout,"EM MPInt.GCD --- depois de Assign(r,raux,rl) \n");
Wr.Flush(Stdio.stdout);

        END; 
        (* For all u,v => GCD(u,v) > 0 *)
        SetSign(r,MP_POS);

Wr.PutText(Stdio.stdout,"EM MPInt.GCD --- depois de SetSign(r,MP_POS) \n");
Wr.Flush(Stdio.stdout);

        RETURN Copy(r)   (* to normalize "r"  *) 
      END 
    END;
  END GCD;

PROCEDURE Sqrt(z : T; VAR r: T) : T RAISES {Thread.Alerted, Wr.Failure} =
  VAR zsign : Sign;
  BEGIN
    zsign:= GetSign(z);
    IF zsign = MP_NEG THEN 
      Wr.PutText(Stdio.stdout,"******* Sqrt of negative number ******\n");
      Wr.Flush(Stdio.stdout);
      WITH c=-1.0d0 DO 
        EVAL Math.sqrt(c);  (* to abort the program *)
        RETURN Create(1)
      END
    ELSIF zsign = MP_ZERO THEN 
      r:=Create(1);
      RETURN Create(1)  
    ELSE 
      WITH
        zl = NumDigits(z),
        zp = IntToNat(z),
        srsize = (zl DIV 2) + (zl MOD 2),
        sr  = Create(srsize),
        srp = IntToNat(sr),
        aux = Create(srsize),
        auxp= IntToNat(aux),
        rsize = MPNat.Sqrt(srp,auxp,zp,zl)       
      DO
        SetSign(sr,MP_POS);  (* Since z > 0 *)
        IF rsize > 0 THEN 
          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(z: T) : BOOLEAN =
  BEGIN
    IF GetSign(z) = MP_NEG THEN 
      RETURN FALSE 
    ELSE 
      WITH 
        zl = NumDigits(z),
        zp = IntToNat(z),
        r  = MPNat.PerfectSquare(zp,zl)
      DO  
        RETURN (r # 0)
      END 
    END 
  END IsPerfectSquare;

PROCEDURE ToString(z : T; base : BaseT:=10) : TEXT = 
CONST dig = ARRAY [0..15] OF TEXT {"0","1","2","3","4","5","6","7","8","9","A","B","C","D","E","F"};
  VAR zaux  : T;
      str   : TEXT;
      zl    : SizeT;
      sarr  : SizeT;
      zsign : Sign;
      i     : CARDINAL;
  BEGIN
    IF base < 2 OR base > 16 THEN RETURN NIL END;
    
    zsign:=GetSign(z);
    zl:=NumDigits(z);

    IF (zsign = MP_ZERO) OR (zl = 0) THEN RETURN "0" END;

    zaux:=Copy(z);

    WITH
      (* computes the size that the string must have 
         MPN requires 3 extras positions                   *)
         
      l2base= Math.log(2.0d0)/Math.log(FLOAT(base,LONGREAL)),
      MaxCharsPerLimb=CEILING(FLOAT(MPNat.BITS_PER_MP_LIMB,LONGREAL)*l2base),
      size_str = zl*MaxCharsPerLimb+3,
      vstr= NEW(UNTRACED REF ARRAY OF CHAR,size_str),
      arr = LOOPHOLE(ADR(vstr[0]),Ctypes.unsigned_char_star)
    DO
      sarr:=MPNat.getstr(arr,base,IntToNat(zaux),zl); 
      str:="";
      FOR i:=0 TO sarr-1 DO 
        str:=str&dig[ORD(vstr[i])]
      END; 

      (* eliminates leading zeros in str *)
      WITH
        strlen = Text.Length(str)
      DO
        i:=0;
        WHILE (i<strlen) AND (Text.GetChar(str,i) = '0') DO i:=i+1 END;
        IF i=strlen THEN 
          str:="0" 
        ELSE
          str:=Text.Sub(str,i,strlen);
        END
      END
    END;
    IF zsign < 0 THEN str:="-"&str END;
    RETURN str
  END ToString;

PROCEDURE FromString(s : TEXT; base : BaseT:=10) : T =
  TYPE pstrT= UNTRACED REF ARRAY OF CHAR;
  VAR z                  : T;
      zl,i,slen, strlen  : CARDINAL;
      str                : TEXT;
      vstr               : pstrT;
      sign               : Sign;
      nr                 : SizeT;
  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 (MPNat.DIGIT_SIZE,LONGREAL),
      zr = nr * Math.log(br) / (Math.log(2.0d0) * dgr) + 1.0d0
    DO
      zl:=CEILING(zr)
    END;

    z:=Create(zl);

    (* allocates an array of char overlapping str *)
    vstr:= NEW(pstrT,strlen);
    FOR i:=0 TO strlen-1 DO 
      vstr[i]:=DigitValueInBase(Text.GetChar(str,i),base)
    END;

    WITH
      arr  = LOOPHOLE(ADR(vstr[0]),Ctypes.unsigned_char_star)
    DO
      nr:=MPNat.setstr(IntToNat(z),arr,strlen,base)
    END; 
    
    (* Set sign of result  *)
    IF Num_IsZero(z)  THEN sign:=MP_ZERO END;
    SetSign(z,sign);

    RETURN 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[FirstPosition(z)]:=num;
        SetSign(z,sign);
        RETURN z
      END
    END;
  END FromInteger;

PROCEDURE ToInteger(z : T) : INTEGER =
  BEGIN
    IF NumDigits(z) > 1 THEN 
      RETURN LAST(INTEGER)
    END;
    
    IF GetSign(z) = MP_NEG THEN 
      RETURN -z[FirstPosition(z)]
    ELSE
      RETURN z[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),
        z  = Create(zl),
        zp = IntToNat(z),
        zs = l DIV zl
      DO
        MPNat.random2(zp,zl);
        IF Num_IsZero(z) THEN 
          SetSign(z,MP_ZERO)
        ELSE
          SetSign(z,zs);
        END;
        RETURN z
      END;
    END 
  END Random;

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