gmp_si.cpp

伯克利做的SFTP安全文件传输协议

// gmp_si.cc
// implementation of Integer using GMP
// copyright SafeTP Development Group, Inc., 2000  Terms of use are as specified in license.txt
// last updated: 9/08/00 23:27

#include <stdio.h>     // FILE, etc.
#include "gmp.h"       // gmp
#include "sint.h"      // Integer
#include <iostream.h>  // cout
#include <stdlib.h>    // rand
#include "nonport.h"   // portableLock, portableUnlock

#ifndef SAFETP
  #include <assert.h>    // assert
  #define xassert assert
#else
  #include "xassert.h"   // xassert
#endif


// the rule for the following code is that every mpz_ function must be
// preceeded by LOCKER.  portableLock() and portableUnlock() are required
// to support nested calls, but nevertheless the code below tries to
// avoid nesting them
//
// UPDATE: it turns out the GMP problem was caused by misconfiguring it
// such that it thought 'alloca' wasn't available, and therefore used a
// non-thread-safe allocator; I've left the LOCKER calls in the code in
// case we need them again for some reason, but the following two lines
// make them do nothing for now
#undef LOCKER
#define LOCKER


// ----------- mapping between 'void *impl' and mpz_t ----------------
// return an mpz_t for modification
inline __mpz_struct *mpzm(Integer &obj)
{
  return (__mpz_struct*)(obj.impl);
}

// mpz_t that won't be modified
inline __mpz_struct *mpzc(Integer const &obj)
{
  return (__mpz_struct*)(obj.impl);
}

// NOTE: since the GMP functions don't include constness,
//       const-correctness is *not* fully automatically
//       enforced within this file

// for accessing the mpz_t stored in 'this'
#define MPZM mpzm(*this)
#define MPZC mpzc(*this)



// ---------------- private helpers ----------------------
void init(Integer &ths)
{
  ((__mpz_struct*&)(ths.impl)) = new __mpz_struct;
  LOCKER
  mpz_init(mpzm(ths));
}


// because mpz/random.c requires a fn called random() to exist,
// and it doesn't on Borland
#ifdef __BORLANDC__
extern "C" long random();
long random()
{
  // rand() returns 16 bits, we want 32
  return (rand() << 16) ^ rand();
}
#endif


// ---------- implementation of public interface ------------
Integer::Integer()
{
  init(*this);
  LOCKER
  mpz_set_ui(MPZM, 0);
}


Integer::Integer(Integer const &obj)
{
  init(*this);
  LOCKER
  mpz_set(MPZM, mpzc(obj));
}


Integer::Integer(long n)
{
  init(*this);
  fromLong(n);
}


Integer::Integer(char const *str, int radix)
{
  init(*this);
  fromString(str, radix);
}


Integer::~Integer()
{
  LOCKER
  mpz_clear(MPZM);
  delete MPZM;
}


bool Integer::isZero() const
{
  LOCKER
  return mpz_sgn(MPZC) == 0;
}

bool Integer::isPositive() const
{
  LOCKER
  return mpz_sgn(MPZC) > 0;
}

bool Integer::isNegative() const
{
  LOCKER
  return mpz_sgn(MPZC) < 0;
}

int Integer::getSign() const
{
  LOCKER
  return mpz_sgn(MPZC);
}


Integer& Integer::operator= (Integer const &obj)
{
  if (this != &obj) {
    LOCKER
    mpz_set(MPZM, mpzc(obj));
  }
  return *this;
}


Integer& Integer::operator= (long n)
{
  fromLong(n);
  return *this;
}


// these values assume 2's complement representation
const long mostNegativeLong = 1 << (sizeof(long)*8 - 1);
const long mostPositiveLong = ~mostNegativeLong;

const Integer mostNegativeLongInt(mostNegativeLong);
const Integer mostPositiveLongInt(mostPositiveLong);

bool Integer::fitsInLong() const
{
  return mostNegativeLongInt <= *this &&
                                *this <= mostPositiveLongInt;
}

long Integer::toLong() const
{
  xassert(fitsInLong());
  LOCKER
  return mpz_get_si(MPZC);
}


void Integer::fromLong(long n)
{
  LOCKER
  mpz_set_si(MPZM, n);
}


