vnl_bignum.cxx

InsightToolkit-1.4.0(有大量的优化算法程序)

  if (b2.is_infinity()) return -1;
  if (b1.count > b2.count) return 1;    // If one has more data than
  if (b2.count > b1.count) return -1;   //   the other, it wins
  Counter i = b1.count;                 // Else same number of elmts
  while (i > 0) {                       // Do lexicographic comparison
    if (b1.data[i - 1] > b2.data[i - 1])
      return 1;
    else if (b1.data[i - 1] < b2.data[i - 1])
      return -1;
    i--;
  }                                     // No data, or all elmts same
  return 0;                             //  so must be equal
}


//: multiply a non-infinite vnl_bignum by a "single digit"
// - Inputs:  vnl_bignum reference, single word multiplier, reference to the product,
//            and index of starting storage location to use in product

void multiply_aux (const vnl_bignum& b, Data d, vnl_bignum& prod, Counter i)
{
  // this function works just like normal multiplication by hand, in that the
  // top number is multiplied by the first digit of the bottom number, then the
  // second digit, and so on.  The only difference is that instead of doing all
  // of the multiplication before adding the rows, addition is done
  // concurrently.
  if (i == 0) {                         // if index is zero
    Counter j = 0;                      //   then zero out all of
    while (j < prod.count)              //   prod's data elements
      prod.data[j++] = 0;
  }
  if (d != 0) {                         // if d == 0, nothing to do
    unsigned long temp;
    Data carry = 0;

    Counter j = 0;
    for (; j < b.count; j++) {
      // for each of b's data elmts, multiply times d and add running product
      temp = (unsigned long)b.data[j] * (unsigned long)d
           + (unsigned long)prod.data[i + j] + carry;
      prod.data[i + j] = Data(temp % 0x10000L); //   store result in product
      carry = Data(temp/0x10000L);              //   keep track of carry
    }
    if (i+j < prod.count)
      prod.data[i + j] = carry;                 // Done.  Store the final carry
  }
}

//: normalize two vnl_bignums
// (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
//  A digit here is one data element in the radix 2**2.)
// - Inputs:  references to two vnl_bignums b1, b2, and their normalized counterparts
// - Outputs: the integral normalization factor used

Data normalize (const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& u, vnl_bignum& v)
{
  Data d = Data(0x10000L/(long(b2.data[b2.count - 1]) + 1)); // Calculate normalization factor.
  u.data = new Data[u.count = b1.count + 1];    // Get data for u (plus extra)
  v.data = new Data[v.count = b2.count];        // Get data for v
  u.data[b1.count] = 0;                         // Set u's leading digit to 0
  multiply_aux(b1,d,u,0);                       // u = b1 * d
  multiply_aux(b2,d,v,0);                       // v = b2 * d
  return d;                                     // return normalization factor
}


//: divide a vnl_bignum by a "single digit"
// (Refer to Knuth, V.2, Section 4.3.2, exercise 16 for details.
//  A digit here is one data element in the radix 2**2.)
// - Inputs:  reference to vnl_bignum dividend, single digit divisor d, vnl_bignum
//            quotient, and single digit remainder r

void divide_aux (const vnl_bignum& b1, Data d, vnl_bignum& q, Data& r)
{
  r = 0;                                        // init remainder to zero
  unsigned long temp;
  for (Counter j = b1.count; j > 0; j--) {
    temp = (unsigned long)r*0x10000L + (unsigned long)b1.data[j - 1]; // get remainder, append next
    if (j < 1 + q.count)
      q.data[j - 1] = Data(temp/d);             //   digit, then divide
    r = Data(temp % d);                         // calculate new remainder
  }
}


//: estimate next dividend
// (Refer to Knuth, V.2, Section 4.3.1, Algorithm D for details.
//  This function estimates how many times v goes into u, starting at u's
//  jth digit.  A digit here is actually a data element, thought of as
//  being in the radix 2**2.)
// - Inputs:  reference to vnl_bignum dividend and divisor and starting digit j
// - Outputs: estimated number of times v goes into u

