/*

  mpi.c

  by Michael J. Fromberger <http://www.dartmouth.edu/~sting/>

  Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved

  Arbitrary precision integer arithmetic library

  modified for use in Meanwhile as a convenience library

*/

#include "mpi.h"

#include <stdlib.h>

#include <string.h>

#include <ctype.h>

#if MP_DEBUG

#include <stdio.h>

#define DIAG(T,V) {fprintf(stderr,T);mw_mp_print(V,stderr);fputc('\n',stderr);}

#else

#define DIAG(T,V)

#endif

/* 

   If MP_LOGTAB is not defined, use the math library to compute the

   logarithms on the fly.  Otherwise, use the static table below.

   Pick which works best for your system.

 */

#if MP_LOGTAB

/* {{{ s_logv_2[] - log table for 2 in various bases */

/*

  A table of the logs of 2 for various bases (the 0 and 1 entries of

  this table are meaningless and should not be referenced).  

  This table is used to compute output lengths for the mw_mp_toradix()

  function.  Since a number n in radix r takes up about log_r(n)

  digits, we estimate the output size by taking the least integer

  greater than log_r(n), where:

  log_r(n) = log_2(n) * log_r(2)

  This table, therefore, is a table of log_r(2) for 2 <= r <= 36,

  which are the output bases supported.  

 */

#include "logtab.h"

/* }}} */

#define LOG_V_2(R)  s_logv_2[(R)]

#else

#include <math.h>

#define LOG_V_2(R)  (log(2.0)/log(R))

#endif

/* Default precision for newly created mw_mp_int's      */

static unsigned int s_mw_mp_defprec = MP_DEFPREC;

/* {{{ Digit arithmetic macros */

/*

  When adding and multiplying digits, the results can be larger than

  can be contained in an mw_mp_digit.  Thus, an mw_mp_word is used.  These

  macros mask off the upper and lower digits of the mw_mp_word (the

  mw_mp_word may be more than 2 mw_mp_digits wide, but we only concern

  ourselves with the low-order 2 mw_mp_digits)

  If your mw_mp_word DOES have more than 2 mw_mp_digits, you need to

  uncomment the first line, and comment out the second.

 */

/* #define  CARRYOUT(W)  (((W)>>DIGIT_BIT)&MP_DIGIT_MAX) */

#define  CARRYOUT(W)  ((W)>>DIGIT_BIT)

#define  ACCUM(W)     ((W)&MP_DIGIT_MAX)

/* }}} */

/* {{{ Comparison constants */

#define  MP_LT       -1

#define  MP_EQ        0

#define  MP_GT        1

/* }}} */

/* {{{ Constant strings */

/* Constant strings returned by mw_mp_strerror() */

static const char *mw_mp_err_string[] = {

  "unknown result code",     /* say what?            */

  "boolean true",            /* MP_OKAY, MP_YES      */

  "boolean false",           /* MP_NO                */

  "out of memory",           /* MP_MEM               */

  "argument out of range",   /* MP_RANGE             */

  "invalid input parameter", /* MP_BADARG            */

  "result is undefined"      /* MP_UNDEF             */

};

/* Value to digit maps for radix conversion   */

/* s_dmap_1 - standard digits and letters */

static const char *s_dmap_1 = 

  "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";

#if 0

/* s_dmap_2 - base64 ordering for digits  */

static const char *s_dmap_2 =

  "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";

#endif

/* }}} */

/* {{{ Static function declarations */

/* 

   If MP_MACRO is false, these will be defined as actual functions;

   otherwise, suitable macro definitions will be used.  This works

   around the fact that ANSI C89 doesn't support an 'inline' keyword

   (although I hear C9x will ... about bloody time).  At present, the

   macro definitions are identical to the function bodies, but they'll

   expand in place, instead of generating a function call.

   I chose these particular functions to be made into macros because

   some profiling showed they are called a lot on a typical workload,

   and yet they are primarily housekeeping.

 */

