mpi.c

Dropbear is an SSH 2 server, designed to be usable in small memory environments. It supports:

  mp_div_d(a, d, q, r)

  Compute the quotient q = a / d and remainder r = a mod d, for a
  single digit d.  Respects the sign of its divisor (single digits are
  unsigned anyway).
 */

mp_err mp_div_d(mp_int *a, mp_digit d, mp_int *q, mp_digit *r)
{
  mp_err   res;
  mp_digit rem;
  int      pow;

  ARGCHK(a != NULL, MP_BADARG);

  if(d == 0)
    return MP_RANGE;

  /* Shortcut for powers of two ... */
  if((pow = s_mp_ispow2d(d)) >= 0) {
    mp_digit  mask;

    mask = (1 << pow) - 1;
    rem = DIGIT(a, 0) & mask;

    if(q) {
      mp_copy(a, q);
      s_mp_div_2d(q, pow);
    }

    if(r)
      *r = rem;

    return MP_OKAY;
  }

  /*
    If the quotient is actually going to be returned, we'll try to
    avoid hitting the memory allocator by copying the dividend into it
    and doing the division there.  This can't be any _worse_ than
    always copying, and will sometimes be better (since it won't make
    another copy)

    If it's not going to be returned, we need to allocate a temporary
    to hold the quotient, which will just be discarded.
   */
  if(q) {
    if((res = mp_copy(a, q)) != MP_OKAY)
      return res;

    res = s_mp_div_d(q, d, &rem);
    if(s_mp_cmp_d(q, 0) == MP_EQ)
      SIGN(q) = MP_ZPOS;

  } else {
    mp_int  qp;

    if((res = mp_init_copy(&qp, a)) != MP_OKAY)
      return res;

    res = s_mp_div_d(&qp, d, &rem);
    if(s_mp_cmp_d(&qp, 0) == 0)
      SIGN(&qp) = MP_ZPOS;

    mp_clear(&qp);
  }

  if(r)
    *r = rem;

  return res;

} /* end mp_div_d() */

/* }}} */

/* {{{ mp_div_2(a, c) */

/*
  mp_div_2(a, c)

  Compute c = a / 2, disregarding the remainder.
 */

mp_err mp_div_2(mp_int *a, mp_int *c)
{
  mp_err  res;

  ARGCHK(a != NULL && c != NULL, MP_BADARG);

  if((res = mp_copy(a, c)) != MP_OKAY)
    return res;

  s_mp_div_2(c);

  return MP_OKAY;

} /* end mp_div_2() */

/* }}} */

/* {{{ mp_expt_d(a, d, b) */

mp_err mp_expt_d(mp_int *a, mp_digit d, mp_int *c)
{
  mp_int   s, x;
  mp_err   res;

  ARGCHK(a != NULL && c != NULL, MP_BADARG);

  if((res = mp_init(&s)) != MP_OKAY)
    return res;
  if((res = mp_init_copy(&x, a)) != MP_OKAY)
    goto X;

  DIGIT(&s, 0) = 1;

  while(d != 0) {
    if(d & 1) {
      if((res = s_mp_mul(&s, &x)) != MP_OKAY)
	goto CLEANUP;
    }

    d >>= 1;

    if((res = s_mp_sqr(&x)) != MP_OKAY)
      goto CLEANUP;
  }

  s_mp_exch(&s, c);

CLEANUP:
  mp_clear(&x);
X:
  mp_clear(&s);

  return res;

} /* end mp_expt_d() */

/* }}} */

/* }}} */

/*------------------------------------------------------------------------*/
/* {{{ Full arithmetic */

/* {{{ mp_abs(a, b) */

/*
  mp_abs(a, b)

  Compute b = |a|.  'a' and 'b' may be identical.
 */

mp_err mp_abs(mp_int *a, mp_int *b)
{
  mp_err   res;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if((res = mp_copy(a, b)) != MP_OKAY)
    return res;

  SIGN(b) = MP_ZPOS;

  return MP_OKAY;

} /* end mp_abs() */

