UNSAFE INTERFACE MPNat;
(* 
   Arithmetic on arbitrary-precision unsigned integers (``bignums'').
   
   This interface provides access to the low-level unsigned
   integer operations in the GNU Multi Precision library, 
   as defined in "gmp.h".
 
   The basic data type is a multiple-precision unsigned integer, or
   ``bignum'', consisting of one or more big digits, the ``bigdigs''.
   The addresses of each bigdig increase??? with the bigdig's 
   multiplier.

   Created Jun 03, 1997 by Marcus Vinicius A. Andrade.
   Last edited on Jun 20, 1998 by J. Stolfi.
*)

IMPORT Ctypes;

TYPE 
  mp_limbT       = Ctypes.unsigned_int;      (* Unsigned big digit. *)
  mp_ptr         = UNTRACED REF mp_limbT;    (* Pointer to LSD of a bignum. *)
  (*
    An "mp_ptr" is a pointer to the least significant bigdig of a bignum. *)

  mp_srcptr      = UNTRACED REF mp_limbT;
  mp_limbSignedT = Ctypes.INT;               (* Signed big digit. *)
  mp_sizeT       = Ctypes.long;              (* Number of bigdigs in a bignum. *)
  mp_expT        = Ctypes.long;              (* Exponent of a  power of two. *)
  ShiftAmount    = Ctypes.unsigned_long_int; (* Number of bits to shift. *)
  
CONST 
  BYTE_SIZE  = 8;
  WORD_SIZE  = BITSIZE(INTEGER);
  DIGIT_SIZE = BITSIZE(mp_limbT);
  MaxLimb    = LAST(mp_limbT);

(* --------------Constants defined in MPN ("gmp-mparam.h") ------------ *)

CONST 
  BITS_PER_MP_LIMB  = 32;
  BYTES_PER_MP_LIMB = BITS_PER_MP_LIMB DIV BYTE_SIZE;
  BITS_PER_LONGINT  = 32;
  BITS_PER_INT      = 32;
  BITS_PER_SHORTINT = 16;
  BITS_PER_CHAR     = 8;
  
(*
  In all the procedures below, the input arguments must have
  at least one bigdig.  Except where noted, the result area
  must either coincide with one or more of the operands,
  or be entirely disjoint from them. 
*)

(* ---------------------- Addition -------------------------- *)

                                 <* EXTERNAL "__mpn_add_n" *> 
PROCEDURE add_n(
    r: mp_ptr; 
    a, b: mp_ptr; 
    size: mp_sizeT;
  ): mp_limbT;  
  (*
    Adds "a[0..size-1]" to "b[0..size-1]", storing the 
    result in "r[0..size-1]". Returns the carry, which is 
    either 0 or 1. *)

                                 <* EXTERNAL "__mpn_add_1" *> 
PROCEDURE add_1(
    r: mp_ptr; 
    a: mp_ptr; size: mp_sizeT; 
    b: mp_limbT;
  ): mp_limbT;
  (*
    Adds the bigdig "b" to the bignum "a[0..size-1]", stores
    the result in "r[0..size-1]". Returns the carry, which is 
    either 0 or 1. *)

                                 <* EXTERNAL "__mpn_add" *>   
PROCEDURE add(
    r: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT; 
    b: mp_ptr; b_size: mp_sizeT;
  ): mp_limbT;
  (*
    Adds "a[0..a_size-1]" to "b[0..b_size-1]", stores the result in
    "r[0..a_size-1]". Returns the carry, which is either 0 or 1.
    Requires "a_size >= b_size". *)

(* ------------------ Subtraction -------------------------- *)

                                 <* EXTERNAL "__mpn_sub_n" *> 
PROCEDURE sub_n(
    r: mp_ptr; 
    a, b: mp_ptr; 
    size: mp_sizeT;
  ): mp_limbT; 
  (*
    Subtracts "b[0..size-1]" from "a[0..size-1]",
    storing the result in "r[0..size-1]". Returns the
    borrow, which is either 0 or 1. *)

                                 <* EXTERNAL "__mpn_sub_1" *> 