#if MP_MACRO == 0

 void     s_mw_mp_setz(mw_mp_digit *dp, mw_mp_size count); /* zero digits           */

 void     s_mw_mp_copy(mw_mp_digit *sp, mw_mp_digit *dp, mw_mp_size count); /* copy    */

 void    *s_mw_mp_alloc(size_t nb, size_t ni);       /* general allocator     */

 void     s_mw_mp_free(void *ptr);                   /* general free function */

#else

 /* Even if these are defined as macros, we need to respect the settings

    of the MP_MEMSET and MP_MEMCPY configuration options...

  */

 #if MP_MEMSET == 0

  #define  s_mw_mp_setz(dp, count) \

       {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=0;}

 #else

  #define  s_mw_mp_setz(dp, count) memset(dp, 0, (count) * sizeof(mw_mp_digit))

 #endif /* MP_MEMSET */

 #if MP_MEMCPY == 0

  #define  s_mw_mp_copy(sp, dp, count) \

       {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=(sp)[ix];}

 #else

  #define  s_mw_mp_copy(sp, dp, count) memcpy(dp, sp, (count) * sizeof(mw_mp_digit))

 #endif /* MP_MEMCPY */

 #define  s_mw_mp_alloc(nb, ni)  calloc(nb, ni)

 #define  s_mw_mp_free(ptr) {if(ptr) free(ptr);}

#endif /* MP_MACRO */

mw_mp_err   s_mw_mp_grow(mw_mp_int *mp, mw_mp_size min);   /* increase allocated size */

mw_mp_err   s_mw_mp_pad(mw_mp_int *mp, mw_mp_size min);    /* left pad with zeroes    */

void     s_mw_mp_clamp(mw_mp_int *mp);               /* clip leading zeroes     */

void     s_mw_mp_exch(mw_mp_int *a, mw_mp_int *b);      /* swap a and b in place   */

mw_mp_err   s_mw_mp_lshd(mw_mp_int *mp, mw_mp_size p);     /* left-shift by p digits  */

void     s_mw_mp_rshd(mw_mp_int *mp, mw_mp_size p);     /* right-shift by p digits */

void     s_mw_mp_div_2d(mw_mp_int *mp, mw_mp_digit d);  /* divide by 2^d in place  */

void     s_mw_mp_mod_2d(mw_mp_int *mp, mw_mp_digit d);  /* modulo 2^d in place     */

mw_mp_err   s_mw_mp_mul_2d(mw_mp_int *mp, mw_mp_digit d);  /* multiply by 2^d in place*/

void     s_mw_mp_div_2(mw_mp_int *mp);               /* divide by 2 in place    */

mw_mp_err   s_mw_mp_mul_2(mw_mp_int *mp);               /* multiply by 2 in place  */

mw_mp_digit s_mw_mp_norm(mw_mp_int *a, mw_mp_int *b);      /* normalize for division  */

mw_mp_err   s_mw_mp_add_d(mw_mp_int *mp, mw_mp_digit d);   /* unsigned digit addition */

mw_mp_err   s_mw_mp_sub_d(mw_mp_int *mp, mw_mp_digit d);   /* unsigned digit subtract */

mw_mp_err   s_mw_mp_mul_d(mw_mp_int *mp, mw_mp_digit d);   /* unsigned digit multiply */

mw_mp_err   s_mw_mp_div_d(mw_mp_int *mp, mw_mp_digit d, mw_mp_digit *r);

		                               /* unsigned digit divide   */

mw_mp_err   s_mw_mp_reduce(mw_mp_int *x, mw_mp_int *m, mw_mp_int *mu);

                                               /* Barrett reduction       */

mw_mp_err   s_mw_mp_add(mw_mp_int *a, mw_mp_int *b);       /* magnitude addition      */

mw_mp_err   s_mw_mp_sub(mw_mp_int *a, mw_mp_int *b);       /* magnitude subtract      */

mw_mp_err   s_mw_mp_mul(mw_mp_int *a, mw_mp_int *b);       /* magnitude multiply      */

#if 0

void     s_mw_mp_kmul(mw_mp_digit *a, mw_mp_digit *b, mw_mp_digit *out, mw_mp_size len);

                                               /* multiply buffers in place */

#endif

#if MP_SQUARE

mw_mp_err   s_mw_mp_sqr(mw_mp_int *a);                  /* magnitude square        */

