Tcl Source Code

      if (c != NULL) {
         mp_zero(c);
      }
      return err;
   }

   /* init our temps */
   if ((err = mp_init_multi(&ta, &tb, &tq, &q, (void *)NULL)) != MP_OKAY) {
      return err;
   }


   mp_set(&tq, 1uL);
   n = mp_count_bits(a) - mp_count_bits(b);
   if ((err = mp_abs(a, &ta)) != MP_OKAY)                         goto LBL_ERR;

      c->sign  = MP_IS_ZERO(c) ? MP_ZPOS : n2;
   }
   if (d != NULL) {
      mp_exch(d, &ta);
      d->sign = MP_IS_ZERO(d) ? MP_ZPOS : n;
   }
LBL_ERR:
   mp_clear_multi(&ta, &tb, &tq, &q, (void *)NULL);
   return err;
}

#else

/* integer signed division.
 * c*b + d == a [e.g. a/b, c=quotient, d=remainder]

      mp_int tmpG, tmpX;
      mp_err err;

      if (!MP_HAS(MP_INVMOD)) {
         return MP_VAL;
      }

      if ((err = mp_init_multi(&tmpG, &tmpX, (void *)NULL)) != MP_OKAY) {
         return err;
      }

      /* first compute 1/G mod P */
      if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
         goto LBL_ERR;
      }

      /* now get |X| */
      if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
         goto LBL_ERR;
      }

      /* and now compute (1/G)**|X| instead of G**X [X < 0] */
      err = mp_exptmod(&tmpG, &tmpX, P, Y);