PROCEDURE sub_1(
    r: mp_ptr; 
    a: mp_ptr; size: mp_sizeT; 
    b: mp_limbT;
  ): mp_limbT;
  (*
    Subtracts the bigdig "b" from the bignum "a[0..size-1]", stores
    the result in "r[0..size-1]". Returns the borrow, which is
    either 0 or 1. *)

                                 <* EXTERNAL "__mpn_sub" *>   
PROCEDURE sub(
    r: mp_ptr;
    a: mp_ptr; a_size: mp_sizeT;
    b: mp_ptr; b_size: mp_sizeT
  ): mp_limbT;
  (*
    Subtracts "b[0..b_size-1]" from "a[0..a_size-1]", stores the result in
    "r[0..a_size-1]". Returns the borrow, which is either 0 or 1.
    Requires "a_size >= b_size". *)

(* ------------------ Multiplication ---------------------- *)

                                 <* EXTERNAL "__mpn_mul_n" *> 
PROCEDURE mul_n(
    r: mp_ptr; 
    a, b: mp_ptr; 
    size: mp_sizeT;
  );
  (*
    Multiplies "a[0..size-1]" by "b[0..size-1]", stores 
    the result in "r[0..2*size-1]". *)

                                 <* EXTERNAL "__mpn_mul_1" *> 
PROCEDURE mul_1(
    r: mp_ptr; 
    a: mp_ptr; size: mp_sizeT; 
    b: mp_limbT;
  ): mp_limbT; 
  (*
    Multiplies the bigdig "b" by the bignum "a[0..size-1]", 
    returns the result in "r[0..size-1]". Returns the 
    spillover, which is a bigdig in "0..MaxLimb-1". *)

                                 <* EXTERNAL "__mpn_addmul_1" *> 
PROCEDURE addmul_1(
    r: mp_ptr; 
    a: mp_ptr; size: mp_sizeT;
    b: mp_limbT;
  ): mp_limbT; 
  (*
    Mutiplies the bigdig "b" by the bignum "a[0..size-1]", 
    adds the product to "r[0..size-1]". Returns the spillover
    plus the addition carry, which is a bigdig in "0..MaxLimb". *)

                                 <* EXTERNAL "__mpn_submul_1" *> 
PROCEDURE submul_1(
    r: mp_ptr; a: mp_ptr; size: mp_sizeT; b: mp_limbT): mp_limbT;
  (*
    Multiplies the bigdig "b" by the bignum "a[0..size-1]", 
    subtracts the product from "r[0..size-1]". Returns the 
    spillover plus the borrow, which is a bigdig in "0..MaxLimb". *)

                                 <* EXTERNAL "__mpn_mul" *>   
PROCEDURE mul(
    r: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT;
    b: mp_ptr; b_size: mp_sizeT;
  ): mp_limbT;
  (*
    Multiplies the bignum "a[0..a_size-1]" by the bignum "b[0..b_size-1]", 
    stores the result in "r[0..a_size+b_size-1]". Returns the 
    most significant bigdig of the product. "r[a_size+b_size-1]". 
    
    Requires "a_size >= b_size". The result area must be disjoint from 
    both operands. *)

(* ------------------ Division ---------------------------- *)

                                 <* EXTERNAL "__mpn_divrem" *> 
PROCEDURE divrem(
    q: mp_ptr; f_size: mp_sizeT;
    a: mp_ptr; a_size: mp_sizeT; 
    b: mp_ptr; b_size: mp_sizeT;
  ): mp_limbT;
  (*
    Divides the bignum "a[0..a_size-1]" by the bignum "b[0..b_size-1]", 
    with "f_size" fractional digits; stores the quotient 
    (minus its most significant limb) in "q[0..q_size+f_size-1]"
    where "q_size = ???". 
    
    Returns the most significant bigdig of the quotient, which
    is either 0 or 1. Leaves the remainder of the division
    in "a[0..a_size-1]".
    
    The bigdigs "q[0..q_size-1]???" are the integer part
    of the quotient (without the MSB), and "q[q_size..q_size+f_size-1]???"
    are the fractional part. 
    
    Requires "a_size >= b_size". The most significant bit of "a" must
    be 1.  If the quotient is not needed, pass "q = a + b_size".
    Aside from this special case, no overlap between arguments is
    permitted. *)

                                 <* EXTERNAL "__mpn_divrem_1" *> 