#else

#define  s_mw_mp_sqr(a) s_mw_mp_mul(a, a)

#endif

mw_mp_err   s_mw_mp_div(mw_mp_int *a, mw_mp_int *b);       /* magnitude divide        */

mw_mp_err   s_mw_mp_2expt(mw_mp_int *a, mw_mp_digit k);    /* a = 2^k                 */

int      s_mw_mp_cmp(mw_mp_int *a, mw_mp_int *b);       /* magnitude comparison    */

int      s_mw_mp_cmw_mp_d(mw_mp_int *a, mw_mp_digit d);    /* magnitude digit compare */

int      s_mw_mp_ispow2(mw_mp_int *v);               /* is v a power of 2?      */

int      s_mw_mp_ispow2d(mw_mp_digit d);             /* is d a power of 2?      */

int      s_mw_mp_tovalue(char ch, int r);          /* convert ch to value    */

char     s_mw_mp_todigit(int val, int r, int low); /* convert val to digit   */

int      s_mw_mp_outlen(int bits, int r);          /* output length in bytes */

/* }}} */

/* {{{ Default precision manipulation */

unsigned int mw_mp_get_prec(void)

{

  return s_mw_mp_defprec;

} /* end mw_mp_get_prec() */

void         mw_mp_set_prec(unsigned int prec)

{

  if(prec == 0)

    s_mw_mp_defprec = MP_DEFPREC;

  else

    s_mw_mp_defprec = prec;

} /* end mw_mp_set_prec() */

/* }}} */

/*------------------------------------------------------------------------*/

/* {{{ mw_mp_init(mp) */

/*

  mw_mp_init(mp)

  Initialize a new zero-valued mw_mp_int.  Returns MP_OKAY if successful,

  MP_MEM if memory could not be allocated for the structure.

 */

mw_mp_err mw_mp_init(mw_mp_int *mp)

{

  return mw_mp_init_size(mp, s_mw_mp_defprec);

} /* end mw_mp_init() */

/* }}} */

/* {{{ mw_mp_init_array(mp[], count) */

mw_mp_err mw_mp_init_array(mw_mp_int mp[], int count)

{

  mw_mp_err  res;

  int     pos;

  ARGCHK(mp !=NULL && count > 0, MP_BADARG);

  for(pos = 0; pos < count; ++pos) {

    if((res = mw_mp_init(&mp[pos])) != MP_OKAY)

      goto CLEANUP;

  }

  return MP_OKAY;

 CLEANUP:

  while(--pos >= 0) 

    mw_mp_clear(&mp[pos]);

  return res;

} /* end mw_mp_init_array() */

/* }}} */

/* {{{ mw_mp_init_size(mp, prec) */

/*

  mw_mp_init_size(mp, prec)

  Initialize a new zero-valued mw_mp_int with at least the given

  precision; returns MP_OKAY if successful, or MP_MEM if memory could

  not be allocated for the structure.

 */

mw_mp_err mw_mp_init_size(mw_mp_int *mp, mw_mp_size prec)

{

  ARGCHK(mp != NULL && prec > 0, MP_BADARG);

  if((DIGITS(mp) = s_mw_mp_alloc(prec, sizeof(mw_mp_digit))) == NULL)

    return MP_MEM;

  SIGN(mp) = MP_ZPOS;

  USED(mp) = 1;

  ALLOC(mp) = prec;

  return MP_OKAY;

} /* end mw_mp_init_size() */

/* }}} */

/* {{{ mw_mp_init_copy(mp, from) */

/*

  mw_mp_init_copy(mp, from)

  Initialize mp as an exact copy of from.  Returns MP_OKAY if

  successful, MP_MEM if memory could not be allocated for the new

  structure.

 */

mw_mp_err mw_mp_init_copy(mw_mp_int *mp, mw_mp_int *from)

{

  ARGCHK(mp != NULL && from != NULL, MP_BADARG);

  if(mp == from)

    return MP_OKAY;

  if((DIGITS(mp) = s_mw_mp_alloc(USED(from), sizeof(mw_mp_digit))) == NULL)

    return MP_MEM;

  s_mw_mp_copy(DIGITS(from), DIGITS(mp), USED(from));

  USED(mp) = USED(from);

  ALLOC(mp) = USED(from);

  SIGN(mp) = SIGN(from);

  return MP_OKAY;

} /* end mw_mp_init_copy() */

