diff --git a/src/addsub/fp_addmod.c b/src/addsub/fp_addmod.c
index fb52a57..f063994 100644
--- a/src/addsub/fp_addmod.c
+++ b/src/addsub/fp_addmod.c
@@ -3,7 +3,7 @@
 #include <tfm_private.h>
 
 /* d = a + b (mod c) */
-int fp_addmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
+int fp_addmod(const fp_int *a, const fp_int *b, const fp_int *c, fp_int *d)
 {
   fp_int tmp;
   fp_zero(&tmp);
diff --git a/src/addsub/fp_submod.c b/src/addsub/fp_submod.c
index f23403c..d00a00d 100644
--- a/src/addsub/fp_submod.c
+++ b/src/addsub/fp_submod.c
@@ -3,7 +3,7 @@
 #include <tfm_private.h>
 
 /* d = a - b (mod c) */
-int fp_submod(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
+int fp_submod(const fp_int *a, const fp_int *b, const fp_int *c, fp_int *d)
 {
   fp_int tmp;
   fp_zero(&tmp);
diff --git a/src/exptmod/fp_exptmod.c b/src/exptmod/fp_exptmod.c
index 3c34262..da3c430 100644
--- a/src/exptmod/fp_exptmod.c
+++ b/src/exptmod/fp_exptmod.c
@@ -74,7 +74,7 @@ static int s_fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
 /* y = g**x (mod b)
  * Some restrictions... x must be positive and < b
  */
-static int s_fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
+static int s_fp_exptmod(const fp_int * G, const fp_int * X, const fp_int * P, fp_int * Y)
 {
   fp_int   M[64], res;
   fp_digit buf, mp;
@@ -232,8 +232,12 @@ static int s_fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
 
 #endif
 
-
-int fp_exptmod(fp_int * G, fp_int * X, fp_int * P, fp_int * Y)
+/*
+ * X should really be const...  however, if it's negative, this function
+ * temporarly changes it to be positive, so we leave it non-const for the
+ * sake of efficiency.
+ */
+int fp_exptmod(const fp_int * G, fp_int * X, const fp_int * P, fp_int * Y)
 {
    fp_int tmp;
    int    err;
diff --git a/src/headers/tfm.h b/src/headers/tfm.h
index 6d6f23a..95811ef 100644
--- a/src/headers/tfm.h
+++ b/src/headers/tfm.h
@@ -415,50 +415,50 @@ int fp_mod_d(const fp_int *a, fp_digit b, fp_digit *c);
 
 /* ---> number theory <--- */
 /* d = a + b (mod c) */
-int fp_addmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
+int fp_addmod(const fp_int *a, const fp_int *b, const fp_int *c, fp_int *d);
 
 /* d = a - b (mod c) */
-int fp_submod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
+int fp_submod(const fp_int *a, const fp_int *b, const fp_int *c, fp_int *d);
 
 /* d = a * b (mod c) */
-int fp_mulmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
+int fp_mulmod(const fp_int *a, const fp_int *b, const fp_int *c, fp_int *d);
 
 /* c = a * a (mod b) */
-int fp_sqrmod(fp_int *a, fp_int *b, fp_int *c);
+int fp_sqrmod(const fp_int *a, const fp_int *b, fp_int *c);
 
 /* c = 1/a (mod b) */
-int fp_invmod(fp_int *a, fp_int *b, fp_int *c);
+int fp_invmod(const fp_int *a, const fp_int *b, fp_int *c);
 
 /* c = (a, b) */
-void fp_gcd(fp_int *a, fp_int *b, fp_int *c);
+void fp_gcd(const fp_int *a, const fp_int *b, fp_int *c);
 
 /* c = [a, b] */
-void fp_lcm(fp_int *a, fp_int *b, fp_int *c);
+void fp_lcm(const fp_int *a, const fp_int *b, fp_int *c);
 
 /* setups the montgomery reduction */
-int fp_montgomery_setup(fp_int *a, fp_digit *mp);
+int fp_montgomery_setup(const fp_int *a, fp_digit *mp);
 
 /* computes a = B**n mod b without division or multiplication useful for
  * normalizing numbers in a Montgomery system.
  */
-void fp_montgomery_calc_normalization(fp_int *a, fp_int *b);
+void fp_montgomery_calc_normalization(fp_int *a, const fp_int *b);
 
 /* computes x/R == x (mod N) via Montgomery Reduction */
-void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp);
+void fp_montgomery_reduce(fp_int *a, const fp_int *m, fp_digit mp);
 
 /* d = a**b (mod c) */
-int fp_exptmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d);
+int fp_exptmod(const fp_int *a, fp_int *b, const fp_int *c, fp_int *d);
 
 /* primality stuff */
 
 /* perform a Miller-Rabin test of a to the base b and store result in "result" */
-void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result);
+void fp_prime_miller_rabin (const fp_int * a, const fp_int * b, int *result);
 
 #define FP_PRIME_SIZE      256
 /* 256 trial divisions + 8 Miller-Rabins, returns FP_YES if probable prime  */
-int fp_isprime(fp_int *a);
+int fp_isprime(const fp_int *a);
 /* extended version of fp_isprime, do 't' Miller-Rabins instead of only 8 */
-int fp_isprime_ex(fp_int *a, int t);
+int fp_isprime_ex(const fp_int *a, int t);
 
 /* Primality generation flags */
 #define TFM_PRIME_BBS      0x0001 /* BBS style prime */
diff --git a/src/mont/fp_montgomery_calc_normalization.c b/src/mont/fp_montgomery_calc_normalization.c
index 6e2e379..b5c1b63 100644
--- a/src/mont/fp_montgomery_calc_normalization.c
+++ b/src/mont/fp_montgomery_calc_normalization.c
@@ -5,7 +5,7 @@
 /* computes a = B**n mod b without division or multiplication useful for
  * normalizing numbers in a Montgomery system.
  */
-void fp_montgomery_calc_normalization(fp_int *a, fp_int *b)
+void fp_montgomery_calc_normalization(fp_int *a, const fp_int *b)
 {
   int     x, bits;
 
diff --git a/src/mont/fp_montgomery_reduce.c b/src/mont/fp_montgomery_reduce.c
index 95d459a..3f9a0fc 100644
--- a/src/mont/fp_montgomery_reduce.c
+++ b/src/mont/fp_montgomery_reduce.c
@@ -456,9 +456,10 @@ asm(                                 \
 #endif
 
 /* computes x/R == x (mod N) via Montgomery Reduction */
-void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
+void fp_montgomery_reduce(fp_int *a, const fp_int *m, fp_digit mp)
 {
-   fp_digit c[FP_SIZE], *_c, *tmpm, mu;
+   const fp_digit *tmpm;
+   fp_digit c[FP_SIZE], *_c, *tmpa, mu;
    int      oldused, x, y, pa;
 
    /* bail if too large */
@@ -519,13 +520,13 @@ void fp_montgomery_reduce(fp_int *a, fp_int *m, fp_digit mp)
 
   /* now copy out */
   _c   = c + pa;
-  tmpm = a->dp;
+  tmpa = a->dp;
   for (x = 0; x < pa+1; x++) {
-     *tmpm++ = *_c++;
+     *tmpa++ = *_c++;
   }
 
   for (; x < oldused; x++)   {
-     *tmpm++ = 0;
+     *tmpa++ = 0;
   }
 
   MONT_FINI;
diff --git a/src/mont/fp_montgomery_setup.c b/src/mont/fp_montgomery_setup.c
index 8b94c9b..4e3740a 100644
--- a/src/mont/fp_montgomery_setup.c
+++ b/src/mont/fp_montgomery_setup.c
@@ -3,7 +3,7 @@
 #include <tfm_private.h>
 
 /* setups the montgomery reduction */
-int fp_montgomery_setup(fp_int *a, fp_digit *rho)
+int fp_montgomery_setup(const fp_int *a, fp_digit *rho)
 {
   fp_digit x, b;
 
diff --git a/src/mul/fp_mulmod.c b/src/mul/fp_mulmod.c
index 44d04be..cb4fe14 100644
--- a/src/mul/fp_mulmod.c
+++ b/src/mul/fp_mulmod.c
@@ -2,7 +2,7 @@
 /* SPDX-License-Identifier: Unlicense */
 #include <tfm_private.h>
 /* d = a * b (mod c) */
-int fp_mulmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
+int fp_mulmod(const fp_int *a, const fp_int *b, const fp_int *c, fp_int *d)
 {
   fp_int tmp;
   fp_zero(&tmp);
diff --git a/src/numtheory/fp_gcd.c b/src/numtheory/fp_gcd.c
index 3bab148..35952ef 100644
--- a/src/numtheory/fp_gcd.c
+++ b/src/numtheory/fp_gcd.c
@@ -3,7 +3,7 @@
 #include <tfm_private.h>
 
 /* c = (a, b) */
-void fp_gcd(fp_int *a, fp_int *b, fp_int *c)
+void fp_gcd(const fp_int *a, const fp_int *b, fp_int *c)
 {
    fp_int u, v, r;
 
diff --git a/src/numtheory/fp_invmod.c b/src/numtheory/fp_invmod.c
index 307f2cb..3001dc0 100644
--- a/src/numtheory/fp_invmod.c
+++ b/src/numtheory/fp_invmod.c
@@ -2,7 +2,7 @@
 /* SPDX-License-Identifier: Unlicense */
 #include <tfm_private.h>
 
-static int s_fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
+static int s_fp_invmod_slow (const fp_int * a, const fp_int * b, fp_int * c)
 {
   fp_int  x, y, u, v, A, B, C, D;
   int     res;
@@ -108,7 +108,7 @@ static int s_fp_invmod_slow (fp_int * a, fp_int * b, fp_int * c)
 }
 
 /* c = 1/a (mod b) for odd b only */
-int fp_invmod(fp_int *a, fp_int *b, fp_int *c)
+int fp_invmod(const fp_int *a, const fp_int *b, fp_int *c)
 {
   fp_int  x, y, u, v, B, D;
   int     neg;
diff --git a/src/numtheory/fp_isprime.c b/src/numtheory/fp_isprime.c
index 837b555..e10c194 100644
--- a/src/numtheory/fp_isprime.c
+++ b/src/numtheory/fp_isprime.c
@@ -2,7 +2,7 @@
 /* SPDX-License-Identifier: Unlicense */
 #include <tfm_private.h>
 
-int fp_isprime(fp_int *a)
+int fp_isprime(const fp_int *a)
 {
   return fp_isprime_ex(a, 8);
 }
diff --git a/src/numtheory/fp_isprime_ex.c b/src/numtheory/fp_isprime_ex.c
index 2b61f8c..b9a47d6 100644
--- a/src/numtheory/fp_isprime_ex.c
+++ b/src/numtheory/fp_isprime_ex.c
@@ -41,7 +41,7 @@ static const fp_digit primes[FP_PRIME_SIZE] = {
   0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653
 };
 
-int fp_isprime_ex(fp_int *a, int t)
+int fp_isprime_ex(const fp_int *a, int t)
 {
    fp_int   b;
    fp_digit d;
diff --git a/src/numtheory/fp_lcm.c b/src/numtheory/fp_lcm.c
index 72d43b5..2b03d41 100644
--- a/src/numtheory/fp_lcm.c
+++ b/src/numtheory/fp_lcm.c
@@ -3,7 +3,7 @@
 #include <tfm_private.h>
 
 /* c = [a, b] */
-void fp_lcm(fp_int *a, fp_int *b, fp_int *c)
+void fp_lcm(const fp_int *a, const fp_int *b, fp_int *c)
 {
    fp_int  t1, t2;
 
diff --git a/src/numtheory/fp_prime_miller_rabin.c b/src/numtheory/fp_prime_miller_rabin.c
index 837d4f3..4c48a70 100644
--- a/src/numtheory/fp_prime_miller_rabin.c
+++ b/src/numtheory/fp_prime_miller_rabin.c
@@ -9,7 +9,7 @@
  * Randomly the chance of error is no more than 1/4 and often
  * very much lower.
  */
-void fp_prime_miller_rabin (fp_int * a, fp_int * b, int *result)
+void fp_prime_miller_rabin (const fp_int * a, const fp_int * b, int *result)
 {
   fp_int  n1, y, r;
   int     s, j;
diff --git a/src/sqr/fp_sqrmod.c b/src/sqr/fp_sqrmod.c
index ebc4e96..7475635 100644
--- a/src/sqr/fp_sqrmod.c
+++ b/src/sqr/fp_sqrmod.c
@@ -3,7 +3,7 @@
 #include <tfm_private.h>
 
 /* c = a * a (mod b) */
-int fp_sqrmod(fp_int *a, fp_int *b, fp_int *c)
+int fp_sqrmod(const fp_int *a, const fp_int *b, fp_int *c)
 {
   fp_int tmp;
   fp_zero(&tmp);