/* }}} */

/* {{{ mp_neg(a, b) */

/*
  mp_neg(a, b)

  Compute b = -a.  'a' and 'b' may be identical.
 */

mp_err mp_neg(mp_int *a, mp_int *b)
{
  mp_err   res;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if((res = mp_copy(a, b)) != MP_OKAY)
    return res;

  if(s_mp_cmp_d(b, 0) == MP_EQ) 
    SIGN(b) = MP_ZPOS;
  else 
    SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG;

  return MP_OKAY;

} /* end mp_neg() */

/* }}} */

/* {{{ mp_add(a, b, c) */

/*
  mp_add(a, b, c)

  Compute c = a + b.  All parameters may be identical.
 */

mp_err mp_add(mp_int *a, mp_int *b, mp_int *c)
{
  mp_err  res;
  int     cmp;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  if(SIGN(a) == SIGN(b)) { /* same sign:  add values, keep sign */

    /* Commutativity of addition lets us do this in either order,
       so we avoid having to use a temporary even if the result 
       is supposed to replace the output
     */
    if(c == b) {
      if((res = s_mp_add(c, a)) != MP_OKAY)
	return res;
    } else {
      if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
	return res;

      if((res = s_mp_add(c, b)) != MP_OKAY) 
	return res;
    }

  } else if((cmp = s_mp_cmp(a, b)) > 0) {  /* different sign: a > b   */

    /* If the output is going to be clobbered, we will use a temporary
       variable; otherwise, we'll do it without touching the memory 
       allocator at all, if possible
     */
    if(c == b) {
      mp_int  tmp;

      if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
	return res;
      if((res = s_mp_sub(&tmp, b)) != MP_OKAY) {
	mp_clear(&tmp);
	return res;
      }

      s_mp_exch(&tmp, c);
      mp_clear(&tmp);

    } else {

      if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
	return res;
      if((res = s_mp_sub(c, b)) != MP_OKAY)
	return res;

    }

  } else if(cmp == 0) {             /* different sign, a == b   */

    mp_zero(c);
    return MP_OKAY;

  } else {                          /* different sign: a < b    */

    /* See above... */
    if(c == a) {
      mp_int  tmp;

      if((res = mp_init_copy(&tmp, b)) != MP_OKAY)
	return res;
      if((res = s_mp_sub(&tmp, a)) != MP_OKAY) {
	mp_clear(&tmp);
	return res;
      }

      s_mp_exch(&tmp, c);
      mp_clear(&tmp);

    } else {

      if(c != b && (res = mp_copy(b, c)) != MP_OKAY)
	return res;
      if((res = s_mp_sub(c, a)) != MP_OKAY)
	return res;

    }
  }

  if(USED(c) == 1 && DIGIT(c, 0) == 0)
    SIGN(c) = MP_ZPOS;

  return MP_OKAY;

} /* end mp_add() */

/* }}} */

/* {{{ mp_sub(a, b, c) */

/*
  mp_sub(a, b, c)

  Compute c = a - b.  All parameters may be identical.
 */