PROCEDURE div_1(
    q: mp_ptr; f_size: mp_sizeT;
    a: mp_ptr; a_size: mp_sizeT; 
    b: mp_limbT;
  ): mp_limbT;
  (*
    Divides the bignum "a[0..a_size-1]" by the bigdig "b",
    with "f_size" fraction bigdigs.  Stores the result in
    "q[0..???-1]" where ???. Returns the remainder, which
    is a bigdig in the range "0..b-1". 
    
    The bigdigs "q[0..q_size-1]???" are the integer part
    of the quotient (without the MSB), and "q[q_size..q_size+f_size-1]???"
    are the fractional part. *)

                                 <* EXTERNAL "__mpn_mod_1" *> 
PROCEDURE mod_1(
    a: mp_ptr; 
    a_size: mp_sizeT; 
    b: mp_limbT;
  ): mp_limbT;
  (*
    Divides the bignum "a[0..a_size-1]" by the bigdig "b", 
    returns the remainder, which is a bigdig in the range "0..b-1". *)

                                 <* EXTERNAL "__mpn_divmod_1" *> 
PROCEDURE divmod_1(
    r: mp_ptr;  
    a: mp_ptr; a_size: mp_sizeT; 
    d: mp_limbT;
  ): mp_limbT;
  (*
    Divides the bignum "a[0..a_size-1]" by the bigdig "d", 
    stores the quotient in "q[0..???-1]" where ???, and
    returns the remainder, which is a bigdig in "[0..b-1]". *)

(* ------------------ Shift ------------------------------- *)

                                 <* EXTERNAL "__mpn_lshift" *> 
PROCEDURE lshift(
    r: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT; 
    d: ShiftAmount;  (* Number of bits to shift, in "1..DIGIT_SIZE-1". *)
  ): mp_limbT;
  (*
    Shifts the bignum "a[0..a_size-1]" left by "d" bits, that is,
    multiplies it by "2^d". Stores the result in "r[0..a_size-1]",
    and returns the spill-over, which is a bigdig.
    
    Overlapping of "a[0..a_size-1]" and "r[0..a_size-1]" is
    allowed in this function, provided "r >= a". *)

                                 <* EXTERNAL "__mpn_rshift" *> 
PROCEDURE rshift(
    r: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT; 
    d: ShiftAmount;
  ): mp_limbT;
  (*
    Shifts the bignum "a[0..a_size-1]"  right by "d" bits, that is,
    divides it it by "2^d". Stores the result in "r[0..a_size-1]",
    and returns the spill-over, which is a bigdig. 
    
    Overlapping of "a[0..a_size-1]" and "r[0..a_size-1]" is
    allowed in this function, provided "r <= a". *)

(* ------------------ Comparison ------------------------- *)

                                 <* EXTERNAL "__mpn_cmp" *> 
PROCEDURE cmp(
    a, b: mp_ptr; size: mp_sizeT;
  ): Ctypes.int;
  (*
    Numerically compares the bignums "a[0..size-1]" AND
    "b[0..size-1]". Returns a positive value if "a > b", 0 of they are equal,
    and a negative value if "a < b". *)
   
(* ------------------ GCD --------------------------------- *)

                                 <* EXTERNAL "__mpn_gcd" *> 
PROCEDURE gcd(
    r: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT; 
    b: mp_ptr; b_size: mp_sizeT;
  ): mp_limbT;
  (*
    Computes the greatest common divisor of the bignums 
    "a[0..a_size-1]" and "b[0..b_size-1]".  Stores the 
    result in "r[0..r_size-1]", where "r_size" is the 
    returned value of the call.  
    
    The bignum "a[0..a_size-1]" must be odd, and "b[0..b_size-1]" must
    have at least as many bits as "a[0..a_size-1]".  Both operands are
    destroyed by this call. *)

                                 <* EXTERNAL "__mpn_gcd_1" *> 
