Blame src/mpi/mpi.c

Packit 16808d
/*
Packit 16808d
  mpi.c
Packit 16808d
  
Packit 16808d
  by Michael J. Fromberger <http://www.dartmouth.edu/~sting/>
Packit 16808d
  Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved
Packit 16808d
  
Packit 16808d
  Arbitrary precision integer arithmetic library
Packit 16808d
Packit 16808d
  modified for use in Meanwhile as a convenience library
Packit 16808d
*/
Packit 16808d
Packit 16808d
#include "mpi.h"
Packit 16808d
#include <stdlib.h>
Packit 16808d
#include <string.h>
Packit 16808d
#include <ctype.h>
Packit 16808d
Packit 16808d
#if MP_DEBUG
Packit 16808d
#include <stdio.h>
Packit 16808d
Packit 16808d
#define DIAG(T,V) {fprintf(stderr,T);mw_mp_print(V,stderr);fputc('\n',stderr);}
Packit 16808d
#else
Packit 16808d
#define DIAG(T,V)
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* 
Packit 16808d
   If MP_LOGTAB is not defined, use the math library to compute the
Packit 16808d
   logarithms on the fly.  Otherwise, use the static table below.
Packit 16808d
   Pick which works best for your system.
Packit 16808d
 */
Packit 16808d
#if MP_LOGTAB
Packit 16808d
Packit 16808d
/* {{{ s_logv_2[] - log table for 2 in various bases */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  A table of the logs of 2 for various bases (the 0 and 1 entries of
Packit 16808d
  this table are meaningless and should not be referenced).  
Packit 16808d
Packit 16808d
  This table is used to compute output lengths for the mw_mp_toradix()
Packit 16808d
  function.  Since a number n in radix r takes up about log_r(n)
Packit 16808d
  digits, we estimate the output size by taking the least integer
Packit 16808d
  greater than log_r(n), where:
Packit 16808d
Packit 16808d
  log_r(n) = log_2(n) * log_r(2)
Packit 16808d
Packit 16808d
  This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
Packit 16808d
  which are the output bases supported.  
Packit 16808d
 */
Packit 16808d
Packit 16808d
#include "logtab.h"
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
#define LOG_V_2(R)  s_logv_2[(R)]
Packit 16808d
Packit 16808d
#else
Packit 16808d
Packit 16808d
#include <math.h>
Packit 16808d
#define LOG_V_2(R)  (log(2.0)/log(R))
Packit 16808d
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* Default precision for newly created mw_mp_int's      */
Packit 16808d
static unsigned int s_mw_mp_defprec = MP_DEFPREC;
Packit 16808d
Packit 16808d
/* {{{ Digit arithmetic macros */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  When adding and multiplying digits, the results can be larger than
Packit 16808d
  can be contained in an mw_mp_digit.  Thus, an mw_mp_word is used.  These
Packit 16808d
  macros mask off the upper and lower digits of the mw_mp_word (the
Packit 16808d
  mw_mp_word may be more than 2 mw_mp_digits wide, but we only concern
Packit 16808d
  ourselves with the low-order 2 mw_mp_digits)
Packit 16808d
Packit 16808d
  If your mw_mp_word DOES have more than 2 mw_mp_digits, you need to
Packit 16808d
  uncomment the first line, and comment out the second.
Packit 16808d
 */
Packit 16808d
Packit 16808d
/* #define  CARRYOUT(W)  (((W)>>DIGIT_BIT)&MP_DIGIT_MAX) */
Packit 16808d
#define  CARRYOUT(W)  ((W)>>DIGIT_BIT)
Packit 16808d
#define  ACCUM(W)     ((W)&MP_DIGIT_MAX)
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Comparison constants */
Packit 16808d
Packit 16808d
#define  MP_LT       -1
Packit 16808d
#define  MP_EQ        0
Packit 16808d
#define  MP_GT        1
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Constant strings */
Packit 16808d
Packit 16808d
/* Constant strings returned by mw_mp_strerror() */
Packit 16808d
static const char *mw_mp_err_string[] = {
Packit 16808d
  "unknown result code",     /* say what?            */
Packit 16808d
  "boolean true",            /* MP_OKAY, MP_YES      */
Packit 16808d
  "boolean false",           /* MP_NO                */
Packit 16808d
  "out of memory",           /* MP_MEM               */
Packit 16808d
  "argument out of range",   /* MP_RANGE             */
Packit 16808d
  "invalid input parameter", /* MP_BADARG            */
Packit 16808d
  "result is undefined"      /* MP_UNDEF             */
Packit 16808d
};
Packit 16808d
Packit 16808d
/* Value to digit maps for radix conversion   */
Packit 16808d
Packit 16808d
/* s_dmap_1 - standard digits and letters */
Packit 16808d
static const char *s_dmap_1 = 
Packit 16808d
  "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
Packit 16808d
Packit 16808d
#if 0
Packit 16808d
/* s_dmap_2 - base64 ordering for digits  */
Packit 16808d
static const char *s_dmap_2 =
Packit 16808d
  "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Static function declarations */
Packit 16808d
Packit 16808d
/* 
Packit 16808d
   If MP_MACRO is false, these will be defined as actual functions;
Packit 16808d
   otherwise, suitable macro definitions will be used.  This works
Packit 16808d
   around the fact that ANSI C89 doesn't support an 'inline' keyword
Packit 16808d
   (although I hear C9x will ... about bloody time).  At present, the
Packit 16808d
   macro definitions are identical to the function bodies, but they'll
Packit 16808d
   expand in place, instead of generating a function call.
Packit 16808d
Packit 16808d
   I chose these particular functions to be made into macros because
Packit 16808d
   some profiling showed they are called a lot on a typical workload,
Packit 16808d
   and yet they are primarily housekeeping.
Packit 16808d
 */
Packit 16808d
#if MP_MACRO == 0
Packit 16808d
 void     s_mw_mp_setz(mw_mp_digit *dp, mw_mp_size count); /* zero digits           */
Packit 16808d
 void     s_mw_mp_copy(mw_mp_digit *sp, mw_mp_digit *dp, mw_mp_size count); /* copy    */
Packit 16808d
 void    *s_mw_mp_alloc(size_t nb, size_t ni);       /* general allocator     */
Packit 16808d
 void     s_mw_mp_free(void *ptr);                   /* general free function */
Packit 16808d
#else
Packit 16808d
Packit 16808d
 /* Even if these are defined as macros, we need to respect the settings
Packit 16808d
    of the MP_MEMSET and MP_MEMCPY configuration options...
Packit 16808d
  */
Packit 16808d
 #if MP_MEMSET == 0
Packit 16808d
  #define  s_mw_mp_setz(dp, count) \
Packit 16808d
       {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=0;}
Packit 16808d
 #else
Packit 16808d
  #define  s_mw_mp_setz(dp, count) memset(dp, 0, (count) * sizeof(mw_mp_digit))
Packit 16808d
 #endif /* MP_MEMSET */
Packit 16808d
Packit 16808d
 #if MP_MEMCPY == 0
Packit 16808d
  #define  s_mw_mp_copy(sp, dp, count) \
Packit 16808d
       {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=(sp)[ix];}
Packit 16808d
 #else
Packit 16808d
  #define  s_mw_mp_copy(sp, dp, count) memcpy(dp, sp, (count) * sizeof(mw_mp_digit))
Packit 16808d
 #endif /* MP_MEMCPY */
Packit 16808d
Packit 16808d
 #define  s_mw_mp_alloc(nb, ni)  calloc(nb, ni)
Packit 16808d
 #define  s_mw_mp_free(ptr) {if(ptr) free(ptr);}
Packit 16808d
#endif /* MP_MACRO */
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_grow(mw_mp_int *mp, mw_mp_size min);   /* increase allocated size */
Packit 16808d
mw_mp_err   s_mw_mp_pad(mw_mp_int *mp, mw_mp_size min);    /* left pad with zeroes    */
Packit 16808d
Packit 16808d
void     s_mw_mp_clamp(mw_mp_int *mp);               /* clip leading zeroes     */
Packit 16808d
Packit 16808d
void     s_mw_mp_exch(mw_mp_int *a, mw_mp_int *b);      /* swap a and b in place   */
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_lshd(mw_mp_int *mp, mw_mp_size p);     /* left-shift by p digits  */
Packit 16808d
void     s_mw_mp_rshd(mw_mp_int *mp, mw_mp_size p);     /* right-shift by p digits */
Packit 16808d
void     s_mw_mp_div_2d(mw_mp_int *mp, mw_mp_digit d);  /* divide by 2^d in place  */
Packit 16808d
void     s_mw_mp_mod_2d(mw_mp_int *mp, mw_mp_digit d);  /* modulo 2^d in place     */
Packit 16808d
mw_mp_err   s_mw_mp_mul_2d(mw_mp_int *mp, mw_mp_digit d);  /* multiply by 2^d in place*/
Packit 16808d
void     s_mw_mp_div_2(mw_mp_int *mp);               /* divide by 2 in place    */
Packit 16808d
mw_mp_err   s_mw_mp_mul_2(mw_mp_int *mp);               /* multiply by 2 in place  */
Packit 16808d
mw_mp_digit s_mw_mp_norm(mw_mp_int *a, mw_mp_int *b);      /* normalize for division  */
Packit 16808d
mw_mp_err   s_mw_mp_add_d(mw_mp_int *mp, mw_mp_digit d);   /* unsigned digit addition */
Packit 16808d
mw_mp_err   s_mw_mp_sub_d(mw_mp_int *mp, mw_mp_digit d);   /* unsigned digit subtract */
Packit 16808d
mw_mp_err   s_mw_mp_mul_d(mw_mp_int *mp, mw_mp_digit d);   /* unsigned digit multiply */
Packit 16808d
mw_mp_err   s_mw_mp_div_d(mw_mp_int *mp, mw_mp_digit d, mw_mp_digit *r);
Packit 16808d
		                               /* unsigned digit divide   */
Packit 16808d
mw_mp_err   s_mw_mp_reduce(mw_mp_int *x, mw_mp_int *m, mw_mp_int *mu);
Packit 16808d
                                               /* Barrett reduction       */
Packit 16808d
mw_mp_err   s_mw_mp_add(mw_mp_int *a, mw_mp_int *b);       /* magnitude addition      */
Packit 16808d
mw_mp_err   s_mw_mp_sub(mw_mp_int *a, mw_mp_int *b);       /* magnitude subtract      */
Packit 16808d
mw_mp_err   s_mw_mp_mul(mw_mp_int *a, mw_mp_int *b);       /* magnitude multiply      */
Packit 16808d
#if 0
Packit 16808d
void     s_mw_mp_kmul(mw_mp_digit *a, mw_mp_digit *b, mw_mp_digit *out, mw_mp_size len);
Packit 16808d
                                               /* multiply buffers in place */
Packit 16808d
#endif
Packit 16808d
#if MP_SQUARE
Packit 16808d
mw_mp_err   s_mw_mp_sqr(mw_mp_int *a);                  /* magnitude square        */
Packit 16808d
#else
Packit 16808d
#define  s_mw_mp_sqr(a) s_mw_mp_mul(a, a)
Packit 16808d
#endif
Packit 16808d
mw_mp_err   s_mw_mp_div(mw_mp_int *a, mw_mp_int *b);       /* magnitude divide        */
Packit 16808d
mw_mp_err   s_mw_mp_2expt(mw_mp_int *a, mw_mp_digit k);    /* a = 2^k                 */
Packit 16808d
int      s_mw_mp_cmp(mw_mp_int *a, mw_mp_int *b);       /* magnitude comparison    */
Packit 16808d
int      s_mw_mp_cmw_mp_d(mw_mp_int *a, mw_mp_digit d);    /* magnitude digit compare */
Packit 16808d
int      s_mw_mp_ispow2(mw_mp_int *v);               /* is v a power of 2?      */
Packit 16808d
int      s_mw_mp_ispow2d(mw_mp_digit d);             /* is d a power of 2?      */
Packit 16808d
Packit 16808d
int      s_mw_mp_tovalue(char ch, int r);          /* convert ch to value    */
Packit 16808d
char     s_mw_mp_todigit(int val, int r, int low); /* convert val to digit   */
Packit 16808d
int      s_mw_mp_outlen(int bits, int r);          /* output length in bytes */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Default precision manipulation */
Packit 16808d
Packit 16808d
unsigned int mw_mp_get_prec(void)
Packit 16808d
{
Packit 16808d
  return s_mw_mp_defprec;
Packit 16808d
Packit 16808d
} /* end mw_mp_get_prec() */
Packit 16808d
Packit 16808d
void         mw_mp_set_prec(unsigned int prec)
Packit 16808d
{
Packit 16808d
  if(prec == 0)
Packit 16808d
    s_mw_mp_defprec = MP_DEFPREC;
Packit 16808d
  else
Packit 16808d
    s_mw_mp_defprec = prec;
Packit 16808d
Packit 16808d
} /* end mw_mp_set_prec() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ mw_mp_init(mp) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_init(mp)
Packit 16808d
Packit 16808d
  Initialize a new zero-valued mw_mp_int.  Returns MP_OKAY if successful,
Packit 16808d
  MP_MEM if memory could not be allocated for the structure.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_init(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  return mw_mp_init_size(mp, s_mw_mp_defprec);
Packit 16808d
Packit 16808d
} /* end mw_mp_init() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_init_array(mp[], count) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_init_array(mw_mp_int mp[], int count)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
  int     pos;
Packit 16808d
Packit 16808d
  ARGCHK(mp !=NULL && count > 0, MP_BADARG);
Packit 16808d
Packit 16808d
  for(pos = 0; pos < count; ++pos) {
Packit 16808d
    if((res = mw_mp_init(&mp[pos])) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  while(--pos >= 0) 
Packit 16808d
    mw_mp_clear(&mp[pos]);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_init_array() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_init_size(mp, prec) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_init_size(mp, prec)
Packit 16808d
Packit 16808d
  Initialize a new zero-valued mw_mp_int with at least the given
Packit 16808d
  precision; returns MP_OKAY if successful, or MP_MEM if memory could
Packit 16808d
  not be allocated for the structure.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_init_size(mw_mp_int *mp, mw_mp_size prec)
Packit 16808d
{
Packit 16808d
  ARGCHK(mp != NULL && prec > 0, MP_BADARG);
Packit 16808d
Packit 16808d
  if((DIGITS(mp) = s_mw_mp_alloc(prec, sizeof(mw_mp_digit))) == NULL)
Packit 16808d
    return MP_MEM;
Packit 16808d
Packit 16808d
  SIGN(mp) = MP_ZPOS;
Packit 16808d
  USED(mp) = 1;
Packit 16808d
  ALLOC(mp) = prec;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_init_size() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_init_copy(mp, from) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_init_copy(mp, from)
Packit 16808d
Packit 16808d
  Initialize mp as an exact copy of from.  Returns MP_OKAY if
Packit 16808d
  successful, MP_MEM if memory could not be allocated for the new
Packit 16808d
  structure.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_init_copy(mw_mp_int *mp, mw_mp_int *from)
Packit 16808d
{
Packit 16808d
  ARGCHK(mp != NULL && from != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(mp == from)
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  if((DIGITS(mp) = s_mw_mp_alloc(USED(from), sizeof(mw_mp_digit))) == NULL)
Packit 16808d
    return MP_MEM;
Packit 16808d
Packit 16808d
  s_mw_mp_copy(DIGITS(from), DIGITS(mp), USED(from));
Packit 16808d
  USED(mp) = USED(from);
Packit 16808d
  ALLOC(mp) = USED(from);
Packit 16808d
  SIGN(mp) = SIGN(from);
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_init_copy() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_copy(from, to) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_copy(from, to)
Packit 16808d
Packit 16808d
  Copies the mw_mp_int 'from' to the mw_mp_int 'to'.  It is presumed that
Packit 16808d
  'to' has already been initialized (if not, use mw_mp_init_copy()
Packit 16808d
  instead). If 'from' and 'to' are identical, nothing happens.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_copy(mw_mp_int *from, mw_mp_int *to)
Packit 16808d
{
Packit 16808d
  ARGCHK(from != NULL && to != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(from == to)
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  { /* copy */
Packit 16808d
    mw_mp_digit   *tmp;
Packit 16808d
Packit 16808d
    /*
Packit 16808d
      If the allocated buffer in 'to' already has enough space to hold
Packit 16808d
      all the used digits of 'from', we'll re-use it to avoid hitting
Packit 16808d
      the memory allocater more than necessary; otherwise, we'd have
Packit 16808d
      to grow anyway, so we just allocate a hunk and make the copy as
Packit 16808d
      usual
Packit 16808d
     */