int Integer::toStringBytes(int radix) const
{
  LOCKER
  return mpz_sizeinbase(MPZC, radix) + 2;
    // +2 for possible minus sign and null terminator
}

void Integer::toString(char *str, int length, int radix) const
{
  xassert(toStringBytes() <= length);
  LOCKER
  mpz_get_str(str, radix, MPZC);
}


void Integer::fromString(char const *str, int radix)
{
  xassert(str != NULL);
  xassert(2 <= radix && radix <= 36);
  LOCKER
  if (mpz_set_str(MPZM, str, radix) != 0) {
    xassert(!"string contains invalid digits");
  }
}


void Integer::toOstream(ostream &os) const
{
  int len = toStringBytes();
  char *temp = new char[len];
  toString(temp, len, 10);
  os << temp;
  delete[] temp;
}


int Integer::compare(Integer const &obj) const
{
  LOCKER
  return mpz_cmp(MPZC, mpzc(obj));
}

int Integer::compareToLong(long n) const
{
  LOCKER
  return mpz_cmp_si(MPZC, n);
}


#define ARITHOP(op, mpz_fn)                                \
  Integer Integer::operator op (Integer const &obj) const  \
  {                                                        \
    Integer ret;                                           \
    LOCKER                                                 \
    mpz_fn(mpzm(ret), MPZC, mpzc(obj));                    \
    return ret;                                            \
  }
ARITHOP(+, mpz_add)
ARITHOP(-, mpz_sub)
ARITHOP(*, mpz_mul)
ARITHOP(/, mpz_tdiv_q)
ARITHOP(%, mpz_tdiv_r)
#undef ARITHOP


#define ARITHEQOP(op, mpz_fn)                         \
  Integer& Integer::operator op (Integer const &obj)  \
  {                                                   \
    LOCKER                                            \
    mpz_fn(MPZM, MPZC, mpzc(obj));                    \
    return *this;                                     \
  }
ARITHEQOP(+=, mpz_add)
ARITHEQOP(-=, mpz_sub)
ARITHEQOP(*=, mpz_mul)
ARITHEQOP(/=, mpz_tdiv_q)
ARITHEQOP(%=, mpz_tdiv_r)
#undef ARITHEQOP


//Integer Integer::exp(Integer const &exponent) const
// GMP doesn't provide a function to implement this..!


void divAndRemainder(
  Integer &quotient, Integer &remainder,
  Integer const &numer, Integer const &denom)
{
  LOCKER
  mpz_tdiv_qr(mpzm(quotient), mpzm(remainder),
              mpzc(numer), mpzc(denom));
}


Integer log(Integer const &base, Integer const &power)
{
  xassert(power >= 1 && base >= 2);

  // I don't have a good way to implement this, but the
  // current callers of this fn are not time-critical...
  // so let's use successive division (horrors)

  // possible optimization: I can special-case all
  // bases that are powers of 2, using mpz_sizeinbase(2)

  Integer ret(0);
  Integer temp(power);
  while (temp > base-1) {
    temp /= base;
    ret += 1;
  }

  return ret;
}


long log_long(long base, Integer const &power)
{
  xassert(power >= 1 && base >= 2);

  // this is the special case I'm interested in
  // optimizing
  if (base == 256) {
    int bits = power.numBits();
    return (bits-1)/8;
      // bits == 8  ->  power < 256  ->   return 0
      // bits == 9  ->  power >= 256  ->   return 1
  }

  return log(base, power).toLong();
}


Integer abs(Integer const &value)
{
  if (value >= 0) {
    return value;
  }
  else {
    return value * -1;
  }
}


bool probablyPrime(Integer const &value, int iters)
{
  LOCKER
  return !!mpz_probab_prime_p(mpzc(value), iters);
}


int Integer::numBits() const
{
  xassert(*this >= 0);

  if (*this == 0) {
    return 1;
  }

  LOCKER
  return mpz_sizeinbase(MPZC, 2);
}


unsigned char Integer::leastSigByte() const
{
  xassert(*this >= 0);

  Integer temp;
  LOCKER
  mpz_fdiv_r_2exp(mpzm(temp), MPZC, 8);
  return (unsigned char)temp.toLong();
}


