bignum.c

来自「开放源码的编译器open watcom 1.6.0版的源代码」· C语言 代码 · 共 1,299 行 · 第 1/4 页

//                   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**sizeof(Data).)
// Inputs:  reference to CoolBignum dividend and divisor and starting digit j
// Outputs:  estimated number of times v goes into u

Data estimate_q_hat (const CoolBignum& u, const CoolBignum& 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 = Data(u0 == v1 ?                       
               radix_s - 1 :                    
               (u0*radix_s + u1)/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*radix_s + u1 - q_hat*v1)*radix_s + 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 * radix_s + 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(radix_s);
    // OLD WAY: if (rhs > rhs * radix_s)// if multiplication causes overflow
    // NEW WAY: see if result won't fit into a long.
    if ( temp_rhs * temp_radix_s > double(MAXLONG) )
      break;                            //   then rhs > lhs, so test fails
    rhs *= radix_s;                     // 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(MAXLONG) ) // 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
}

// Original version.
// Data estimate_q_hat (const CoolBignum& u, const CoolBignum& 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 = Data(u0 == v1 ?                    
//             radix_s - 1 :                    
//             (u0*radix_s + u1)/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*radix_s + u1 - q_hat*v1)*radix_s + 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 * radix_s + u1;         // Calculate part of right-hand side 
//     rhs -= (q_hat * v1);             // Now subtract off part
//     if (rhs > rhs * radix_s)         // if multiplication causes overflow
//       break;                         //   then rhs > lhs, so test fails
//     rhs *= radix_s;                  // No overflow:  ok to multiply
//     if (rhs > rhs + u2)                      // if addition yields overflow
//       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
// }



// multiply_subtract -- 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**sizeof(Data).)
// Inputs:  reference to CoolBignum dividend, divisor, estimated result, and index
//          into jth digit of dividend
// Outputs:  true number of times v goes into u

Data multiply_subtract (CoolBignum& u, const CoolBignum& 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
  CoolBignum rslt;                                      // create a temporary CoolBignum
  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;
  for (i = 0; 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]
           + (unsigned long)radix_s - borrow;
    diff -= (unsigned long)Data(prod);          //   form proper digit of u
    rslt.data[i] = Data(diff);                  //   save the result
    borrow = (diff/radix_s == 0);               //   keep track of any borrows
    carry = Data(prod/radix_s);                 //   keep track of carries
  }
  tmpcnt = u.count - v.count - 1 - j + i;
  diff = (unsigned long)u.data[tmpcnt]          //  special case for the last
         + (unsigned long)radix_s - borrow;     //     digit
  diff -= carry;
  rslt.data[i] = Data(diff);
  borrow = (diff/radix_s == 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/radix_s);
      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 - 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 CoolBignum dividend b1, divisor b2, quotient q, and
//          remainder r.
// Outputs:  none

void divide (const CoolBignum& b1, const CoolBignum& b2, CoolBignum& q, CoolBignum& r) {
  q = r = zero_s;
  if (b1 == zero_s)                  // 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 = one_s;                       //   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
      CoolBignum 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
}



// operator/ -- overload division operator for CoolBignums
// Inputs:  references to CoolBignum numerator and denominator
// Outputs:  CoolBignum quotient

CoolBignumE operator/ (const CoolBignum& b1, const CoolBignum& b2) {
  CoolBignum q,r;                               // Quotient and remainder
  if (b2.count == 0) {                          // If b2 is zero
    if (b1.count == 0) {                        // If b1 == zero
      q.state = N_DIVIDE_BY_ZERO;               // Both num & den are zero
      b1.divide_by_zero("operator/ ()");        // Raise the exception
    }
    else if (b1.sign > 0) {                     // If b1 > zero
      q.state = N_PLUS_INFINITY;                // Positive infinity
      b1.plus_infinity("operator/ ()");         // Raise the exception
    }
    else if (b1.sign < 0) {                     // If b1 < zero
      q.state = N_MINUS_INFINITY;               // Negative infinity
      b1.minus_infinity("operator/ ()");        // Raise the exception
    }
  }
  else {
    divide(b1,b2,q,r);                          // Call divide fn
  }
  CoolBignumE& result = *((CoolBignumE*) &q);   // same physical object
  return result;                                // shallow swap on return
}



// operator% -- overload modulus operator for CoolBignums
// Inputs:  references to CoolBignum number and modulus
// Outputs:  CoolBignum modulus result

CoolBignumE operator% (const CoolBignum& b1, const CoolBignum& b2) {
  CoolBignum q,r;                             // Temporary quotient and remainder
  if (b2.count == 0) {                // if b2 == 0
    r.state = N_DIVIDE_BY_ZERO;       //   implies division by zero
    b2.divide_by_zero("operator% ()"); //  raise an exception
  }
  else {                              // otherwise,
    divide(b1,b2,q,r);                //   divide b1 by b2 and save remainder
  }                                  
  CoolBignumE& result = *((CoolBignumE*) &r);   // same physical object
  return result;                                // shallow swap on return
}



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

CoolBignumE left_shift (const CoolBignum& b1, long 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.
  CoolBignum rslt;                                      // result of shift
  rslt.sign = b1.sign;                          // result follows sign of input
  Counter growth = Counter(l / data_bits_s);    // # of words rslt will grow by
  Data shift = Data(l % data_bits_s);           // amount to actually shift
  Data rshift = data_bits_s - shift;            // amount to shift next word by
  Data carry =                                  // value that will be shifted
    b1.data[b1.count - 1] >> (data_bits_s - 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);
  }
  CoolBignumE& result = *((CoolBignumE*) &rslt);// same physical object
  return result;                                // shallow swap on return
}



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

CoolBignumE right_shift (const CoolBignum& b1, long l) {
  CoolBignum rslt;                                      // result of shift
  Counter shrinkage = Counter(l / data_bits_s); // # of words rslt will shrink
  Data shift = Data(l % data_bits_s);           // 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 = data_bits_s - 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);
    }
  }
  CoolBignumE& result = *((CoolBignumE*) &rslt);// same physical object
  return result;                                // shallow swap on return
}