/* }}} */

/* {{{ mw_mp_copy(from, to) */

/*

  mw_mp_copy(from, to)

  Copies the mw_mp_int 'from' to the mw_mp_int 'to'.  It is presumed that

  'to' has already been initialized (if not, use mw_mp_init_copy()

  instead). If 'from' and 'to' are identical, nothing happens.

 */

mw_mp_err mw_mp_copy(mw_mp_int *from, mw_mp_int *to)

{

  ARGCHK(from != NULL && to != NULL, MP_BADARG);

  if(from == to)

    return MP_OKAY;

  { /* copy */

    mw_mp_digit   *tmp;

    /*

      If the allocated buffer in 'to' already has enough space to hold

      all the used digits of 'from', we'll re-use it to avoid hitting

      the memory allocater more than necessary; otherwise, we'd have

      to grow anyway, so we just allocate a hunk and make the copy as

      usual

     */

    if(ALLOC(to) >= USED(from)) {

      s_mw_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));

      s_mw_mp_copy(DIGITS(from), DIGITS(to), USED(from));

    } else {

      if((tmp = s_mw_mp_alloc(USED(from), sizeof(mw_mp_digit))) == NULL)

	return MP_MEM;

      s_mw_mp_copy(DIGITS(from), tmp, USED(from));

      if(DIGITS(to) != NULL) {

#if MP_CRYPTO

	s_mw_mp_setz(DIGITS(to), ALLOC(to));

#endif

	s_mw_mp_free(DIGITS(to));

      }

      DIGITS(to) = tmp;

      ALLOC(to) = USED(from);

    }

    /* Copy the precision and sign from the original */

    USED(to) = USED(from);

    SIGN(to) = SIGN(from);

  } /* end copy */

  return MP_OKAY;

} /* end mw_mp_copy() */

/* }}} */

/* {{{ mw_mp_exch(mp1, mp2) */

/*

  mw_mp_exch(mp1, mp2)

  Exchange mp1 and mp2 without allocating any intermediate memory

  (well, unless you count the stack space needed for this call and the

  locals it creates...).  This cannot fail.

 */

void mw_mp_exch(mw_mp_int *mp1, mw_mp_int *mp2)

{

#if MP_ARGCHK == 2

  assert(mp1 != NULL && mp2 != NULL);

#else

  if(mp1 == NULL || mp2 == NULL)

    return;

#endif

  s_mw_mp_exch(mp1, mp2);

} /* end mw_mp_exch() */

/* }}} */

/* {{{ mw_mp_clear(mp) */

/*

  mw_mp_clear(mp)

  Release the storage used by an mw_mp_int, and void its fields so that

  if someone calls mw_mp_clear() again for the same int later, we won't

  get tollchocked.

 */

void   mw_mp_clear(mw_mp_int *mp)

{

  if(mp == NULL)

    return;

  if(DIGITS(mp) != NULL) {

#if MP_CRYPTO

    s_mw_mp_setz(DIGITS(mp), ALLOC(mp));

#endif

    s_mw_mp_free(DIGITS(mp));

    DIGITS(mp) = NULL;

  }

  USED(mp) = 0;

  ALLOC(mp) = 0;

} /* end mw_mp_clear() */

/* }}} */

/* {{{ mw_mp_clear_array(mp[], count) */

void   mw_mp_clear_array(mw_mp_int mp[], int count)

{

  ARGCHK(mp != NULL && count > 0, MP_BADARG);

  while(--count >= 0) 

    mw_mp_clear(&mp[count]);

} /* end mw_mp_clear_array() */

/* }}} */

/* {{{ mw_mp_zero(mp) */

/*

  mw_mp_zero(mp) 

  Set mp to zero.  Does not change the allocated size of the structure,

  and therefore cannot fail (except on a bad argument, which we ignore)

 */

void   mw_mp_zero(mw_mp_int *mp)

