aboutsummaryrefslogtreecommitdiffstats
path: root/src/divide.cpp
blob: c7f3e1b424e18f4faf4dbc9722ac346bb6330dbf (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
/*************************************************
* Division Algorithm Source File                 *
* (C) 1999-2008 The Botan Project                *
*************************************************/

#include <botan/numthry.h>
#include <botan/mp_core.h>

namespace Botan {

namespace {

/*************************************************
* Handle signed operands, if necessary           *
*************************************************/
void sign_fixup(const BigInt& x, const BigInt& y, BigInt& q, BigInt& r)
   {
   if(x.sign() == BigInt::Negative)
      {
      q.flip_sign();
      if(r.is_nonzero()) { --q; r = y.abs() - r; }
      }
   if(y.sign() == BigInt::Negative)
      q.flip_sign();
   }

}

/*************************************************
* Solve x = q * y + r                            *
*************************************************/
void divide(const BigInt& x, const BigInt& y_arg, BigInt& q, BigInt& r)
   {
   if(y_arg.is_zero())
      throw BigInt::DivideByZero();

   BigInt y = y_arg;
   const u32bit y_words = y.sig_words();
   r = x;

   r.set_sign(BigInt::Positive);
   y.set_sign(BigInt::Positive);

   s32bit compare = r.cmp(y);

   if(compare < 0)
      q = 0;
   else if(compare ==  0)
      {
      q = 1;
      r = 0;
      }
   else
      {
      u32bit shifts = 0;
      word y_top = y[y.sig_words()-1];
      while(y_top < MP_WORD_TOP_BIT) { y_top <<= 1; ++shifts; }
      y <<= shifts;
      r <<= shifts;

      const u32bit n = r.sig_words() - 1, t = y_words - 1;

      q.get_reg().create(n - t + 1);
      if(n <= t)
         {
         while(r > y) { r -= y; ++q; }
         r >>= shifts;
         sign_fixup(x, y_arg, q, r);
         return;
         }

      BigInt temp = y << (MP_WORD_BITS * (n-t));

      while(r >= temp) { r -= temp; ++q[n-t]; }

      for(u32bit j = n; j != t; --j)
         {
         const word x_j0  = r.word_at(j);
         const word x_j1 = r.word_at(j-1);
         const word y_t  = y.word_at(t);

         if(x_j0 == y_t)
            q[j-t-1] = MP_WORD_MAX;
         else
            q[j-t-1] = bigint_divop(x_j0, x_j1, y_t);

         while(bigint_divcore(q[j-t-1], y_t, y.word_at(t-1),
                              x_j0, x_j1, r.word_at(j-2)))
            --q[j-t-1];

         r -= (q[j-t-1] * y) << (MP_WORD_BITS * (j-t-1));
         if(r.is_negative())
            {
            r += y << (MP_WORD_BITS * (j-t-1));
            --q[j-t-1];
            }
         }
      r >>= shifts;
      }

   sign_fixup(x, y_arg, q, r);
   }

}