Another round of libtommath const'ification. To be submitted to the libtommath folks

38 files changed, 105 insertions, 108 deletions

diff --git a/generic/tclTomMath.h b/generic/tclTomMath.h
index 14cf5b6..1f22d6f 100644
--- a/generic/tclTomMath.h
+++ b/generic/tclTomMath.h
@@ -497,42 +497,42 @@ int mp_mod_d(const mp_int *a, mp_digit b, mp_digit *c);
 
 /* d = a + b (mod c) */
 /*
-int mp_addmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_addmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 */
 
 /* d = a - b (mod c) */
 /*
-int mp_submod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_submod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 */
 
 /* d = a * b (mod c) */
 /*
-int mp_mulmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_mulmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 */
 
 /* c = a * a (mod b) */
 /*
-int mp_sqrmod(mp_int *a, mp_int *b, mp_int *c);
+int mp_sqrmod(const mp_int *a, const mp_int *b, mp_int *c);
 */
 
 /* c = 1/a (mod b) */
 /*
-int mp_invmod(mp_int *a, mp_int *b, mp_int *c);
+int mp_invmod(const mp_int *a, const mp_int *b, mp_int *c);
 */
 
 /* c = (a, b) */
 /*
-int mp_gcd(mp_int *a, mp_int *b, mp_int *c);
+int mp_gcd(const mp_int *a, const mp_int *b, mp_int *c);
 */
 
 /* produces value such that U1*a + U2*b = U3 */
 /*
-int mp_exteuclid(mp_int *a, mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3);
+int mp_exteuclid(const mp_int *a, const mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3);
 */
 
 /* c = [a, b] or (a*b)/(a, b) */
 /*
-int mp_lcm(mp_int *a, mp_int *b, mp_int *c);
+int mp_lcm(const mp_int *a, const mp_int *b, mp_int *c);
 */
 
 /* finds one of the b'th root of a, such that |c|**b <= |a|
@@ -540,10 +540,10 @@ int mp_lcm(mp_int *a, mp_int *b, mp_int *c);
  * returns error if a < 0 and b is even
  */
 /*
-int mp_n_root(mp_int *a, mp_digit b, mp_int *c);
+int mp_n_root(const mp_int *a, mp_digit b, mp_int *c);
 */
 /*
-int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast);
+int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast);
 */
 
 /* special sqrt algo */
@@ -553,22 +553,22 @@ int mp_sqrt(const mp_int *arg, mp_int *ret);
 
 /* special sqrt (mod prime) */
 /*
-int mp_sqrtmod_prime(mp_int *arg, mp_int *prime, mp_int *ret);
+int mp_sqrtmod_prime(const mp_int *arg, const mp_int *prime, mp_int *ret);
 */
 
 /* is number a square? */
 /*
-int mp_is_square(mp_int *arg, int *ret);
+int mp_is_square(const mp_int *arg, int *ret);
 */
 
 /* computes the jacobi c = (a | n) (or Legendre if b is prime)  */
 /*
-int mp_jacobi(mp_int *a, mp_int *n, int *c);
+int mp_jacobi(const mp_int *a, const mp_int *n, int *c);
 */
 
 /* used to setup the Barrett reduction for a given modulus b */
 /*
-int mp_reduce_setup(mp_int *a, mp_int *b);
+int mp_reduce_setup(mp_int *a, const mp_int *b);
 */
 
 /* Barrett Reduction, computes a (mod b) with a precomputed value c
@@ -577,74 +577,74 @@ int mp_reduce_setup(mp_int *a, mp_int *b);
  * compute the reduction as -1 * mp_reduce(mp_abs(a)) [pseudo code].
  */
 /*
-int mp_reduce(mp_int *a, mp_int *b, mp_int *c);
+int mp_reduce(mp_int *a, const mp_int *b, mp_int *c);
 */
 
 /* setups the montgomery reduction */
 /*
-int mp_montgomery_setup(mp_int *a, mp_digit *mp);
+int mp_montgomery_setup(const mp_int *a, mp_digit *mp);
 */
 
 /* computes a = B**n mod b without division or multiplication useful for
  * normalizing numbers in a Montgomery system.
  */
 /*
-int mp_montgomery_calc_normalization(mp_int *a, mp_int *b);
+int mp_montgomery_calc_normalization(mp_int *a, const mp_int *b);
 */
 
 /* computes x/R == x (mod N) via Montgomery Reduction */
 /*
-int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
+int mp_montgomery_reduce(mp_int *a, const mp_int *m, mp_digit mp);
 */
 
 /* returns 1 if a is a valid DR modulus */
 /*
-int mp_dr_is_modulus(mp_int *a);
+int mp_dr_is_modulus(const mp_int *a);
 */
 
 /* sets the value of "d" required for mp_dr_reduce */
 /*
-void mp_dr_setup(mp_int *a, mp_digit *d);
+void mp_dr_setup(const mp_int *a, mp_digit *d);
 */
 
 /* reduces a modulo b using the Diminished Radix method */
 /*
-int mp_dr_reduce(mp_int *a, mp_int *b, mp_digit mp);
+int mp_dr_reduce(mp_int *a, const mp_int *b, mp_digit mp);
 */
 
 /* returns true if a can be reduced with mp_reduce_2k */
 /*
-int mp_reduce_is_2k(mp_int *a);
+int mp_reduce_is_2k(const mp_int *a);
 */
 
 /* determines k value for 2k reduction */
 /*
-int mp_reduce_2k_setup(mp_int *a, mp_digit *d);
+int mp_reduce_2k_setup(const mp_int *a, mp_digit *d);
 */
 
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
 /*
-int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d);
+int mp_reduce_2k(mp_int *a, const mp_int *n, mp_digit d);
 */
 
 /* returns true if a can be reduced with mp_reduce_2k_l */
 /*
-int mp_reduce_is_2k_l(mp_int *a);
+int mp_reduce_is_2k_l(const mp_int *a);
 */
 
 /* determines k value for 2k reduction */
 /*
-int mp_reduce_2k_setup_l(mp_int *a, mp_int *d);
+int mp_reduce_2k_setup_l(const mp_int *a, mp_int *d);
 */
 
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
 /*
-int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d);
+int mp_reduce_2k_l(mp_int *a, const mp_int *n, mp_int *d);
 */
 
 /* d = a**b (mod c) */
 /*
-int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_exptmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 */
 
 /* ---> Primes <--- */
