kc3-lang/libtommath

diff --git a/bn_fast_mp_invmod.c b/bn_fast_mp_invmod.c
index 7771136..08389dd 100644
--- a/bn_fast_mp_invmod.c
+++ b/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/bn_fast_mp_montgomery_reduce.c b/bn_fast_mp_montgomery_reduce.c
index f2c38bf..54d9b0a 100644
--- a/bn_fast_mp_montgomery_reduce.c
+++ b/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/bn_fast_s_mp_mul_digs.c b/bn_fast_s_mp_mul_digs.c
index 763dbb1..558d151 100644
--- a/bn_fast_s_mp_mul_digs.c
+++ b/bn_fast_s_mp_mul_digs.c
@@ -31,7 +31,7 @@
  * Based on Algorithm 14.12 on pp.595 of HAC.
  *
  */
-int fast_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
+int fast_s_mp_mul_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
 {
    int     olduse, res, pa, ix, iz;
    mp_digit W[MP_WARRAY];
diff --git a/bn_fast_s_mp_mul_high_digs.c b/bn_fast_s_mp_mul_high_digs.c
index 588d80b..8b662ed 100644
--- a/bn_fast_s_mp_mul_high_digs.c
+++ b/bn_fast_s_mp_mul_high_digs.c
@@ -24,7 +24,7 @@
  *
  * Based on Algorithm 14.12 on pp.595 of HAC.
  */
-int fast_s_mp_mul_high_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
+int fast_s_mp_mul_high_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
 {
    int     olduse, res, pa, ix, iz;
    mp_digit W[MP_WARRAY];
diff --git a/bn_fast_s_mp_sqr.c b/bn_fast_s_mp_sqr.c
index ceed82b..161f785 100644
--- a/bn_fast_s_mp_sqr.c
+++ b/bn_fast_s_mp_sqr.c
@@ -25,7 +25,7 @@
 After that loop you do the squares and add them in.
 */
 
-int fast_s_mp_sqr(mp_int *a, mp_int *b)
+int fast_s_mp_sqr(const mp_int *a, mp_int *b)
 {
    int       olduse, res, pa, ix, iz;
    mp_digit   W[MP_WARRAY], *tmpx;
diff --git a/bn_mp_abs.c b/bn_mp_abs.c
index d5fc012..9b6bcec 100644
--- a/bn_mp_abs.c
+++ b/bn_mp_abs.c
@@ -19,7 +19,7 @@
  *
  * Simple function copies the input and fixes the sign to positive
  */
-int mp_abs(mp_int *a, mp_int *b)
+int mp_abs(const mp_int *a, mp_int *b)
 {
    int     res;
 
diff --git a/bn_mp_add.c b/bn_mp_add.c
index 4df4c81..d31d5a0 100644
--- a/bn_mp_add.c
+++ b/bn_mp_add.c
@@ -16,7 +16,7 @@
  */
 
 /* high level addition (handles signs) */
-int mp_add(mp_int *a, mp_int *b, mp_int *c)
+int mp_add(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     sa, sb, res;
 
diff --git a/bn_mp_add_d.c b/bn_mp_add_d.c
index 1e6ff63..e5ede1f 100644
--- a/bn_mp_add_d.c
+++ b/bn_mp_add_d.c
@@ -16,7 +16,7 @@
  */
 
 /* single digit addition */
-int mp_add_d(mp_int *a, mp_digit b, mp_int *c)
+int mp_add_d(const mp_int *a, mp_digit b, mp_int *c)
 {
    int     res, ix, oldused;
    mp_digit *tmpa, *tmpc, mu;
@@ -30,14 +30,15 @@ int mp_add_d(mp_int *a, mp_digit b, mp_int *c)
 
    /* if a is negative and |a| >= b, call c = |a| - b */
    if ((a->sign == MP_NEG) && ((a->used > 1) || (a->dp[0] >= b))) {
+      mp_int a_ = *a;
       /* temporarily fix sign of a */
-      a->sign = MP_ZPOS;
+      a_.sign = MP_ZPOS;
 
       /* c = |a| - b */
-      res = mp_sub_d(a, b, c);
+      res = mp_sub_d(&a_, b, c);
 
       /* fix sign  */
-      a->sign = c->sign = MP_NEG;
+      c->sign = MP_NEG;
 
       /* clamp */
       mp_clamp(c);
diff --git a/bn_mp_addmod.c b/bn_mp_addmod.c
index 229a716..0d612c3 100644
--- a/bn_mp_addmod.c
+++ b/bn_mp_addmod.c
@@ -16,7 +16,7 @@
  */
 
 /* 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)
 {
    int     res;
    mp_int  t;
diff --git a/bn_mp_and.c b/bn_mp_and.c
index 2f1472a..09ff772 100644
--- a/bn_mp_and.c
+++ b/bn_mp_and.c
@@ -16,10 +16,11 @@
  */
 
 /* AND two ints together */
-int mp_and(mp_int *a, mp_int *b, mp_int *c)
+int mp_and(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     res, ix, px;
-   mp_int  t, *x;
+   mp_int  t;
+   const mp_int *x;
 
    if (a->used > b->used) {
       if ((res = mp_init_copy(&t, a)) != MP_OKAY) {
diff --git a/bn_mp_div.c b/bn_mp_div.c
index fdb3453..dbfdc03 100644
--- a/bn_mp_div.c
+++ b/bn_mp_div.c
@@ -18,7 +18,7 @@
 #ifdef BN_MP_DIV_SMALL
 
 /* slower bit-bang division... also smaller */
-int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
+int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
 {
    mp_int ta, tb, tq, q;
    int    res, n, n2;
@@ -100,7 +100,7 @@ LBL_ERR:
  * The overall algorithm is as described as
  * 14.20 from HAC but fixed to treat these cases.
 */
-int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d)
+int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
 {
    mp_int  q, x, y, t1, t2;
    int     res, n, t, i, norm, neg;
diff --git a/bn_mp_div_2.c b/bn_mp_div_2.c
index b9d5339..edc8982 100644
--- a/bn_mp_div_2.c
+++ b/bn_mp_div_2.c
@@ -16,7 +16,7 @@
  */
 
 /* b = a/2 */
-int mp_div_2(mp_int *a, mp_int *b)
+int mp_div_2(const mp_int *a, mp_int *b)
 {
    int     x, res, oldused;
 
diff --git a/bn_mp_div_3.c b/bn_mp_div_3.c
index c3a023a..9cc8caa 100644
--- a/bn_mp_div_3.c
+++ b/bn_mp_div_3.c
@@ -16,7 +16,7 @@
  */
 
 /* divide by three (based on routine from MPI and the GMP manual) */
-int mp_div_3(mp_int *a, mp_int *c, mp_digit *d)
+int mp_div_3(const mp_int *a, mp_int *c, mp_digit *d)
 {
    mp_int   q;
    mp_word  w, t;
diff --git a/bn_mp_div_d.c b/bn_mp_div_d.c
index 141db1d..db4a0a2 100644
--- a/bn_mp_div_d.c
+++ b/bn_mp_div_d.c
@@ -34,7 +34,7 @@ static int s_is_power_of_two(mp_digit b, int *p)
 }
 
 /* single digit division (based on routine from MPI) */
-int mp_div_d(mp_int *a, mp_digit b, mp_int *c, mp_digit *d)
+int mp_div_d(const mp_int *a, mp_digit b, mp_int *c, mp_digit *d)
 {
    mp_int  q;
    mp_word w;
diff --git a/bn_mp_dr_is_modulus.c b/bn_mp_dr_is_modulus.c
index 4631daa..bf4ed8b 100644
--- a/bn_mp_dr_is_modulus.c
+++ b/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/bn_mp_dr_reduce.c b/bn_mp_dr_reduce.c
index 25079be..1ccb669 100644
--- a/bn_mp_dr_reduce.c
+++ b/bn_mp_dr_reduce.c
@@ -29,7 +29,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)
+int mp_dr_reduce(mp_int *x, const mp_int *n, mp_digit k)
 {
    int      err, i, m;
    mp_word  r;
diff --git a/bn_mp_dr_setup.c b/bn_mp_dr_setup.c
index 97f31ba..af0e213 100644
--- a/bn_mp_dr_setup.c
+++ b/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/bn_mp_expt_d.c b/bn_mp_expt_d.c
index 38bf679..f5ce3c1 100644
--- a/bn_mp_expt_d.c
+++ b/bn_mp_expt_d.c
@@ -16,7 +16,7 @@
  */
 
 /* wrapper function for mp_expt_d_ex() */
-int mp_expt_d(mp_int *a, mp_digit b, mp_int *c)
+int mp_expt_d(const mp_int *a, mp_digit b, mp_int *c)
 {
    return mp_expt_d_ex(a, b, c, 0);
 }
diff --git a/bn_mp_expt_d_ex.c b/bn_mp_expt_d_ex.c
index bece77b..99319a5 100644
--- a/bn_mp_expt_d_ex.c
+++ b/bn_mp_expt_d_ex.c
@@ -16,7 +16,7 @@
  */
 
 /* calculate c = a**b  using a square-multiply algorithm */
-int mp_expt_d_ex(mp_int *a, mp_digit b, mp_int *c, int fast)
+int mp_expt_d_ex(const mp_int *a, mp_digit b, mp_int *c, int fast)
 {
    int     res;
    unsigned int x;
diff --git a/bn_mp_exptmod.c b/bn_mp_exptmod.c
index c4f392b..934fd25 100644
--- a/bn_mp_exptmod.c
+++ b/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/bn_mp_exptmod_fast.c b/bn_mp_exptmod_fast.c
index 38e0265..4a188d0 100644
--- a/bn_mp_exptmod_fast.c
+++ b/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/bn_mp_exteuclid.c b/bn_mp_exteuclid.c
index 98eef76..08e5ff2 100644
--- a/bn_mp_exteuclid.c
+++ b/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/bn_mp_fwrite.c b/bn_mp_fwrite.c
index 3641823..829dd4a 100644
--- a/bn_mp_fwrite.c
+++ b/bn_mp_fwrite.c
@@ -16,7 +16,7 @@
  */
 
 #ifndef LTM_NO_FILE
-int mp_fwrite(mp_int *a, int radix, FILE *stream)
+int mp_fwrite(const mp_int *a, int radix, FILE *stream)
 {
    char *buf;
    int err, len, x;
diff --git a/bn_mp_gcd.c b/bn_mp_gcd.c
index 18f6dc3..f5aa78b 100644
--- a/bn_mp_gcd.c
+++ b/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/bn_mp_get_int.c b/bn_mp_get_int.c
index a3d1602..f4a347f 100644
--- a/bn_mp_get_int.c
+++ b/bn_mp_get_int.c
@@ -16,7 +16,7 @@
  */
 
 /* get the lower 32-bits of an mp_int */
-unsigned long mp_get_int(mp_int *a)
+unsigned long mp_get_int(const mp_int *a)
 {
    int i;
    mp_min_u32 res;
diff --git a/bn_mp_get_long.c b/bn_mp_get_long.c
index 053930c..3fc7c35 100644
--- a/bn_mp_get_long.c
+++ b/bn_mp_get_long.c
@@ -16,7 +16,7 @@
  */
 
 /* get the lower unsigned long of an mp_int, platform dependent */
-unsigned long mp_get_long(mp_int *a)
+unsigned long mp_get_long(const mp_int *a)
 {
    int i;
    unsigned long res;
diff --git a/bn_mp_get_long_long.c b/bn_mp_get_long_long.c
index 131571a..838c3c3 100644
--- a/bn_mp_get_long_long.c
+++ b/bn_mp_get_long_long.c
@@ -16,7 +16,7 @@
  */
 
 /* get the lower unsigned long long of an mp_int, platform dependent */
-unsigned long long mp_get_long_long(mp_int *a)
+unsigned long long mp_get_long_long(const mp_int *a)
 {
    int i;
    unsigned long long res;
diff --git a/bn_mp_invmod.c b/bn_mp_invmod.c
index b70fe18..525493a 100644
--- a/bn_mp_invmod.c
+++ b/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/bn_mp_invmod_slow.c b/bn_mp_invmod_slow.c
index 2bdd2b1..2bb5196 100644
--- a/bn_mp_invmod_slow.c
+++ b/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/bn_mp_is_square.c b/bn_mp_is_square.c
index 303fab6..dd5150e 100644
--- a/bn_mp_is_square.c
+++ b/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/bn_mp_jacobi.c b/bn_mp_jacobi.c
index 8981393..c314c82 100644
--- a/bn_mp_jacobi.c
+++ b/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/bn_mp_karatsuba_mul.c b/bn_mp_karatsuba_mul.c
index 353c37c..1a84211 100644
--- a/bn_mp_karatsuba_mul.c
+++ b/bn_mp_karatsuba_mul.c
@@ -44,7 +44,7 @@
  * Generally though the overhead of this method doesn't pay off
  * until a certain size (N ~ 80) is reached.
  */
-int mp_karatsuba_mul(mp_int *a, mp_int *b, mp_int *c)
+int mp_karatsuba_mul(const mp_int *a, const mp_int *b, mp_int *c)
 {
    mp_int  x0, x1, y0, y1, t1, x0y0, x1y1;
    int     B, err;
diff --git a/bn_mp_karatsuba_sqr.c b/bn_mp_karatsuba_sqr.c
index fe39a33..c566b06 100644
--- a/bn_mp_karatsuba_sqr.c
+++ b/bn_mp_karatsuba_sqr.c
@@ -22,7 +22,7 @@
  * is essentially the same algorithm but merely
  * tuned to perform recursive squarings.
  */
-int mp_karatsuba_sqr(mp_int *a, mp_int *b)
+int mp_karatsuba_sqr(const mp_int *a, mp_int *b)
 {
    mp_int  x0, x1, t1, t2, x0x0, x1x1;
    int     B, err;
diff --git a/bn_mp_lcm.c b/bn_mp_lcm.c
index dc661f3..24b621c 100644
--- a/bn_mp_lcm.c
+++ b/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/bn_mp_mod.c b/bn_mp_mod.c
index 267688c..64e73ea 100644
--- a/bn_mp_mod.c
+++ b/bn_mp_mod.c
@@ -16,7 +16,7 @@
  */
 
 /* c = a mod b, 0 <= c < b if b > 0, b < c <= 0 if b < 0 */
-int mp_mod(mp_int *a, mp_int *b, mp_int *c)
+int mp_mod(const mp_int *a, const mp_int *b, mp_int *c)
 {
    mp_int  t;
    int     res;
diff --git a/bn_mp_mod_d.c b/bn_mp_mod_d.c
index ff77346..9a24e78 100644
--- a/bn_mp_mod_d.c
+++ b/bn_mp_mod_d.c
@@ -15,7 +15,7 @@
  * Tom St Denis, tstdenis82@gmail.com, http://libtom.org
  */
 
-int mp_mod_d(mp_int *a, mp_digit b, mp_digit *c)
+int mp_mod_d(const mp_int *a, mp_digit b, mp_digit *c)
 {
    return mp_div_d(a, b, NULL, c);
 }
diff --git a/bn_mp_montgomery_calc_normalization.c b/bn_mp_montgomery_calc_normalization.c
index 2d95140..f2b0856 100644
--- a/bn_mp_montgomery_calc_normalization.c
+++ b/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/bn_mp_montgomery_reduce.c b/bn_mp_montgomery_reduce.c
index 1909997..a38173e 100644
--- a/bn_mp_montgomery_reduce.c
+++ b/bn_mp_montgomery_reduce.c
@@ -16,7 +16,7 @@
  */
 
 /* computes xR**-1 == x (mod N) via Montgomery Reduction */
-int mp_montgomery_reduce(mp_int *x, mp_int *n, mp_digit rho)
+int mp_montgomery_reduce(mp_int *x, const mp_int *n, mp_digit rho)
 {
    int     ix, res, digs;
    mp_digit mu;
diff --git a/bn_mp_montgomery_setup.c b/bn_mp_montgomery_setup.c
index 46f560d..685ba51 100644
--- a/bn_mp_montgomery_setup.c
+++ b/bn_mp_montgomery_setup.c
@@ -16,7 +16,7 @@
  */
 
 /* setups the montgomery reduction stuff */
-int mp_montgomery_setup(mp_int *n, mp_digit *rho)
+int mp_montgomery_setup(const mp_int *n, mp_digit *rho)
 {
    mp_digit x, b;
 
diff --git a/bn_mp_mul.c b/bn_mp_mul.c
index 315a520..71d523d 100644
--- a/bn_mp_mul.c
+++ b/bn_mp_mul.c
@@ -16,7 +16,7 @@
  */
 
 /* high level multiplication (handles sign) */
-int mp_mul(mp_int *a, mp_int *b, mp_int *c)
+int mp_mul(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     res, neg;
    neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
diff --git a/bn_mp_mul_2.c b/bn_mp_mul_2.c
index 33bdc4c..1744681 100644
--- a/bn_mp_mul_2.c
+++ b/bn_mp_mul_2.c
@@ -16,7 +16,7 @@
  */
 
 /* b = a*2 */
-int mp_mul_2(mp_int *a, mp_int *b)
+int mp_mul_2(const mp_int *a, mp_int *b)
 {
    int     x, res, oldused;
 
diff --git a/bn_mp_mul_d.c b/bn_mp_mul_d.c
index 1aa448c..0f6d03e 100644
--- a/bn_mp_mul_d.c
+++ b/bn_mp_mul_d.c
@@ -16,7 +16,7 @@
  */
 
 /* multiply by a digit */
-int mp_mul_d(mp_int *a, mp_digit b, mp_int *c)
+int mp_mul_d(const mp_int *a, mp_digit b, mp_int *c)
 {
    mp_digit u, *tmpa, *tmpc;
    mp_word  r;
diff --git a/bn_mp_mulmod.c b/bn_mp_mulmod.c
index b1e6a33..aeee4ee 100644
--- a/bn_mp_mulmod.c
+++ b/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/bn_mp_n_root.c b/bn_mp_n_root.c
index 8211c0a..a09804f 100644
--- a/bn_mp_n_root.c
+++ b/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/bn_mp_n_root_ex.c b/bn_mp_n_root_ex.c
index 9546745..ca50649 100644
--- a/bn_mp_n_root_ex.c
+++ b/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/bn_mp_or.c b/bn_mp_or.c
index 2318c79..f411509 100644
--- a/bn_mp_or.c
+++ b/bn_mp_or.c
@@ -16,10 +16,11 @@
  */
 
 /* OR two ints together */
-int mp_or(mp_int *a, mp_int *b, mp_int *c)
+int mp_or(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     res, ix, px;
-   mp_int  t, *x;
+   mp_int  t;
+   const mp_int *x;
 
    if (a->used > b->used) {
       if ((res = mp_init_copy(&t, a)) != MP_OKAY) {
diff --git a/bn_mp_prime_fermat.c b/bn_mp_prime_fermat.c
index 37ad6ec..9c15435 100644
--- a/bn_mp_prime_fermat.c
+++ b/bn_mp_prime_fermat.c
@@ -23,7 +23,7 @@
  *
  * Sets result to 1 if the congruence holds, or zero otherwise.
  */
-int mp_prime_fermat(mp_int *a, mp_int *b, int *result)
+int mp_prime_fermat(const mp_int *a, const mp_int *b, int *result)
 {
    mp_int  t;
    int     err;
diff --git a/bn_mp_prime_is_divisible.c b/bn_mp_prime_is_divisible.c
index 92af330..c1e1158 100644
--- a/bn_mp_prime_is_divisible.c
+++ b/bn_mp_prime_is_divisible.c
@@ -20,7 +20,7 @@
  *
  * sets result to 0 if not, 1 if yes
  */
-int mp_prime_is_divisible(mp_int *a, int *result)
+int mp_prime_is_divisible(const mp_int *a, int *result)
 {
    int     err, ix;
    mp_digit res;
diff --git a/bn_mp_prime_is_prime.c b/bn_mp_prime_is_prime.c
index 20a7d1f..e97712d 100644
--- a/bn_mp_prime_is_prime.c
+++ b/bn_mp_prime_is_prime.c
@@ -22,7 +22,7 @@
  *
  * Sets result to 1 if probably prime, 0 otherwise
  */
-int mp_prime_is_prime(mp_int *a, int t, int *result)
+int mp_prime_is_prime(const mp_int *a, int t, int *result)
 {
    mp_int  b;
    int     ix, err, res;
diff --git a/bn_mp_prime_miller_rabin.c b/bn_mp_prime_miller_rabin.c
index 917dc01..5de5f05 100644
--- a/bn_mp_prime_miller_rabin.c
+++ b/bn_mp_prime_miller_rabin.c
@@ -22,7 +22,7 @@
  * Randomly the chance of error is no more than 1/4 and often
  * very much lower.
  */
-int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result)
+int mp_prime_miller_rabin(const mp_int *a, const mp_int *b, int *result)
 {
    mp_int  n1, y, r;
    int     s, j, err;
diff --git a/bn_mp_reduce.c b/bn_mp_reduce.c
index a2b9bf7..6665acb 100644
--- a/bn_mp_reduce.c
+++ b/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/bn_mp_reduce_2k.c b/bn_mp_reduce_2k.c
index 6bc96d1..2922cad 100644
--- a/bn_mp_reduce_2k.c
+++ b/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/bn_mp_reduce_2k_l.c b/bn_mp_reduce_2k_l.c
index 8e6eeb0..3b23a37 100644
--- a/bn_mp_reduce_2k_l.c
+++ b/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/bn_mp_reduce_2k_setup.c b/bn_mp_reduce_2k_setup.c
index 802a5ba..e6ae839 100644
--- a/bn_mp_reduce_2k_setup.c
+++ b/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/bn_mp_reduce_2k_setup_l.c b/bn_mp_reduce_2k_setup_l.c
index 34367ed..af81b5b 100644
--- a/bn_mp_reduce_2k_setup_l.c
+++ b/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/bn_mp_reduce_is_2k.c b/bn_mp_reduce_is_2k.c
index c733ca9..932521e 100644
--- a/bn_mp_reduce_is_2k.c
+++ b/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/bn_mp_reduce_is_2k_l.c b/bn_mp_reduce_is_2k_l.c
index d4804d5..22c7582 100644
--- a/bn_mp_reduce_is_2k_l.c
+++ b/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/bn_mp_reduce_setup.c b/bn_mp_reduce_setup.c
index 00ff61c..70e193a 100644
--- a/bn_mp_reduce_setup.c
+++ b/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/bn_mp_signed_bin_size.c b/bn_mp_signed_bin_size.c
index 082aeca..1fdfd85 100644
--- a/bn_mp_signed_bin_size.c
+++ b/bn_mp_signed_bin_size.c
@@ -16,7 +16,7 @@
  */
 
 /* get the size for an signed equivalent */
-int mp_signed_bin_size(mp_int *a)
+int mp_signed_bin_size(const mp_int *a)
 {
    return 1 + mp_unsigned_bin_size(a);
 }
diff --git a/bn_mp_sqr.c b/bn_mp_sqr.c
index e2e8641..2b71097 100644
--- a/bn_mp_sqr.c
+++ b/bn_mp_sqr.c
@@ -16,7 +16,7 @@
  */
 
 /* computes b = a*a */
-int mp_sqr(mp_int *a, mp_int *b)
+int mp_sqr(const mp_int *a, mp_int *b)
 {
    int     res;
 
diff --git a/bn_mp_sqrmod.c b/bn_mp_sqrmod.c
index 96a7574..c3c7ec9 100644
--- a/bn_mp_sqrmod.c
+++ b/bn_mp_sqrmod.c
@@ -16,7 +16,7 @@
  */
 
 /* 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)
 {
    int     res;
    mp_int  t;
diff --git a/bn_mp_sqrt.c b/bn_mp_sqrt.c
index 95a6892..d70c523 100644
--- a/bn_mp_sqrt.c
+++ b/bn_mp_sqrt.c
@@ -16,7 +16,7 @@
  */
 
 /* this function is less generic than mp_n_root, simpler and faster */
-int mp_sqrt(mp_int *arg, mp_int *ret)
+int mp_sqrt(const mp_int *arg, mp_int *ret)
 {
    int res;
    mp_int t1, t2;
diff --git a/bn_mp_sqrtmod_prime.c b/bn_mp_sqrtmod_prime.c
index 12b427c..261723e 100644
--- a/bn_mp_sqrtmod_prime.c
+++ b/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/bn_mp_sub.c b/bn_mp_sub.c
index 75c7c2d..19cb65e 100644
--- a/bn_mp_sub.c
+++ b/bn_mp_sub.c
@@ -16,7 +16,7 @@
  */
 
 /* high level subtraction (handles signs) */
-int mp_sub(mp_int *a, mp_int *b, mp_int *c)
+int mp_sub(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     sa, sb, res;
 
diff --git a/bn_mp_sub_d.c b/bn_mp_sub_d.c
index 7016abc..4d66a90 100644
--- a/bn_mp_sub_d.c
+++ b/bn_mp_sub_d.c
@@ -16,7 +16,7 @@
  */
 
 /* single digit subtraction */
-int mp_sub_d(mp_int *a, mp_digit b, mp_int *c)
+int mp_sub_d(const mp_int *a, mp_digit b, mp_int *c)
 {
    mp_digit *tmpa, *tmpc, mu;
    int       res, ix, oldused;
@@ -32,9 +32,10 @@ int mp_sub_d(mp_int *a, mp_digit b, mp_int *c)
     * addition [with fudged signs]
     */
    if (a->sign == MP_NEG) {
-      a->sign = MP_ZPOS;
-      res     = mp_add_d(a, b, c);
-      a->sign = c->sign = MP_NEG;
+      mp_int a_ = *a;
+      a_.sign = MP_ZPOS;
+      res     = mp_add_d(&a_, b, c);
+      c->sign = MP_NEG;
 
       /* clamp */
       mp_clamp(c);
diff --git a/bn_mp_submod.c b/bn_mp_submod.c
index 510fb19..c4db397 100644
--- a/bn_mp_submod.c
+++ b/bn_mp_submod.c
@@ -16,7 +16,7 @@
  */
 
 /* 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)
 {
    int     res;
    mp_int  t;
diff --git a/bn_mp_to_signed_bin.c b/bn_mp_to_signed_bin.c
index 615bf32..4d4be88 100644
--- a/bn_mp_to_signed_bin.c
+++ b/bn_mp_to_signed_bin.c
@@ -16,7 +16,7 @@
  */
 
 /* store in signed [big endian] format */
-int mp_to_signed_bin(mp_int *a, unsigned char *b)
+int mp_to_signed_bin(const mp_int *a, unsigned char *b)
 {
    int     res;
 
diff --git a/bn_mp_to_signed_bin_n.c b/bn_mp_to_signed_bin_n.c
index 501f849..1447624 100644
--- a/bn_mp_to_signed_bin_n.c
+++ b/bn_mp_to_signed_bin_n.c
@@ -16,7 +16,7 @@
  */
 
 /* store in signed [big endian] format */
-int mp_to_signed_bin_n(mp_int *a, unsigned char *b, unsigned long *outlen)
+int mp_to_signed_bin_n(const mp_int *a, unsigned char *b, unsigned long *outlen)
 {
    if (*outlen < (unsigned long)mp_signed_bin_size(a)) {
       return MP_VAL;
diff --git a/bn_mp_to_unsigned_bin.c b/bn_mp_to_unsigned_bin.c
index e103383..9339cce 100644
--- a/bn_mp_to_unsigned_bin.c
+++ b/bn_mp_to_unsigned_bin.c
@@ -16,7 +16,7 @@
  */
 
 /* store in unsigned [big endian] format */
-int mp_to_unsigned_bin(mp_int *a, unsigned char *b)
+int mp_to_unsigned_bin(const mp_int *a, unsigned char *b)
 {
    int     x, res;
    mp_int  t;
diff --git a/bn_mp_to_unsigned_bin_n.c b/bn_mp_to_unsigned_bin_n.c
index 5ee28f1..707dc82 100644
--- a/bn_mp_to_unsigned_bin_n.c
+++ b/bn_mp_to_unsigned_bin_n.c
@@ -16,7 +16,7 @@
  */
 
 /* store in unsigned [big endian] format */
-int mp_to_unsigned_bin_n(mp_int *a, unsigned char *b, unsigned long *outlen)
+int mp_to_unsigned_bin_n(const mp_int *a, unsigned char *b, unsigned long *outlen)
 {
    if (*outlen < (unsigned long)mp_unsigned_bin_size(a)) {
       return MP_VAL;
diff --git a/bn_mp_toom_mul.c b/bn_mp_toom_mul.c
index 8b771bc..3554ea8 100644
--- a/bn_mp_toom_mul.c
+++ b/bn_mp_toom_mul.c
@@ -22,7 +22,7 @@
  * only particularly useful on VERY large inputs
  * (we're talking 1000s of digits here...).
 */
-int mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
+int mp_toom_mul(const mp_int *a, const mp_int *b, mp_int *c)
 {
    mp_int w0, w1, w2, w3, w4, tmp1, tmp2, a0, a1, a2, b0, b1, b2;
    int res, B;
diff --git a/bn_mp_toom_sqr.c b/bn_mp_toom_sqr.c
index 5e1e452..b985435 100644
--- a/bn_mp_toom_sqr.c
+++ b/bn_mp_toom_sqr.c
@@ -16,7 +16,7 @@
  */
 
 /* squaring using Toom-Cook 3-way algorithm */
-int mp_toom_sqr(mp_int *a, mp_int *b)
+int mp_toom_sqr(const mp_int *a, mp_int *b)
 {
    mp_int w0, w1, w2, w3, w4, tmp1, a0, a1, a2;
    int res, B;
diff --git a/bn_mp_toradix.c b/bn_mp_toradix.c
index 5016219..7dd6e4f 100644
--- a/bn_mp_toradix.c
+++ b/bn_mp_toradix.c
@@ -16,7 +16,7 @@
  */
 
 /* stores a bignum as a ASCII string in a given radix (2..64) */
-int mp_toradix(mp_int *a, char *str, int radix)
+int mp_toradix(const mp_int *a, char *str, int radix)
 {
    int     res, digs;
    mp_int  t;
diff --git a/bn_mp_toradix_n.c b/bn_mp_toradix_n.c
index 287d3f8..ef885fc 100644
--- a/bn_mp_toradix_n.c
+++ b/bn_mp_toradix_n.c
@@ -19,7 +19,7 @@
  *
  * Stores upto maxlen-1 chars and always a NULL byte
  */
-int mp_toradix_n(mp_int *a, char *str, int radix, int maxlen)
+int mp_toradix_n(const mp_int *a, char *str, int radix, int maxlen)
 {
    int     res, digs;
    mp_int  t;
diff --git a/bn_mp_unsigned_bin_size.c b/bn_mp_unsigned_bin_size.c
index 174cdc1..04107fe 100644
--- a/bn_mp_unsigned_bin_size.c
+++ b/bn_mp_unsigned_bin_size.c
@@ -16,7 +16,7 @@
  */
 
 /* get the size for an unsigned equivalent */
-int mp_unsigned_bin_size(mp_int *a)
+int mp_unsigned_bin_size(const mp_int *a)
 {
    int     size = mp_count_bits(a);
    return (size / 8) + (((size & 7) != 0) ? 1 : 0);
diff --git a/bn_mp_xor.c b/bn_mp_xor.c
index 4224a7f..9ebc53a 100644
--- a/bn_mp_xor.c
+++ b/bn_mp_xor.c
@@ -16,10 +16,11 @@
  */
 
 /* XOR two ints together */
-int mp_xor(mp_int *a, mp_int *b, mp_int *c)
+int mp_xor(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     res, ix, px;
-   mp_int  t, *x;
+   mp_int  t;
+   const mp_int *x;
 
    if (a->used > b->used) {
       if ((res = mp_init_copy(&t, a)) != MP_OKAY) {
diff --git a/bn_s_mp_add.c b/bn_s_mp_add.c
index 6ba65da..2046722 100644
--- a/bn_s_mp_add.c
+++ b/bn_s_mp_add.c
@@ -16,9 +16,9 @@
  */
 
 /* low level addition, based on HAC pp.594, Algorithm 14.7 */
-int s_mp_add(mp_int *a, mp_int *b, mp_int *c)
+int s_mp_add(const mp_int *a, const mp_int *b, mp_int *c)
 {
-   mp_int *x;
+   const mp_int *x;
    int     olduse, res, min, max;
 
    /* find sizes, we let |a| <= |b| which means we have to sort
diff --git a/bn_s_mp_exptmod.c b/bn_s_mp_exptmod.c
index f8e6b3b..0d0145d 100644
--- a/bn_s_mp_exptmod.c
+++ b/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/bn_s_mp_mul_digs.c b/bn_s_mp_mul_digs.c
index c72c2a8..af13a02 100644
--- a/bn_s_mp_mul_digs.c
+++ b/bn_s_mp_mul_digs.c
@@ -19,7 +19,7 @@
  * HAC pp. 595, Algorithm 14.12  Modified so you can control how
  * many digits of output are created.
  */
-int s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
+int s_mp_mul_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
 {
    mp_int  t;
    int     res, pa, pb, ix, iy;
diff --git a/bn_s_mp_mul_high_digs.c b/bn_s_mp_mul_high_digs.c
index 0e12f3b..37c108e 100644
--- a/bn_s_mp_mul_high_digs.c
+++ b/bn_s_mp_mul_high_digs.c
@@ -18,7 +18,7 @@
 /* multiplies |a| * |b| and does not compute the lower digs digits
  * [meant to get the higher part of the product]
  */
-int s_mp_mul_high_digs(mp_int *a, mp_int *b, mp_int *c, int digs)
+int s_mp_mul_high_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs)
 {
    mp_int  t;
    int     res, pa, pb, ix, iy;
diff --git a/bn_s_mp_sqr.c b/bn_s_mp_sqr.c
index 2910f24..aae06eb 100644
--- a/bn_s_mp_sqr.c
+++ b/bn_s_mp_sqr.c
@@ -16,7 +16,7 @@
  */
 
 /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
-int s_mp_sqr(mp_int *a, mp_int *b)
+int s_mp_sqr(const mp_int *a, mp_int *b)
 {
    mp_int  t;
    int     res, ix, iy, pa;
diff --git a/bn_s_mp_sub.c b/bn_s_mp_sub.c
index f61aa83..52b8096 100644
--- a/bn_s_mp_sub.c
+++ b/bn_s_mp_sub.c
@@ -16,7 +16,7 @@
  */
 
 /* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
-int s_mp_sub(mp_int *a, mp_int *b, mp_int *c)
+int s_mp_sub(const mp_int *a, const mp_int *b, mp_int *c)
 {
    int     olduse, res, min, max;
 
diff --git a/tommath.h b/tommath.h
index 3f49767..fb5fa0b 100644
--- a/tommath.h
+++ b/tommath.h
@@ -223,13 +223,13 @@ int mp_set_long(mp_int *a, unsigned long b);
 int mp_set_long_long(mp_int *a, unsigned long long b);
 
 /* get a 32-bit value */
-unsigned long mp_get_int(mp_int *a);
+unsigned long mp_get_int(const mp_int *a);
 
 /* get a platform dependent unsigned long value */
-unsigned long mp_get_long(mp_int *a);
+unsigned long mp_get_long(const mp_int *a);
 
 /* get a platform dependent unsigned long long value */
-unsigned long long mp_get_long_long(mp_int *a);
+unsigned long long mp_get_long_long(const mp_int *a);
 
 /* initialize and set a digit */
 int mp_init_set(mp_int *a, mp_digit b);
@@ -264,13 +264,13 @@ int mp_lshd(mp_int *a, int b);
 int mp_div_2d(const mp_int *a, int b, mp_int *c, mp_int *d);
 
 /* b = a/2 */
-int mp_div_2(mp_int *a, mp_int *b);
+int mp_div_2(const mp_int *a, mp_int *b);
 
 /* c = a * 2**b, implemented as c = a << b */
 int mp_mul_2d(const mp_int *a, int b, mp_int *c);
 
 /* b = a*2 */
-int mp_mul_2(mp_int *a, mp_int *b);
+int mp_mul_2(const mp_int *a, mp_int *b);
 
 /* c = a mod 2**b */
 int mp_mod_2d(const mp_int *a, int b, mp_int *c);
@@ -288,13 +288,13 @@ int mp_rand(mp_int *a, int digits);
 
 /* ---> binary operations <--- */
 /* c = a XOR b  */
-int mp_xor(mp_int *a, mp_int *b, mp_int *c);
+int mp_xor(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = a OR b */
-int mp_or(mp_int *a, mp_int *b, mp_int *c);
+int mp_or(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = a AND b */
-int mp_and(mp_int *a, mp_int *b, mp_int *c);
+int mp_and(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* ---> Basic arithmetic <--- */
 
@@ -302,7 +302,7 @@ int mp_and(mp_int *a, mp_int *b, mp_int *c);
 int mp_neg(const mp_int *a, mp_int *b);
 
 /* b = |a| */
-int mp_abs(mp_int *a, mp_int *b);
+int mp_abs(const mp_int *a, mp_int *b);
 
 /* compare a to b */
 int mp_cmp(const mp_int *a, const mp_int *b);
@@ -311,22 +311,22 @@ int mp_cmp(const mp_int *a, const mp_int *b);
 int mp_cmp_mag(const mp_int *a, const mp_int *b);
 
 /* c = a + b */
-int mp_add(mp_int *a, mp_int *b, mp_int *c);
+int mp_add(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = a - b */
-int mp_sub(mp_int *a, mp_int *b, mp_int *c);
+int mp_sub(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* c = a * b */
-int mp_mul(mp_int *a, mp_int *b, mp_int *c);
+int mp_mul(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* b = a*a  */
-int mp_sqr(mp_int *a, mp_int *b);
+int mp_sqr(const mp_int *a, mp_int *b);
 
 /* a/b => cb + d == a */
-int mp_div(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
+int mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d);
 
 /* c = a mod b, 0 <= c < b  */
-int mp_mod(mp_int *a, mp_int *b, mp_int *c);
+int mp_mod(const mp_int *a, const mp_int *b, mp_int *c);
 
 /* ---> single digit functions <--- */
 
@@ -334,122 +334,122 @@ int mp_mod(mp_int *a, mp_int *b, mp_int *c);
 int mp_cmp_d(const mp_int *a, mp_digit b);
 
 /* c = a + b */
-int mp_add_d(mp_int *a, mp_digit b, mp_int *c);
+int mp_add_d(const mp_int *a, mp_digit b, mp_int *c);
 
 /* c = a - b */
-int mp_sub_d(mp_int *a, mp_digit b, mp_int *c);
+int mp_sub_d(const mp_int *a, mp_digit b, mp_int *c);
 
 /* c = a * b */
-int mp_mul_d(mp_int *a, mp_digit b, mp_int *c);
+int mp_mul_d(const mp_int *a, mp_digit b, mp_int *c);
 
 /* a/b => cb + d == a */
-int mp_div_d(mp_int *a, mp_digit b, mp_int *c, mp_digit *d);
+int mp_div_d(const mp_int *a, mp_digit b, mp_int *c, mp_digit *d);
 
 /* a/3 => 3c + d == a */
-int mp_div_3(mp_int *a, mp_int *c, mp_digit *d);
+int mp_div_3(const mp_int *a, mp_int *c, mp_digit *d);
 
 /* c = a**b */
-int mp_expt_d(mp_int *a, mp_digit b, mp_int *c);
-int mp_expt_d_ex(mp_int *a, mp_digit b, mp_int *c, int fast);
+int mp_expt_d(const mp_int *a, mp_digit b, mp_int *c);
+int mp_expt_d_ex(const mp_int *a, mp_digit b, mp_int *c, int fast);
 
 /* c = a mod b, 0 <= c < b  */
-int mp_mod_d(mp_int *a, mp_digit b, mp_digit *c);
+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(mp_int *arg, mp_int *ret);
+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 <--- */
 
@@ -464,17 +464,17 @@ int mp_exptmod(mp_int *a, mp_int *b, mp_int *c, mp_int *d);
 extern const mp_digit ltm_prime_tab[PRIME_SIZE];
 
 /* result=1 if a is divisible by one of the first PRIME_SIZE primes */
-int mp_prime_is_divisible(mp_int *a, int *result);
+int mp_prime_is_divisible(const mp_int *a, int *result);
 
 /* performs one Fermat test of "a" using base "b".
  * Sets result to 0 if composite or 1 if probable prime
  */
-int mp_prime_fermat(mp_int *a, mp_int *b, int *result);
+int mp_prime_fermat(const mp_int *a, const mp_int *b, int *result);
 
 /* performs one Miller-Rabin test of "a" using base "b".
  * Sets result to 0 if composite or 1 if probable prime
  */
-int mp_prime_miller_rabin(mp_int *a, mp_int *b, int *result);
+int mp_prime_miller_rabin(const mp_int *a, const mp_int *b, int *result);
 
 /* This gives [for a given bit size] the number of trials required
  * such that Miller-Rabin gives a prob of failure lower than 2^-96
@@ -488,7 +488,7 @@ int mp_prime_rabin_miller_trials(int size);
  *
  * Sets result to 1 if probably prime, 0 otherwise
  */
-int mp_prime_is_prime(mp_int *a, int t, int *result);
+int mp_prime_is_prime(const mp_int *a, int t, int *result);
 
 /* finds the next prime after the number "a" using "t" trials
  * of Miller-Rabin.
@@ -526,24 +526,24 @@ int mp_prime_random_ex(mp_int *a, int t, int size, int flags, ltm_prime_callback
 /* ---> radix conversion <--- */
 int mp_count_bits(const mp_int *a);
 
-int mp_unsigned_bin_size(mp_int *a);
+int mp_unsigned_bin_size(const mp_int *a);
 int mp_read_unsigned_bin(mp_int *a, const unsigned char *b, int c);
-int mp_to_unsigned_bin(mp_int *a, unsigned char *b);
-int mp_to_unsigned_bin_n(mp_int *a, unsigned char *b, unsigned long *outlen);
+int mp_to_unsigned_bin(const mp_int *a, unsigned char *b);
+int mp_to_unsigned_bin_n(const mp_int *a, unsigned char *b, unsigned long *outlen);
 
-int mp_signed_bin_size(mp_int *a);
+int mp_signed_bin_size(const mp_int *a);
 int mp_read_signed_bin(mp_int *a, const unsigned char *b, int c);
-int mp_to_signed_bin(mp_int *a,  unsigned char *b);
-int mp_to_signed_bin_n(mp_int *a, unsigned char *b, unsigned long *outlen);
+int mp_to_signed_bin(const mp_int *a,  unsigned char *b);
+int mp_to_signed_bin_n(const mp_int *a, unsigned char *b, unsigned long *outlen);
 
 int mp_read_radix(mp_int *a, const char *str, int radix);
-int mp_toradix(mp_int *a, char *str, int radix);
-int mp_toradix_n(mp_int *a, char *str, int radix, int maxlen);
+int mp_toradix(const mp_int *a, char *str, int radix);
+int mp_toradix_n(const mp_int *a, char *str, int radix, int maxlen);
 int mp_radix_size(const mp_int *a, int radix, int *size);
 
 #ifndef LTM_NO_FILE
 int mp_fread(mp_int *a, int radix, FILE *stream);
-int mp_fwrite(mp_int *a, int radix, FILE *stream);
+int mp_fwrite(const mp_int *a, int radix, FILE *stream);
 #endif
 
 #define mp_read_raw(mp, str, len) mp_read_signed_bin((mp), (str), (len))
diff --git a/tommath_private.h b/tommath_private.h
index ee39694..087ddcd 100644
--- a/tommath_private.h
+++ b/tommath_private.h
@@ -55,24 +55,24 @@ extern void XFREE(void *p);
 #endif
 
 /* lowlevel functions, do not call! */
-int s_mp_add(mp_int *a, mp_int *b, mp_int *c);
-int s_mp_sub(mp_int *a, mp_int *b, mp_int *c);
+int s_mp_add(const mp_int *a, const mp_int *b, mp_int *c);
+int s_mp_sub(const mp_int *a, const mp_int *b, mp_int *c);
 #define s_mp_mul(a, b, c) s_mp_mul_digs(a, b, c, (a)->used + (b)->used + 1)
-int fast_s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs);
-int s_mp_mul_digs(mp_int *a, mp_int *b, mp_int *c, int digs);
-int fast_s_mp_mul_high_digs(mp_int *a, mp_int *b, mp_int *c, int digs);
-int s_mp_mul_high_digs(mp_int *a, mp_int *b, mp_int *c, int digs);
-int fast_s_mp_sqr(mp_int *a, mp_int *b);
-int s_mp_sqr(mp_int *a, mp_int *b);
-int mp_karatsuba_mul(mp_int *a, mp_int *b, mp_int *c);
-int mp_toom_mul(mp_int *a, mp_int *b, mp_int *c);
-int mp_karatsuba_sqr(mp_int *a, mp_int *b);
-int mp_toom_sqr(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_s_mp_mul_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs);
+int s_mp_mul_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs);
+int fast_s_mp_mul_high_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs);
+int s_mp_mul_high_digs(const mp_int *a, const mp_int *b, mp_int *c, int digs);
+int fast_s_mp_sqr(const mp_int *a, mp_int *b);
+int s_mp_sqr(const mp_int *a, mp_int *b);
+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(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;