Integer Integer::operator >> (unsigned bits) const
{
  xassert(*this >= 0);
  Integer ret;
  LOCKER
  mpz_tdiv_q_2exp(mpzm(ret), MPZC, bits);
  return ret;
}

Integer& Integer::operator >>= (unsigned bits)
{
  xassert(*this >= 0);
  LOCKER
  mpz_tdiv_q_2exp(MPZM, MPZC, bits);
  return *this;
}


Integer Integer::operator << (unsigned bits) const
{
  xassert(*this >= 0);
  Integer ret;
  LOCKER
  mpz_mul_2exp(mpzm(ret), MPZC, bits);
  return ret;
}

Integer& Integer::operator <<= (unsigned bits)
{
  xassert(*this >= 0);
  LOCKER
  mpz_mul_2exp(MPZM, MPZC, bits);
  return *this;
}


#define MODARITHFN(name, op)                                               \
  Integer Integer::name(Integer const &obj, Integer const &modulus) const  \
  {                                                                        \
    return (*this op obj) % modulus;                                       \
  }
MODARITHFN(plusMod, +)
MODARITHFN(minusMod, -)
MODARITHFN(timesMod, *)
#undef MODARITHFN


Integer Integer::divMod(Integer const &denom, Integer const &modulus) const
{
  return timesMod(denom.inverseMod(modulus), modulus);
}


Integer Integer::expMod(Integer const &exponent, Integer const &modulus) const
{
  Integer ret;
  LOCKER
  mpz_powm(mpzm(ret), MPZC, mpzc(exponent), mpzc(modulus));
  return ret;
}


Integer Integer::inverseMod(Integer const &modulus) const
{
  Integer ret;
  //cout << "inverting " << *this
  //     << " mod " << modulus << endl;
  LOCKER
  if (!mpz_invert(mpzm(ret), MPZC, mpzc(modulus))) {
    //cout << "NO INVERSE!\n";
    ret = 0;
  }
  else {
    //cout << "answer: " << ret << endl;
  }

  // apparently, mpz_invert sometimes returns negative
  // answers, which need to be adjusted positive
  // (this is a known bug in GMP 2.0.2)
  if (ret < 0) {
    ret += modulus;
  }

  return ret;
}


Integer gcd(Integer const &a, Integer const &b)
{
  Integer ret;
  LOCKER
  mpz_gcd(mpzm(ret), mpzc(a), mpzc(b));
  return ret;
}


#define divRoundUp(numer, denom) (numer+denom-1)/denom

inline int nthBit(unsigned char const *bits, int which)
{
  return 1 & (bits[which/8] >> (which % 8));
}

Integer random(Integer::RandomSource &src,
               Integer const &top)
{
  xassert(top >= 1);
  int bits;
  {
    LOCKER
    bits = mpz_sizeinbase(mpzc(top), 2);
    // lock is released here
  }
  int bytes = divRoundUp(bits, 8);
  xassert(bytes > 0);

  unsigned char *randBits = new unsigned char[bytes];
  Integer ret;

  //cout << "max: " << top << endl;

  int iters = 0;
  do {
    // count iterations, to guard against infinite loop
    xassert(++iters < 200);    // no way this fails unless there's a bug

    // grab some more random bits
    src.getRandomBytes(randBits, bytes);

    // start with a zero value
    ret = 0;

    // set bits, starting at top, so we allocate full width at first
    for (int i=bits-1; i>=0; i--) {
      // the mpz bits are ordered from lsb (at index 0)
      // to msb (at index bits-1).. though I guess that
      // knowledge is not actually important here..
      if (nthBit(randBits, i)) {
        LOCKER
        mpz_setbit(mpzm(ret), i);
      }
      else {
        LOCKER
        mpz_clrbit(mpzm(ret), i);
      }
    }

    //cout << "generated: " << ret << endl;

    // if we generated a number that is too big, try again
  } while (ret >= top);

  delete[] randBits;
  return ret;
}


void extendedEuclidean(Integer &d, Integer &x, Integer &y,
                       Integer const &a, Integer const &b)
{
  LOCKER
  mpz_gcdext(mpzm(d), mpzm(x), mpzm(y), mpzc(a), mpzc(b));
}