Packit 16808d
    if(ALLOC(to) >= USED(from)) {
Packit 16808d
      s_mw_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));
Packit 16808d
      s_mw_mp_copy(DIGITS(from), DIGITS(to), USED(from));
Packit 16808d
      
Packit 16808d
    } else {
Packit 16808d
      if((tmp = s_mw_mp_alloc(USED(from), sizeof(mw_mp_digit))) == NULL)
Packit 16808d
	return MP_MEM;
Packit 16808d
Packit 16808d
      s_mw_mp_copy(DIGITS(from), tmp, USED(from));
Packit 16808d
Packit 16808d
      if(DIGITS(to) != NULL) {
Packit 16808d
#if MP_CRYPTO
Packit 16808d
	s_mw_mp_setz(DIGITS(to), ALLOC(to));
Packit 16808d
#endif
Packit 16808d
	s_mw_mp_free(DIGITS(to));
Packit 16808d
      }
Packit 16808d
Packit 16808d
      DIGITS(to) = tmp;
Packit 16808d
      ALLOC(to) = USED(from);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    /* Copy the precision and sign from the original */
Packit 16808d
    USED(to) = USED(from);
Packit 16808d
    SIGN(to) = SIGN(from);
Packit 16808d
  } /* end copy */
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_copy() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_exch(mp1, mp2) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_exch(mp1, mp2)
Packit 16808d
Packit 16808d
  Exchange mp1 and mp2 without allocating any intermediate memory
Packit 16808d
  (well, unless you count the stack space needed for this call and the
Packit 16808d
  locals it creates...).  This cannot fail.
Packit 16808d
 */
Packit 16808d
Packit 16808d
void mw_mp_exch(mw_mp_int *mp1, mw_mp_int *mp2)
Packit 16808d
{
Packit 16808d
#if MP_ARGCHK == 2
Packit 16808d
  assert(mp1 != NULL && mp2 != NULL);
Packit 16808d
#else
Packit 16808d
  if(mp1 == NULL || mp2 == NULL)
Packit 16808d
    return;
Packit 16808d
#endif
Packit 16808d
Packit 16808d
  s_mw_mp_exch(mp1, mp2);
Packit 16808d
Packit 16808d
} /* end mw_mp_exch() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_clear(mp) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_clear(mp)
Packit 16808d
Packit 16808d
  Release the storage used by an mw_mp_int, and void its fields so that
Packit 16808d
  if someone calls mw_mp_clear() again for the same int later, we won't
Packit 16808d
  get tollchocked.
Packit 16808d
 */