mp_err mp_sub(mp_int *a, mp_int *b, mp_int *c)
{
  mp_err  res;
  int     cmp;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  if(SIGN(a) != SIGN(b)) {
    if(c == a) {
      if((res = s_mp_add(c, b)) != MP_OKAY)
	return res;
    } else {
      if(c != b && ((res = mp_copy(b, c)) != MP_OKAY))
	return res;
      if((res = s_mp_add(c, a)) != MP_OKAY)
	return res;
      SIGN(c) = SIGN(a);
    }

  } else if((cmp = s_mp_cmp(a, b)) > 0) { /* Same sign, a > b */
    if(c == b) {
      mp_int  tmp;

      if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
	return res;
      if((res = s_mp_sub(&tmp, b)) != MP_OKAY) {
	mp_clear(&tmp);
	return res;
      }
      s_mp_exch(&tmp, c);
      mp_clear(&tmp);

    } else {
      if(c != a && ((res = mp_copy(a, c)) != MP_OKAY))
	return res;

      if((res = s_mp_sub(c, b)) != MP_OKAY)
	return res;
    }

  } else if(cmp == 0) {  /* Same sign, equal magnitude */
    mp_zero(c);
    return MP_OKAY;

  } else {               /* Same sign, b > a */
    if(c == a) {
      mp_int  tmp;

      if((res = mp_init_copy(&tmp, b)) != MP_OKAY)
	return res;

      if((res = s_mp_sub(&tmp, a)) != MP_OKAY) {
	mp_clear(&tmp);
	return res;
      }
      s_mp_exch(&tmp, c);
      mp_clear(&tmp);

    } else {
      if(c != b && ((res = mp_copy(b, c)) != MP_OKAY)) 
	return res;

      if((res = s_mp_sub(c, a)) != MP_OKAY)
	return res;
    }

    SIGN(c) = !SIGN(b);
  }

  if(USED(c) == 1 && DIGIT(c, 0) == 0)
    SIGN(c) = MP_ZPOS;

  return MP_OKAY;

} /* end mp_sub() */

/* }}} */

/* {{{ mp_mul(a, b, c) */

/*
  mp_mul(a, b, c)

  Compute c = a * b.  All parameters may be identical.
 */

mp_err mp_mul(mp_int *a, mp_int *b, mp_int *c)
{
  mp_err   res;
  mp_sign  sgn;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  sgn = (SIGN(a) == SIGN(b)) ? MP_ZPOS : MP_NEG;

  if(c == b) {
    if((res = s_mp_mul(c, a)) != MP_OKAY)
      return res;

  } else {
    if((res = mp_copy(a, c)) != MP_OKAY)
      return res;

    if((res = s_mp_mul(c, b)) != MP_OKAY)
      return res;
  }
  
  if(sgn == MP_ZPOS || s_mp_cmp_d(c, 0) == MP_EQ)
    SIGN(c) = MP_ZPOS;
  else
    SIGN(c) = sgn;
  
  return MP_OKAY;

} /* end mp_mul() */

/* }}} */

/* {{{ mp_mul_2d(a, d, c) */

/*
  mp_mul_2d(a, d, c)

  Compute c = a * 2^d.  a may be the same as c.
 */

mp_err mp_mul_2d(mp_int *a, mp_digit d, mp_int *c)
{
  mp_err   res;

  ARGCHK(a != NULL && c != NULL, MP_BADARG);

  if((res = mp_copy(a, c)) != MP_OKAY)
    return res;

  if(d == 0)
    return MP_OKAY;

  return s_mp_mul_2d(c, d);

} /* end mp_mul() */

/* }}} */

/* {{{ mp_sqr(a, b) */

#if MP_SQUARE
mp_err mp_sqr(mp_int *a, mp_int *b)
{
  mp_err   res;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if((res = mp_copy(a, b)) != MP_OKAY)
    return res;

  if((res = s_mp_sqr(b)) != MP_OKAY)
    return res;

  SIGN(b) = MP_ZPOS;

  return MP_OKAY;

} /* end mp_sqr() */
#endif

/* }}} */

/* {{{ mp_div(a, b, q, r) */

/*
  mp_div(a, b, q, r)

  Compute q = a / b and r = a mod b.  Input parameters may be re-used
  as output parameters.  If q or r is NULL, that portion of the
  computation will be discarded (although it will still be computed)

  Pay no attention to the hacker behind the curtain.
 */

