Constify all the number theory functions

This covers:

- fp_addmod()
- fp_submod()
- fp_mulmod()
- fp_sqrmod()
- fp_invmod()
- fp_gcd()
- fp_lcm()
- fp_montgomery_setup()
- fp_montgomery_calc_normalization()
- fp_montgomery_reduce()
- fp_exptmod()
- fp_prime_miller_rabin()
- fp_isprime()
- fp_isprime_ex()

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);

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);

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;

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 */

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;
 

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;

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;
 

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);

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;
 

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;

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);
 }

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;

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;
 

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;

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);