LibTomMath Documentation

This example demonstrates using a mix of Libtommath functions to determine if a number $n$ is a perfect power ($n = root^{exponent}$). It uses mp_prime_next_prime to iterate through possible prime exponents and mp_expt_n for exponentiation.

static mp_err mp_is_perfect_power(const mp_int *n, bool *result,
                                    mp_int *rootout, mp_int *exponent)
{
   mp_int root, power, prime, max;
   int highbit, lowbit, p;
   mp_err err = MP_OKAY;

   *result = false;

   if (mp_cmp_d(n, 4) == MP_LT) {
      err = MP_VAL;
      goto LTM_OUT;
   }

   highbit = mp_count_bits(n) - 1;

   if (MP_IS_POWER_OF_TWO(n)) {
      *result = true;
      if (exponent != NULL) {
         if ((err = mp_set_l(exponent, (long)highbit)) != MP_OKAY) return err;
      }
      if (rootout != NULL) {
         if ((err = mp_set_l(rootout, 2l)) != MP_OKAY) return err;
      }
      return err;
   }

   if ((err = mp_init_multi(&root, &power, &prime, &max, NULL)) != MP_OKAY) return err;
   if ((err = mp_set_l(&max, (long)highbit)) != MP_OKAY) goto LTM_ERR;

   mp_set(&prime, 2u);

   while (mp_cmp(&prime, &max) != MP_GT) {
      if ((err = mp_prime_next_prime(&prime, -1, false)) != MP_OKAY) goto LTM_ERR;
      p = (int)mp_get_l(&prime);
      if ((err = mp_root_n(n, p, &root)) != MP_OKAY) goto LTM_ERR;
      if ((err = mp_expt_n(&root, p, &power)) != MP_OKAY) goto LTM_ERR;

      if (mp_cmp(n, &power) == MP_EQ) {
         *result = true;
         if (rootout != NULL) mp_exch(&root, rootout);
         if (exponent != NULL) {
            if ((e = mp_set_l(exponent, (long)p)) != MP_OKAY) goto LTM_ERR;
         }
      }
   }

LTM_OUT:
   if (rootout != NULL) mp_set(rootout, 0u);
   if (exponent != NULL) mp_set(exponent, 0u);
LTM_ERR:
   mp_clear_multi(&root, &power, &prime, NULL);
   return err;
}

For high-performance modular arithmetic, use Montgomery reduction. This requires a pre-computation step.

1. Setup: Use mp_montgomery_setup(const mp_int *a, mp_digit *mp) to compute the pre-computation value mp for an odd modulus a.
2. Normalization: Use mp_montgomery_calc_normalization(&R, b) to compute $R = r^n$.
3. Reduction: Use mp_montgomery_reduce(mp_int *a, const mp_int *m, mp_digit mp) to reduce a in place modulo m using the pre-computed mp.

Example workflow for computing $a^3 ext{ mod } b$:

/* ... initialization and setup ... */

/* normalize 'a' so now a is equal to aR */
mp_mulmod(&a, &R, &b, &a);

/* square a to get c = a^2R^2 */
mp_sqr(&a, &c);

/* reduce 'c' back down to c = a^2R^2 * R^-1 == a^2R */
mp_montgomery_reduce(&c, &b, mp);

/* multiply a to get c = a^3R^2 */
mp_mul(&a, &c, &c);

/* reduce 'c' back down to c = a^3R^2 * R^-1 == a^3R */
mp_montgomery_reduce(&c, &b, mp);

/* reduce 'c' again to get c = a^3R * R^-1 == a^3 */
mp_montgomery_reduce(&c, &b, mp);

/* c now equals a^3 mod b */

int main(void)
{
   mp_int   a, b, c, R;
   mp_digit mp;
   mp_err      result;

/* initialize a,b to desired values,
    * mp_init R, c and set c to 1....
    */

/* get normalization */
   if ((result = mp_montgomery_calc_normalization(&R, b)) != MP_OKAY) {
      printf("Error getting norm.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }
   /* get mp value */
   if ((result = mp_montgomery_setup(&c, &mp)) != MP_OKAY) {
      printf("Error setting up montgomery.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }

/* normalize `a' so now a is equal to aR */
   if ((result = mp_mulmod(&a, &R, &b, &a)) != MP_OKAY) {
      printf("Error computing aR.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }

/* square a to get c = a^2R^2 */
   if ((result = mp_sqr(&a, &c)) != MP_OKAY) {
      printf("Error squaring.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }
   /* now reduce `c' back down to c = a^2R^2 * R^-1 == a^2R */
   if ((result = mp_montgomery_reduce(&c, &b, mp)) != MP_OKAY) {
      printf("Error reducing.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }

/* multiply a to get c = a^3R^2 */
   if ((result = mp_mul(&a, &c, &c)) != MP_OKAY) {
      printf("Error reducing.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }

/* now reduce `c' back down to c = a^3R^2 * R^-1 == a^3R */
   if ((result = mp_montgomery_reduce(&c, &b, mp)) != MP_OKAY) {
      printf("Error reducing.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }

/* now reduce (again) `c' back down to c = a^3R * R^-1 == a^3 */
   if ((result = mp_montgomery_reduce(&c, &b, mp)) != MP_OKAY) {
      printf("Error reducing.  %s",
             mp_error_to_string(result));
      return EXIT_FAILURE;
   }

/* c now equals a^3 mod b */

return EXIT_SUCCESS;
}