mp_err mp_div(mp_int *a, mp_int *b, mp_int *q, mp_int *r)
{
  mp_err   res;
  mp_int   qtmp, rtmp;
  int      cmp;

  ARGCHK(a != NULL && b != NULL, MP_BADARG);

  if(mp_cmp_z(b) == MP_EQ)
    return MP_RANGE;

  /* If a <= b, we can compute the solution without division, and
     avoid any memory allocation
   */
  if((cmp = s_mp_cmp(a, b)) < 0) {
    if(r) {
      if((res = mp_copy(a, r)) != MP_OKAY)
	return res;
    }

    if(q) 
      mp_zero(q);

    return MP_OKAY;

  } else if(cmp == 0) {

    /* Set quotient to 1, with appropriate sign */
    if(q) {
      int qneg = (SIGN(a) != SIGN(b));

      mp_set(q, 1);
      if(qneg)
	SIGN(q) = MP_NEG;
    }

    if(r)
      mp_zero(r);

    return MP_OKAY;
  }

  /* If we get here, it means we actually have to do some division */

  /* Set up some temporaries... */
  if((res = mp_init_copy(&qtmp, a)) != MP_OKAY)
    return res;
  if((res = mp_init_copy(&rtmp, b)) != MP_OKAY)
    goto CLEANUP;

  if((res = s_mp_div(&qtmp, &rtmp)) != MP_OKAY)
    goto CLEANUP;

  /* Compute the signs for the output  */
  SIGN(&rtmp) = SIGN(a); /* Sr = Sa              */
  if(SIGN(a) == SIGN(b))
    SIGN(&qtmp) = MP_ZPOS;  /* Sq = MP_ZPOS if Sa = Sb */
  else
    SIGN(&qtmp) = MP_NEG;   /* Sq = MP_NEG if Sa != Sb */

  if(s_mp_cmp_d(&qtmp, 0) == MP_EQ)
    SIGN(&qtmp) = MP_ZPOS;
  if(s_mp_cmp_d(&rtmp, 0) == MP_EQ)
    SIGN(&rtmp) = MP_ZPOS;

  /* Copy output, if it is needed      */
  if(q) 
    s_mp_exch(&qtmp, q);

  if(r) 
    s_mp_exch(&rtmp, r);

CLEANUP:
  mp_clear(&rtmp);
  mp_clear(&qtmp);

  return res;

} /* end mp_div() */

/* }}} */

/* {{{ mp_div_2d(a, d, q, r) */

mp_err mp_div_2d(mp_int *a, mp_digit d, mp_int *q, mp_int *r)
{
  mp_err  res;

  ARGCHK(a != NULL, MP_BADARG);

  if(q) {
    if((res = mp_copy(a, q)) != MP_OKAY)
      return res;

    s_mp_div_2d(q, d);
  }

  if(r) {
    if((res = mp_copy(a, r)) != MP_OKAY)
      return res;

    s_mp_mod_2d(r, d);
  }

  return MP_OKAY;

} /* end mp_div_2d() */

/* }}} */

/* {{{ mp_expt(a, b, c) */

/*
  mp_expt(a, b, c)

  Compute c = a ** b, that is, raise a to the b power.  Uses a
  standard iterative square-and-multiply technique.
 */

mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c)
{
  mp_int   s, x;
  mp_err   res;
  mp_digit d;
  int      dig, bit;

  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);

  if(mp_cmp_z(b) < 0)
    return MP_RANGE;

  if((res = mp_init(&s)) != MP_OKAY)
    return res;

  mp_set(&s, 1);

  if((res = mp_init_copy(&x, a)) != MP_OKAY)
    goto X;

  /* Loop over low-order digits in ascending order */
  for(dig = 0; dig < (USED(b) - 1); dig++) {
    d = DIGIT(b, dig);

    /* Loop over bits of each non-maximal digit */
    for(bit = 0; bit < DIGIT_BIT; bit++) {
      if(d & 1) {
	if((res = s_mp_mul(&s, &x)) != MP_OKAY) 
	  goto CLEANUP;
      }

      d >>= 1;
      
      if((res = s_mp_sqr(&x)) != MP_OKAY)
	goto CLEANUP;
    }
  }

  /* Consider now the last digit... */
  d = DIGIT(b, dig);

  while(d) {
    if(d & 1) {
      if((res = s_mp_mul(&s, &x)) != MP_OKAY)
	goto CLEANUP;
    }