PROCEDURE gcd_1(
    a: mp_ptr; a_size: mp_sizeT; 
    b: mp_limbT;
  ): mp_limbT;
  (*
    Returns the greatest common divisor of the bignum 
    "a[0..a_size-1]" and the (non-zero) bigdig "b".  
    The result lies in the ramge "0..b". *)

                                 <* EXTERNAL "__mpn_gcdext" *> 
PROCEDURE gcdext(
    r1, r2: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT; 
    b: mp_ptr; b_size: mp_sizeT;
  ): mp_sizeT;
  (*
    Computes the greatest common divisor of the bignums 
    "a[0..a_size-1]" and "b[0..b_size-1]".  Stores ???
    in "r1[0..???-1]" and ??? in "r2[0..??-1]", where ???.
    Returns ???. *)

(* ------------------ Sqrt -------------------------------- *)

                                 <* EXTERNAL "__mpn_sqrtrem" *> 
PROCEDURE sqrtrem(
    r1, r2: mp_ptr; 
    a: mp_ptr; a_size: mp_sizeT;
  ): mp_sizeT;
  (*
    Computes the square root of the bignum 
    "a[0..a_size-1]", rounded downwards.  Stores ???
    in "r1[0..???-1]" and ??? in "r2[0..??-1]", where ???.
    Returns ???. *)
   
(* ------------------ Input/Output ------------------------ *)

                                 <* EXTERNAL "__mpn_get_str" *> 
PROCEDURE get_str(
    str: Ctypes.unsigned_char_star; 
    base: Ctypes.int; 
    s: mp_ptr; size: mp_sizeT;
  ): mp_sizeT;
  (*
    Converts a C character string "str"
    to a bignum.  Each character 
    of "str" must be an ASCII digit in the given "base" which must be ???.
    The string should be terminated 
    by the NUL character. The result is returned in "s[0..size-1]"
    where size is ???. Returns ???. *)

                                 <* EXTERNAL "__mpn_set_str" *> 
PROCEDURE set_str(
    r: mp_ptr; 
    str: Ctypes.unsigned_char_star; strsize: mp_sizeT; 
    base: Ctypes.int;
  ): mp_sizeT;
  (*
    Converts a bignum "r[0..???-1]" into a C character string,
    its representation in the given "base" (which must be ???).  
    The string "str" is assumed to have enough space ????
    and will be terminated by a NUL character.  Returns ???. *)

                                 <* EXTERNAL "__mpn_scan0" *> 
PROCEDURE scan0(
    r: mp_ptr; 
    bit: Ctypes.unsigned_long_int;
  ): Ctypes.unsigned_long_int; 
  (*
    ??? *)
    
                                 <* EXTERNAL "__mpn_scan1" *> 
PROCEDURE scan1(
    r: mp_ptr; 
    bit: Ctypes.unsigned_long_int;
  ): Ctypes.unsigned_long_int; 
  (*
    ??? *)

                                 <* EXTERNAL "__mpn_random2" *> 
PROCEDURE random2(
    r: mp_ptr; 
    rsize: mp_sizeT;
  );
  (*
    Fills "r[0..rsize-1]" with a random unsigned??? number, uniformly
    distributed between 0??? and the maximum representable number. *)

                                 <* EXTERNAL "__mpn_popcount" *> 
PROCEDURE popcount(
    r: mp_ptr; 
    size: Ctypes.unsigned_long_int;
  ): Ctypes.unsigned_long_int;
  (*
    ??? *)

                                 <* EXTERNAL "__mpn_hamdist" *> 
PROCEDURE hamdist(
    r, s: mp_ptr; 
    size: Ctypes.unsigned_long_INT;
  ): Ctypes.unsigned_long_int;
  (*
    Computes the Hamming distance between two unsigned???
    bignums "r[0..size-1]" and "s[0..size-1]". *)

(* ------------------ Perfect Square --------------------- *)

                                 <* EXTERNAL "__mpn_perfect_square_p" *> 
PROCEDURE perfect_square_p(
    s: mp_ptr; size: mp_sizeT;
  ): Ctypes.int;
  (*
    Returns 1 if the bignum "s[0..size-1]" is a perfect square, 
    and 0 otherwise. *)

END MPNat.