Data estimate_q_hat (const vnl_bignum& u, const vnl_bignum& v, Counter j)
{
  Data q_hat,
       v1 = v.data[v.count - 1],                // localize frequent data
       v2 = v.data[v.count - 2],
       u0 = u.data[u.count - 1 - j],
       u1 = u.data[u.count - 2 - j],
       u2 = u.data[u.count - 3 - j];

  // Initial Knuth estimate, usually correct
  q_hat = (u0 == v1 ? Data(0xffff) : Data(long(u0 * 0x10000L + u1) / long(v1)));

  // high speed test to determine most of the cases in which initial
  // estimate is too large.  Eliminates most cases in which q_hat is one too
  // large.  Eliminates all cases in which q_hat is two too large.  The test
  // looks hairy because we have to watch out for overflow.  In the book, this
  // test is the simple inequality:
  //     v2*q_hat > (u0*0x10000L + u1 - q_hat*v1)*0x10000L + u2.
  // If the inequality is true, decrease q_hat by 1.  If inequality is still
  // true, decrease q_hat again.
  unsigned long lhs, rhs;               // lhs, rhs of Knuth inequality
  for (Counter i = 0; i < 2; i++) {     // loop at most twice
    lhs = v2 * q_hat;                   // Calculate left-hand side of ineq.
    rhs = u0 * 0x10000L + u1;// Calculate part of right-hand side
    rhs -= (q_hat * v1);                // Now subtract off part

    // DML:  My attempt to fix the overflow testing bug..
    double temp_rhs = double(rhs);
    double temp_radix_s = double(0x10000);
    // OLD WAY: if (rhs > rhs * 0x10000)// if multiplication causes overflow
    // NEW WAY: see if result won't fit into a long.
    if ( temp_rhs * temp_radix_s > double(0x7fffffffL) )
      break;                            //   then rhs > lhs, so test fails
    rhs *= 0x10000L;                    // No overflow:  ok to multiply

    temp_rhs = double(rhs);
    double temp_u2 = double(u2);
    // OLD WAY: if (rhs > rhs + u2)     // if addition yields overflow
    if ( temp_rhs + temp_u2 > double(0x7fffffffL) ) // NEW WAY.
      break;                            //   then rhs > lhs, so test fails
    rhs += u2;                          // No overflow: ok to add.
    if (lhs <= rhs)                     // if lhs <= rhs
      break;                            //   test fails
    q_hat--;                            // Test passes:  decrement q_hat
  }                                     // Loop again
  return q_hat;                         // Return estimate
}


//: calculate u - v*q_hat
// (Refer to Knuth, V. 2, Section 4.3.1, Algorithm D for details.
//  A digit here is a data element, thought of as being in the radix 2**2.)
// - Inputs:  reference to vnl_bignum dividend, divisor, estimated result, and index
//            into jth digit of dividend
// - Outputs: true number of times v goes into u

Data multiply_subtract (vnl_bignum& u, const vnl_bignum& v, Data q_hat, Counter j)
{
  // At this point it has been estimated that v goes into the jth and higher
  // digits of u about q_hat times, and in fact that q_hat is either the
  // correct number of times or one too large.

  if (q_hat == 0) return q_hat;         // if q_hat 0, nothing to do
  vnl_bignum rslt;                      // create a temporary vnl_bignum
  Counter tmpcnt;
  rslt.data =                           // allocate data for it
     new Data[rslt.count = v.count + 1];

  // simultaneous computation of u - v*q_hat
  unsigned long prod, diff;
  Data carry = 0, borrow = 0;
  Counter i = 0;
  for (; i < v.count; i++) {
    // for each digit of v, multiply it by q_hat and subtract the result
    prod = (unsigned long)v.data[i] * (unsigned long)q_hat + carry;
    diff = (unsigned long)u.data[u.count - v.count - 1 - j + i] + 0x10000L - borrow;
    diff -= (unsigned long)Data(prod);  //   form proper digit of u
    rslt.data[i] = Data(diff);          //   save the result
    borrow = (diff/0x10000L == 0);      //   keep track of any borrows
    carry = Data(prod/0x10000L);        //   keep track of carries
  }
  tmpcnt = u.count - v.count - 1 - j + i;
  diff = (unsigned long)u.data[tmpcnt]  //  special case for the last
         + 0x10000L - borrow;           //  digit
  diff -= carry;
  rslt.data[i] = Data(diff);
  borrow = (diff/0x10000L == 0);

  // A leftover borrow indicates that u - v*q_hat is negative, i.e., that
  // q_hat was one too large.  So to get correct result, decrement q_hat and
  // add back one multiple of v
  if (borrow) {
    q_hat--;
    carry = 0;
    unsigned long sum;
    for (i = 0; i < v.count; i++) {
      sum = (unsigned long)rslt.data[i] + (unsigned long)v.data[i] + carry;
      carry = Data(sum/0x10000L);
      u.data[u.count - v.count - 1 - j + i] = Data(sum);
    }
    u.data[u.count - v.count - 1 - j + i] = rslt.data[i] + carry;
  }
  else {                                // otherwise, result is ok
    for (i = 0; i < rslt.count; i++)    // store result back into u
      u.data[u.count - v.count - 1 - j + i] = rslt.data[i];
  }
  return q_hat;                         // return corrected q_hat
}


//: divide b2 into b1, getting quotient q and remainder r.
// (Refer to Knuth, V.2, Seciton 4.3.1, Algorithm D for details.
//  This function implements Algorithm D.)
// - Inputs:  references to a vnl_bignum dividend b1, divisor b2, quotient q, and
//            remainder r.