LBL_ERR:
      mp_clear_multi(&tmpG, &tmpX, (void *)NULL);
      return err;
   }

   /* modified diminished radix reduction */
   if (MP_HAS(MP_REDUCE_IS_2K_L) && MP_HAS(MP_REDUCE_2K_L) && MP_HAS(S_MP_EXPTMOD) &&
       (mp_reduce_is_2k_l(P) == MP_YES)) {
      return s_mp_exptmod(G, X, P, Y, 1);

#include "tommath_private.h"
#ifdef BN_MP_EXTEUCLID_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

/* Extended euclidean algorithm of (a, b) produces
   a*u1 + b*u2 = u3
 */
mp_err 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;
   mp_err err;

   if ((err = mp_init_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* initialize, (u1,u2,u3) = (1,0,a) */
   mp_set(&u1, 1uL);
   if ((err = mp_copy(a, &u3)) != MP_OKAY)                        goto LBL_ERR;

   }
   if (U3 != NULL) {
      mp_exch(U3, &u3);
   }

   err = MP_OKAY;
LBL_ERR:
   mp_clear_multi(&u1, &u2, &u3, &v1, &v2, &v3, &t1, &t2, &t3, &q, &tmp, (void *)NULL);
   return err;
}
#endif

#include "tommath_private.h"
#ifdef BN_MP_LCM_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

/* computes least common multiple as |a*b|/(a, b) */
mp_err mp_lcm(const mp_int *a, const mp_int *b, mp_int *c)
{
   mp_err  err;
   mp_int  t1, t2;


   if ((err = mp_init_multi(&t1, &t2, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* t1 = get the GCD of the two inputs */
   if ((err = mp_gcd(a, b, &t1)) != MP_OKAY) {
      goto LBL_T;
   }

      err = mp_mul(a, &t2, c);
   }

   /* fix the sign to positive */
   c->sign = MP_ZPOS;

LBL_T:
   mp_clear_multi(&t1, &t2, (void *)NULL);
   return err;
}
#endif

   mp_sign neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;

   if (MP_HAS(S_MP_BALANCE_MUL) &&
       /* Check sizes. The smaller one needs to be larger than the Karatsuba cut-off.
        * The bigger one needs to be at least about one MP_KARATSUBA_MUL_CUTOFF bigger
        * to make some sense, but it depends on architecture, OS, position of the
        * stars... so YMMV.
        * Using it to cut the input into slices small enough for s_mp_mul_digs_fast
        * was actually slower on the author's machine, but YMMV.
        */
       (min_len >= MP_KARATSUBA_MUL_CUTOFF) &&
       ((max_len / 2) >= MP_KARATSUBA_MUL_CUTOFF) &&
       /* Not much effect was observed below a ratio of 1:2, but again: YMMV. */
       (max_len >= (2 * min_len))) {
      err = s_mp_balance_mul(a,b,c);

   mp_int T1z, T2z, Np1z, sz, tz;

   int a, ap2, length, i, j;
   mp_err err;

   *result = MP_NO;

   if ((err = mp_init_multi(&T1z, &T2z, &Np1z, &sz, &tz, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   for (a = 0; a < LTM_FROBENIUS_UNDERWOOD_A; a++) {
      /* TODO: That's ugly! No, really, it is! */
      if ((a==2) || (a==4) || (a==7) || (a==8) || (a==10) ||
          (a==14) || (a==18) || (a==23) || (a==26) || (a==28)) {

   mp_set_u32(&T1z, (uint32_t)((2 * a) + 5));
   if ((err = mp_mod(&T1z, N, &T1z)) != MP_OKAY)                  goto LBL_FU_ERR;
   if (MP_IS_ZERO(&sz) && (mp_cmp(&tz, &T1z) == MP_EQ)) {
      *result = MP_YES;
   }

LBL_FU_ERR:
   mp_clear_multi(&tz, &sz, &Np1z, &T2z, &T1z, (void *)NULL);
   return err;
}

#endif
#endif

   such that Jacobi(D,N) = -1 (Selfridge's algorithm). Theory
   indicates that, if N is not a perfect square, D will "nearly
   always" be "small." Just in case, an overflow trap for D is
   included.
   */

   if ((err = mp_init_multi(&Dz, &gcd, &Np1, &Uz, &Vz, &U2mz, &V2mz, &Qmz, &Q2mz, &Qkdz, &T1z, &T2z, &T3z, &T4z, &Q2kdz,
                            (void *)NULL)) != MP_OKAY) {
      return err;
   }

   D = 5;
   sign = 1;

   for (;;) {

      if (r < (s - 1)) {
         if ((err = mp_sqr(&Qkdz, &Qkdz)) != MP_OKAY)             goto LBL_LS_ERR;
         if ((err = mp_mod(&Qkdz, a, &Qkdz)) != MP_OKAY)          goto LBL_LS_ERR;
         if ((err = mp_mul_2(&Qkdz, &Q2kdz)) != MP_OKAY)          goto LBL_LS_ERR;
      }
   }
LBL_LS_ERR:
   mp_clear_multi(&Q2kdz, &T4z, &T3z, &T2z, &T1z, &Qkdz, &Q2mz, &Qmz, &V2mz, &U2mz, &Vz, &Uz, &Np1, &gcd, &Dz, (void *)NULL);
   return err;
}
#endif
#endif
#endif

   }

   /* input must be positive if b is even */
   if (((b & 1) == 0) && mp_isneg(a)) {
      return MP_VAL;
   }

   if ((err = mp_init_multi(&t1, &t2, &t3, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* if a is negative fudge the sign but keep track */
   a_ = *a;
   a_.sign = MP_ZPOS;

   /* set the result */
   mp_exch(&t1, c);

   /* set the sign of the result */
   c->sign = a->sign;

LBL_ERR:
   mp_clear_multi(&t1, &t2, &t3, (void *)NULL);
   return err;
}

#endif

#include "tommath_private.h"
#ifdef BN_MP_SQRT_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

#ifndef NO_FLOATING_POINT
#include <float.h>
#include <math.h>
#if (MP_DIGIT_BIT != 28) || (FLT_RADIX != 2) || (DBL_MANT_DIG != 53) || (DBL_MAX_EXP != 1024)
#define NO_FLOATING_POINT
#endif
#endif

/* this function is less generic than mp_n_root, simpler and faster */
mp_err mp_sqrt(const mp_int *arg, mp_int *ret)
{
   mp_err err;
   mp_int t1, t2;
#ifndef NO_FLOATING_POINT
   int i, j, k;
   volatile double d;
   mp_digit dig;
#endif

   /* must be positive */
   if (arg->sign == MP_NEG) {
      return MP_VAL;
   }

   /* easy out */
   if (MP_IS_ZERO(arg)) {
      mp_zero(ret);
      return MP_OKAY;
   }

#ifndef NO_FLOATING_POINT

   i = (arg->used / 2) - 1;
   j = 2 * i;
   if ((err = mp_init_size(&t1, i+2)) != MP_OKAY) {
      return err;
   }

   if ((err = mp_init(&t2)) != MP_OKAY) {
      goto E2;
   }

   for (k = 0; k < i; ++k) {
      t1.dp[k] = (mp_digit) 0;
   }

   /* Estimate the square root using the hardware floating point unit. */

   d = 0.0;
   for (k = arg->used-1; k >= j; --k) {
      d = ldexp(d, MP_DIGIT_BIT) + (double)(arg->dp[k]);
   }

   /*
    * At this point, d is the nearest floating point number to the most
    * significant 1 or 2 mp_digits of arg. Extract its square root.
    */

   d = sqrt(d);

   /* dig is the most significant mp_digit of the square root */

   dig = (mp_digit) ldexp(d, -MP_DIGIT_BIT);

   /*
    * If the most significant digit is nonzero, find the next digit down
    * by subtracting MP_DIGIT_BIT times thie most significant digit.
    * Subtract one from the result so that our initial estimate is always
    * low.
    */

   if (dig) {
      t1.used = i+2;
      d -= ldexp((double) dig, MP_DIGIT_BIT);
      if (d >= 1.0) {
         t1.dp[i+1] = dig;
         t1.dp[i] = ((mp_digit) d) - 1;
      } else {
         t1.dp[i+1] = dig-1;
         t1.dp[i] = MP_DIGIT_MAX;
      }
   } else {
      t1.used = i+1;
      t1.dp[i] = ((mp_digit) d) - 1;
   }

#else

   if ((err = mp_init_copy(&t1, arg)) != MP_OKAY) {
      return err;
   }

   if ((err = mp_init(&t2)) != MP_OKAY) {
      goto E2;
   }

   /* First approx. (not very bad for large arg) */
   mp_rshd(&t1, t1.used/2);

#endif

   /* t1 > 0  */
   if ((err = mp_div(arg, &t1, &t2, NULL)) != MP_OKAY) {
      goto E1;
   }
   if ((err = mp_add(&t1, &t2, &t1)) != MP_OKAY) {
      goto E1;

      mp_zero(ret);
      return MP_OKAY;
   }
   if (mp_cmp_d(prime, 2uL) == MP_EQ)                            return MP_VAL; /* prime must be odd */
   if ((err = mp_kronecker(n, prime, &legendre)) != MP_OKAY)        return err;
   if (legendre == -1)                                           return MP_VAL; /* quadratic non-residue mod prime */

   if ((err = mp_init_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* SPECIAL CASE: if prime mod 4 == 3
    * compute directly: err = n^(prime+1)/4 mod prime
    * Handbook of Applied Cryptography algorithm 3.36
    */

      if ((err = mp_mulmod(&T, &C, prime, &T)) != MP_OKAY)        goto cleanup;
      /* T = (T * C) mod prime */
      mp_set(&M, i);
      /* M = i */
   }

cleanup:
   mp_clear_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, (void *)NULL);
   return err;
}

#endif

   nblocks = MP_MAX(a->used, b->used) / MP_MIN(a->used, b->used);
   bsize = MP_MIN(a->used, b->used) ;

   if ((err = mp_init_size(&a0, bsize + 2)) != MP_OKAY) {
      return err;
   }
   if ((err = mp_init_multi(&tmp, &r, (void *)NULL)) != MP_OKAY) {
      mp_clear(&a0);
      return err;
   }

   /* Make sure that A is the larger one*/
   if (len_a < len_b) {
      B = *a;

      if ((err = mp_add(&r, &tmp, &r)) != MP_OKAY) {
         goto LBL_ERR;
      }
   }

   mp_exch(&r,c);
LBL_ERR:
   mp_clear_multi(&a0, &tmp, &r, (void *)NULL);
   return err;
}
#endif

   /* 2. [modified] b must be odd   */
   if (MP_IS_EVEN(b)) {
      return MP_VAL;
   }

   /* init all our temps */
   if ((err = mp_init_multi(&x, &y, &u, &v, &B, &D, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* x == modulus, y == value to invert */
   if ((err = mp_copy(b, &x)) != MP_OKAY)                         goto LBL_ERR;

   /* we need y = |a| */

   }

   mp_exch(&D, c);
   c->sign = neg;
   err = MP_OKAY;

LBL_ERR:
   mp_clear_multi(&x, &y, &u, &v, &B, &D, (void *)NULL);
   return err;
}
#endif

   /* b cannot be negative */
   if ((b->sign == MP_NEG) || MP_IS_ZERO(b)) {
      return MP_VAL;
   }

   /* init temps */
   if ((err = mp_init_multi(&x, &y, &u, &v,
                            &A, &B, &C, &D, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* x = a, y = b */
   if ((err = mp_mod(a, b, &x)) != MP_OKAY)                       goto LBL_ERR;
   if ((err = mp_copy(b, &y)) != MP_OKAY)                         goto LBL_ERR;

      if ((err = mp_sub(&C, b, &C)) != MP_OKAY)                   goto LBL_ERR;
   }

   /* C is now the inverse */
   mp_exch(&C, c);
   err = MP_OKAY;
LBL_ERR:
   mp_clear_multi(&x, &y, &u, &v, &A, &B, &C, &D, (void *)NULL);
   return err;
}
#endif

   if ((cmp == MP_LT) || (cmp == MP_EQ)) {
      *c = cmp == MP_EQ;
      return MP_OKAY;
   }

   if ((err =
           mp_init_multi(&bracket_low, &bracket_high,
                         &bracket_mid, &t, &bi_base, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   low = 0;
   mp_set(&bracket_low, 1uL);
   high = 1;

      }
   }

   *c = (mp_cmp(&bracket_high, a) == MP_EQ) ? high : low;

LBL_END:
   mp_clear_multi(&bracket_low, &bracket_high, &bracket_mid,
                  &t, &bi_base, (void *)NULL);
   return err;
}


#endif

#include "tommath_private.h"
#ifdef BN_S_MP_MUL_HIGH_DIGS_FAST_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

/* this is a modified version of s_mp_mul_digs_fast  that only produces
 * output digits *above* digs.  See the comments for s_mp_mul_digs_fast
 * to see how it works.
 *
 * This is used in the Barrett reduction since for one of the multiplications
 * only the higher digits were needed.  This essentially halves the work.
 *
 * Based on Algorithm 14.12 on pp.595 of HAC.
 */

   jenkins_x.c = jenkins_x.d + e;
   jenkins_x.d = e + jenkins_x.a;
   return jenkins_x.d;
}

void s_mp_rand_jenkins_init(uint64_t seed)
{
   int i;
   jenkins_x.a = 0xf1ea5eedULL;
   jenkins_x.b = jenkins_x.c = jenkins_x.d = seed;
   for (i = 0; i < 20; ++i) {
      (void)s_rand_jenkins_val();
   }
}

mp_err s_mp_rand_jenkins(void *p, size_t n)
{
   char *q = (char *)p;

mp_err s_mp_toom_mul(const mp_int *a, const mp_int *b, mp_int *c)
{
   mp_int S1, S2, T1, a0, a1, a2, b0, b1, b2;
   int B, count;
   mp_err err;

   /* init temps */
   if ((err = mp_init_multi(&S1, &S2, &T1, (void *)NULL)) != MP_OKAY) {
      return err;
   }

   /* B */
   B = MP_MIN(a->used, b->used) / 3;

   /** a = a2 * x^2 + a1 * x + a0; */

LBL_ERRb0:
   mp_clear(&a2);
LBL_ERRa2:
   mp_clear(&a1);
LBL_ERRa1:
   mp_clear(&a0);
LBL_ERRa0:
   mp_clear_multi(&S1, &S2, &T1, (void *)NULL);
   return err;
}

#endif

          to other functions like mp_invmod, mp_div, etc...
       -- Sped up mp_exptmod_fast by using new code to find R mod m [e.g. B^n mod m]
       -- minor fixes

Jan 17th, 2003
v0.12  -- re-wrote the majority of the makefile so its more portable and will
          install via "make install" on most *nix platforms
       -- Re-packaged all the source as separate files.  Means the library a single
          file packagage any more.  Instead of just adding "bn.c" you have to add
          libtommath.a
       -- Renamed "bn.h" to "tommath.h"
       -- Changes to the manual to reflect all of this
       -- Used GNU Indent to clean up the source

Jan 15th, 2003