{

  if(mp == NULL)

    return;

  s_mw_mp_setz(DIGITS(mp), ALLOC(mp));

  USED(mp) = 1;

  SIGN(mp) = MP_ZPOS;

} /* end mw_mp_zero() */

/* }}} */

/* {{{ mw_mp_set(mp, d) */

void   mw_mp_set(mw_mp_int *mp, mw_mp_digit d)

{

  if(mp == NULL)

    return;

  mw_mp_zero(mp);

  DIGIT(mp, 0) = d;

} /* end mw_mp_set() */

/* }}} */

/* {{{ mw_mp_set_int(mp, z) */

mw_mp_err mw_mp_set_int(mw_mp_int *mp, long z)

{

  int            ix;

  unsigned long  v = abs(z);

  mw_mp_err         res;

  ARGCHK(mp != NULL, MP_BADARG);

  mw_mp_zero(mp);

  if(z == 0)

    return MP_OKAY;  /* shortcut for zero */

  for(ix = sizeof(long) - 1; ix >= 0; ix--) {

    if((res = s_mw_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)

      return res;

    res = s_mw_mp_add_d(mp, 

		     (mw_mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));

    if(res != MP_OKAY)

      return res;

  }

  if(z < 0)

    SIGN(mp) = MP_NEG;

  return MP_OKAY;

} /* end mw_mp_set_int() */

/* }}} */

/*------------------------------------------------------------------------*/

/* {{{ Digit arithmetic */

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

/*

  mw_mp_add_d(a, d, b)

  Compute the sum b = a + d, for a single digit d.  Respects the sign of

  its primary addend (single digits are unsigned anyway).

 */

mw_mp_err mw_mp_add_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *b)

{

  mw_mp_err   res = MP_OKAY;

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

  if((res = mw_mp_copy(a, b)) != MP_OKAY)

    return res;

  if(SIGN(b) == MP_ZPOS) {

    res = s_mw_mp_add_d(b, d);

  } else if(s_mw_mp_cmw_mp_d(b, d) >= 0) {

    res = s_mw_mp_sub_d(b, d);

  } else {

    SIGN(b) = MP_ZPOS;

    DIGIT(b, 0) = d - DIGIT(b, 0);

  }

  return res;

} /* end mw_mp_add_d() */

/* }}} */

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

/*

  mw_mp_sub_d(a, d, b)

  Compute the difference b = a - d, for a single digit d.  Respects the

  sign of its subtrahend (single digits are unsigned anyway).

 */

mw_mp_err mw_mp_sub_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *b)

{

  mw_mp_err   res;

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

  if((res = mw_mp_copy(a, b)) != MP_OKAY)

    return res;

  if(SIGN(b) == MP_NEG) {

    if((res = s_mw_mp_add_d(b, d)) != MP_OKAY)

      return res;

  } else if(s_mw_mp_cmw_mp_d(b, d) >= 0) {

    if((res = s_mw_mp_sub_d(b, d)) != MP_OKAY)

      return res;

  } else {

    mw_mp_neg(b, b);

    DIGIT(b, 0) = d - DIGIT(b, 0);

    SIGN(b) = MP_NEG;

  }

  if(s_mw_mp_cmw_mp_d(b, 0) == 0)

    SIGN(b) = MP_ZPOS;

  return MP_OKAY;

} /* end mw_mp_sub_d() */

/* }}} */

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

/*

  mw_mp_mul_d(a, d, b)

  Compute the product b = a * d, for a single digit d.  Respects the sign

  of its multiplicand (single digits are unsigned anyway)

 */

mw_mp_err mw_mp_mul_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *b)

{

  mw_mp_err  res;

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

  if(d == 0) {

    mw_mp_zero(b);

    return MP_OKAY;

  }

  if((res = mw_mp_copy(a, b)) != MP_OKAY)

    return res;

  res = s_mw_mp_mul_d(b, d);

  return res;

} /* end mw_mp_mul_d() */

/* }}} */

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

mw_mp_err mw_mp_mul_2(mw_mp_int *a, mw_mp_int *c)