void divide (const vnl_bignum& b1, const vnl_bignum& b2, vnl_bignum& q, vnl_bignum& r)
{
  // Note that q or r may *not* be identical to either b1 or b2 !
  assert(&b1 != &q && &b2 != &q && &b1 != &r && &b2 != &r);
  q = r = 0L;
  if (b1 == 0L)                      // If divisor is zero
    return;                          //   return zero quotient and remainder
  int mag = magnitude_cmp(b1,b2);    // Compare magnitudes
  if (mag < 0)                       // if abs(b1) < abs(b2)
    r = b1;                          //   return zero quotient, b1 remainder
  else if (mag == 0)                 // if abs(b1) == abs(b2)
    q = 1L;                          //   quotient is 1, remainder is 0
  else {                             // otherwise abs(b1) > abs(b2), so divide
    q.data = new Data[q.count = b1.count - b2.count + 1]; // Allocate quotient
    r.data = new Data[r.count = b2.count];                // Allocate remainder
    if (b2.count == 1) {                        // Single digit divisor?
      divide_aux(b1,b2.data[0],q,r.data[0]);    // Do single digit divide
    }
    else {                                      // Else full-blown divide
      vnl_bignum u,v;
      Data d = normalize(b1,b2,u,v);            // Set u = b1/d, v = b2/d
      Data q_hat;                               // Multiplier
      Counter j = 0;
      while (j <= b1.count - b2.count) {        // Main division loop
        q_hat = estimate_q_hat(u,v,j);          // Estimate # times v divides
        q.data[q.count - 1 - j] =               // Do division, get true answ.
          multiply_subtract(u,v,q_hat,j);
        j++;
      }
      static Data dufus;                // dummy variable
      divide_aux(u,d,r,dufus);          // Unnormalize u for remainder
    }
    q.trim();                           // Trim leading zeros of quot.
    r.trim();                           // Trim leading zeros of rem.
  }
  q.sign = r.sign = b1.sign * b2.sign;  // Calculate signs
}


//: left shift (arithmetic) non-infinite vnl_bignum by positive number.
// - Inputs:  reference to vnl_bignum, positive shift value
// - Outputs: vnl_bignum, multiplied by the corresponding power of two

vnl_bignum left_shift (const vnl_bignum& b1, int l)
{
  // to carry out this arithmetic left shift, we cheat.  Instead of physically
  // shifting the data array l bits to the left, we shift just enough to get
  // the correct word alignment, and then pad the array on the right with as
  // many zeros as we need.
  vnl_bignum rslt;                      // result of shift
  rslt.sign = b1.sign;                  // result follows sign of input
  Counter growth = Counter(l / 16);     // # of words rslt will grow by
  Data shift = Data(l % 16);            // amount to actually shift
  Data rshift = 16 - shift;             // amount to shift next word by
  Data carry =                          // value that will be shifted
    b1.data[b1.count - 1] >> (16 - shift); // out end of current array
  rslt.data =                              // allocate new data array
    new Data[rslt.count = b1.count + growth + (carry != 0 ? 1 : 0)];
  Counter i = 0;
  while (i < growth)                            // zero out padded elements
    rslt.data[i++] = 0;
  rslt.data[i++] = b1.data[0] << shift;         // shift first non-zero element
  while (i < rslt.count - 1) {                  // for remaining data words
    rslt.data[i] = (b1.data[i - growth] << shift) + // shift current data word
      (b1.data[i - 1 - growth] >> rshift);          // propagate adjacent
    i++;                                        // carry into current word
  }
  if (i < rslt.count) {
    if (carry)                                  // if last word had overflow
      rslt.data[i] = carry;                     //   store it new data
    else                                        // otherwise,
      rslt.data[i] = (b1.data[i - growth] << shift) // do like the rest
                   + (b1.data[i - 1 - growth] >> rshift);
  }
  vnl_bignum& result = *((vnl_bignum*) &rslt);// same physical object
  return result;                              // shallow swap on return
}


//: right shift (arithmetic) non-infinite vnl_bignum by positive number.
// - Inputs:  reference to vnl_bignum, positive shift value
// - Outputs: vnl_bignum, divided by the corresponding power of two

vnl_bignum right_shift (const vnl_bignum& b1, int l)
{
  vnl_bignum rslt;                              // result of shift
  Counter shrinkage = Counter(l / 16);          // # of words rslt will shrink
  Data shift = Data(l % 16);                    // amount to actually shift
  Data dregs = (b1.data[b1.count-1] >> shift);  // high end data to save
  if (shrinkage + (dregs == 0) < b1.count) {    // if not all data shifted out
    rslt.sign = b1.sign;                        // rslt follows sign of input
                                                // allocate new data
    rslt.data = new Data[rslt.count = b1.count - shrinkage - (dregs == 0 ? 1 : 0)];
    Data lshift = 16 - shift;                   // amount to shift high word
    Counter i = 0;
    while (i < rslt.count - 1) {                // shift current word
      rslt.data[i] = (b1.data[i + shrinkage] >> shift) + // propagate adjacent
        (b1.data[i + shrinkage + 1] << lshift); //   word into current word
      i++;
    }
    if (dregs)                                  // don't lose dregs
      rslt.data[i] = dregs;
    else {
      rslt.data[i] = (b1.data[i + shrinkage] >> shift) +
        (b1.data[i + shrinkage + 1] << lshift);
    }
  }
  vnl_bignum& result = *((vnl_bignum*) &rslt);  // same physical object
  return result;                                // shallow swap on return
}