Botan  1.11.15
src/lib/math/bigint/divide.cpp
Go to the documentation of this file.
00001 /*
00002 * Division Algorithm
00003 * (C) 1999-2007,2012 Jack Lloyd
00004 *
00005 * Botan is released under the Simplified BSD License (see license.txt)
00006 */
00007 
00008 #include <botan/divide.h>
00009 #include <botan/internal/mp_core.h>
00010 #include <botan/internal/mp_madd.h>
00011 
00012 namespace Botan {
00013 
00014 namespace {
00015 
00016 /*
00017 * Handle signed operands, if necessary
00018 */
00019 void sign_fixup(const BigInt& x, const BigInt& y, BigInt& q, BigInt& r)
00020    {
00021    if(x.sign() == BigInt::Negative)
00022       {
00023       q.flip_sign();
00024       if(r.is_nonzero()) { --q; r = y.abs() - r; }
00025       }
00026    if(y.sign() == BigInt::Negative)
00027       q.flip_sign();
00028    }
00029 
00030 bool division_check(word q, word y2, word y1,
00031                     word x3, word x2, word x1)
00032    {
00033    // Compute (y3,y2,y1) = (y2,y1) * q
00034 
00035    word y3 = 0;
00036    y1 = word_madd2(q, y1, &y3);
00037    y2 = word_madd2(q, y2, &y3);
00038 
00039    // Return (y3,y2,y1) >? (x3,x2,x1)
00040 
00041    if(y3 > x3) return true;
00042    if(y3 < x3) return false;
00043 
00044    if(y2 > x2) return true;
00045    if(y2 < x2) return false;
00046 
00047    if(y1 > x1) return true;
00048    if(y1 < x1) return false;
00049 
00050    return false;
00051    }
00052 
00053 }
00054 
00055 /*
00056 * Solve x = q * y + r
00057 */
00058 void divide(const BigInt& x, const BigInt& y_arg, BigInt& q, BigInt& r)
00059    {
00060    if(y_arg.is_zero())
00061       throw BigInt::DivideByZero();
00062 
00063    BigInt y = y_arg;
00064    const size_t y_words = y.sig_words();
00065 
00066    r = x;
00067    q = 0;
00068 
00069    r.set_sign(BigInt::Positive);
00070    y.set_sign(BigInt::Positive);
00071 
00072    s32bit compare = r.cmp(y);
00073 
00074    if(compare == 0)
00075       {
00076       q = 1;
00077       r = 0;
00078       }
00079    else if(compare > 0)
00080       {
00081       size_t shifts = 0;
00082       word y_top = y.word_at(y.sig_words()-1);
00083       while(y_top < MP_WORD_TOP_BIT) { y_top <<= 1; ++shifts; }
00084       y <<= shifts;
00085       r <<= shifts;
00086 
00087       const size_t n = r.sig_words() - 1, t = y_words - 1;
00088 
00089       if(n < t)
00090          throw Internal_Error("BigInt division word sizes");
00091 
00092       q.grow_to(n - t + 1);
00093 
00094       word* q_words = q.mutable_data();
00095 
00096       if(n <= t)
00097          {
00098          while(r > y) { r -= y; ++q; }
00099          r >>= shifts;
00100          sign_fixup(x, y_arg, q, r);
00101          return;
00102          }
00103 
00104       BigInt temp = y << (MP_WORD_BITS * (n-t));
00105 
00106       while(r >= temp) { r -= temp; q_words[n-t] += 1; }
00107 
00108       for(size_t j = n; j != t; --j)
00109          {
00110          const word x_j0  = r.word_at(j);
00111          const word x_j1 = r.word_at(j-1);
00112          const word y_t  = y.word_at(t);
00113 
00114          if(x_j0 == y_t)
00115             q_words[j-t-1] = MP_WORD_MAX;
00116          else
00117             q_words[j-t-1] = bigint_divop(x_j0, x_j1, y_t);
00118 
00119          while(division_check(q_words[j-t-1],
00120                               y_t, y.word_at(t-1),
00121                               x_j0, x_j1, r.word_at(j-2)))
00122             {
00123             q_words[j-t-1] -= 1;
00124             }
00125 
00126          r -= (q_words[j-t-1] * y) << (MP_WORD_BITS * (j-t-1));
00127 
00128          if(r.is_negative())
00129             {
00130             r += y << (MP_WORD_BITS * (j-t-1));
00131             q_words[j-t-1] -= 1;
00132             }
00133          }
00134       r >>= shifts;
00135       }
00136 
00137    sign_fixup(x, y_arg, q, r);
00138    }
00139 
00140 }