{

  mw_mp_err  res;

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

  if((res = mw_mp_copy(a, c)) != MP_OKAY)

    return res;

  return s_mw_mp_mul_2(c);

} /* end mw_mp_mul_2() */

/* }}} */

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

/*

  mw_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).

 */

mw_mp_err mw_mp_div_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *q, mw_mp_digit *r)

{

  mw_mp_err   res;

  mw_mp_digit rem;

  int      pow;

  ARGCHK(a != NULL, MP_BADARG);

  if(d == 0)

    return MP_RANGE;

  /* Shortcut for powers of two ... */

  if((pow = s_mw_mp_ispow2d(d)) >= 0) {

    mw_mp_digit  mask;

    mask = (1 << pow) - 1;

    rem = DIGIT(a, 0) & mask;

    if(q) {

      mw_mp_copy(a, q);

      s_mw_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 = mw_mp_copy(a, q)) != MP_OKAY)

      return res;

    res = s_mw_mp_div_d(q, d, &rem;;

    if(s_mw_mp_cmw_mp_d(q, 0) == MP_EQ)

      SIGN(q) = MP_ZPOS;

  } else {

    mw_mp_int  qp;

    if((res = mw_mp_init_copy(&qp, a)) != MP_OKAY)

      return res;

    res = s_mw_mp_div_d(&qp, d, &rem;;

    if(s_mw_mp_cmw_mp_d(&qp, 0) == 0)

      SIGN(&qp) = MP_ZPOS;

    mw_mp_clear(&qp);

  }

  if(r)

    *r = rem;

  return res;

} /* end mw_mp_div_d() */

/* }}} */

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

/*

  mw_mp_div_2(a, c)

  Compute c = a / 2, disregarding the remainder.

 */

mw_mp_err mw_mp_div_2(mw_mp_int *a, mw_mp_int *c)

{

  mw_mp_err  res;

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

  if((res = mw_mp_copy(a, c)) != MP_OKAY)

    return res;

  s_mw_mp_div_2(c);

  return MP_OKAY;

} /* end mw_mp_div_2() */

/* }}} */

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

mw_mp_err mw_mp_expt_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *c)

{

  mw_mp_int   s, x;

  mw_mp_err   res;

  mw_mp_sign  cs = MP_ZPOS;

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

  if((res = mw_mp_init(&s)) != MP_OKAY)

    return res;

  if((res = mw_mp_init_copy(&x, a)) != MP_OKAY)

    goto X;

  DIGIT(&s, 0) = 1;

  if((d % 2) == 1)

    cs = SIGN(a);

  while(d != 0) {

    if(d & 1) {

      if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY)

	goto CLEANUP;

    }

    d >>= 1;

    if((res = s_mw_mp_sqr(&x)) != MP_OKAY)

      goto CLEANUP;

  }

  SIGN(&s) = cs;

  s_mw_mp_exch(&s, c);

CLEANUP:

  mw_mp_clear(&x);

X:

  mw_mp_clear(&s);

  return res;

} /* end mw_mp_expt_d() */

/* }}} */

/* }}} */

/*------------------------------------------------------------------------*/

/* {{{ Full arithmetic */

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

/*

  mw_mp_abs(a, b)

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

 */

mw_mp_err mw_mp_abs(mw_mp_int *a, mw_mp_int *b)

{

  mw_mp_err   res;

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

  if((res = mw_mp_copy(a, b)) != MP_OKAY)

    return res;

  SIGN(b) = MP_ZPOS;

  return MP_OKAY;

} /* end mw_mp_abs() */

/* }}} */

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

/*

  mw_mp_neg(a, b)

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

 */

mw_mp_err mw_mp_neg(mw_mp_int *a, mw_mp_int *b)

{

  mw_mp_err   res;

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

  if((res = mw_mp_copy(a, b)) != MP_OKAY)

    return res;

  if(s_mw_mp_cmw_mp_d(b, 0) == MP_EQ) 

    SIGN(b) = MP_ZPOS;

  else 

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

  return MP_OKAY;

} /* end mw_mp_neg() */

/* }}} */

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

/*

  mw_mp_add(a, b, c)

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

 */

mw_mp_err mw_mp_add(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)

{

  mw_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