diff --git a/libtommath/bn_fast_mp_invmod.c b/libtommath/bn_fast_mp_invmod.c
index 7771136..08389dd 100644
--- a/libtommath/bn_fast_mp_invmod.c
+++ b/libtommath/bn_fast_mp_invmod.c
@@ -21,7 +21,7 @@
  * Based on slow invmod except this is optimized for the case where b is
  * odd as per HAC Note 14.64 on pp. 610
  */
-int fast_mp_invmod(mp_int *a, mp_int *b, mp_int *c)
+int fast_mp_invmod(const mp_int *a, const mp_int *b, mp_int *c)
 {
    mp_int  x, y, u, v, B, D;
    int     res, neg;
diff --git a/libtommath/bn_fast_mp_montgomery_reduce.c b/libtommath/bn_fast_mp_montgomery_reduce.c
index f2c38bf..54d9b0a 100644
--- a/libtommath/bn_fast_mp_montgomery_reduce.c
+++ b/libtommath/bn_fast_mp_montgomery_reduce.c
@@ -23,7 +23,7 @@
  *
  * Based on Algorithm 14.32 on pp.601 of HAC.
 */
-int fast_mp_montgomery_reduce(mp_int *x, mp_int *n, mp_digit rho)
+int fast_mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho)
 {
    int     ix, res, olduse;
    mp_word W[MP_WARRAY];
diff --git a/libtommath/bn_mp_abs.c b/libtommath/bn_mp_abs.c
index 343a102..03904d2 100644
--- a/libtommath/bn_mp_abs.c
+++ b/libtommath/bn_mp_abs.c
@@ -20,7 +20,7 @@
  * Simple function copies the input and fixes the sign to positive
  */
 int
-mp_abs(mp_int *a, mp_int *b)
+mp_abs(const mp_int *a, mp_int *b)
 {
    int     res;
 
diff --git a/libtommath/bn_mp_addmod.c b/libtommath/bn_mp_addmod.c
index 0521974..5aee233 100644
--- a/libtommath/bn_mp_addmod.c
+++ b/libtommath/bn_mp_addmod.c
@@ -17,7 +17,7 @@
 
 /* d = a + b (mod c) */
 int
-mp_addmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
+mp_addmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d)
 {
    int     res;
    mp_int  t;
diff --git a/libtommath/bn_mp_dr_is_modulus.c b/libtommath/bn_mp_dr_is_modulus.c
index 4631daa..bf4ed8b 100644
--- a/libtommath/bn_mp_dr_is_modulus.c
+++ b/libtommath/bn_mp_dr_is_modulus.c
@@ -16,7 +16,7 @@
  */
 
 /* determines if a number is a valid DR modulus */
-int mp_dr_is_modulus(mp_int *a)
+int mp_dr_is_modulus(const mp_int *a)
 {
    int ix;
 
diff --git a/libtommath/bn_mp_dr_reduce.c b/libtommath/bn_mp_dr_reduce.c
index e90dc55..17ab289 100644
--- a/libtommath/bn_mp_dr_reduce.c
+++ b/libtommath/bn_mp_dr_reduce.c
@@ -30,7 +30,7 @@
  * Input x must be in the range 0 <= x <= (n-1)**2
  */
 int
-mp_dr_reduce(mp_int *x, mp_int *n, mp_digit k)
+mp_dr_reduce(mp_int *x, const mp_int *n, mp_digit k)
 {
    int      err, i, m;
    mp_word  r;
diff --git a/libtommath/bn_mp_dr_setup.c b/libtommath/bn_mp_dr_setup.c
index b49c81a..5dff995 100644
--- a/libtommath/bn_mp_dr_setup.c
+++ b/libtommath/bn_mp_dr_setup.c
@@ -16,7 +16,7 @@
  */
 
 /* determines the setup value */
-void mp_dr_setup(mp_int *a, mp_digit *d)
+void mp_dr_setup(const mp_int *a, mp_digit *d)
 {
    /* the casts are required if DIGIT_BIT is one less than
     * the number of bits in a mp_digit [e.g. DIGIT_BIT==31]
diff --git a/libtommath/bn_mp_exptmod.c b/libtommath/bn_mp_exptmod.c
index c4f392b..934fd25 100644
--- a/libtommath/bn_mp_exptmod.c
+++ b/libtommath/bn_mp_exptmod.c
@@ -21,7 +21,7 @@
  * embedded in the normal function but that wasted alot of stack space
  * for nothing (since 99% of the time the Montgomery code would be called)
  */
-int mp_exptmod(mp_int *G, mp_int *X, mp_int *P, mp_int *Y)
+int mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y)
 {
    int dr;
 
diff --git a/libtommath/bn_mp_exptmod_fast.c b/libtommath/bn_mp_exptmod_fast.c
index 7278b9f..08c6bc3 100644
--- a/libtommath/bn_mp_exptmod_fast.c
+++ b/libtommath/bn_mp_exptmod_fast.c
@@ -29,7 +29,7 @@
 #define TAB_SIZE 256
 #endif
 
-int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode)
+int mp_exptmod_fast(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
 {
    mp_int  M[TAB_SIZE], res;
    mp_digit buf, mp;
@@ -39,7 +39,7 @@ int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode)
     * one of many reduction algorithms without modding the guts of
     * the code with if statements everywhere.
     */
-   int (*redux)(mp_int *,mp_int *,mp_digit);
+   int (*redux)(mp_int *,const mp_int *,mp_digit);
 
    /* find window size */
    x = mp_count_bits(X);
diff --git a/libtommath/bn_mp_exteuclid.c b/libtommath/bn_mp_exteuclid.c
index 419146b..adaea86 100644
--- a/libtommath/bn_mp_exteuclid.c
+++ b/libtommath/bn_mp_exteuclid.c
@@ -18,7 +18,7 @@
 /* Extended euclidean algorithm of (a, b) produces
    a*u1 + b*u2 = u3
  */
-int mp_exteuclid(mp_int *a, mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3)
+int mp_exteuclid(const mp_int *a, const mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3)
 {
    mp_int u1,u2,u3,v1,v2,v3,t1,t2,t3,q,tmp;
    int err;
diff --git a/libtommath/bn_mp_gcd.c b/libtommath/bn_mp_gcd.c
index 18f6dc3..f5aa78b 100644
--- a/libtommath/bn_mp_gcd.c
+++ b/libtommath/bn_mp_gcd.c
@@ -16,7 +16,7 @@
  */
 
 /* Greatest Common Divisor using the binary method */
-int mp_gcd(mp_int *a, mp_int *b, mp_int *c)
+int mp_gcd(const mp_int *a, const mp_int *b, mp_int *c)
 {
    mp_int  u, v;
    int     k, u_lsb, v_lsb, res;
diff --git a/libtommath/bn_mp_invmod.c b/libtommath/bn_mp_invmod.c
index b70fe18..525493a 100644
--- a/libtommath/bn_mp_invmod.c
+++ b/libtommath/bn_mp_invmod.c
@@ -16,7 +16,7 @@
  */
 
 /* hac 14.61, pp608 */
-int mp_invmod(mp_int *a, mp_int *b, mp_int *c)
+int mp_invmod(const mp_int *a, const mp_int *b, mp_int *c)
 {
    /* b cannot be negative */
    if ((b->sign == MP_NEG) || (mp_iszero(b) == MP_YES)) {
diff --git a/libtommath/bn_mp_invmod_slow.c b/libtommath/bn_mp_invmod_slow.c
index 2bdd2b1..2bb5196 100644
--- a/libtommath/bn_mp_invmod_slow.c
+++ b/libtommath/bn_mp_invmod_slow.c
@@ -16,7 +16,7 @@
  */
 
 /* hac 14.61, pp608 */
-int mp_invmod_slow(mp_int *a, mp_int *b, mp_int *c)
+int mp_invmod_slow(const mp_int *a, const mp_int *b, mp_int *c)
 {
    mp_int  x, y, u, v, A, B, C, D;
    int     res;
diff --git a/libtommath/bn_mp_is_square.c b/libtommath/bn_mp_is_square.c
index 84229bf..4d8612f 100644
--- a/libtommath/bn_mp_is_square.c
+++ b/libtommath/bn_mp_is_square.c
@@ -38,7 +38,7 @@ static const char rem_105[105] = {
 };
 
 /* Store non-zero to ret if arg is square, and zero if not */
-int mp_is_square(mp_int *arg,int *ret)
+int mp_is_square(const mp_int *arg,int *ret)
 {
    int           res;
    mp_digit      c;
diff --git a/libtommath/bn_mp_jacobi.c b/libtommath/bn_mp_jacobi.c
index 8981393..c314c82 100644
--- a/libtommath/bn_mp_jacobi.c
+++ b/libtommath/bn_mp_jacobi.c
@@ -20,7 +20,7 @@
  * HAC is wrong here, as the special case of (0 | 1) is not
  * handled correctly.
  */
-int mp_jacobi(mp_int *a, mp_int *n, int *c)
+int mp_jacobi(const mp_int *a, const mp_int *n, int *c)
 {
    mp_int  a1, p1;
    int     k, s, r, res;
diff --git a/libtommath/bn_mp_lcm.c b/libtommath/bn_mp_lcm.c
index dc661f3..24b621c 100644
--- a/libtommath/bn_mp_lcm.c
+++ b/libtommath/bn_mp_lcm.c
@@ -16,7 +16,7 @@
  */
 
 /* computes least common multiple as |a*b|/(a, b) */
-int mp_lcm(mp_int *a, mp_int *b, mp_int *c)
+int mp_lcm(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     res;
    mp_int  t1, t2;
diff --git a/libtommath/bn_mp_mod_d.c b/libtommath/bn_mp_mod_d.c
index 6afe4f1..5217aa4 100644
--- a/libtommath/bn_mp_mod_d.c
+++ b/libtommath/bn_mp_mod_d.c
@@ -16,7 +16,7 @@
  */
 
 int
-mp_mod_d(mp_int *a, mp_digit b, mp_digit *c)
+mp_mod_d(const mp_int *a, mp_digit b, mp_digit *c)
 {
    return mp_div_d(a, b, NULL, c);
 }
diff --git a/libtommath/bn_mp_montgomery_calc_normalization.c b/libtommath/bn_mp_montgomery_calc_normalization.c
index 2d95140..f2b0856 100644
--- a/libtommath/bn_mp_montgomery_calc_normalization.c
+++ b/libtommath/bn_mp_montgomery_calc_normalization.c
@@ -21,7 +21,7 @@
  * The method is slightly modified to shift B unconditionally upto just under
  * the leading bit of b.  This saves alot of multiple precision shifting.
  */
-int mp_montgomery_calc_normalization(mp_int *a, mp_int *b)
+int mp_montgomery_calc_normalization(mp_int *a, const mp_int *b)
 {
    int     x, bits, res;
 
diff --git a/libtommath/bn_mp_montgomery_reduce.c b/libtommath/bn_mp_montgomery_reduce.c
index 1ee69b2..e9f6c1c 100644
--- a/libtommath/bn_mp_montgomery_reduce.c
+++ b/libtommath/bn_mp_montgomery_reduce.c
@@ -17,7 +17,7 @@
 
 /* computes xR**-1 == x (mod N) via Montgomery Reduction */
 int
-mp_montgomery_reduce(mp_int *x, mp_int *n, mp_digit rho)
+mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho)
 {
    int     ix, res, digs;
    mp_digit mu;
diff --git a/libtommath/bn_mp_montgomery_setup.c b/libtommath/bn_mp_montgomery_setup.c
index fac71a6..37069c3 100644
--- a/libtommath/bn_mp_montgomery_setup.c
+++ b/libtommath/bn_mp_montgomery_setup.c
@@ -17,7 +17,7 @@
 
 /* setups the montgomery reduction stuff */
 int
-mp_montgomery_setup(mp_int *n, mp_digit *rho)
+mp_montgomery_setup(const mp_int *n, mp_digit *rho)
 {
    mp_digit x, b;
 
diff --git a/libtommath/bn_mp_mulmod.c b/libtommath/bn_mp_mulmod.c
index b1e6a33..aeee4ee 100644
--- a/libtommath/bn_mp_mulmod.c
+++ b/libtommath/bn_mp_mulmod.c
@@ -16,7 +16,7 @@
  */
 
 /* d = a * b (mod c) */
-int mp_mulmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
+int mp_mulmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d)
 {
    int     res;
    mp_int  t;
diff --git a/libtommath/bn_mp_n_root.c b/libtommath/bn_mp_n_root.c
index 8211c0a..a09804f 100644
--- a/libtommath/bn_mp_n_root.c
+++ b/libtommath/bn_mp_n_root.c
@@ -18,7 +18,7 @@
 /* wrapper function for mp_n_root_ex()
  * computes c = (a)**(1/b) such that (c)**b <= a and (c+1)**b > a
  */
-int mp_n_root(mp_int *a, mp_digit b, mp_int *c)
+int mp_n_root(const mp_int *a, mp_digit b, mp_int *c)
 {
    return mp_n_root_ex(a, b, c, 0);
 }
diff --git a/libtommath/bn_mp_n_root_ex.c b/libtommath/bn_mp_n_root_ex.c
index 9546745..ca50649 100644
--- a/libtommath/bn_mp_n_root_ex.c
+++ b/libtommath/bn_mp_n_root_ex.c
@@ -25,10 +25,10 @@
  * each step involves a fair bit.  This is not meant to
  * find huge roots [square and cube, etc].
  */
-int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
+int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
 {
-   mp_int  t1, t2, t3;
-   int     res, neg;
+   mp_int  t1, t2, t3, a_;
+   int     res;
 
    /* input must be positive if b is even */
    if (((b & 1) == 0) && (a->sign == MP_NEG)) {
@@ -48,8 +48,8 @@ int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
    }
 
    /* if a is negative fudge the sign but keep track */
-   neg     = a->sign;
-   a->sign = MP_ZPOS;
+   a_ = *a;
+   a_.sign = MP_ZPOS;
 
    /* t2 = 2 */
    mp_set(&t2, 2);
@@ -74,7 +74,7 @@ int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
       }
 
       /* t2 = t1**b - a */
-      if ((res = mp_sub(&t2, a, &t2)) != MP_OKAY) {
+      if ((res = mp_sub(&t2, &a_, &t2)) != MP_OKAY) {
          goto LBL_T3;
       }
 
@@ -100,7 +100,7 @@ int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
          goto LBL_T3;
       }
 
-      if (mp_cmp(&t2, a) == MP_GT) {
+      if (mp_cmp(&t2, &a_) == MP_GT) {
          if ((res = mp_sub_d(&t1, 1, &t1)) != MP_OKAY) {
             goto LBL_T3;
          }
@@ -109,14 +109,11 @@ int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
       }
    }
 
-   /* reset the sign of a first */
-   a->sign = neg;
-
    /* set the result */
    mp_exch(&t1, c);
 
    /* set the sign of the result */
-   c->sign = neg;
+   c->sign = a->sign;
 
    res = MP_OKAY;
 
diff --git a/libtommath/bn_mp_reduce.c b/libtommath/bn_mp_reduce.c
index a2b9bf7..6665acb 100644
--- a/libtommath/bn_mp_reduce.c
+++ b/libtommath/bn_mp_reduce.c
@@ -19,7 +19,7 @@
  * precomputed via mp_reduce_setup.
  * From HAC pp.604 Algorithm 14.42
  */
-int mp_reduce(mp_int *x, mp_int *m, mp_int *mu)
+int mp_reduce(mp_int *x, const mp_int *m, mp_int *mu)
 {
    mp_int  q;
    int     res, um = m->used;
diff --git a/libtommath/bn_mp_reduce_2k.c b/libtommath/bn_mp_reduce_2k.c
index 6bc96d1..2922cad 100644
--- a/libtommath/bn_mp_reduce_2k.c
+++ b/libtommath/bn_mp_reduce_2k.c
@@ -16,7 +16,7 @@
  */
 
 /* reduces a modulo n where n is of the form 2**p - d */
-int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d)
+int mp_reduce_2k(mp_int *a, const mp_int *n, mp_digit d)
 {
    mp_int q;
    int    p, res;
diff --git a/libtommath/bn_mp_reduce_2k_l.c b/libtommath/bn_mp_reduce_2k_l.c
index 8e6eeb0..3b23a37 100644
--- a/libtommath/bn_mp_reduce_2k_l.c
+++ b/libtommath/bn_mp_reduce_2k_l.c
@@ -19,7 +19,7 @@
    This differs from reduce_2k since "d" can be larger
    than a single digit.
 */
-int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
+int mp_reduce_2k_l(mp_int *a, const mp_int *n, mp_int *d)
 {
    mp_int q;
    int    p, res;
diff --git a/libtommath/bn_mp_reduce_2k_setup.c b/libtommath/bn_mp_reduce_2k_setup.c
index 802a5ba..e6ae839 100644
--- a/libtommath/bn_mp_reduce_2k_setup.c
+++ b/libtommath/bn_mp_reduce_2k_setup.c
@@ -16,7 +16,7 @@
  */
 
 /* determines the setup value */
-int mp_reduce_2k_setup(mp_int *a, mp_digit *d)
+int mp_reduce_2k_setup(const mp_int *a, mp_digit *d)
 {
    int res, p;
    mp_int tmp;
diff --git a/libtommath/bn_mp_reduce_2k_setup_l.c b/libtommath/bn_mp_reduce_2k_setup_l.c
index 34367ed..af81b5b 100644
--- a/libtommath/bn_mp_reduce_2k_setup_l.c
+++ b/libtommath/bn_mp_reduce_2k_setup_l.c
@@ -16,7 +16,7 @@
  */
 
 /* determines the setup value */
-int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
+int mp_reduce_2k_setup_l(const mp_int *a, mp_int *d)
 {
    int    res;
    mp_int tmp;
diff --git a/libtommath/bn_mp_reduce_is_2k.c b/libtommath/bn_mp_reduce_is_2k.c
index c733ca9..932521e 100644
--- a/libtommath/bn_mp_reduce_is_2k.c
+++ b/libtommath/bn_mp_reduce_is_2k.c
@@ -16,7 +16,7 @@
  */
 
 /* determines if mp_reduce_2k can be used */
-int mp_reduce_is_2k(mp_int *a)
+int mp_reduce_is_2k(const mp_int *a)
 {
    int ix, iy, iw;
    mp_digit iz;
diff --git a/libtommath/bn_mp_reduce_is_2k_l.c b/libtommath/bn_mp_reduce_is_2k_l.c
index d4804d5..22c7582 100644
--- a/libtommath/bn_mp_reduce_is_2k_l.c
+++ b/libtommath/bn_mp_reduce_is_2k_l.c
@@ -16,7 +16,7 @@
  */
 
 /* determines if reduce_2k_l can be used */
-int mp_reduce_is_2k_l(mp_int *a)
+int mp_reduce_is_2k_l(const mp_int *a)
 {
    int ix, iy;
 
diff --git a/libtommath/bn_mp_reduce_setup.c b/libtommath/bn_mp_reduce_setup.c
index 00ff61c..70e193a 100644
--- a/libtommath/bn_mp_reduce_setup.c
+++ b/libtommath/bn_mp_reduce_setup.c
@@ -18,7 +18,7 @@
 /* pre-calculate the value required for Barrett reduction
  * For a given modulus "b" it calulates the value required in "a"
  */
-int mp_reduce_setup(mp_int *a, mp_int *b)
+int mp_reduce_setup(mp_int *a, const mp_int *b)
 {
    int     res;
 
diff --git a/libtommath/bn_mp_sqrmod.c b/libtommath/bn_mp_sqrmod.c
index ebb1b53..b8265fe 100644
--- a/libtommath/bn_mp_sqrmod.c
+++ b/libtommath/bn_mp_sqrmod.c
@@ -17,7 +17,7 @@
 
 /* c = a * a (mod b) */
 int
-mp_sqrmod(mp_int *a, mp_int *b, mp_int *c)
+mp_sqrmod(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     res;
    mp_int  t;
diff --git a/libtommath/bn_mp_sqrtmod_prime.c b/libtommath/bn_mp_sqrtmod_prime.c
index 12b427c..261723e 100644
--- a/libtommath/bn_mp_sqrtmod_prime.c
+++ b/libtommath/bn_mp_sqrtmod_prime.c
@@ -15,7 +15,7 @@
  *
  */
 
-int mp_sqrtmod_prime(mp_int *n, mp_int *prime, mp_int *ret)
+int mp_sqrtmod_prime(const mp_int *n, const mp_int *prime, mp_int *ret)
 {
    int res, legendre;
    mp_int t1, C, Q, S, Z, M, T, R, two;
diff --git a/libtommath/bn_mp_submod.c b/libtommath/bn_mp_submod.c
index 96a3a1b..9a45f6d 100644
--- a/libtommath/bn_mp_submod.c
+++ b/libtommath/bn_mp_submod.c
@@ -17,7 +17,7 @@
 
 /* d = a - b (mod c) */
 int
-mp_submod(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
+mp_submod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d)
 {
    int     res;
    mp_int  t;
diff --git a/libtommath/bn_s_mp_exptmod.c b/libtommath/bn_s_mp_exptmod.c
index bd37169..e73c12e 100644
--- a/libtommath/bn_s_mp_exptmod.c
+++ b/libtommath/bn_s_mp_exptmod.c
@@ -20,12 +20,12 @@
 #define TAB_SIZE 256
 #endif
 
-int s_mp_exptmod(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode)
+int s_mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode)
 {
    mp_int  M[TAB_SIZE], res, mu;
    mp_digit buf;
    int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
-   int (*redux)(mp_int *,mp_int *,mp_int *);
+   int (*redux)(mp_int *,const mp_int *,mp_int *);
 
    /* find window size */
    x = mp_count_bits(X);
diff --git a/libtommath/tommath.h b/libtommath/tommath.h
index 9d811fa..513b5b9 100644
--- a/libtommath/tommath.h
+++ b/libtommath/tommath.h
@@ -358,98 +358,98 @@ int mp_mod_d(const mp_int *a, mp_digit b, mp_digit *c);
 /* ---> number theory <--- */
 
 /* d = a + b (mod c) */
-int mp_addmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_addmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* d = a - b (mod c) */
-int mp_submod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_submod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* d = a * b (mod c) */
-int mp_mulmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_mulmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* c = a * a (mod b) */
-int mp_sqrmod(mp_int *a, mp_int *b, mp_int *c);
+int mp_sqrmod(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = 1/a (mod b) */
-int mp_invmod(mp_int *a, mp_int *b, mp_int *c);
+int mp_invmod(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = (a, b) */
-int mp_gcd(mp_int *a, mp_int *b, mp_int *c);
+int mp_gcd(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* produces value such that U1*a + U2*b = U3 */
-int mp_exteuclid(mp_int *a, mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3);
+int mp_exteuclid(const mp_int *a, const mp_int *b, mp_int *U1, mp_int *U2, mp_int *U3);
 
 /* c = [a, b] or (a*b)/(a, b) */
-int mp_lcm(mp_int *a, mp_int *b, mp_int *c);
+int mp_lcm(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* finds one of the b'th root of a, such that |c|**b <= |a|
  *
  * returns error if a < 0 and b is even
  */
-int mp_n_root(mp_int *a, mp_digit b, mp_int *c);
-int mp_n_root_ex(mp_int *a, mp_digit b, mp_int *c, int fast);
+int mp_n_root(const mp_int *a, mp_digit b, mp_int *c);
+int mp_n_root_ex(const mp_int *a, mp_digit b, mp_int *c, int fast);
 
 /* special sqrt algo */
 int mp_sqrt(const mp_int *arg, mp_int *ret);
 
 /* special sqrt (mod prime) */
-int mp_sqrtmod_prime(mp_int *arg, mp_int *prime, mp_int *ret);
+int mp_sqrtmod_prime(const mp_int *arg, const mp_int *prime, mp_int *ret);
 
 /* is number a square? */
-int mp_is_square(mp_int *arg, int *ret);
+int mp_is_square(const mp_int *arg, int *ret);
 
 /* computes the jacobi c = (a | n) (or Legendre if b is prime)  */
-int mp_jacobi(mp_int *a, mp_int *n, int *c);
+int mp_jacobi(const mp_int *a, const mp_int *n, int *c);
 
 /* used to setup the Barrett reduction for a given modulus b */
-int mp_reduce_setup(mp_int *a, mp_int *b);
+int mp_reduce_setup(mp_int *a, const mp_int *b);
 
 /* Barrett Reduction, computes a (mod b) with a precomputed value c
  *
  * Assumes that 0 < a <= b*b, note if 0 > a > -(b*b) then you can merely
  * compute the reduction as -1 * mp_reduce(mp_abs(a)) [pseudo code].
  */
-int mp_reduce(mp_int *a, mp_int *b, mp_int *c);
+int mp_reduce(mp_int *a, const mp_int *b, mp_int *c);
 
 /* setups the montgomery reduction */
-int mp_montgomery_setup(mp_int *a, mp_digit *mp);
+int mp_montgomery_setup(const mp_int *a, mp_digit *mp);
 
 /* computes a = B**n mod b without division or multiplication useful for
  * normalizing numbers in a Montgomery system.
  */
-int mp_montgomery_calc_normalization(mp_int *a, mp_int *b);
+int mp_montgomery_calc_normalization(mp_int *a, const mp_int *b);
 
 /* computes x/R == x (mod N) via Montgomery Reduction */
-int mp_montgomery_reduce(mp_int *a, mp_int *m, mp_digit mp);
+int mp_montgomery_reduce(mp_int *a, const mp_int *m, mp_digit mp);
 
 /* returns 1 if a is a valid DR modulus */
-int mp_dr_is_modulus(mp_int *a);
+int mp_dr_is_modulus(const mp_int *a);
 
 /* sets the value of "d" required for mp_dr_reduce */
-void mp_dr_setup(mp_int *a, mp_digit *d);
+void mp_dr_setup(const mp_int *a, mp_digit *d);
 
 /* reduces a modulo b using the Diminished Radix method */
-int mp_dr_reduce(mp_int *a, mp_int *b, mp_digit mp);
+int mp_dr_reduce(mp_int *a, const mp_int *b, mp_digit mp);
 
 /* returns true if a can be reduced with mp_reduce_2k */
-int mp_reduce_is_2k(mp_int *a);
+int mp_reduce_is_2k(const mp_int *a);
 
 /* determines k value for 2k reduction */
-int mp_reduce_2k_setup(mp_int *a, mp_digit *d);
+int mp_reduce_2k_setup(const mp_int *a, mp_digit *d);
 
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
-int mp_reduce_2k(mp_int *a, mp_int *n, mp_digit d);
+int mp_reduce_2k(mp_int *a, const mp_int *n, mp_digit d);
 
 /* returns true if a can be reduced with mp_reduce_2k_l */
-int mp_reduce_is_2k_l(mp_int *a);
+int mp_reduce_is_2k_l(const mp_int *a);
 
 /* determines k value for 2k reduction */
-int mp_reduce_2k_setup_l(mp_int *a, mp_int *d);
+int mp_reduce_2k_setup_l(const mp_int *a, mp_int *d);
 
 /* reduces a modulo b where b is of the form 2**p - k [0 <= a] */
-int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d);
+int mp_reduce_2k_l(mp_int *a, const mp_int *n, mp_int *d);
 
 /* d = a**b (mod c) */
-int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_exptmod(const mp_int *a, const mp_int *b, const mp_int *c, mp_int *d);
 
 /* ---> Primes <--- */
 
diff --git a/libtommath/tommath_private.h b/libtommath/tommath_private.h
index 7f33fab..58846bf 100644
--- a/libtommath/tommath_private.h
+++ b/libtommath/tommath_private.h
@@ -70,11 +70,11 @@ int mp_karatsuba_mul(const mp_int *a, const mp_int *b, mp_int *c);
 int mp_toom_mul(const mp_int *a, const mp_int *b, mp_int *c);
 int mp_karatsuba_sqr(const mp_int *a, mp_int *b);
 int mp_toom_sqr(const mp_int *a, mp_int *b);
-int fast_mp_invmod(mp_int *a, mp_int *b, mp_int *c);
-int mp_invmod_slow(mp_int *a, mp_int *b, mp_int *c);
-int fast_mp_montgomery_reduce(mp_int *x, mp_int *n, mp_digit rho);
-int mp_exptmod_fast(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode);
-int s_mp_exptmod(mp_int *G, mp_int *X, mp_int *P, mp_int *Y, int redmode);
+int fast_mp_invmod(const mp_int *a, const mp_int *b, mp_int *c);
+int mp_invmod_slow(const mp_int *a, const mp_int *b, mp_int *c);
+int fast_mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho);
+int mp_exptmod_fast(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode);
+int s_mp_exptmod(const mp_int *G, const mp_int *X, const mp_int *P, mp_int *Y, int redmode);
 void bn_reverse(unsigned char *s, int len);
 
 extern const char *mp_s_rmap;