Packit 16808d
Packit 16808d
void   mw_mp_clear(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  if(mp == NULL)
Packit 16808d
    return;
Packit 16808d
Packit 16808d
  if(DIGITS(mp) != NULL) {
Packit 16808d
#if MP_CRYPTO
Packit 16808d
    s_mw_mp_setz(DIGITS(mp), ALLOC(mp));
Packit 16808d
#endif
Packit 16808d
    s_mw_mp_free(DIGITS(mp));
Packit 16808d
    DIGITS(mp) = NULL;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  USED(mp) = 0;
Packit 16808d
  ALLOC(mp) = 0;
Packit 16808d
Packit 16808d
} /* end mw_mp_clear() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_clear_array(mp[], count) */
Packit 16808d
Packit 16808d
void   mw_mp_clear_array(mw_mp_int mp[], int count)
Packit 16808d
{
Packit 16808d
  ARGCHK(mp != NULL && count > 0, MP_BADARG);
Packit 16808d
Packit 16808d
  while(--count >= 0) 
Packit 16808d
    mw_mp_clear(&mp[count]);
Packit 16808d
Packit 16808d
} /* end mw_mp_clear_array() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_zero(mp) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_zero(mp) 
Packit 16808d
Packit 16808d
  Set mp to zero.  Does not change the allocated size of the structure,
Packit 16808d
  and therefore cannot fail (except on a bad argument, which we ignore)
Packit 16808d
 */
Packit 16808d
void   mw_mp_zero(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  if(mp == NULL)
Packit 16808d
    return;
Packit 16808d
Packit 16808d
  s_mw_mp_setz(DIGITS(mp), ALLOC(mp));
Packit 16808d
  USED(mp) = 1;
Packit 16808d
  SIGN(mp) = MP_ZPOS;
Packit 16808d
Packit 16808d
} /* end mw_mp_zero() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_set(mp, d) */
Packit 16808d
Packit 16808d
void   mw_mp_set(mw_mp_int *mp, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  if(mp == NULL)
Packit 16808d
    return;
Packit 16808d
Packit 16808d
  mw_mp_zero(mp);
Packit 16808d
  DIGIT(mp, 0) = d;
Packit 16808d
Packit 16808d
} /* end mw_mp_set() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_set_int(mp, z) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_set_int(mw_mp_int *mp, long z)
Packit 16808d
{
Packit 16808d
  int            ix;
Packit 16808d
  unsigned long  v = abs(z);
Packit 16808d
  mw_mp_err         res;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  mw_mp_zero(mp);
Packit 16808d
  if(z == 0)
Packit 16808d
    return MP_OKAY;  /* shortcut for zero */
Packit 16808d
Packit 16808d
  for(ix = sizeof(long) - 1; ix >= 0; ix--) {
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    res = s_mw_mp_add_d(mp, 
Packit 16808d
		     (mw_mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));
Packit 16808d
    if(res != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(z < 0)
Packit 16808d
    SIGN(mp) = MP_NEG;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_set_int() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ Digit arithmetic */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_add_d(a, d, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_add_d(a, d, b)
Packit 16808d
Packit 16808d
  Compute the sum b = a + d, for a single digit d.  Respects the sign of
Packit 16808d
  its primary addend (single digits are unsigned anyway).
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_add_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res = MP_OKAY;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  if(SIGN(b) == MP_ZPOS) {
Packit 16808d
    res = s_mw_mp_add_d(b, d);
Packit 16808d
  } else if(s_mw_mp_cmw_mp_d(b, d) >= 0) {
Packit 16808d
    res = s_mw_mp_sub_d(b, d);
Packit 16808d
  } else {
Packit 16808d
    SIGN(b) = MP_ZPOS;
Packit 16808d
Packit 16808d
    DIGIT(b, 0) = d - DIGIT(b, 0);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_add_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_sub_d(a, d, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_sub_d(a, d, b)
Packit 16808d
Packit 16808d
  Compute the difference b = a - d, for a single digit d.  Respects the
Packit 16808d
  sign of its subtrahend (single digits are unsigned anyway).
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_sub_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  if(SIGN(b) == MP_NEG) {
Packit 16808d
    if((res = s_mw_mp_add_d(b, d)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
  } else if(s_mw_mp_cmw_mp_d(b, d) >= 0) {
Packit 16808d
    if((res = s_mw_mp_sub_d(b, d)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
  } else {
Packit 16808d
    mw_mp_neg(b, b);
Packit 16808d
Packit 16808d
    DIGIT(b, 0) = d - DIGIT(b, 0);
Packit 16808d
    SIGN(b) = MP_NEG;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(s_mw_mp_cmw_mp_d(b, 0) == 0)
Packit 16808d
    SIGN(b) = MP_ZPOS;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_sub_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mul_d(a, d, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_mul_d(a, d, b)
Packit 16808d
Packit 16808d
  Compute the product b = a * d, for a single digit d.  Respects the sign
Packit 16808d
  of its multiplicand (single digits are unsigned anyway)
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_mul_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(d == 0) {
Packit 16808d
    mw_mp_zero(b);
Packit 16808d
    return MP_OKAY;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  res = s_mw_mp_mul_d(b, d);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_mul_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mul_2(a, c) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_mul_2(mw_mp_int *a, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  return s_mw_mp_mul_2(c);
Packit 16808d
Packit 16808d
} /* end mw_mp_mul_2() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_div_d(a, d, q, r) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_div_d(a, d, q, r)
Packit 16808d
Packit 16808d
  Compute the quotient q = a / d and remainder r = a mod d, for a
Packit 16808d
  single digit d.  Respects the sign of its divisor (single digits are
Packit 16808d
  unsigned anyway).
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_div_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *q, mw_mp_digit *r)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_digit rem;
Packit 16808d
  int      pow;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(d == 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /* Shortcut for powers of two ... */
Packit 16808d
  if((pow = s_mw_mp_ispow2d(d)) >= 0) {
Packit 16808d
    mw_mp_digit  mask;
Packit 16808d
Packit 16808d
    mask = (1 << pow) - 1;
Packit 16808d
    rem = DIGIT(a, 0) & mask;
Packit 16808d
Packit 16808d
    if(q) {
Packit 16808d
      mw_mp_copy(a, q);
Packit 16808d
      s_mw_mp_div_2d(q, pow);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    if(r)
Packit 16808d
      *r = rem;
Packit 16808d
Packit 16808d
    return MP_OKAY;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /*
Packit 16808d
    If the quotient is actually going to be returned, we'll try to
Packit 16808d
    avoid hitting the memory allocator by copying the dividend into it
Packit 16808d
    and doing the division there.  This can't be any _worse_ than
Packit 16808d
    always copying, and will sometimes be better (since it won't make
Packit 16808d
    another copy)
Packit 16808d
Packit 16808d
    If it's not going to be returned, we need to allocate a temporary
Packit 16808d
    to hold the quotient, which will just be discarded.
Packit 16808d
   */
Packit 16808d
  if(q) {
Packit 16808d
    if((res = mw_mp_copy(a, q)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    res = s_mw_mp_div_d(q, d, &rem;;
Packit 16808d
    if(s_mw_mp_cmw_mp_d(q, 0) == MP_EQ)
Packit 16808d
      SIGN(q) = MP_ZPOS;
Packit 16808d
Packit 16808d
  } else {
Packit 16808d
    mw_mp_int  qp;
Packit 16808d
Packit 16808d
    if((res = mw_mp_init_copy(&qp, a)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    res = s_mw_mp_div_d(&qp, d, &rem;;
Packit 16808d
    if(s_mw_mp_cmw_mp_d(&qp, 0) == 0)
Packit 16808d
      SIGN(&qp) = MP_ZPOS;
Packit 16808d
Packit 16808d
    mw_mp_clear(&qp);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(r)
Packit 16808d
    *r = rem;
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_div_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_div_2(a, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_div_2(a, c)
Packit 16808d
Packit 16808d
  Compute c = a / 2, disregarding the remainder.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_div_2(mw_mp_int *a, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  s_mw_mp_div_2(c);
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_div_2() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_expt_d(a, d, b) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_expt_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_int   s, x;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_sign  cs = MP_ZPOS;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_init(&s)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init_copy(&x, a)) != MP_OKAY)
Packit 16808d
    goto X;
Packit 16808d
Packit 16808d
  DIGIT(&s, 0) = 1;
Packit 16808d
Packit 16808d
  if((d % 2) == 1)
Packit 16808d
    cs = SIGN(a);
Packit 16808d
Packit 16808d
  while(d != 0) {
Packit 16808d
    if(d & 1) {
Packit 16808d
      if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    d >>= 1;
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_sqr(&x)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  SIGN(&s) = cs;
Packit 16808d
Packit 16808d
  s_mw_mp_exch(&s, c);
Packit 16808d
Packit 16808d
CLEANUP:
Packit 16808d
  mw_mp_clear(&x);
Packit 16808d
X:
Packit 16808d
  mw_mp_clear(&s);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_expt_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ Full arithmetic */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_abs(a, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_abs(a, b)
Packit 16808d
Packit 16808d
  Compute b = |a|.  'a' and 'b' may be identical.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_abs(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  SIGN(b) = MP_ZPOS;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_abs() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_neg(a, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_neg(a, b)
Packit 16808d
Packit 16808d
  Compute b = -a.  'a' and 'b' may be identical.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_neg(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  if(s_mw_mp_cmw_mp_d(b, 0) == MP_EQ) 
Packit 16808d
    SIGN(b) = MP_ZPOS;
Packit 16808d
  else 
Packit 16808d
    SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_neg() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_add(a, b, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_add(a, b, c)
Packit 16808d
Packit 16808d
  Compute c = a + b.  All parameters may be identical.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_add(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
  int     cmp;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(SIGN(a) == SIGN(b)) { /* same sign:  add values, keep sign */
Packit 16808d
Packit 16808d
    /* Commutativity of addition lets us do this in either order,
Packit 16808d
       so we avoid having to use a temporary even if the result 
Packit 16808d
       is supposed to replace the output
Packit 16808d
     */
Packit 16808d
    if(c == b) {
Packit 16808d
      if((res = s_mw_mp_add(c, a)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
    } else {
Packit 16808d
      if(c != a && (res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
      if((res = s_mw_mp_add(c, b)) != MP_OKAY) 
Packit 16808d
	return res;
Packit 16808d
    }
Packit 16808d
Packit 16808d
  } else if((cmp = s_mw_mp_cmp(a, b)) > 0) {  /* different sign: a > b   */
Packit 16808d
Packit 16808d
    /* If the output is going to be clobbered, we will use a temporary
Packit 16808d
       variable; otherwise, we'll do it without touching the memory 
Packit 16808d
       allocator at all, if possible
Packit 16808d
     */
Packit 16808d
    if(c == b) {
Packit 16808d
      mw_mp_int  tmp;
Packit 16808d
Packit 16808d
      if((res = mw_mp_init_copy(&tmp, a)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      if((res = s_mw_mp_sub(&tmp, b)) != MP_OKAY) {
Packit 16808d
	mw_mp_clear(&tmp);
Packit 16808d
	return res;
Packit 16808d
      }
Packit 16808d
Packit 16808d
      s_mw_mp_exch(&tmp, c);
Packit 16808d
      mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
    } else {
Packit 16808d
Packit 16808d
      if(c != a && (res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      if((res = s_mw_mp_sub(c, b)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
    }
Packit 16808d
Packit 16808d
  } else if(cmp == 0) {             /* different sign, a == b   */
Packit 16808d
Packit 16808d
    mw_mp_zero(c);
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  } else {                          /* different sign: a < b    */
Packit 16808d
Packit 16808d
    /* See above... */
Packit 16808d
    if(c == a) {
Packit 16808d
      mw_mp_int  tmp;
Packit 16808d
Packit 16808d
      if((res = mw_mp_init_copy(&tmp, b)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      if((res = s_mw_mp_sub(&tmp, a)) != MP_OKAY) {
Packit 16808d
	mw_mp_clear(&tmp);
Packit 16808d
	return res;
Packit 16808d
      }
Packit 16808d
Packit 16808d
      s_mw_mp_exch(&tmp, c);
Packit 16808d
      mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
    } else {
Packit 16808d
Packit 16808d
      if(c != b && (res = mw_mp_copy(b, c)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      if((res = s_mw_mp_sub(c, a)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
    }
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(USED(c) == 1 && DIGIT(c, 0) == 0)
Packit 16808d
    SIGN(c) = MP_ZPOS;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_add() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_sub(a, b, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_sub(a, b, c)
Packit 16808d
Packit 16808d
  Compute c = a - b.  All parameters may be identical.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_sub(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
  int     cmp;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(SIGN(a) != SIGN(b)) {
Packit 16808d
    if(c == a) {
Packit 16808d
      if((res = s_mw_mp_add(c, b)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
    } else {
Packit 16808d
      if(c != b && ((res = mw_mp_copy(b, c)) != MP_OKAY))
Packit 16808d
	return res;
Packit 16808d
      if((res = s_mw_mp_add(c, a)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      SIGN(c) = SIGN(a);
Packit 16808d
    }
Packit 16808d
Packit 16808d
  } else if((cmp = s_mw_mp_cmp(a, b)) > 0) { /* Same sign, a > b */
Packit 16808d
    if(c == b) {
Packit 16808d
      mw_mp_int  tmp;
Packit 16808d
Packit 16808d
      if((res = mw_mp_init_copy(&tmp, a)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      if((res = s_mw_mp_sub(&tmp, b)) != MP_OKAY) {
Packit 16808d
	mw_mp_clear(&tmp);
Packit 16808d
	return res;
Packit 16808d
      }
Packit 16808d
      s_mw_mp_exch(&tmp, c);
Packit 16808d
      mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
    } else {
Packit 16808d
      if(c != a && ((res = mw_mp_copy(a, c)) != MP_OKAY))
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
      if((res = s_mw_mp_sub(c, b)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
    }
Packit 16808d
Packit 16808d
  } else if(cmp == 0) {  /* Same sign, equal magnitude */
Packit 16808d
    mw_mp_zero(c);
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  } else {               /* Same sign, b > a */
Packit 16808d
    if(c == a) {
Packit 16808d
      mw_mp_int  tmp;
Packit 16808d
Packit 16808d
      if((res = mw_mp_init_copy(&tmp, b)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
      if((res = s_mw_mp_sub(&tmp, a)) != MP_OKAY) {
Packit 16808d
	mw_mp_clear(&tmp);
Packit 16808d
	return res;
Packit 16808d
      }
Packit 16808d
      s_mw_mp_exch(&tmp, c);
Packit 16808d
      mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
    } else {
Packit 16808d
      if(c != b && ((res = mw_mp_copy(b, c)) != MP_OKAY)) 
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
      if((res = s_mw_mp_sub(c, a)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    SIGN(c) = !SIGN(b);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(USED(c) == 1 && DIGIT(c, 0) == 0)
Packit 16808d
    SIGN(c) = MP_ZPOS;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_sub() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mul(a, b, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_mul(a, b, c)
Packit 16808d
Packit 16808d
  Compute c = a * b.  All parameters may be identical.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_mul(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_sign  sgn;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  sgn = (SIGN(a) == SIGN(b)) ? MP_ZPOS : MP_NEG;
Packit 16808d
Packit 16808d
  if(c == b) {
Packit 16808d
    if((res = s_mw_mp_mul(c, a)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
  } else {
Packit 16808d
    if((res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_mul(c, b)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
  }
Packit 16808d
  
Packit 16808d
  if(sgn == MP_ZPOS || s_mw_mp_cmw_mp_d(c, 0) == MP_EQ)
Packit 16808d
    SIGN(c) = MP_ZPOS;
Packit 16808d
  else
Packit 16808d
    SIGN(c) = sgn;
Packit 16808d
  
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_mul() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mul_2d(a, d, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_mul_2d(a, d, c)
Packit 16808d
Packit 16808d
  Compute c = a * 2^d.  a may be the same as c.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_mul_2d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  if(d == 0)
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  return s_mw_mp_mul_2d(c, d);
Packit 16808d
Packit 16808d
} /* end mw_mp_mul() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_sqr(a, b) */
Packit 16808d
Packit 16808d
#if MP_SQUARE
Packit 16808d
mw_mp_err mw_mp_sqr(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_copy(a, b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  if((res = s_mw_mp_sqr(b)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  SIGN(b) = MP_ZPOS;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_sqr() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_div(a, b, q, r) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_div(a, b, q, r)
Packit 16808d
Packit 16808d
  Compute q = a / b and r = a mod b.  Input parameters may be re-used
Packit 16808d
  as output parameters.  If q or r is NULL, that portion of the
Packit 16808d
  computation will be discarded (although it will still be computed)
Packit 16808d
Packit 16808d
  Pay no attention to the hacker behind the curtain.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_div(mw_mp_int *a, mw_mp_int *b, mw_mp_int *q, mw_mp_int *r)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_int   qtmp, rtmp;
Packit 16808d
  int      cmp;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(b) == MP_EQ)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /* If a <= b, we can compute the solution without division, and
Packit 16808d
     avoid any memory allocation
Packit 16808d
   */
Packit 16808d
  if((cmp = s_mw_mp_cmp(a, b)) < 0) {
Packit 16808d
    if(r) {
Packit 16808d
      if((res = mw_mp_copy(a, r)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    if(q) 
Packit 16808d
      mw_mp_zero(q);
Packit 16808d
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  } else if(cmp == 0) {
Packit 16808d
Packit 16808d
    /* Set quotient to 1, with appropriate sign */
Packit 16808d
    if(q) {
Packit 16808d
      int qneg = (SIGN(a) != SIGN(b));
Packit 16808d
Packit 16808d
      mw_mp_set(q, 1);
Packit 16808d
      if(qneg)
Packit 16808d
	SIGN(q) = MP_NEG;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    if(r)
Packit 16808d
      mw_mp_zero(r);
Packit 16808d
Packit 16808d
    return MP_OKAY;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* If we get here, it means we actually have to do some division */
Packit 16808d
Packit 16808d
  /* Set up some temporaries... */
Packit 16808d
  if((res = mw_mp_init_copy(&qtmp, a)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init_copy(&rtmp, b)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
Packit 16808d
  if((res = s_mw_mp_div(&qtmp, &rtmp)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
Packit 16808d
  /* Compute the signs for the output  */
Packit 16808d
  SIGN(&rtmp) = SIGN(a); /* Sr = Sa              */
Packit 16808d
  if(SIGN(a) == SIGN(b))
Packit 16808d
    SIGN(&qtmp) = MP_ZPOS;  /* Sq = MP_ZPOS if Sa = Sb */
Packit 16808d
  else
Packit 16808d
    SIGN(&qtmp) = MP_NEG;   /* Sq = MP_NEG if Sa != Sb */
Packit 16808d
Packit 16808d
  if(s_mw_mp_cmw_mp_d(&qtmp, 0) == MP_EQ)
Packit 16808d
    SIGN(&qtmp) = MP_ZPOS;
Packit 16808d
  if(s_mw_mp_cmw_mp_d(&rtmp, 0) == MP_EQ)
Packit 16808d
    SIGN(&rtmp) = MP_ZPOS;
Packit 16808d
Packit 16808d
  /* Copy output, if it is needed      */
Packit 16808d
  if(q) 
Packit 16808d
    s_mw_mp_exch(&qtmp, q);
Packit 16808d
Packit 16808d
  if(r) 
Packit 16808d
    s_mw_mp_exch(&rtmp, r);
Packit 16808d
Packit 16808d
CLEANUP:
Packit 16808d
  mw_mp_clear(&rtmp);
Packit 16808d
  mw_mp_clear(&qtmp);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_div() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_div_2d(a, d, q, r) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_div_2d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *q, mw_mp_int *r)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(q) {
Packit 16808d
    if((res = mw_mp_copy(a, q)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    s_mw_mp_div_2d(q, d);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(r) {
Packit 16808d
    if((res = mw_mp_copy(a, r)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    s_mw_mp_mod_2d(r, d);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_div_2d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_expt(a, b, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_expt(a, b, c)
Packit 16808d
Packit 16808d
  Compute c = a ** b, that is, raise a to the b power.  Uses a
Packit 16808d
  standard iterative square-and-multiply technique.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_expt(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_int   s, x;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_digit d;
Packit 16808d
  int      dig, bit;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(b) < 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  if((res = mw_mp_init(&s)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  mw_mp_set(&s, 1);
Packit 16808d
Packit 16808d
  if((res = mw_mp_init_copy(&x, a)) != MP_OKAY)
Packit 16808d
    goto X;
Packit 16808d
Packit 16808d
  /* Loop over low-order digits in ascending order */
Packit 16808d
  for(dig = 0; dig < (USED(b) - 1); dig++) {
Packit 16808d
    d = DIGIT(b, dig);
Packit 16808d
Packit 16808d
    /* Loop over bits of each non-maximal digit */
Packit 16808d
    for(bit = 0; bit < DIGIT_BIT; bit++) {
Packit 16808d
      if(d & 1) {
Packit 16808d
	if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY) 
Packit 16808d
	  goto CLEANUP;
Packit 16808d
      }
Packit 16808d
Packit 16808d
      d >>= 1;
Packit 16808d
      
Packit 16808d
      if((res = s_mw_mp_sqr(&x)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
    }
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Consider now the last digit... */
Packit 16808d
  d = DIGIT(b, dig);
Packit 16808d
Packit 16808d
  while(d) {
Packit 16808d
    if(d & 1) {
Packit 16808d
      if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    d >>= 1;
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_sqr(&x)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
  
Packit 16808d
  if(mw_mp_iseven(b))
Packit 16808d
    SIGN(&s) = SIGN(a);
Packit 16808d
Packit 16808d
  res = mw_mp_copy(&s, c);
Packit 16808d
Packit 16808d
CLEANUP:
Packit 16808d
  mw_mp_clear(&x);
Packit 16808d
X:
Packit 16808d
  mw_mp_clear(&s);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_expt() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_2expt(a, k) */
Packit 16808d
Packit 16808d
/* Compute a = 2^k */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_2expt(mw_mp_int *a, mw_mp_digit k)
Packit 16808d
{
Packit 16808d
  ARGCHK(a != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  return s_mw_mp_2expt(a, k);
Packit 16808d
Packit 16808d
} /* end mw_mp_2expt() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mod(a, m, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_mod(a, m, c)
Packit 16808d
Packit 16808d
  Compute c = a (mod m).  Result will always be 0 <= c < m.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_mod(mw_mp_int *a, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
  int     mag;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(SIGN(m) == MP_NEG)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /*
Packit 16808d
     If |a| > m, we need to divide to get the remainder and take the
Packit 16808d
     absolute value.  
Packit 16808d
Packit 16808d
     If |a| < m, we don't need to do any division, just copy and adjust
Packit 16808d
     the sign (if a is negative).
Packit 16808d
Packit 16808d
     If |a| == m, we can simply set the result to zero.
Packit 16808d
Packit 16808d
     This order is intended to minimize the average path length of the
Packit 16808d
     comparison chain on common workloads -- the most frequent cases are
Packit 16808d
     that |a| != m, so we do those first.
Packit 16808d
   */
Packit 16808d
  if((mag = s_mw_mp_cmp(a, m)) > 0) {
Packit 16808d
    if((res = mw_mp_div(a, m, NULL, c)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    
Packit 16808d
    if(SIGN(c) == MP_NEG) {
Packit 16808d
      if((res = mw_mp_add(c, m, c)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
    }
Packit 16808d
Packit 16808d
  } else if(mag < 0) {
Packit 16808d
    if((res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    if(mw_mp_cmw_mp_z(a) < 0) {
Packit 16808d
      if((res = mw_mp_add(c, m, c)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
Packit 16808d
    }
Packit 16808d
    
Packit 16808d
  } else {
Packit 16808d
    mw_mp_zero(c);
Packit 16808d
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_mod() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mod_d(a, d, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_mod_d(a, d, c)
Packit 16808d
Packit 16808d
  Compute c = a (mod d).  Result will always be 0 <= c < d
Packit 16808d
 */
Packit 16808d
mw_mp_err mw_mp_mod_d(mw_mp_int *a, mw_mp_digit d, mw_mp_digit *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_digit rem;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(s_mw_mp_cmw_mp_d(a, d) > 0) {
Packit 16808d
    if((res = mw_mp_div_d(a, d, NULL, &rem)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
  } else {
Packit 16808d
    if(SIGN(a) == MP_NEG)
Packit 16808d
      rem = d - DIGIT(a, 0);
Packit 16808d
    else
Packit 16808d
      rem = DIGIT(a, 0);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(c)
Packit 16808d
    *c = rem;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_mod_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_sqrt(a, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_sqrt(a, b)
Packit 16808d
Packit 16808d
  Compute the integer square root of a, and store the result in b.
Packit 16808d
  Uses an integer-arithmetic version of Newton's iterative linear
Packit 16808d
  approximation technique to determine this value; the result has the
Packit 16808d
  following two properties:
Packit 16808d
Packit 16808d
     b^2 <= a
Packit 16808d
     (b+1)^2 >= a
Packit 16808d
Packit 16808d
  It is a range error to pass a negative value.
Packit 16808d
 */
Packit 16808d
mw_mp_err mw_mp_sqrt(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_int   x, t;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  /* Cannot take square root of a negative value */
Packit 16808d
  if(SIGN(a) == MP_NEG)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /* Special cases for zero and one, trivial     */
Packit 16808d
  if(mw_mp_cmw_mp_d(a, 0) == MP_EQ || mw_mp_cmw_mp_d(a, 1) == MP_EQ) 
Packit 16808d
    return mw_mp_copy(a, b);
Packit 16808d
    
Packit 16808d
  /* Initialize the temporaries we'll use below  */
Packit 16808d
  if((res = mw_mp_init_size(&t, USED(a))) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  /* Compute an initial guess for the iteration as a itself */
Packit 16808d
  if((res = mw_mp_init_copy(&x, a)) != MP_OKAY)
Packit 16808d
    goto X;
Packit 16808d
Packit 16808d
  for(;;) {
Packit 16808d
    /* t = (x * x) - a */
Packit 16808d
    mw_mp_copy(&x, &t);      /* can't fail, t is big enough for original x */
Packit 16808d
    if((res = mw_mp_sqr(&t, &t)) != MP_OKAY ||
Packit 16808d
       (res = mw_mp_sub(&t, a, &t)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
Packit 16808d
    /* t = t / 2x       */
Packit 16808d
    s_mw_mp_mul_2(&x);
Packit 16808d
    if((res = mw_mp_div(&t, &x, &t, NULL)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
    s_mw_mp_div_2(&x);
Packit 16808d
Packit 16808d
    /* Terminate the loop, if the quotient is zero */
Packit 16808d
    if(mw_mp_cmw_mp_z(&t) == MP_EQ)
Packit 16808d
      break;
Packit 16808d
Packit 16808d
    /* x = x - t       */
Packit 16808d
    if((res = mw_mp_sub(&x, &t, &x)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Copy result to output parameter */
Packit 16808d
  mw_mp_sub_d(&x, 1, &x);
Packit 16808d
  s_mw_mp_exch(&x, b);
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  mw_mp_clear(&x);
Packit 16808d
 X:
Packit 16808d
  mw_mp_clear(&t); 
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_sqrt() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ Modular arithmetic */
Packit 16808d
Packit 16808d
#if MP_MODARITH
Packit 16808d
/* {{{ mw_mp_addmod(a, b, m, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_addmod(a, b, m, c)
Packit 16808d
Packit 16808d
  Compute c = (a + b) mod m
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_addmod(mw_mp_int *a, mw_mp_int *b, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_add(a, b, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_mod(c, m, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
}
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_submod(a, b, m, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_submod(a, b, m, c)
Packit 16808d
Packit 16808d
  Compute c = (a - b) mod m
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_submod(mw_mp_int *a, mw_mp_int *b, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_sub(a, b, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_mod(c, m, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
}
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_mulmod(a, b, m, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_mulmod(a, b, m, c)
Packit 16808d
Packit 16808d
  Compute c = (a * b) mod m
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_mulmod(mw_mp_int *a, mw_mp_int *b, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_mul(a, b, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_mod(c, m, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
}
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_sqrmod(a, m, c) */
Packit 16808d
Packit 16808d
#if MP_SQUARE
Packit 16808d
mw_mp_err mw_mp_sqrmod(mw_mp_int *a, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_sqr(a, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_mod(c, m, c)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_sqrmod() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_exptmod(a, b, m, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_exptmod(a, b, m, c)
Packit 16808d
Packit 16808d
  Compute c = (a ** b) mod m.  Uses a standard square-and-multiply
Packit 16808d
  method with modular reductions at each step. (This is basically the
Packit 16808d
  same code as mw_mp_expt(), except for the addition of the reductions)
Packit 16808d
  
Packit 16808d
  The modular reductions are done using Barrett's algorithm (see
Packit 16808d
  s_mw_mp_reduce() below for details)
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_exptmod(mw_mp_int *a, mw_mp_int *b, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_int   s, x, mu;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_digit d, *db = DIGITS(b);
Packit 16808d
  mw_mp_size  ub = USED(b);
Packit 16808d
  int      dig, bit;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(b) < 0 || mw_mp_cmw_mp_z(m) <= 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  if((res = mw_mp_init(&s)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init_copy(&x, a)) != MP_OKAY)
Packit 16808d
    goto X;
Packit 16808d
  if((res = mw_mp_mod(&x, m, &x)) != MP_OKAY ||
Packit 16808d
     (res = mw_mp_init(&mu)) != MP_OKAY)
Packit 16808d
    goto MU;
Packit 16808d
Packit 16808d
  mw_mp_set(&s, 1);
Packit 16808d
Packit 16808d
  /* mu = b^2k / m */
Packit 16808d
  s_mw_mp_add_d(&mu, 1); 
Packit 16808d
  s_mw_mp_lshd(&mu, 2 * USED(m));
Packit 16808d
  if((res = mw_mp_div(&mu, m, &mu, NULL)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
Packit 16808d
  /* Loop over digits of b in ascending order, except highest order */
Packit 16808d
  for(dig = 0; dig < (ub - 1); dig++) {
Packit 16808d
    d = *db++;
Packit 16808d
Packit 16808d
    /* Loop over the bits of the lower-order digits */
Packit 16808d
    for(bit = 0; bit < DIGIT_BIT; bit++) {
Packit 16808d
      if(d & 1) {
Packit 16808d
	if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY)
Packit 16808d
	  goto CLEANUP;
Packit 16808d
	if((res = s_mw_mp_reduce(&s, m, &mu)) != MP_OKAY)
Packit 16808d
	  goto CLEANUP;
Packit 16808d
      }
Packit 16808d
Packit 16808d
      d >>= 1;
Packit 16808d
Packit 16808d
      if((res = s_mw_mp_sqr(&x)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
      if((res = s_mw_mp_reduce(&x, m, &mu)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
    }
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Now do the last digit... */
Packit 16808d
  d = *db;
Packit 16808d
Packit 16808d
  while(d) {
Packit 16808d
    if(d & 1) {
Packit 16808d
      if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
      if((res = s_mw_mp_reduce(&s, m, &mu)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    d >>= 1;
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_sqr(&x)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
    if((res = s_mw_mp_reduce(&x, m, &mu)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_exch(&s, c);
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  mw_mp_clear(&mu);
Packit 16808d
 MU:
Packit 16808d
  mw_mp_clear(&x);
Packit 16808d
 X:
Packit 16808d
  mw_mp_clear(&s);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_exptmod() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_exptmod_d(a, d, m, c) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_exptmod_d(mw_mp_int *a, mw_mp_digit d, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_int   s, x;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_init(&s)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init_copy(&x, a)) != MP_OKAY)
Packit 16808d
    goto X;
Packit 16808d
Packit 16808d
  mw_mp_set(&s, 1);
Packit 16808d
Packit 16808d
  while(d != 0) {
Packit 16808d
    if(d & 1) {
Packit 16808d
      if((res = s_mw_mp_mul(&s, &x)) != MP_OKAY ||
Packit 16808d
	 (res = mw_mp_mod(&s, m, &s)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    d /= 2;
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_sqr(&x)) != MP_OKAY ||
Packit 16808d
       (res = mw_mp_mod(&x, m, &x)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_exch(&s, c);
Packit 16808d
Packit 16808d
CLEANUP:
Packit 16808d
  mw_mp_clear(&x);
Packit 16808d
X:
Packit 16808d
  mw_mp_clear(&s);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_exptmod_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
#endif /* if MP_MODARITH */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ Comparison functions */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_cmw_mp_z(a) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_cmw_mp_z(a)
Packit 16808d
Packit 16808d
  Compare a <=> 0.  Returns <0 if a<0, 0 if a=0, >0 if a>0.
Packit 16808d
 */
Packit 16808d
Packit 16808d
int    mw_mp_cmw_mp_z(mw_mp_int *a)
Packit 16808d
{
Packit 16808d
  if(SIGN(a) == MP_NEG)
Packit 16808d
    return MP_LT;
Packit 16808d
  else if(USED(a) == 1 && DIGIT(a, 0) == 0)
Packit 16808d
    return MP_EQ;
Packit 16808d
  else
Packit 16808d
    return MP_GT;
Packit 16808d
Packit 16808d
} /* end mw_mp_cmw_mp_z() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_cmw_mp_d(a, d) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_cmw_mp_d(a, d)
Packit 16808d
Packit 16808d
  Compare a <=> d.  Returns <0 if a<d, 0 if a=d, >0 if a>d
Packit 16808d
 */
Packit 16808d
Packit 16808d
int    mw_mp_cmw_mp_d(mw_mp_int *a, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  ARGCHK(a != NULL, MP_EQ);
Packit 16808d
Packit 16808d
  if(SIGN(a) == MP_NEG)
Packit 16808d
    return MP_LT;
Packit 16808d
Packit 16808d
  return s_mw_mp_cmw_mp_d(a, d);
Packit 16808d
Packit 16808d
} /* end mw_mp_cmw_mp_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_cmp(a, b) */
Packit 16808d
Packit 16808d
int    mw_mp_cmp(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_EQ);
Packit 16808d
Packit 16808d
  if(SIGN(a) == SIGN(b)) {
Packit 16808d
    int  mag;
Packit 16808d
Packit 16808d
    if((mag = s_mw_mp_cmp(a, b)) == MP_EQ)
Packit 16808d
      return MP_EQ;
Packit 16808d
Packit 16808d
    if(SIGN(a) == MP_ZPOS)
Packit 16808d
      return mag;
Packit 16808d
    else
Packit 16808d
      return -mag;
Packit 16808d
Packit 16808d
  } else if(SIGN(a) == MP_ZPOS) {
Packit 16808d
    return MP_GT;
Packit 16808d
  } else {
Packit 16808d
    return MP_LT;
Packit 16808d
  }
Packit 16808d
Packit 16808d
} /* end mw_mp_cmp() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_cmw_mp_mag(a, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_cmw_mp_mag(a, b)
Packit 16808d
Packit 16808d
  Compares |a| <=> |b|, and returns an appropriate comparison result
Packit 16808d
 */
Packit 16808d
Packit 16808d
int    mw_mp_cmw_mp_mag(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  ARGCHK(a != NULL && b != NULL, MP_EQ);
Packit 16808d
Packit 16808d
  return s_mw_mp_cmp(a, b);
Packit 16808d
Packit 16808d
} /* end mw_mp_cmw_mp_mag() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_cmw_mp_int(a, z) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  This just converts z to an mw_mp_int, and uses the existing comparison
Packit 16808d
  routines.  This is sort of inefficient, but it's not clear to me how
Packit 16808d
  frequently this wil get used anyway.  For small positive constants,
Packit 16808d
  you can always use mw_mp_cmw_mp_d(), and for zero, there is mw_mp_cmw_mp_z().
Packit 16808d
 */
Packit 16808d
int    mw_mp_cmw_mp_int(mw_mp_int *a, long z)
Packit 16808d
{
Packit 16808d
  mw_mp_int  tmp;
Packit 16808d
  int     out;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL, MP_EQ);
Packit 16808d
  
Packit 16808d
  mw_mp_init(&tmp); mw_mp_set_int(&tmp, z);
Packit 16808d
  out = mw_mp_cmp(a, &tmp);
Packit 16808d
  mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
  return out;
Packit 16808d
Packit 16808d
} /* end mw_mp_cmw_mp_int() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_isodd(a) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_isodd(a)
Packit 16808d
Packit 16808d
  Returns a true (non-zero) value if a is odd, false (zero) otherwise.
Packit 16808d
 */
Packit 16808d
int    mw_mp_isodd(mw_mp_int *a)
Packit 16808d
{
Packit 16808d
  ARGCHK(a != NULL, 0);
Packit 16808d
Packit 16808d
  return (DIGIT(a, 0) & 1);
Packit 16808d
Packit 16808d
} /* end mw_mp_isodd() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_iseven(a) */
Packit 16808d
Packit 16808d
int    mw_mp_iseven(mw_mp_int *a)
Packit 16808d
{
Packit 16808d
  return !mw_mp_isodd(a);
Packit 16808d
Packit 16808d
} /* end mw_mp_iseven() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ Number theoretic functions */
Packit 16808d
Packit 16808d
#if MP_NUMTH
Packit 16808d
/* {{{ mw_mp_gcd(a, b, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Like the old mw_mp_gcd() function, except computes the GCD using the
Packit 16808d
  binary algorithm due to Josef Stein in 1961 (via Knuth).
Packit 16808d
 */
Packit 16808d
mw_mp_err mw_mp_gcd(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_int   u, v, t;
Packit 16808d
  mw_mp_size  k = 0;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(a) == MP_EQ && mw_mp_cmw_mp_z(b) == MP_EQ)
Packit 16808d
      return MP_RANGE;
Packit 16808d
  if(mw_mp_cmw_mp_z(a) == MP_EQ) {
Packit 16808d
    if((res = mw_mp_copy(b, c)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    SIGN(c) = MP_ZPOS; return MP_OKAY;
Packit 16808d
  } else if(mw_mp_cmw_mp_z(b) == MP_EQ) {
Packit 16808d
    if((res = mw_mp_copy(a, c)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    SIGN(c) = MP_ZPOS; return MP_OKAY;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if((res = mw_mp_init(&t)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init_copy(&u, a)) != MP_OKAY)
Packit 16808d
    goto U;
Packit 16808d
  if((res = mw_mp_init_copy(&v, b)) != MP_OKAY)
Packit 16808d
    goto V;
Packit 16808d
Packit 16808d
  SIGN(&u) = MP_ZPOS;
Packit 16808d
  SIGN(&v) = MP_ZPOS;
Packit 16808d
Packit 16808d
  /* Divide out common factors of 2 until at least 1 of a, b is even */
Packit 16808d
  while(mw_mp_iseven(&u) && mw_mp_iseven(&v)) {
Packit 16808d
    s_mw_mp_div_2(&u);
Packit 16808d
    s_mw_mp_div_2(&v);
Packit 16808d
    ++k;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Initialize t */
Packit 16808d
  if(mw_mp_isodd(&u)) {
Packit 16808d
    if((res = mw_mp_copy(&v, &t)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
    
Packit 16808d
    /* t = -v */
Packit 16808d
    if(SIGN(&v) == MP_ZPOS)
Packit 16808d
      SIGN(&t) = MP_NEG;
Packit 16808d
    else
Packit 16808d
      SIGN(&t) = MP_ZPOS;
Packit 16808d
    
Packit 16808d
  } else {
Packit 16808d
    if((res = mw_mp_copy(&u, &t)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
Packit 16808d
  }
Packit 16808d
Packit 16808d
  for(;;) {
Packit 16808d
    while(mw_mp_iseven(&t)) {
Packit 16808d
      s_mw_mp_div_2(&t);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    if(mw_mp_cmw_mp_z(&t) == MP_GT) {
Packit 16808d
      if((res = mw_mp_copy(&t, &u)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
Packit 16808d
    } else {
Packit 16808d
      if((res = mw_mp_copy(&t, &v)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
Packit 16808d
      /* v = -t */
Packit 16808d
      if(SIGN(&t) == MP_ZPOS)
Packit 16808d
	SIGN(&v) = MP_NEG;
Packit 16808d
      else
Packit 16808d
	SIGN(&v) = MP_ZPOS;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    if((res = mw_mp_sub(&u, &v, &t)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
Packit 16808d
    if(s_mw_mp_cmw_mp_d(&t, 0) == MP_EQ)
Packit 16808d
      break;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_2expt(&v, k);       /* v = 2^k   */
Packit 16808d
  res = mw_mp_mul(&u, &v, c); /* c = u * v */
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  mw_mp_clear(&v);
Packit 16808d
 V:
Packit 16808d
  mw_mp_clear(&u);
Packit 16808d
 U:
Packit 16808d
  mw_mp_clear(&t);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_bgcd() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_lcm(a, b, c) */
Packit 16808d
Packit 16808d
/* We compute the least common multiple using the rule:
Packit 16808d
Packit 16808d
   ab = [a, b](a, b)
Packit 16808d
Packit 16808d
   ... by computing the product, and dividing out the gcd.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_lcm(mw_mp_int *a, mw_mp_int *b, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_int  gcd, prod;
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  /* Set up temporaries */
Packit 16808d
  if((res = mw_mp_init(&gcd)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init(&prod)) != MP_OKAY)
Packit 16808d
    goto GCD;
Packit 16808d
Packit 16808d
  if((res = mw_mp_mul(a, b, &prod)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
  if((res = mw_mp_gcd(a, b, &gcd)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
Packit 16808d
  res = mw_mp_div(&prod, &gcd, c, NULL);
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  mw_mp_clear(&prod);
Packit 16808d
 GCD:
Packit 16808d
  mw_mp_clear(&gcd;;
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_lcm() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_xgcd(a, b, g, x, y) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_xgcd(a, b, g, x, y)
Packit 16808d
Packit 16808d
  Compute g = (a, b) and values x and y satisfying Bezout's identity
Packit 16808d
  (that is, ax + by = g).  This uses the extended binary GCD algorithm
Packit 16808d
  based on the Stein algorithm used for mw_mp_gcd()
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_xgcd(mw_mp_int *a, mw_mp_int *b, mw_mp_int *g, mw_mp_int *x, mw_mp_int *y)
Packit 16808d
{
Packit 16808d
  mw_mp_int   gx, xc, yc, u, v, A, B, C, D;
Packit 16808d
  mw_mp_int  *clean[9];
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  int      last = -1;
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(b) == 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /* Initialize all these variables we need */
Packit 16808d
  if((res = mw_mp_init(&u)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &u;
Packit 16808d
  if((res = mw_mp_init(&v)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &v;
Packit 16808d
  if((res = mw_mp_init(&gx)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &gx;
Packit 16808d
  if((res = mw_mp_init(&A)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &A;
Packit 16808d
  if((res = mw_mp_init(&B)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &B;
Packit 16808d
  if((res = mw_mp_init(&C)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &C;
Packit 16808d
  if((res = mw_mp_init(&D)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &D;
Packit 16808d
  if((res = mw_mp_init_copy(&xc, a)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &xc;
Packit 16808d
  mw_mp_abs(&xc, &xc);
Packit 16808d
  if((res = mw_mp_init_copy(&yc, b)) != MP_OKAY) goto CLEANUP;
Packit 16808d
  clean[++last] = &yc;
Packit 16808d
  mw_mp_abs(&yc, &yc);
Packit 16808d
Packit 16808d
  mw_mp_set(&gx, 1);
Packit 16808d
Packit 16808d
  /* Divide by two until at least one of them is even */
Packit 16808d
  while(mw_mp_iseven(&xc) && mw_mp_iseven(&yc)) {
Packit 16808d
    s_mw_mp_div_2(&xc);
Packit 16808d
    s_mw_mp_div_2(&yc);
Packit 16808d
    if((res = s_mw_mp_mul_2(&gx)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  mw_mp_copy(&xc, &u);
Packit 16808d
  mw_mp_copy(&yc, &v);
Packit 16808d
  mw_mp_set(&A, 1); mw_mp_set(&D, 1);
Packit 16808d
Packit 16808d
  /* Loop through binary GCD algorithm */
Packit 16808d
  for(;;) {
Packit 16808d
    while(mw_mp_iseven(&u)) {
Packit 16808d
      s_mw_mp_div_2(&u);
Packit 16808d
Packit 16808d
      if(mw_mp_iseven(&A) && mw_mp_iseven(&B)) {
Packit 16808d
	s_mw_mp_div_2(&A); s_mw_mp_div_2(&B);
Packit 16808d
      } else {
Packit 16808d
	if((res = mw_mp_add(&A, &yc, &A)) != MP_OKAY) goto CLEANUP;
Packit 16808d
	s_mw_mp_div_2(&A);
Packit 16808d
	if((res = mw_mp_sub(&B, &xc, &B)) != MP_OKAY) goto CLEANUP;
Packit 16808d
	s_mw_mp_div_2(&B);
Packit 16808d
      }
Packit 16808d
    }
Packit 16808d
Packit 16808d
    while(mw_mp_iseven(&v)) {
Packit 16808d
      s_mw_mp_div_2(&v);
Packit 16808d
Packit 16808d
      if(mw_mp_iseven(&C) && mw_mp_iseven(&D)) {
Packit 16808d
	s_mw_mp_div_2(&C); s_mw_mp_div_2(&D);
Packit 16808d
      } else {
Packit 16808d
	if((res = mw_mp_add(&C, &yc, &C)) != MP_OKAY) goto CLEANUP;
Packit 16808d
	s_mw_mp_div_2(&C);
Packit 16808d
	if((res = mw_mp_sub(&D, &xc, &D)) != MP_OKAY) goto CLEANUP;
Packit 16808d
	s_mw_mp_div_2(&D);
Packit 16808d
      }
Packit 16808d
    }
Packit 16808d
Packit 16808d
    if(mw_mp_cmp(&u, &v) >= 0) {
Packit 16808d
      if((res = mw_mp_sub(&u, &v, &u)) != MP_OKAY) goto CLEANUP;
Packit 16808d
      if((res = mw_mp_sub(&A, &C, &A)) != MP_OKAY) goto CLEANUP;
Packit 16808d
      if((res = mw_mp_sub(&B, &D, &B)) != MP_OKAY) goto CLEANUP;
Packit 16808d
Packit 16808d
    } else {
Packit 16808d
      if((res = mw_mp_sub(&v, &u, &v)) != MP_OKAY) goto CLEANUP;
Packit 16808d
      if((res = mw_mp_sub(&C, &A, &C)) != MP_OKAY) goto CLEANUP;
Packit 16808d
      if((res = mw_mp_sub(&D, &B, &D)) != MP_OKAY) goto CLEANUP;
Packit 16808d
Packit 16808d
    }
Packit 16808d
Packit 16808d
    /* If we're done, copy results to output */
Packit 16808d
    if(mw_mp_cmw_mp_z(&u) == 0) {
Packit 16808d
      if(x)
Packit 16808d
	if((res = mw_mp_copy(&C, x)) != MP_OKAY) goto CLEANUP;
Packit 16808d
Packit 16808d
      if(y)
Packit 16808d
	if((res = mw_mp_copy(&D, y)) != MP_OKAY) goto CLEANUP;
Packit 16808d
      
Packit 16808d
      if(g)
Packit 16808d
	if((res = mw_mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP;
Packit 16808d
Packit 16808d
      break;
Packit 16808d
    }
Packit 16808d
  }
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  while(last >= 0)
Packit 16808d
    mw_mp_clear(clean[last--]);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_xgcd() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_invmod(a, m, c) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_invmod(a, m, c)
Packit 16808d
Packit 16808d
  Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
Packit 16808d
  This is equivalent to the question of whether (a, m) = 1.  If not,
Packit 16808d
  MP_UNDEF is returned, and there is no inverse.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_invmod(mw_mp_int *a, mw_mp_int *m, mw_mp_int *c)
Packit 16808d
{
Packit 16808d
  mw_mp_int  g, x;
Packit 16808d
  mw_mp_sign sa;
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(a && m && c, MP_BADARG);
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(a) == 0 || mw_mp_cmw_mp_z(m) == 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  sa = SIGN(a);
Packit 16808d
Packit 16808d
  if((res = mw_mp_init(&g)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  if((res = mw_mp_init(&x)) != MP_OKAY)
Packit 16808d
    goto X;
Packit 16808d
Packit 16808d
  if((res = mw_mp_xgcd(a, m, &g, &x, NULL)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_d(&g, 1) != MP_EQ) {
Packit 16808d
    res = MP_UNDEF;
Packit 16808d
    goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  res = mw_mp_mod(&x, m, c);
Packit 16808d
  SIGN(c) = sa;
Packit 16808d
Packit 16808d
CLEANUP:
Packit 16808d
  mw_mp_clear(&x);
Packit 16808d
X:
Packit 16808d
  mw_mp_clear(&g);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_invmod() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
#endif /* if MP_NUMTH */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ mw_mp_print(mp, ofp) */
Packit 16808d
Packit 16808d
#if MP_IOFUNC
Packit 16808d
/*
Packit 16808d
  mw_mp_print(mp, ofp)
Packit 16808d
Packit 16808d
  Print a textual representation of the given mw_mp_int on the output
Packit 16808d
  stream 'ofp'.  Output is generated using the internal radix.
Packit 16808d
 */
Packit 16808d
Packit 16808d
void   mw_mp_print(mw_mp_int *mp, FILE *ofp)
Packit 16808d
{
Packit 16808d
  int   ix;
Packit 16808d
Packit 16808d
  if(mp == NULL || ofp == NULL)
Packit 16808d
    return;
Packit 16808d
Packit 16808d
  fputc((SIGN(mp) == MP_NEG) ? '-' : '+', ofp);
Packit 16808d
Packit 16808d
  for(ix = USED(mp) - 1; ix >= 0; ix--) {
Packit 16808d
    fprintf(ofp, DIGIT_FMT, DIGIT(mp, ix));
Packit 16808d
  }
Packit 16808d
Packit 16808d
} /* end mw_mp_print() */
Packit 16808d
Packit 16808d
#endif /* if MP_IOFUNC */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* {{{ More I/O Functions */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_read_signed_bin(mp, str, len) */
Packit 16808d
Packit 16808d
/* 
Packit 16808d
   mw_mp_read_signed_bin(mp, str, len)
Packit 16808d
Packit 16808d
   Read in a raw value (base 256) into the given mw_mp_int
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err  mw_mp_read_signed_bin(mw_mp_int *mp, unsigned char *str, int len)
Packit 16808d
{
Packit 16808d
  mw_mp_err         res;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
Packit 16808d
Packit 16808d
  if((res = mw_mp_read_unsigned_bin(mp, str + 1, len - 1)) == MP_OKAY) {
Packit 16808d
    /* Get sign from first byte */
Packit 16808d
    if(str[0])
Packit 16808d
      SIGN(mp) = MP_NEG;
Packit 16808d
    else
Packit 16808d
      SIGN(mp) = MP_ZPOS;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end mw_mp_read_signed_bin() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_signed_bin_size(mp) */
Packit 16808d
Packit 16808d
int    mw_mp_signed_bin_size(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  ARGCHK(mp != NULL, 0);
Packit 16808d
Packit 16808d
  return mw_mp_unsigned_bin_size(mp) + 1;
Packit 16808d
Packit 16808d
} /* end mw_mp_signed_bin_size() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_to_signed_bin(mp, str) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_to_signed_bin(mw_mp_int *mp, unsigned char *str)
Packit 16808d
{
Packit 16808d
  ARGCHK(mp != NULL && str != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  /* Caller responsible for allocating enough memory (use mw_mp_raw_size(mp)) */
Packit 16808d
  str[0] = (char)SIGN(mp);
Packit 16808d
Packit 16808d
  return mw_mp_to_unsigned_bin(mp, str + 1);
Packit 16808d
Packit 16808d
} /* end mw_mp_to_signed_bin() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_read_unsigned_bin(mp, str, len) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_read_unsigned_bin(mp, str, len)
Packit 16808d
Packit 16808d
  Read in an unsigned value (base 256) into the given mw_mp_int
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err  mw_mp_read_unsigned_bin(mw_mp_int *mp, unsigned char *str, int len)
Packit 16808d
{
Packit 16808d
  int     ix;
Packit 16808d
  mw_mp_err  res;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
Packit 16808d
Packit 16808d
  mw_mp_zero(mp);
Packit 16808d
Packit 16808d
  for(ix = 0; ix < len; ix++) {
Packit 16808d
    if((res = s_mw_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    if((res = mw_mp_add_d(mp, str[ix], mp)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
  }
Packit 16808d
  
Packit 16808d
  return MP_OKAY;
Packit 16808d
  
Packit 16808d
} /* end mw_mp_read_unsigned_bin() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_unsigned_bin_size(mp) */
Packit 16808d
Packit 16808d
int     mw_mp_unsigned_bin_size(mw_mp_int *mp) 
Packit 16808d
{
Packit 16808d
  mw_mp_digit   topdig;
Packit 16808d
  int        count;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL, 0);
Packit 16808d
Packit 16808d
  /* Special case for the value zero */
Packit 16808d
  if(USED(mp) == 1 && DIGIT(mp, 0) == 0)
Packit 16808d
    return 1;
Packit 16808d
Packit 16808d
  count = (USED(mp) - 1) * sizeof(mw_mp_digit);
Packit 16808d
  topdig = DIGIT(mp, USED(mp) - 1);
Packit 16808d
Packit 16808d
  while(topdig != 0) {
Packit 16808d
    ++count;
Packit 16808d
    topdig >>= CHAR_BIT;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return count;
Packit 16808d
Packit 16808d
} /* end mw_mp_unsigned_bin_size() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_to_unsigned_bin(mp, str) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_to_unsigned_bin(mw_mp_int *mp, unsigned char *str)
Packit 16808d
{
Packit 16808d
  mw_mp_digit      *dp, *end, d;
Packit 16808d
  unsigned char *spos;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL && str != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  dp = DIGITS(mp);
Packit 16808d
  end = dp + USED(mp) - 1;
Packit 16808d
  spos = str;
Packit 16808d
Packit 16808d
  /* Special case for zero, quick test */
Packit 16808d
  if(dp == end && *dp == 0) {
Packit 16808d
    *str = '\0';
Packit 16808d
    return MP_OKAY;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Generate digits in reverse order */
Packit 16808d
  while(dp < end) {
Packit 16808d
    int      ix;
Packit 16808d
Packit 16808d
    d = *dp;
Packit 16808d
    for(ix = 0; ix < sizeof(mw_mp_digit); ++ix) {
Packit 16808d
      *spos = d & UCHAR_MAX;
Packit 16808d
      d >>= CHAR_BIT;
Packit 16808d
      ++spos;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    ++dp;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Now handle last digit specially, high order zeroes are not written */
Packit 16808d
  d = *end;
Packit 16808d
  while(d != 0) {
Packit 16808d
    *spos = d & UCHAR_MAX;
Packit 16808d
    d >>= CHAR_BIT;
Packit 16808d
    ++spos;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Reverse everything to get digits in the correct order */
Packit 16808d
  while(--spos > str) {
Packit 16808d
    unsigned char t = *str;
Packit 16808d
    *str = *spos;
Packit 16808d
    *spos = t;
Packit 16808d
Packit 16808d
    ++str;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_to_unsigned_bin() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_count_bits(mp) */
Packit 16808d
Packit 16808d
int    mw_mp_count_bits(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  int      len;
Packit 16808d
  mw_mp_digit d;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL, MP_BADARG);
Packit 16808d
Packit 16808d
  len = DIGIT_BIT * (USED(mp) - 1);
Packit 16808d
  d = DIGIT(mp, USED(mp) - 1);
Packit 16808d
Packit 16808d
  while(d != 0) {
Packit 16808d
    ++len;
Packit 16808d
    d >>= 1;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return len;
Packit 16808d
  
Packit 16808d
} /* end mw_mp_count_bits() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_read_radix(mp, str, radix) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_read_radix(mp, str, radix)
Packit 16808d
Packit 16808d
  Read an integer from the given string, and set mp to the resulting
Packit 16808d
  value.  The input is presumed to be in base 10.  Leading non-digit
Packit 16808d
  characters are ignored, and the function reads until a non-digit
Packit 16808d
  character or the end of the string.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err  mw_mp_read_radix(mw_mp_int *mp, unsigned char *str, int radix)
Packit 16808d
{
Packit 16808d
  int     ix = 0, val = 0;
Packit 16808d
  mw_mp_err  res;
Packit 16808d
  mw_mp_sign sig = MP_ZPOS;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX, 
Packit 16808d
	 MP_BADARG);
Packit 16808d
Packit 16808d
  mw_mp_zero(mp);
Packit 16808d
Packit 16808d
  /* Skip leading non-digit characters until a digit or '-' or '+' */
Packit 16808d
  while(str[ix] && 
Packit 16808d
	(s_mw_mp_tovalue(str[ix], radix) < 0) && 
Packit 16808d
	str[ix] != '-' &&
Packit 16808d
	str[ix] != '+') {
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(str[ix] == '-') {
Packit 16808d
    sig = MP_NEG;
Packit 16808d
    ++ix;
Packit 16808d
  } else if(str[ix] == '+') {
Packit 16808d
    sig = MP_ZPOS; /* this is the default anyway... */
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  while((val = s_mw_mp_tovalue(str[ix], radix)) >= 0) {
Packit 16808d
    if((res = s_mw_mp_mul_d(mp, radix)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    if((res = s_mw_mp_add_d(mp, val)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(s_mw_mp_cmw_mp_d(mp, 0) == MP_EQ)
Packit 16808d
    SIGN(mp) = MP_ZPOS;
Packit 16808d
  else
Packit 16808d
    SIGN(mp) = sig;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_read_radix() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_radix_size(mp, radix) */
Packit 16808d
Packit 16808d
int    mw_mp_radix_size(mw_mp_int *mp, int radix)
Packit 16808d
{
Packit 16808d
  int  len;
Packit 16808d
  ARGCHK(mp != NULL, 0);
Packit 16808d
Packit 16808d
  len = s_mw_mp_outlen(mw_mp_count_bits(mp), radix) + 1; /* for NUL terminator */
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(mp) < 0)
Packit 16808d
    ++len; /* for sign */
Packit 16808d
Packit 16808d
  return len;
Packit 16808d
Packit 16808d
} /* end mw_mp_radix_size() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_value_radix_size(num, qty, radix) */
Packit 16808d
Packit 16808d
/* num = number of digits
Packit 16808d
   qty = number of bits per digit
Packit 16808d
   radix = target base
Packit 16808d
   
Packit 16808d
   Return the number of digits in the specified radix that would be
Packit 16808d
   needed to express 'num' digits of 'qty' bits each.
Packit 16808d
 */
Packit 16808d
int    mw_mp_value_radix_size(int num, int qty, int radix)
Packit 16808d
{
Packit 16808d
  ARGCHK(num >= 0 && qty > 0 && radix >= 2 && radix <= MAX_RADIX, 0);
Packit 16808d
Packit 16808d
  return s_mw_mp_outlen(num * qty, radix);
Packit 16808d
Packit 16808d
} /* end mw_mp_value_radix_size() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_toradix(mp, str, radix) */
Packit 16808d
Packit 16808d
mw_mp_err mw_mp_toradix(mw_mp_int *mp, unsigned char *str, int radix)
Packit 16808d
{
Packit 16808d
  int  ix, pos = 0;
Packit 16808d
Packit 16808d
  ARGCHK(mp != NULL && str != NULL, MP_BADARG);
Packit 16808d
  ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(mp) == MP_EQ) {
Packit 16808d
    str[0] = '0';
Packit 16808d
    str[1] = '\0';
Packit 16808d
  } else {
Packit 16808d
    mw_mp_err   res;
Packit 16808d
    mw_mp_int   tmp;
Packit 16808d
    mw_mp_sign  sgn;
Packit 16808d
    mw_mp_digit rem, rdx = (mw_mp_digit)radix;
Packit 16808d
    char     ch;
Packit 16808d
Packit 16808d
    if((res = mw_mp_init_copy(&tmp, mp)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    /* Save sign for later, and take absolute value */
Packit 16808d
    sgn = SIGN(&tmp); SIGN(&tmp) = MP_ZPOS;
Packit 16808d
Packit 16808d
    /* Generate output digits in reverse order      */
Packit 16808d
    while(mw_mp_cmw_mp_z(&tmp) != 0) {
Packit 16808d
      if((res = s_mw_mp_div_d(&tmp, rdx, &rem)) != MP_OKAY) {
Packit 16808d
	mw_mp_clear(&tmp);
Packit 16808d
	return res;
Packit 16808d
      }
Packit 16808d
Packit 16808d
      /* Generate digits, use capital letters */
Packit 16808d
      ch = s_mw_mp_todigit(rem, radix, 0);
Packit 16808d
Packit 16808d
      str[pos++] = ch;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    /* Add - sign if original value was negative */
Packit 16808d
    if(sgn == MP_NEG)
Packit 16808d
      str[pos++] = '-';
Packit 16808d
Packit 16808d
    /* Add trailing NUL to end the string        */
Packit 16808d
    str[pos--] = '\0';
Packit 16808d
Packit 16808d
    /* Reverse the digits and sign indicator     */
Packit 16808d
    ix = 0;
Packit 16808d
    while(ix < pos) {
Packit 16808d
      char tmp = str[ix];
Packit 16808d
Packit 16808d
      str[ix] = str[pos];
Packit 16808d
      str[pos] = tmp;
Packit 16808d
      ++ix;
Packit 16808d
      --pos;
Packit 16808d
    }
Packit 16808d
    
Packit 16808d
    mw_mp_clear(&tmp);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end mw_mp_toradix() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_char2value(ch, r) */
Packit 16808d
Packit 16808d
int    mw_mp_char2value(char ch, int r)
Packit 16808d
{
Packit 16808d
  return s_mw_mp_tovalue(ch, r);
Packit 16808d
Packit 16808d
} /* end mw_mp_tovalue() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ mw_mp_strerror(ec) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  mw_mp_strerror(ec)
Packit 16808d
Packit 16808d
  Return a string describing the meaning of error code 'ec'.  The
Packit 16808d
  string returned is allocated in static memory, so the caller should
Packit 16808d
  not attempt to modify or free the memory associated with this
Packit 16808d
  string.
Packit 16808d
 */
Packit 16808d
const char  *mw_mp_strerror(mw_mp_err ec)
Packit 16808d
{
Packit 16808d
  int   aec = (ec < 0) ? -ec : ec;
Packit 16808d
Packit 16808d
  /* Code values are negative, so the senses of these comparisons
Packit 16808d
     are accurate */
Packit 16808d
  if(ec < MP_LAST_CODE || ec > MP_OKAY) {
Packit 16808d
    return mw_mp_err_string[0];  /* unknown error code */
Packit 16808d
  } else {
Packit 16808d
    return mw_mp_err_string[aec + 1];
Packit 16808d
  }
Packit 16808d
Packit 16808d
} /* end mw_mp_strerror() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*========================================================================*/
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* Static function definitions (internal use only)                        */
Packit 16808d
Packit 16808d
/* {{{ Memory management */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_grow(mp, min) */
Packit 16808d
Packit 16808d
/* Make sure there are at least 'min' digits allocated to mp              */
Packit 16808d
mw_mp_err   s_mw_mp_grow(mw_mp_int *mp, mw_mp_size min)
Packit 16808d
{
Packit 16808d
  if(min > ALLOC(mp)) {
Packit 16808d
    mw_mp_digit   *tmp;
Packit 16808d
Packit 16808d
    /* Set min to next nearest default precision block size */
Packit 16808d
    min = ((min + (s_mw_mp_defprec - 1)) / s_mw_mp_defprec) * s_mw_mp_defprec;
Packit 16808d
Packit 16808d
    if((tmp = s_mw_mp_alloc(min, sizeof(mw_mp_digit))) == NULL)
Packit 16808d
      return MP_MEM;
Packit 16808d
Packit 16808d
    s_mw_mp_copy(DIGITS(mp), tmp, USED(mp));
Packit 16808d
Packit 16808d
#if MP_CRYPTO
Packit 16808d
    s_mw_mp_setz(DIGITS(mp), ALLOC(mp));
Packit 16808d
#endif
Packit 16808d
    s_mw_mp_free(DIGITS(mp));
Packit 16808d
    DIGITS(mp) = tmp;
Packit 16808d
    ALLOC(mp) = min;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_grow() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_pad(mp, min) */
Packit 16808d
Packit 16808d
/* Make sure the used size of mp is at least 'min', growing if needed     */
Packit 16808d
mw_mp_err   s_mw_mp_pad(mw_mp_int *mp, mw_mp_size min)
Packit 16808d
{
Packit 16808d
  if(min > USED(mp)) {
Packit 16808d
    mw_mp_err  res;
Packit 16808d
Packit 16808d
    /* Make sure there is room to increase precision  */
Packit 16808d
    if(min > ALLOC(mp) && (res = s_mw_mp_grow(mp, min)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    /* Increase precision; should already be 0-filled */
Packit 16808d
    USED(mp) = min;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_pad() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_setz(dp, count) */
Packit 16808d
Packit 16808d
#if MP_MACRO == 0
Packit 16808d
/* Set 'count' digits pointed to by dp to be zeroes                       */
Packit 16808d
void s_mw_mp_setz(mw_mp_digit *dp, mw_mp_size count)
Packit 16808d
{
Packit 16808d
#if MP_MEMSET == 0
Packit 16808d
  int  ix;
Packit 16808d
Packit 16808d
  for(ix = 0; ix < count; ix++)
Packit 16808d
    dp[ix] = 0;
Packit 16808d
#else
Packit 16808d
  memset(dp, 0, count * sizeof(mw_mp_digit));
Packit 16808d
#endif
Packit 16808d
Packit 16808d
} /* end s_mw_mp_setz() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_copy(sp, dp, count) */
Packit 16808d
Packit 16808d
#if MP_MACRO == 0
Packit 16808d
/* Copy 'count' digits from sp to dp                                      */
Packit 16808d
void s_mw_mp_copy(mw_mp_digit *sp, mw_mp_digit *dp, mw_mp_size count)
Packit 16808d
{
Packit 16808d
#if MP_MEMCPY == 0
Packit 16808d
  int  ix;
Packit 16808d
Packit 16808d
  for(ix = 0; ix < count; ix++)
Packit 16808d
    dp[ix] = sp[ix];
Packit 16808d
#else
Packit 16808d
  memcpy(dp, sp, count * sizeof(mw_mp_digit));
Packit 16808d
#endif
Packit 16808d
Packit 16808d
} /* end s_mw_mp_copy() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_alloc(nb, ni) */
Packit 16808d
Packit 16808d
#if MP_MACRO == 0
Packit 16808d
/* Allocate ni records of nb bytes each, and return a pointer to that     */
Packit 16808d
void    *s_mw_mp_alloc(size_t nb, size_t ni)
Packit 16808d
{
Packit 16808d
  return calloc(nb, ni);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_alloc() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_free(ptr) */
Packit 16808d
Packit 16808d
#if MP_MACRO == 0
Packit 16808d
/* Free the memory pointed to by ptr                                      */
Packit 16808d
void     s_mw_mp_free(void *ptr)
Packit 16808d
{
Packit 16808d
  if(ptr)
Packit 16808d
    free(ptr);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_free() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_clamp(mp) */
Packit 16808d
Packit 16808d
/* Remove leading zeroes from the given value                             */
Packit 16808d
void     s_mw_mp_clamp(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  mw_mp_size   du = USED(mp);
Packit 16808d
  mw_mp_digit *zp = DIGITS(mp) + du - 1;
Packit 16808d
Packit 16808d
  while(du > 1 && !*zp--)
Packit 16808d
    --du;
Packit 16808d
Packit 16808d
  if(du == 1 && *zp == 0)
Packit 16808d
    SIGN(mp) = MP_ZPOS;
Packit 16808d
Packit 16808d
  USED(mp) = du;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_clamp() */
Packit 16808d
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_exch(a, b) */
Packit 16808d
Packit 16808d
/* Exchange the data for a and b; (b, a) = (a, b)                         */
Packit 16808d
void     s_mw_mp_exch(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_int   tmp;
Packit 16808d
Packit 16808d
  tmp = *a;
Packit 16808d
  *a = *b;
Packit 16808d
  *b = tmp;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_exch() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Arithmetic helpers */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_lshd(mp, p) */
Packit 16808d
Packit 16808d
/* 
Packit 16808d
   Shift mp leftward by p digits, growing if needed, and zero-filling
Packit 16808d
   the in-shifted digits at the right end.  This is a convenient
Packit 16808d
   alternative to multiplication by powers of the radix
Packit 16808d
 */   
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_lshd(mw_mp_int *mp, mw_mp_size p)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_size  pos;
Packit 16808d
  mw_mp_digit *dp;
Packit 16808d
  int     ix;
Packit 16808d
Packit 16808d
  if(p == 0)
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
  if((res = s_mw_mp_pad(mp, USED(mp) + p)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  pos = USED(mp) - 1;
Packit 16808d
  dp = DIGITS(mp);
Packit 16808d
Packit 16808d
  /* Shift all the significant figures over as needed */
Packit 16808d
  for(ix = pos - p; ix >= 0; ix--) 
Packit 16808d
    dp[ix + p] = dp[ix];
Packit 16808d
Packit 16808d
  /* Fill the bottom digits with zeroes */
Packit 16808d
  for(ix = 0; ix < p; ix++)
Packit 16808d
    dp[ix] = 0;
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_lshd() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_rshd(mp, p) */
Packit 16808d
Packit 16808d
/* 
Packit 16808d
   Shift mp rightward by p digits.  Maintains the invariant that
Packit 16808d
   digits above the precision are all zero.  Digits shifted off the
Packit 16808d
   end are lost.  Cannot fail.
Packit 16808d
 */
Packit 16808d
Packit 16808d
void     s_mw_mp_rshd(mw_mp_int *mp, mw_mp_size p)
Packit 16808d
{
Packit 16808d
  mw_mp_size  ix;
Packit 16808d
  mw_mp_digit *dp;
Packit 16808d
Packit 16808d
  if(p == 0)
Packit 16808d
    return;
Packit 16808d
Packit 16808d
  /* Shortcut when all digits are to be shifted off */
Packit 16808d
  if(p >= USED(mp)) {
Packit 16808d
    s_mw_mp_setz(DIGITS(mp), ALLOC(mp));
Packit 16808d
    USED(mp) = 1;
Packit 16808d
    SIGN(mp) = MP_ZPOS;
Packit 16808d
    return;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Shift all the significant figures over as needed */
Packit 16808d
  dp = DIGITS(mp);
Packit 16808d
  for(ix = p; ix < USED(mp); ix++)
Packit 16808d
    dp[ix - p] = dp[ix];
Packit 16808d
Packit 16808d
  /* Fill the top digits with zeroes */
Packit 16808d
  ix -= p;
Packit 16808d
  while(ix < USED(mp))
Packit 16808d
    dp[ix++] = 0;
Packit 16808d
Packit 16808d
  /* Strip off any leading zeroes    */
Packit 16808d
  s_mw_mp_clamp(mp);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_rshd() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_div_2(mp) */
Packit 16808d
Packit 16808d
/* Divide by two -- take advantage of radix properties to do it fast      */
Packit 16808d
void     s_mw_mp_div_2(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  s_mw_mp_div_2d(mp, 1);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_div_2() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_mul_2(mp) */
Packit 16808d
Packit 16808d
mw_mp_err s_mw_mp_mul_2(mw_mp_int *mp)
Packit 16808d
{
Packit 16808d
  int      ix;
Packit 16808d
  mw_mp_digit kin = 0, kout, *dp = DIGITS(mp);
Packit 16808d
  mw_mp_err   res;
Packit 16808d
Packit 16808d
  /* Shift digits leftward by 1 bit */
Packit 16808d
  for(ix = 0; ix < USED(mp); ix++) {
Packit 16808d
    kout = (dp[ix] >> (DIGIT_BIT - 1)) & 1;
Packit 16808d
    dp[ix] = (dp[ix] << 1) | kin;
Packit 16808d
Packit 16808d
    kin = kout;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Deal with rollover from last digit */
Packit 16808d
  if(kin) {
Packit 16808d
    if(ix >= ALLOC(mp)) {
Packit 16808d
      if((res = s_mw_mp_grow(mp, ALLOC(mp) + 1)) != MP_OKAY)
Packit 16808d
	return res;
Packit 16808d
      dp = DIGITS(mp);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    dp[ix] = kin;
Packit 16808d
    USED(mp) += 1;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_mul_2() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_mod_2d(mp, d) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Remainder the integer by 2^d, where d is a number of bits.  This
Packit 16808d
  amounts to a bitwise AND of the value, and does not require the full
Packit 16808d
  division code
Packit 16808d
 */
Packit 16808d
void     s_mw_mp_mod_2d(mw_mp_int *mp, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  unsigned int  ndig = (d / DIGIT_BIT), nbit = (d % DIGIT_BIT);
Packit 16808d
  unsigned int  ix;
Packit 16808d
  mw_mp_digit      dmask, *dp = DIGITS(mp);
Packit 16808d
Packit 16808d
  if(ndig >= USED(mp))
Packit 16808d
    return;
Packit 16808d
Packit 16808d
  /* Flush all the bits above 2^d in its digit */
Packit 16808d
  dmask = (1 << nbit) - 1;
Packit 16808d
  dp[ndig] &= dmask;
Packit 16808d
Packit 16808d
  /* Flush all digits above the one with 2^d in it */
Packit 16808d
  for(ix = ndig + 1; ix < USED(mp); ix++)
Packit 16808d
    dp[ix] = 0;
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(mp);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_mod_2d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_mul_2d(mp, d) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Multiply by the integer 2^d, where d is a number of bits.  This
Packit 16808d
  amounts to a bitwise shift of the value, and does not require the
Packit 16808d
  full multiplication code.
Packit 16808d
 */
Packit 16808d
mw_mp_err    s_mw_mp_mul_2d(mw_mp_int *mp, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_digit save, next, mask, *dp;
Packit 16808d
  mw_mp_size  used;
Packit 16808d
  int      ix;
Packit 16808d
Packit 16808d
  if((res = s_mw_mp_lshd(mp, d / DIGIT_BIT)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  dp = DIGITS(mp); used = USED(mp);
Packit 16808d
  d %= DIGIT_BIT;
Packit 16808d
Packit 16808d
  mask = (1 << d) - 1;
Packit 16808d
Packit 16808d
  /* If the shift requires another digit, make sure we've got one to
Packit 16808d
     work with */
Packit 16808d
  if((dp[used - 1] >> (DIGIT_BIT - d)) & mask) {
Packit 16808d
    if((res = s_mw_mp_grow(mp, used + 1)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    dp = DIGITS(mp);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Do the shifting... */
Packit 16808d
  save = 0;
Packit 16808d
  for(ix = 0; ix < used; ix++) {
Packit 16808d
    next = (dp[ix] >> (DIGIT_BIT - d)) & mask;
Packit 16808d
    dp[ix] = (dp[ix] << d) | save;
Packit 16808d
    save = next;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* If, at this point, we have a nonzero carryout into the next
Packit 16808d
     digit, we'll increase the size by one digit, and store it...
Packit 16808d
   */
Packit 16808d
  if(save) {
Packit 16808d
    dp[used] = save;
Packit 16808d
    USED(mp) += 1;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(mp);
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_mul_2d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_div_2d(mp, d) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Divide the integer by 2^d, where d is a number of bits.  This
Packit 16808d
  amounts to a bitwise shift of the value, and does not require the
Packit 16808d
  full division code (used in Barrett reduction, see below)
Packit 16808d
 */
Packit 16808d
void     s_mw_mp_div_2d(mw_mp_int *mp, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  int       ix;
Packit 16808d
  mw_mp_digit  save, next, mask, *dp = DIGITS(mp);
Packit 16808d
Packit 16808d
  s_mw_mp_rshd(mp, d / DIGIT_BIT);
Packit 16808d
  d %= DIGIT_BIT;
Packit 16808d
Packit 16808d
  mask = (1 << d) - 1;
Packit 16808d
Packit 16808d
  save = 0;
Packit 16808d
  for(ix = USED(mp) - 1; ix >= 0; ix--) {
Packit 16808d
    next = dp[ix] & mask;
Packit 16808d
    dp[ix] = (dp[ix] >> d) | (save << (DIGIT_BIT - d));
Packit 16808d
    save = next;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(mp);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_div_2d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_norm(a, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  s_mw_mp_norm(a, b)
Packit 16808d
Packit 16808d
  Normalize a and b for division, where b is the divisor.  In order
Packit 16808d
  that we might make good guesses for quotient digits, we want the
Packit 16808d
  leading digit of b to be at least half the radix, which we
Packit 16808d
  accomplish by multiplying a and b by a constant.  This constant is
Packit 16808d
  returned (so that it can be divided back out of the remainder at the
Packit 16808d
  end of the division process).
Packit 16808d
Packit 16808d
  We multiply by the smallest power of 2 that gives us a leading digit
Packit 16808d
  at least half the radix.  By choosing a power of 2, we simplify the 
Packit 16808d
  multiplication and division steps to simple shifts.
Packit 16808d
 */
Packit 16808d
mw_mp_digit s_mw_mp_norm(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_digit  t, d = 0;
Packit 16808d
Packit 16808d
  t = DIGIT(b, USED(b) - 1);
Packit 16808d
  while(t < (RADIX / 2)) {
Packit 16808d
    t <<= 1;
Packit 16808d
    ++d;
Packit 16808d
  }
Packit 16808d
    
Packit 16808d
  if(d != 0) {
Packit 16808d
    s_mw_mp_mul_2d(a, d);
Packit 16808d
    s_mw_mp_mul_2d(b, d);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return d;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_norm() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Primitive digit arithmetic */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_add_d(mp, d) */
Packit 16808d
Packit 16808d
/* Add d to |mp| in place                                                 */
Packit 16808d
mw_mp_err   s_mw_mp_add_d(mw_mp_int *mp, mw_mp_digit d)    /* unsigned digit addition */
Packit 16808d
{
Packit 16808d
  mw_mp_word   w, k = 0;
Packit 16808d
  mw_mp_size   ix = 1, used = USED(mp);
Packit 16808d
  mw_mp_digit *dp = DIGITS(mp);
Packit 16808d
Packit 16808d
  w = dp[0] + d;
Packit 16808d
  dp[0] = ACCUM(w);
Packit 16808d
  k = CARRYOUT(w);
Packit 16808d
Packit 16808d
  while(ix < used && k) {
Packit 16808d
    w = dp[ix] + k;
Packit 16808d
    dp[ix] = ACCUM(w);
Packit 16808d
    k = CARRYOUT(w);
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(k != 0) {
Packit 16808d
    mw_mp_err  res;
Packit 16808d
Packit 16808d
    if((res = s_mw_mp_pad(mp, USED(mp) + 1)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    DIGIT(mp, ix) = k;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_add_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_sub_d(mp, d) */
Packit 16808d
Packit 16808d
/* Subtract d from |mp| in place, assumes |mp| > d                        */
Packit 16808d
mw_mp_err   s_mw_mp_sub_d(mw_mp_int *mp, mw_mp_digit d)    /* unsigned digit subtract */
Packit 16808d
{
Packit 16808d
  mw_mp_word   w, b = 0;
Packit 16808d
  mw_mp_size   ix = 1, used = USED(mp);
Packit 16808d
  mw_mp_digit *dp = DIGITS(mp);
Packit 16808d
Packit 16808d
  /* Compute initial subtraction    */
Packit 16808d
  w = (RADIX + dp[0]) - d;
Packit 16808d
  b = CARRYOUT(w) ? 0 : 1;
Packit 16808d
  dp[0] = ACCUM(w);
Packit 16808d
Packit 16808d
  /* Propagate borrows leftward     */
Packit 16808d
  while(b && ix < used) {
Packit 16808d
    w = (RADIX + dp[ix]) - b;
Packit 16808d
    b = CARRYOUT(w) ? 0 : 1;
Packit 16808d
    dp[ix] = ACCUM(w);
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Remove leading zeroes          */
Packit 16808d
  s_mw_mp_clamp(mp);
Packit 16808d
Packit 16808d
  /* If we have a borrow out, it's a violation of the input invariant */
Packit 16808d
  if(b)
Packit 16808d
    return MP_RANGE;
Packit 16808d
  else
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_sub_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_mul_d(a, d) */
Packit 16808d
Packit 16808d
/* Compute a = a * d, single digit multiplication                         */
Packit 16808d
mw_mp_err   s_mw_mp_mul_d(mw_mp_int *a, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  mw_mp_word w, k = 0;
Packit 16808d
  mw_mp_size ix, max;
Packit 16808d
  mw_mp_err  res;
Packit 16808d
  mw_mp_digit *dp = DIGITS(a);
Packit 16808d
Packit 16808d
  /*
Packit 16808d
    Single-digit multiplication will increase the precision of the
Packit 16808d
    output by at most one digit.  However, we can detect when this
Packit 16808d
    will happen -- if the high-order digit of a, times d, gives a
Packit 16808d
    two-digit result, then the precision of the result will increase;
Packit 16808d
    otherwise it won't.  We use this fact to avoid calling s_mw_mp_pad()
Packit 16808d
    unless absolutely necessary.
Packit 16808d
   */
Packit 16808d
  max = USED(a);
Packit 16808d
  w = dp[max - 1] * d;
Packit 16808d
  if(CARRYOUT(w) != 0) {
Packit 16808d
    if((res = s_mw_mp_pad(a, max + 1)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
    dp = DIGITS(a);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  for(ix = 0; ix < max; ix++) {
Packit 16808d
    w = (dp[ix] * d) + k;
Packit 16808d
    dp[ix] = ACCUM(w);
Packit 16808d
    k = CARRYOUT(w);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* If there is a precision increase, take care of it here; the above
Packit 16808d
     test guarantees we have enough storage to do this safely.
Packit 16808d
   */
Packit 16808d
  if(k) {
Packit 16808d
    dp[max] = k; 
Packit 16808d
    USED(a) = max + 1;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(a);
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
  
Packit 16808d
} /* end s_mw_mp_mul_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_div_d(mp, d, r) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  s_mw_mp_div_d(mp, d, r)
Packit 16808d
Packit 16808d
  Compute the quotient mp = mp / d and remainder r = mp mod d, for a
Packit 16808d
  single digit d.  If r is null, the remainder will be discarded.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_div_d(mw_mp_int *mp, mw_mp_digit d, mw_mp_digit *r)
Packit 16808d
{
Packit 16808d
  mw_mp_word   w = 0, t;
Packit 16808d
  mw_mp_int    quot;
Packit 16808d
  mw_mp_err    res;
Packit 16808d
  mw_mp_digit *dp = DIGITS(mp), *qp;
Packit 16808d
  int       ix;
Packit 16808d
Packit 16808d
  if(d == 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /* Make room for the quotient */
Packit 16808d
  if((res = mw_mp_init_size(&quot, USED(mp))) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  USED(&quot) = USED(mp); /* so clamping will work below */
Packit 16808d
  qp = DIGITS(");
Packit 16808d
Packit 16808d
  /* Divide without subtraction */
Packit 16808d
  for(ix = USED(mp) - 1; ix >= 0; ix--) {
Packit 16808d
    w = (w << DIGIT_BIT) | dp[ix];
Packit 16808d
Packit 16808d
    if(w >= d) {
Packit 16808d
      t = w / d;
Packit 16808d
      w = w % d;
Packit 16808d
    } else {
Packit 16808d
      t = 0;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    qp[ix] = t;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Deliver the remainder, if desired */
Packit 16808d
  if(r)
Packit 16808d
    *r = w;
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(");
Packit 16808d
  mw_mp_exch(&quot, mp);
Packit 16808d
  mw_mp_clear(");
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_div_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Primitive full arithmetic */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_add(a, b) */
Packit 16808d
Packit 16808d
/* Compute a = |a| + |b|                                                  */
Packit 16808d
mw_mp_err   s_mw_mp_add(mw_mp_int *a, mw_mp_int *b)        /* magnitude addition      */
Packit 16808d
{
Packit 16808d
  mw_mp_word   w = 0;
Packit 16808d
  mw_mp_digit *pa, *pb;
Packit 16808d
  mw_mp_size   ix, used = USED(b);
Packit 16808d
  mw_mp_err    res;
Packit 16808d
Packit 16808d
  /* Make sure a has enough precision for the output value */
Packit 16808d
  if((used > USED(a)) && (res = s_mw_mp_pad(a, used)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  /*
Packit 16808d
    Add up all digits up to the precision of b.  If b had initially
Packit 16808d
    the same precision as a, or greater, we took care of it by the
Packit 16808d
    padding step above, so there is no problem.  If b had initially
Packit 16808d
    less precision, we'll have to make sure the carry out is duly
Packit 16808d
    propagated upward among the higher-order digits of the sum.
Packit 16808d
   */
Packit 16808d
  pa = DIGITS(a);
Packit 16808d
  pb = DIGITS(b);
Packit 16808d
  for(ix = 0; ix < used; ++ix) {
Packit 16808d
    w += *pa + *pb++;
Packit 16808d
    *pa++ = ACCUM(w);
Packit 16808d
    w = CARRYOUT(w);
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* If we run out of 'b' digits before we're actually done, make
Packit 16808d
     sure the carries get propagated upward...  
Packit 16808d
   */
Packit 16808d
  used = USED(a);
Packit 16808d
  while(w && ix < used) {
Packit 16808d
    w += *pa;
Packit 16808d
    *pa++ = ACCUM(w);
Packit 16808d
    w = CARRYOUT(w);
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* If there's an overall carry out, increase precision and include
Packit 16808d
     it.  We could have done this initially, but why touch the memory
Packit 16808d
     allocator unless we're sure we have to?
Packit 16808d
   */
Packit 16808d
  if(w) {
Packit 16808d
    if((res = s_mw_mp_pad(a, used + 1)) != MP_OKAY)
Packit 16808d
      return res;
Packit 16808d
Packit 16808d
    DIGIT(a, ix) = w;  /* pa may not be valid after s_mw_mp_pad() call */
Packit 16808d
  }
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_add() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_sub(a, b) */
Packit 16808d
Packit 16808d
/* Compute a = |a| - |b|, assumes |a| >= |b|                              */
Packit 16808d
mw_mp_err   s_mw_mp_sub(mw_mp_int *a, mw_mp_int *b)        /* magnitude subtract      */
Packit 16808d
{
Packit 16808d
  mw_mp_word   w = 0;
Packit 16808d
  mw_mp_digit *pa, *pb;
Packit 16808d
  mw_mp_size   ix, used = USED(b);
Packit 16808d
Packit 16808d
  /*
Packit 16808d
    Subtract and propagate borrow.  Up to the precision of b, this
Packit 16808d
    accounts for the digits of b; after that, we just make sure the
Packit 16808d
    carries get to the right place.  This saves having to pad b out to
Packit 16808d
    the precision of a just to make the loops work right...
Packit 16808d
   */
Packit 16808d
  pa = DIGITS(a);
Packit 16808d
  pb = DIGITS(b);
Packit 16808d
Packit 16808d
  for(ix = 0; ix < used; ++ix) {
Packit 16808d
    w = (RADIX + *pa) - w - *pb++;
Packit 16808d
    *pa++ = ACCUM(w);
Packit 16808d
    w = CARRYOUT(w) ? 0 : 1;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  used = USED(a);
Packit 16808d
  while(ix < used) {
Packit 16808d
    w = RADIX + *pa - w;
Packit 16808d
    *pa++ = ACCUM(w);
Packit 16808d
    w = CARRYOUT(w) ? 0 : 1;
Packit 16808d
    ++ix;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Clobber any leading zeroes we created    */
Packit 16808d
  s_mw_mp_clamp(a);
Packit 16808d
Packit 16808d
  /* 
Packit 16808d
     If there was a borrow out, then |b| > |a| in violation
Packit 16808d
     of our input invariant.  We've already done the work,
Packit 16808d
     but we'll at least complain about it...
Packit 16808d
   */
Packit 16808d
  if(w)
Packit 16808d
    return MP_RANGE;
Packit 16808d
  else
Packit 16808d
    return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_sub() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_mul(a, b) */
Packit 16808d
Packit 16808d
/* Compute a = |a| * |b|                                                  */
Packit 16808d
mw_mp_err   s_mw_mp_mul(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_word   w, k = 0;
Packit 16808d
  mw_mp_int    tmp;
Packit 16808d
  mw_mp_err    res;
Packit 16808d
  mw_mp_size   ix, jx, ua = USED(a), ub = USED(b);
Packit 16808d
  mw_mp_digit *pa, *pb, *pt, *pbt;
Packit 16808d
Packit 16808d
  if((res = mw_mp_init_size(&tmp, ua + ub)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  /* This has the effect of left-padding with zeroes... */
Packit 16808d
  USED(&tmp) = ua + ub;
Packit 16808d
Packit 16808d
  /* We're going to need the base value each iteration */
Packit 16808d
  pbt = DIGITS(&tmp);
Packit 16808d
Packit 16808d
  /* Outer loop:  Digits of b */
Packit 16808d
Packit 16808d
  pb = DIGITS(b);
Packit 16808d
  for(ix = 0; ix < ub; ++ix, ++pb) {
Packit 16808d
    if(*pb == 0) 
Packit 16808d
      continue;
Packit 16808d
Packit 16808d
    /* Inner product:  Digits of a */
Packit 16808d
    pa = DIGITS(a);
Packit 16808d
    for(jx = 0; jx < ua; ++jx, ++pa) {
Packit 16808d
      pt = pbt + ix + jx;
Packit 16808d
      w = *pb * *pa + k + *pt;
Packit 16808d
      *pt = ACCUM(w);
Packit 16808d
      k = CARRYOUT(w);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    pbt[ix + jx] = k;
Packit 16808d
    k = 0;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(&tmp);
Packit 16808d
  s_mw_mp_exch(&tmp, a);
Packit 16808d
Packit 16808d
  mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_mul() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_kmul(a, b, out, len) */
Packit 16808d
Packit 16808d
#if 0
Packit 16808d
void   s_mw_mp_kmul(mw_mp_digit *a, mw_mp_digit *b, mw_mp_digit *out, mw_mp_size len)
Packit 16808d
{
Packit 16808d
  mw_mp_word   w, k = 0;
Packit 16808d
  mw_mp_size   ix, jx;
Packit 16808d
  mw_mp_digit *pa, *pt;
Packit 16808d
Packit 16808d
  for(ix = 0; ix < len; ++ix, ++b) {
Packit 16808d
    if(*b == 0)
Packit 16808d
      continue;
Packit 16808d
    
Packit 16808d
    pa = a;
Packit 16808d
    for(jx = 0; jx < len; ++jx, ++pa) {
Packit 16808d
      pt = out + ix + jx;
Packit 16808d
      w = *b * *pa + k + *pt;
Packit 16808d
      *pt = ACCUM(w);
Packit 16808d
      k = CARRYOUT(w);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    out[ix + jx] = k;
Packit 16808d
    k = 0;
Packit 16808d
  }
Packit 16808d
Packit 16808d
} /* end s_mw_mp_kmul() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_sqr(a) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Computes the square of a, in place.  This can be done more
Packit 16808d
  efficiently than a general multiplication, because many of the
Packit 16808d
  computation steps are redundant when squaring.  The inner product
Packit 16808d
  step is a bit more complicated, but we save a fair number of
Packit 16808d
  iterations of the multiplication loop.
Packit 16808d
 */
Packit 16808d
#if MP_SQUARE
Packit 16808d
mw_mp_err   s_mw_mp_sqr(mw_mp_int *a)
Packit 16808d
{
Packit 16808d
  mw_mp_word  w, k = 0;
Packit 16808d
  mw_mp_int   tmp;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_size  ix, jx, kx, used = USED(a);
Packit 16808d
  mw_mp_digit *pa1, *pa2, *pt, *pbt;
Packit 16808d
Packit 16808d
  if((res = mw_mp_init_size(&tmp, 2 * used)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  /* Left-pad with zeroes */
Packit 16808d
  USED(&tmp) = 2 * used;
Packit 16808d
Packit 16808d
  /* We need the base value each time through the loop */
Packit 16808d
  pbt = DIGITS(&tmp);
Packit 16808d
Packit 16808d
  pa1 = DIGITS(a);
Packit 16808d
  for(ix = 0; ix < used; ++ix, ++pa1) {
Packit 16808d
    if(*pa1 == 0)
Packit 16808d
      continue;
Packit 16808d
Packit 16808d
    w = DIGIT(&tmp, ix + ix) + (*pa1 * *pa1);
Packit 16808d
Packit 16808d
    pbt[ix + ix] = ACCUM(w);
Packit 16808d
    k = CARRYOUT(w);
Packit 16808d
Packit 16808d
    /*
Packit 16808d
      The inner product is computed as:
Packit 16808d
Packit 16808d
         (C, S) = t[i,j] + 2 a[i] a[j] + C
Packit 16808d
Packit 16808d
      This can overflow what can be represented in an mw_mp_word, and
Packit 16808d
      since C arithmetic does not provide any way to check for
Packit 16808d
      overflow, we have to check explicitly for overflow conditions
Packit 16808d
      before they happen.
Packit 16808d
     */
Packit 16808d
    for(jx = ix + 1, pa2 = DIGITS(a) + jx; jx < used; ++jx, ++pa2) {
Packit 16808d
      mw_mp_word  u = 0, v;
Packit 16808d
      
Packit 16808d
      /* Store this in a temporary to avoid indirections later */
Packit 16808d
      pt = pbt + ix + jx;
Packit 16808d
Packit 16808d
      /* Compute the multiplicative step */
Packit 16808d
      w = *pa1 * *pa2;
Packit 16808d
Packit 16808d
      /* If w is more than half MP_WORD_MAX, the doubling will
Packit 16808d
	 overflow, and we need to record a carry out into the next
Packit 16808d
	 word */
Packit 16808d
      u = (w >> (MP_WORD_BIT - 1)) & 1;
Packit 16808d
Packit 16808d
      /* Double what we've got, overflow will be ignored as defined
Packit 16808d
	 for C arithmetic (we've already noted if it is to occur)
Packit 16808d
       */
Packit 16808d
      w *= 2;
Packit 16808d
Packit 16808d
      /* Compute the additive step */
Packit 16808d
      v = *pt + k;
Packit 16808d
Packit 16808d
      /* If we do not already have an overflow carry, check to see
Packit 16808d
	 if the addition will cause one, and set the carry out if so 
Packit 16808d
       */
Packit 16808d
      u |= ((MP_WORD_MAX - v) < w);
Packit 16808d
Packit 16808d
      /* Add in the rest, again ignoring overflow */
Packit 16808d
      w += v;
Packit 16808d
Packit 16808d
      /* Set the i,j digit of the output */
Packit 16808d
      *pt = ACCUM(w);
Packit 16808d
Packit 16808d
      /* Save carry information for the next iteration of the loop.
Packit 16808d
	 This is why k must be an mw_mp_word, instead of an mw_mp_digit */
Packit 16808d
      k = CARRYOUT(w) | (u << DIGIT_BIT);
Packit 16808d
Packit 16808d
    } /* for(jx ...) */
Packit 16808d
Packit 16808d
    /* Set the last digit in the cycle and reset the carry */
Packit 16808d
    k = DIGIT(&tmp, ix + jx) + k;
Packit 16808d
    pbt[ix + jx] = ACCUM(k);
Packit 16808d
    k = CARRYOUT(k);
Packit 16808d
Packit 16808d
    /* If we are carrying out, propagate the carry to the next digit
Packit 16808d
       in the output.  This may cascade, so we have to be somewhat
Packit 16808d
       circumspect -- but we will have enough precision in the output
Packit 16808d
       that we won't overflow 
Packit 16808d
     */
Packit 16808d
    kx = 1;
Packit 16808d
    while(k) {
Packit 16808d
      k = pbt[ix + jx + kx] + 1;
Packit 16808d
      pbt[ix + jx + kx] = ACCUM(k);
Packit 16808d
      k = CARRYOUT(k);
Packit 16808d
      ++kx;
Packit 16808d
    }
Packit 16808d
  } /* for(ix ...) */
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(&tmp);
Packit 16808d
  s_mw_mp_exch(&tmp, a);
Packit 16808d
Packit 16808d
  mw_mp_clear(&tmp);
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_sqr() */
Packit 16808d
#endif
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_div(a, b) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  s_mw_mp_div(a, b)
Packit 16808d
Packit 16808d
  Compute a = a / b and b = a mod b.  Assumes b > a.
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_div(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_int   quot, rem, t;
Packit 16808d
  mw_mp_word  q;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_digit d;
Packit 16808d
  int      ix;
Packit 16808d
Packit 16808d
  if(mw_mp_cmw_mp_z(b) == 0)
Packit 16808d
    return MP_RANGE;
Packit 16808d
Packit 16808d
  /* Shortcut if b is power of two */
Packit 16808d
  if((ix = s_mw_mp_ispow2(b)) >= 0) {
Packit 16808d
    mw_mp_copy(a, b);  /* need this for remainder */
Packit 16808d
    s_mw_mp_div_2d(a, (mw_mp_digit)ix);
Packit 16808d
    s_mw_mp_mod_2d(b, (mw_mp_digit)ix);
Packit 16808d
Packit 16808d
    return MP_OKAY;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Allocate space to store the quotient */
Packit 16808d
  if((res = mw_mp_init_size(&quot, USED(a))) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  /* A working temporary for division     */
Packit 16808d
  if((res = mw_mp_init_size(&t, USED(a))) != MP_OKAY)
Packit 16808d
    goto T;
Packit 16808d
Packit 16808d
  /* Allocate space for the remainder     */
Packit 16808d
  if((res = mw_mp_init_size(&rem, USED(a))) != MP_OKAY)
Packit 16808d
    goto REM;
Packit 16808d
Packit 16808d
  /* Normalize to optimize guessing       */
Packit 16808d
  d = s_mw_mp_norm(a, b);
Packit 16808d
Packit 16808d
  /* Perform the division itself...woo!   */
Packit 16808d
  ix = USED(a) - 1;
Packit 16808d
Packit 16808d
  while(ix >= 0) {
Packit 16808d
    /* Find a partial substring of a which is at least b */
Packit 16808d
    while(s_mw_mp_cmp(&rem, b) < 0 && ix >= 0) {
Packit 16808d
      if((res = s_mw_mp_lshd(&rem, 1)) != MP_OKAY) 
Packit 16808d
	goto CLEANUP;
Packit 16808d
Packit 16808d
      if((res = s_mw_mp_lshd(&quot, 1)) != MP_OKAY)
Packit 16808d
	goto CLEANUP;
Packit 16808d
Packit 16808d
      DIGIT(&rem, 0) = DIGIT(a, ix);
Packit 16808d
      s_mw_mp_clamp(&rem;;
Packit 16808d
      --ix;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    /* If we didn't find one, we're finished dividing    */
Packit 16808d
    if(s_mw_mp_cmp(&rem, b) < 0) 
Packit 16808d
      break;    
Packit 16808d
Packit 16808d
    /* Compute a guess for the next quotient digit       */
Packit 16808d
    q = DIGIT(&rem, USED(&rem) - 1);
Packit 16808d
    if(q <= DIGIT(b, USED(b) - 1) && USED(&rem) > 1)
Packit 16808d
      q = (q << DIGIT_BIT) | DIGIT(&rem, USED(&rem) - 2);
Packit 16808d
Packit 16808d
    q /= DIGIT(b, USED(b) - 1);
Packit 16808d
Packit 16808d
    /* The guess can be as much as RADIX + 1 */
Packit 16808d
    if(q >= RADIX)
Packit 16808d
      q = RADIX - 1;
Packit 16808d
Packit 16808d
    /* See what that multiplies out to                   */
Packit 16808d
    mw_mp_copy(b, &t);
Packit 16808d
    if((res = s_mw_mp_mul_d(&t, q)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
Packit 16808d
    /* 
Packit 16808d
       If it's too big, back it off.  We should not have to do this
Packit 16808d
       more than once, or, in rare cases, twice.  Knuth describes a
Packit 16808d
       method by which this could be reduced to a maximum of once, but
Packit 16808d
       I didn't implement that here.
Packit 16808d
     */
Packit 16808d
    while(s_mw_mp_cmp(&t, &rem) > 0) {
Packit 16808d
      --q;
Packit 16808d
      s_mw_mp_sub(&t, b);
Packit 16808d
    }
Packit 16808d
Packit 16808d
    /* At this point, q should be the right next digit   */
Packit 16808d
    if((res = s_mw_mp_sub(&rem, &t)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
Packit 16808d
    /*
Packit 16808d
      Include the digit in the quotient.  We allocated enough memory
Packit 16808d
      for any quotient we could ever possibly get, so we should not
Packit 16808d
      have to check for failures here
Packit 16808d
     */
Packit 16808d
    DIGIT(&quot, 0) = q;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Denormalize remainder                */
Packit 16808d
  if(d != 0) 
Packit 16808d
    s_mw_mp_div_2d(&rem, d);
Packit 16808d
Packit 16808d
  s_mw_mp_clamp(");
Packit 16808d
  s_mw_mp_clamp(&rem;;
Packit 16808d
Packit 16808d
  /* Copy quotient back to output         */
Packit 16808d
  s_mw_mp_exch(&quot, a);
Packit 16808d
  
Packit 16808d
  /* Copy remainder back to output        */
Packit 16808d
  s_mw_mp_exch(&rem, b);
Packit 16808d
Packit 16808d
CLEANUP:
Packit 16808d
  mw_mp_clear(&rem;;
Packit 16808d
REM:
Packit 16808d
  mw_mp_clear(&t);
Packit 16808d
T:
Packit 16808d
  mw_mp_clear(");
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_div() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_2expt(a, k) */
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_2expt(mw_mp_int *a, mw_mp_digit k)
Packit 16808d
{
Packit 16808d
  mw_mp_err    res;
Packit 16808d
  mw_mp_size   dig, bit;
Packit 16808d
Packit 16808d
  dig = k / DIGIT_BIT;
Packit 16808d
  bit = k % DIGIT_BIT;
Packit 16808d
Packit 16808d
  mw_mp_zero(a);
Packit 16808d
  if((res = s_mw_mp_pad(a, dig + 1)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
  
Packit 16808d
  DIGIT(a, dig) |= (1 << bit);
Packit 16808d
Packit 16808d
  return MP_OKAY;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_2expt() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_reduce(x, m, mu) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Compute Barrett reduction, x (mod m), given a precomputed value for
Packit 16808d
  mu = b^2k / m, where b = RADIX and k = #digits(m).  This should be
Packit 16808d
  faster than straight division, when many reductions by the same
Packit 16808d
  value of m are required (such as in modular exponentiation).  This
Packit 16808d
  can nearly halve the time required to do modular exponentiation,
Packit 16808d
  as compared to using the full integer divide to reduce.
Packit 16808d
Packit 16808d
  This algorithm was derived from the _Handbook of Applied
Packit 16808d
  Cryptography_ by Menezes, Oorschot and VanStone, Ch. 14,
Packit 16808d
  pp. 603-604.  
Packit 16808d
 */
Packit 16808d
Packit 16808d
mw_mp_err   s_mw_mp_reduce(mw_mp_int *x, mw_mp_int *m, mw_mp_int *mu)
Packit 16808d
{
Packit 16808d
  mw_mp_int   q;
Packit 16808d
  mw_mp_err   res;
Packit 16808d
  mw_mp_size  um = USED(m);
Packit 16808d
Packit 16808d
  if((res = mw_mp_init_copy(&q, x)) != MP_OKAY)
Packit 16808d
    return res;
Packit 16808d
Packit 16808d
  s_mw_mp_rshd(&q, um - 1);       /* q1 = x / b^(k-1)  */
Packit 16808d
  s_mw_mp_mul(&q, mu);            /* q2 = q1 * mu      */
Packit 16808d
  s_mw_mp_rshd(&q, um + 1);       /* q3 = q2 / b^(k+1) */
Packit 16808d
Packit 16808d
  /* x = x mod b^(k+1), quick (no division) */
Packit 16808d
  s_mw_mp_mod_2d(x, DIGIT_BIT * (um + 1));
Packit 16808d
Packit 16808d
  /* q = q * m mod b^(k+1), quick (no division) */
Packit 16808d
  s_mw_mp_mul(&q, m);
Packit 16808d
  s_mw_mp_mod_2d(&q, DIGIT_BIT * (um + 1));
Packit 16808d
Packit 16808d
  /* x = x - q */
Packit 16808d
  if((res = mw_mp_sub(x, &q, x)) != MP_OKAY)
Packit 16808d
    goto CLEANUP;
Packit 16808d
Packit 16808d
  /* If x < 0, add b^(k+1) to it */
Packit 16808d
  if(mw_mp_cmw_mp_z(x) < 0) {
Packit 16808d
    mw_mp_set(&q, 1);
Packit 16808d
    if((res = s_mw_mp_lshd(&q, um + 1)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
    if((res = mw_mp_add(x, &q, x)) != MP_OKAY)
Packit 16808d
      goto CLEANUP;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  /* Back off if it's too big */
Packit 16808d
  while(mw_mp_cmp(x, m) >= 0) {
Packit 16808d
    if((res = s_mw_mp_sub(x, m)) != MP_OKAY)
Packit 16808d
      break;
Packit 16808d
  }
Packit 16808d
Packit 16808d
 CLEANUP:
Packit 16808d
  mw_mp_clear(&q);
Packit 16808d
Packit 16808d
  return res;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_reduce() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Primitive comparisons */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_cmp(a, b) */
Packit 16808d
Packit 16808d
/* Compare |a| <=> |b|, return 0 if equal, <0 if a<b, >0 if a>b           */
Packit 16808d
int      s_mw_mp_cmp(mw_mp_int *a, mw_mp_int *b)
Packit 16808d
{
Packit 16808d
  mw_mp_size   ua = USED(a), ub = USED(b);
Packit 16808d
Packit 16808d
  if(ua > ub)
Packit 16808d
    return MP_GT;
Packit 16808d
  else if(ua < ub)
Packit 16808d
    return MP_LT;
Packit 16808d
  else {
Packit 16808d
    int      ix = ua - 1;
Packit 16808d
    mw_mp_digit *ap = DIGITS(a) + ix, *bp = DIGITS(b) + ix;
Packit 16808d
Packit 16808d
    while(ix >= 0) {
Packit 16808d
      if(*ap > *bp)
Packit 16808d
	return MP_GT;
Packit 16808d
      else if(*ap < *bp)
Packit 16808d
	return MP_LT;
Packit 16808d
Packit 16808d
      --ap; --bp; --ix;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    return MP_EQ;
Packit 16808d
  }
Packit 16808d
Packit 16808d
} /* end s_mw_mp_cmp() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_cmw_mp_d(a, d) */
Packit 16808d
Packit 16808d
/* Compare |a| <=> d, return 0 if equal, <0 if a<d, >0 if a>d             */
Packit 16808d
int      s_mw_mp_cmw_mp_d(mw_mp_int *a, mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  mw_mp_size  ua = USED(a);
Packit 16808d
  mw_mp_digit *ap = DIGITS(a);
Packit 16808d
Packit 16808d
  if(ua > 1)
Packit 16808d
    return MP_GT;
Packit 16808d
Packit 16808d
  if(*ap < d) 
Packit 16808d
    return MP_LT;
Packit 16808d
  else if(*ap > d)
Packit 16808d
    return MP_GT;
Packit 16808d
  else
Packit 16808d
    return MP_EQ;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_cmw_mp_d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_ispow2(v) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Returns -1 if the value is not a power of two; otherwise, it returns
Packit 16808d
  k such that v = 2^k, i.e. lg(v).
Packit 16808d
 */
Packit 16808d
int      s_mw_mp_ispow2(mw_mp_int *v)
Packit 16808d
{
Packit 16808d
  mw_mp_digit d, *dp;
Packit 16808d
  mw_mp_size  uv = USED(v);
Packit 16808d
  int      extra = 0, ix;
Packit 16808d
Packit 16808d
  d = DIGIT(v, uv - 1); /* most significant digit of v */
Packit 16808d
Packit 16808d
  while(d && ((d & 1) == 0)) {
Packit 16808d
    d >>= 1;
Packit 16808d
    ++extra;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(d == 1) {
Packit 16808d
    ix = uv - 2;
Packit 16808d
    dp = DIGITS(v) + ix;
Packit 16808d
Packit 16808d
    while(ix >= 0) {
Packit 16808d
      if(*dp)
Packit 16808d
	return -1; /* not a power of two */
Packit 16808d
Packit 16808d
      --dp; --ix;
Packit 16808d
    }
Packit 16808d
Packit 16808d
    return ((uv - 1) * DIGIT_BIT) + extra;
Packit 16808d
  } 
Packit 16808d
Packit 16808d
  return -1;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_ispow2() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_ispow2d(d) */
Packit 16808d
Packit 16808d
int      s_mw_mp_ispow2d(mw_mp_digit d)
Packit 16808d
{
Packit 16808d
  int   pow = 0;
Packit 16808d
Packit 16808d
  while((d & 1) == 0) {
Packit 16808d
    ++pow; d >>= 1;
Packit 16808d
  }
Packit 16808d
Packit 16808d
  if(d == 1)
Packit 16808d
    return pow;
Packit 16808d
Packit 16808d
  return -1;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_ispow2d() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ Primitive I/O helpers */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_tovalue(ch, r) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Convert the given character to its digit value, in the given radix.
Packit 16808d
  If the given character is not understood in the given radix, -1 is
Packit 16808d
  returned.  Otherwise the digit's numeric value is returned.
Packit 16808d
Packit 16808d
  The results will be odd if you use a radix < 2 or > 62, you are
Packit 16808d
  expected to know what you're up to.
Packit 16808d
 */
Packit 16808d
int      s_mw_mp_tovalue(char ch, int r)
Packit 16808d
{
Packit 16808d
  int    val, xch;
Packit 16808d
  
Packit 16808d
  if(r > 36)
Packit 16808d
    xch = ch;
Packit 16808d
  else
Packit 16808d
    xch = toupper(ch);
Packit 16808d
Packit 16808d
  if(isdigit(xch))
Packit 16808d
    val = xch - '0';
Packit 16808d
  else if(isupper(xch))
Packit 16808d
    val = xch - 'A' + 10;
Packit 16808d
  else if(islower(xch))
Packit 16808d
    val = xch - 'a' + 36;
Packit 16808d
  else if(xch == '+')
Packit 16808d
    val = 62;
Packit 16808d
  else if(xch == '/')
Packit 16808d
    val = 63;
Packit 16808d
  else 
Packit 16808d
    return -1;
Packit 16808d
Packit 16808d
  if(val < 0 || val >= r)
Packit 16808d
    return -1;
Packit 16808d
Packit 16808d
  return val;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_tovalue() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_todigit(val, r, low) */
Packit 16808d
Packit 16808d
/*
Packit 16808d
  Convert val to a radix-r digit, if possible.  If val is out of range
Packit 16808d
  for r, returns zero.  Otherwise, returns an ASCII character denoting
Packit 16808d
  the value in the given radix.
Packit 16808d
Packit 16808d
  The results may be odd if you use a radix < 2 or > 64, you are
Packit 16808d
  expected to know what you're doing.
Packit 16808d
 */
Packit 16808d
  
Packit 16808d
char     s_mw_mp_todigit(int val, int r, int low)
Packit 16808d
{
Packit 16808d
  char   ch;
Packit 16808d
Packit 16808d
  if(val < 0 || val >= r)
Packit 16808d
    return 0;
Packit 16808d
Packit 16808d
  ch = s_dmap_1[val];
Packit 16808d
Packit 16808d
  if(r <= 36 && low)
Packit 16808d
    ch = tolower(ch);
Packit 16808d
Packit 16808d
  return ch;
Packit 16808d
Packit 16808d
} /* end s_mw_mp_todigit() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* {{{ s_mw_mp_outlen(bits, radix) */
Packit 16808d
Packit 16808d
/* 
Packit 16808d
   Return an estimate for how long a string is needed to hold a radix
Packit 16808d
   r representation of a number with 'bits' significant bits.
Packit 16808d
Packit 16808d
   Does not include space for a sign or a NUL terminator.
Packit 16808d
 */
Packit 16808d
int      s_mw_mp_outlen(int bits, int r)
Packit 16808d
{
Packit 16808d
  return (int)((double)bits * LOG_V_2(r) + 0.5);
Packit 16808d
Packit 16808d
} /* end s_mw_mp_outlen() */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/* }}} */
Packit 16808d
Packit 16808d
/*------------------------------------------------------------------------*/
Packit 16808d
/* HERE THERE BE DRAGONS                                                  */